
        format none;

        function xdot = f (x, t)
            # Zmena hodnoty x za jednotku casu
            xdot = zeros (58,1);
            xdot(1) = 0; # fixed, cytosol 
xdot(2) = 0; # fixed, k4 
xdot(3) = 0; # fixed, k3 
xdot(4) = 0; # fixed, k2 
xdot(5) = 0; # fixed, k5 
xdot(6) = 0; # fixed, k1 
xdot(7) = 0; # fixed, k6 
xdot(8) = 0; # fixed, k7 
xdot(9) = 0; # fixed, k8 
xdot(10) = 0; # fixed, k9 
xdot(11) = 0; # fixed, k10 
xdot(12) = 0; # fixed, k11 
xdot(13) = 0; # fixed, k12 
xdot(14) = 0; # fixed, k21 
xdot(15) = 0; # fixed, k22 
xdot(16) = 0; # fixed, k23 
xdot(17) = 0; # fixed, k24 
xdot(18) = 0; # fixed, ka1 
xdot(19) = 0; # fixed, ka2 
xdot(20) = 0; # fixed, kcat1 
xdot(21) = 0; # fixed, kcat2 
xdot(22) = 0; # fixed, kcat3 
xdot(23) = 0; # fixed, kcat4 
xdot(24) = 0; # fixed, kcat_b1 
xdot(25) = 0; # fixed, kcat_b2 
xdot(26) = 0; # fixed, kmdeg1 
xdot(27) = 0; # fixed, kmdeg2 
xdot(28) = 0; # fixed, kpdeg1 
xdot(29) = 0; # fixed, kpdeg2 
xdot(30) = 0; # fixed, kpdeg3 
xdot(31) = 0; # fixed, kpdeg4 
xdot(32) = 0; # fixed, ktl1 
xdot(33) = 0; # fixed, ktl2 
xdot(34) = 0; # fixed, kbts1 
xdot(35) = 0; # fixed, kbts2 
xdot(36) = 0; # fixed, Km1 
xdot(37) = 0; # fixed, Km2 
xdot(38) = 0; # fixed, Km3 
xdot(39) = 0; # fixed, Km4 
xdot(40) = 0; # fixed, Km_b1 
xdot(41) = 0; # fixed, Km_b2 
xdot(42) = 0; # fixed, RNAP 
xdot(43) = 0; # fixed, L_ 
xdot(57) = x(57) = x(47)/(x(47)+x(48)); # assignment, KaiC_rel  FROM INIT EXPRESSION!!! 
xdot(58) = x(58) = x(48)/(x(47)+x(48)); # assignment, PKaiC_rel  FROM INIT EXPRESSION!!! 
xdot(44) =  - (x(14)*x(44))  - (x(6)*x(44))  + (x(4)*x(47)*x(47)*x(47)*x(47)*x(47)*x(47))  - ((x(20)*x(52)*x(44))/(x(36) + x(44)))  + (x(15)*x(45))  + ((x(21)*x(56)*x(45))/(x(37) + x(45))) ; # reaction, KaiC6 
xdot(45) =  + (x(17)*x(46))  + ((x(23)*x(56)*x(46))/(x(39) + x(46)))  + 0 * (x(18)*(x(34)*x(42))/(1+x(34)*x(42))*(x(46)/x(45))*x(43))  + 0 * ((x(24)*x(45)*x(50))/(x(40) + x(50)))  + 0 * (x(19)*(x(35)*x(42))/(1+x(35)*x(42))*(x(46)/x(45))*x(43))  + (x(14)*x(44))  + ((x(20)*x(52)*x(44))/(x(36) + x(44)))  - (x(15)*x(45))  - (x(16)*x(45))  - ((x(21)*x(56)*x(45))/(x(37) + x(45)))  - (x(3)*x(45))  + (x(2)*x(47)*x(47)*x(47)*x(48)*x(48)*x(48))  - ((x(22)*x(52)*x(45))/(x(38) + x(45))) ; # reaction, PPKaiC6 
xdot(46) =  - (x(17)*x(46))  - ((x(23)*x(56)*x(46))/(x(39) + x(46)))  - (x(5)*x(46))  + (x(7)*x(48)*x(48)*x(48)*x(48)*x(48)*x(48))  + 0 * (x(18)*(x(34)*x(42))/(1+x(34)*x(42))*(x(46)/x(45))*x(43))  + 0 * ((x(25)*x(46)*x(49))/(x(41) + x(49)))  + 0 * (x(19)*(x(35)*x(42))/(1+x(35)*x(42))*(x(46)/x(45))*x(43))  + (x(16)*x(45))  + ((x(22)*x(52)*x(45))/(x(38) + x(45))) ; # reaction, CPKaiC6 
xdot(47) =  + (x(33)*x(43)*x(55))  +  6 * (x(6)*x(44))  -  6 * (x(4)*x(47)*x(47)*x(47)*x(47)*x(47)*x(47))  - (x(30)*x(47)*x(43))  +  3 * (x(3)*x(45))  -  3 * (x(2)*x(47)*x(47)*x(47)*x(48)*x(48)*x(48)) ; # reaction, KaiC 
xdot(48) =  +  6 * (x(5)*x(46))  -  6 * (x(7)*x(48)*x(48)*x(48)*x(48)*x(48)*x(48))  +  3 * (x(3)*x(45))  -  3 * (x(2)*x(47)*x(47)*x(47)*x(48)*x(48)*x(48))  - (x(31)*x(48)*x(43)) ; # reaction, PKaiC 
xdot(49) =  - ((x(25)*x(46)*x(49))/(x(41) + x(49)))  + ((x(24)*x(45)*x(50))/(x(40) + x(50))) ; # reaction, KaiB4i 
xdot(50) =  + (x(13)*x(56))  - (x(12)*1*x(50)*x(52))  + ((x(25)*x(46)*x(49))/(x(41) + x(49)))  - ((x(24)*x(45)*x(50))/(x(40) + x(50)))  - (x(11)*x(50))  + (x(10)*x(51)*x(51)*x(51)*x(51)) ; # reaction, KaiB4 
xdot(51) =  + (x(33)*x(43)*x(55))  - (x(29)*x(51)*x(43))  +  4 * (x(11)*x(50))  -  4 * (x(10)*x(51)*x(51)*x(51)*x(51)) ; # reaction, KaiB 
xdot(52) =  + (x(13)*x(56))  - (x(12)*1*x(50)*x(52))  - (x(8)*x(52))  + (x(9)*x(53)*x(53))  + 0 * ((x(20)*x(52)*x(44))/(x(36) + x(44)))  + 0 * ((x(22)*x(52)*x(45))/(x(38) + x(45))) ; # reaction, KaiA2 
xdot(53) =  +  2 * (x(8)*x(52))  -  2 * (x(9)*x(53)*x(53))  - (x(28)*x(53)*x(43))  + (x(32)*x(43)*x(54)) ; # reaction, KaiA 
xdot(54) =  + (x(18)*(x(34)*x(42))/(1+x(34)*x(42))*(x(46)/x(45))*x(43))  + 0 * (x(32)*x(43)*x(54))  - (x(26)*x(54)) ; # reaction, kaiA_mRNA 
xdot(55) =  + 0 * (x(33)*x(43)*x(55))  - (x(27)*x(55))  + (x(19)*(x(35)*x(42))/(1+x(35)*x(42))*(x(46)/x(45))*x(43)) ; # reaction, kaiBC_mRNA 
xdot(56) =  + 0 * ((x(23)*x(56)*x(46))/(x(39) + x(46)))  - (x(13)*x(56))  + (x(12)*1*x(50)*x(52))  + 0 * ((x(21)*x(56)*x(45))/(x(37) + x(45))) ; # reaction, KaiA2B4 
xdot(57) = 0; # end
xdot(58) = 0; # end

        endfunction

        # time scale
        t = linspace(0, 120, 60);

        # default values
        x0 = zeros (58,1);
        x0(1) = 1; # fixed compartment, cytosol
x0(2) = 1.019e-14; # fixed parameter, k4
x0(3) = 0.0001615; # fixed parameter, k3
x0(4) = 2.039e-16; # fixed parameter, k2
x0(5) = 0.162; # fixed parameter, k5
x0(6) = 1.615; # fixed parameter, k1
x0(7) = 1.019e-10; # fixed parameter, k6
x0(8) = 0.162; # fixed parameter, k7
x0(9) = 0.268; # fixed parameter, k8
x0(10) = 7.393e-17; # fixed parameter, k9
x0(11) = 0.0001615; # fixed parameter, k10
x0(12) = 0.0008756; # fixed parameter, k11
x0(13) = 0.08788; # fixed parameter, k12
x0(14) = 1.079e-08; # fixed parameter, k21
x0(15) = 1.079e-05; # fixed parameter, k22
x0(16) = 1.079e-06; # fixed parameter, k23
x0(17) = 1.079e-08; # fixed parameter, k24
x0(18) = 10170000; # fixed parameter, ka1
x0(19) = 64580000; # fixed parameter, ka2
x0(20) = 0.539; # fixed parameter, kcat1
x0(21) = 0.539; # fixed parameter, kcat2
x0(22) = 1.079; # fixed parameter, kcat3
x0(23) = 0.89; # fixed parameter, kcat4
x0(24) = 2.423; # fixed parameter, kcat_b1
x0(25) = 0.346; # fixed parameter, kcat_b2
x0(26) = 0.133; # fixed parameter, kmdeg1
x0(27) = 0.178; # fixed parameter, kmdeg2
x0(28) = 0.008; # fixed parameter, kpdeg1
x0(29) = 0.49; # fixed parameter, kpdeg2
x0(30) = 1.3; # fixed parameter, kpdeg3
x0(31) = 0.2; # fixed parameter, kpdeg4
x0(32) = 0.008239; # fixed parameter, ktl1
x0(33) = 170.1; # fixed parameter, ktl2
x0(34) = 3.657e-12; # fixed parameter, kbts1
x0(35) = 1e-12; # fixed parameter, kbts2
x0(36) = 602; # fixed parameter, Km1
x0(37) = 602; # fixed parameter, Km2
x0(38) = 0.602; # fixed parameter, Km3
x0(39) = 0.602; # fixed parameter, Km4
x0(40) = 0.602; # fixed parameter, Km_b1
x0(41) = 66.75; # fixed parameter, Km_b2
x0(42) = 5000; # fixed parameter, RNAP
x0(43) = 1; # fixed parameter, L_
x0(44) = 139.22; # fixed specie, KaiC6
x0(45) = 779.158; # fixed specie, PPKaiC6
x0(46) = 1229.563; # fixed specie, CPKaiC6
x0(47) = 932.446; # fixed specie, KaiC
x0(48) = 110.829; # fixed specie, PKaiC
x0(49) = 3251.369; # fixed specie, KaiB4i
x0(50) = 0.173; # fixed specie, KaiB4
x0(51) = 1130; # fixed specie, KaiB
x0(52) = 166.559; # fixed specie, KaiA2
x0(53) = 9.998; # fixed specie, KaiA
x0(54) = 2.856; # fixed specie, kaiA_mRNA
x0(55) = 2.865; # fixed specie, kaiBC_mRNA
x0(56) = 51.022; # fixed specie, KaiA2B4
x0(57) = x0(47)/(x0(47)+x0(48)); # assignment parameter, KaiC_rel
x0(58) = x0(48)/(x0(47)+x0(48)); # assignment parameter, PKaiC_rel


        #lsode_options ( 'absolute tolerance', max (50*eps, 0.5e-28) );
        lsode_options ( 'relative tolerance', max (50*eps, 0.5e-28) ); 

        y = lsode ("f", x0, t);
        
        for i = 1:length(t)
            y(i, 57) = y(i, 47)/(y(i, 47)+y(i, 48)); # assignment, KaiC_rel  FROM INIT EXPRESSION!!! 
y(i, 58) = y(i, 48)/(y(i, 47)+y(i, 48)); # assignment, PKaiC_rel  FROM INIT EXPRESSION!!! 

        endfor
        
        [t',y]
        