how to model kinetic reaction system

Research area: Chemical kinetics

problem description: So this is parallel kinetics. I want to determine rate of reaction (as dCdt) and constant rate (as k) from experimental data that I have, time (as t) and concentration (as C). There are 6 variable called as mag,dag,tag,gly,ffa, and air.

so there are 6 constant rate and rate of reaction that I need to know, using ode23 or ode45. and the equations are described below:

dCdt(mag)=k1*Ctag*Cgly+2*k2*Cdag*Cgly-k3*Ctag*Cmag+k5*Cdag*Cair+k6*Cgly*Cffa;
dCdt(dag)=k1*Ctag*Cgly-k2*Cdag*Cgly+2*k3*Ctag*Cmag+k4*Ctag*Cair-k5*Cdag*Cair;
dCdt(tag)=-k1*Ctag*Cgly-k3*Ctag*Cmag-k4*Ctag*Cair;
dCdt(ffa)= k4*Ctag*Cair+k5*Cdag*Cair-k6*Cgly*Cffa;
dCdt(gly)=-k1*Ctag*Cgly-k2*Cdag*Cgly-k6*Cgly*Cffa;
dCdt(air)=-k4*Ctag*Cair-k5*Cdag*Cair+k6*Cgly*Cffa;

all I want to do is fit my data to this function to solve for k1,k2,k3,k4,k5,k6 and find each dCdt. I've tried make a code by using fminsearch and ode23, but I dont know how to define each concentration for ode so as finding each k. Can anybody help me with the code?

Thanks in advance.

function constantrate
clc
Tspan = [0.1 0.1]
constant=fminsearch(@rate, Tspan);

function constantrate = rate (k)
Data = [
60 90.5 5.2 1.1 10.5 2 3
120 90.6 5.1 1 10.1 1.8 3.3
180 90.7 5.1 1.4 12.5 1.9 4
240 90.7 5.4 2.1 10 1.7 4.5
];

t=Data(:,1); %time
Cmag=Data(:,2); %concentration of mag
Cdag=Data(:,3); %concentration of dag
Ctag=Data(:,4); %concentration of tag
Cgly=Data(:,5); %concentration of gly
Cair=Data(:,6); %concentration of air
Cffa=Data(:,7); %concentration of ffa

[tint, Cint]=ode23(@reaction,t,Cmag,[],k);
constantrate=sum((Cmag-Cint).^2);
plot (tint, Cint, 'r','o',t,Cmag,'b')
legend('calculation','experimental')
xlabel('time (min)')
ylabel('concentration (M)')

function dCdt=reaction(Ctag,Cmag,Cdag,Cgly,Cair,Cffa,k1,k2,k3,k4,k5,k6)
dCdt(mag)=k1*Ctag*Cgly+2*k2*Cdag*Cgly-k3*Ctag*Cmag+k5*Cdag*Cair+k6*Cgly*Cffa;
dCdt(dag)=k1*Ctag*Cgly-k2*Cdag*Cgly+2*k3*Ctag*Cmag+k4*Ctag*Cair-k5*Cdag*Cair;
dCdt(tag)=-k1*Ctag*Cgly-k3*Ctag*Cmag-k4*Ctag*Cair;
dCdt(ffa)= k4*Ctag*Cair+k5*Cdag*Cair-k6*Cgly*Cffa;
dCdt(gly)=-k1*Ctag*Cgly-k2*Cdag*Cgly-k6*Cgly*Cffa;
dCdt(air)=-k4*Ctag*Cair-k5*Cdag*Cair+k6*Cgly*Cffa;

Answer 1

[Star Strider]

There are several examples on fitting differential equations to data in Answers. See: ODE and Data fitting and Parameter Estimation for a System of Differential Equations. There are many others.

[Anastasya Silitonga]

Thank you so much for your reply and reference. But I'm not sure with the result, herewith the code that I use and the plot. I dont understand why all the theta seems not integrated, since I want to know the theta I choose number guess as it explain as theta0.. but as I choose theta0 with different guess, the result is changing.. it also happen for the concentration, the data in plot seems doesnt suit the data I input, do you have an idea to deal with that?

function anastasya
clc

function C=kinetics(theta, t)
c0=[1;0;0;0;0;0];
[~, Cv]=ode45(@DifEq,t,c0);

%theta as constant rate or k
%mag = 1, dag=2, tag=3, ffa=4, gly=5, air=6
function dC=DifEq(~,c)
dCdt=zeros(4,1);
dCdt(1)=theta(1)*c(3)*c(5)+2*theta(2)*c(2)*c(5)-theta(3)*c(3)*c(1)+theta(5)*c(2)*c(6)+theta(6)*c(5)*c(4);
dCdt(2)=theta(1)*c(3)*c(5)-theta(2)*c(2)*c(5)+2*theta(3)*c(3)*c(1)+theta(4)*c(3)*c(6)-theta(5)*c(2)*c(6);
dCdt(3)=-theta(1)*c(3)*c(5)-theta(3)*c(3)*c(1)-theta(4)*c(3)*c(6);
dCdt(4)= theta(4)*c(3)*c(6)+theta(5)*c(2)*c(6)-theta(6)*c(5)*c(4);
dCdt(5)=-theta(1)*c(3)*c(5)-theta(2)*c(2)*c(5)-theta(6)*c(5)*c(4);
dCdt(6)=-theta(4)*c(3)*c(6)-theta(5)*c(2)*c(6)+theta(6)*c(5)*c(4);
dC=dCdt;
end

C=Cv;
end


t=[
0.5
1
1.5
2
2.5
3
3.5
4
];

%row: mag = 1, dag=2, tag=3, ffa=4, gly=5, air=6
c=[
0.026 0.037 1.003 0.016 0.718 0.054
0.013 0.045 1.003 0.016 0.718 0.054
0.039 0.134 0.919 0.049 0.761 0.052
0.078 0.238 0.825 0.066 0.811 0.049
0.259 0.416 0.606 0.131 0.926 0.042
0.621 0.461 0.386 0.262 1.041 0.035
1.112 0.416 0.198 0.328 1.139 0.029
1.293 0.401 0.178 0.197 1.150 0.028
];
theta0=[0.1;1;1;1;1;1];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)', 'C_6(t)','Location','N')

end

[Star Strider]

My pleasure!

I let the ga (genetic algorithm) function estimate these parameters. Some fit reasonably well and some were significantly off.

The best estimate appears to be:

	Rate Constants:
		Theta(1) =  0.00276
		Theta(2) =  1.05513
		Theta(3) =  2.79222
		Theta(4) =  0.24273
		Theta(5) =  4.25839
		Theta(6) =  0.06654

Initial Conditions for compartments 1-6 (in order):

		Theta(7) =  0.00528
		Theta(8) =  0.00931
		Theta(9) =  0.60556
		Theta(10) =  0.00284
		Theta(11) =  1.01283
		Theta(12) =  0.42944

using this code:

function C=kinetics(theta, t)
% c0=[1;0;0;0;0;0];
c0 = theta(7:end);
[~, Cv]=ode45(@DifEq,t,c0);
%theta as constant rate or k
%mag = 1, dag=2, tag=3, ffa=4, gly=5, air=6
function dC=DifEq(~,c)
dCdt=zeros(6,1);
dCdt(1)=theta(1)*c(3)*c(5)+2*theta(2)*c(2)*c(5)-theta(3)*c(3)*c(1)+theta(5)*c(2)*c(6)+theta(6)*c(5)*c(4);
dCdt(2)=theta(1)*c(3)*c(5)-theta(2)*c(2)*c(5)+2*theta(3)*c(3)*c(1)+theta(4)*c(3)*c(6)-theta(5)*c(2)*c(6);
dCdt(3)=-theta(1)*c(3)*c(5)-theta(3)*c(3)*c(1)-theta(4)*c(3)*c(6);
dCdt(4)= theta(4)*c(3)*c(6)+theta(5)*c(2)*c(6)-theta(6)*c(5)*c(4);
dCdt(5)=-theta(1)*c(3)*c(5)-theta(2)*c(2)*c(5)-theta(6)*c(5)*c(4);
dCdt(6)=-theta(4)*c(3)*c(6)-theta(5)*c(2)*c(6)+theta(6)*c(5)*c(4);
dC=dCdt;
end
C=Cv;
end


t=[0.5
1
1.5
2
2.5
3
3.5
4];
%row: mag = 1, dag=2, tag=3, ffa=4, gly=5, air=6
c=[0.026 0.037 1.003 0.016 0.718 0.054
0.013 0.045 1.003 0.016 0.718 0.054
0.039 0.134 0.919 0.049 0.761 0.052
0.078 0.238 0.825 0.066 0.811 0.049
0.259 0.416 0.606 0.131 0.926 0.042
0.621 0.461 0.386 0.262 1.041 0.035
1.112 0.416 0.198 0.328 1.139 0.029
1.293 0.401 0.178 0.197 1.150 0.028];


ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 12;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-6, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
% [B,fval,exitflag,output,population,scores] = ga(ftns, 53, [],[],[],[],zeros(53,1),Inf(53,1),[],[],opts)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
GA_hms = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSSSSS');
fprintf('\nElapsed Time: %s\n', GA_hms)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

% tv = linspace(min(t), max(t));
tv = linspace(0, max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
lgdc = sprintfc('C_%d(t)',1:size(c,2));
legend(hlp, lgdc, 'Location','N')

There could be a problem in the way you coded the kinetic equations, so please check those, particularly to be certain they are in the correct order. (I would generally expect a much better fit.) It could also be that they are coded correctly, and simply do not reflect the actual dynamics of the system you are modeling.

[Anastasya Silitonga]

Thanks a lot Star! I'm sorry to disturb you again, but I had tried to use that code, and why I still get an error ya?

This is what command window say:

Error using optimoptions (line 114)
Invalid solver specified. Provide a solver name or handle (such as 'fmincon' or @fminunc).
Type DOC OPTIMOPTIONS for a list of solvers.

Error in anastasya (line 46)
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-6,
'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

[Star Strider]

You need to have the Global Optimization Toolbox to use the ga function. I am running 2020a, and the code I posted (for reference) ran without error.

My objective otherwise was to present a set ot parameters (and slightly revised code) that you can use as initial estimates with fminsearch or other optimization functions to estimate the parameters. Those parameters should give a reasonable fit.

I still recommend that you be certain that the order and coding of the equations in ‘kinetics’ are correct. I generally expect a much better fit.

[Anastasya Silitonga]

Thanks again! this is new for me and I learn a lot! Do you mind if I ask again? With these equations, is it possible to find the best order and coefficient of rate at the same time? I understand that it could happen by using fminsearch and ode23 for single equation, but I don't understand since it contains so much parallel equations.

[Star Strider]

My pleasure!

This finds the best coefficients of the differential equation that minimise the norm or least-square difference between the differential equation and the data. I do not know the system you are modeling (it appears to be lipolysis or something similar), so I do not know what the parameters mean. As I remember (from undergraduate physical chemistry back when log tables and slide rules were all the rage, and computers were behemoths that ate punchcards), the order would be given by the exponents of the reactions in kinetic equations, and I do not see any exponents here (so I assume everything is first-order in the respective reactants).

So if you want to optimise or estimate the order and rate, and the differential equations are set up to do that, optimising the fit will estimate the parameters you want.

[Daniel Zhou]

Hi, Star. I'm new to Matlab. And I just copied and pasted the code you posted above. Why is Matlab saying"Error: File: kinetics.m Line: 23 Column: 1

This statement is not inside any function.

(It follows the END that terminates the definition of the function "kinetics".)". Could you please help me with this?

[Ryan Patrick]

Hi, Zhou, I'm a new starter to Matlab too. Thanks to Stars' code firstly. I suggest you could extract line 23 till the end line, and put them into a new .m file ( in the same folder), I encountered the same problem as you and solved it by this method. better solutions are welcome.

[Star Strider]

It is necessary to put them inside their own function, as I did in Parameter Estimation for a System of Differential Equations with ‘Igor_Moura.m’. That is simply the most efficient way to code them so that everything works across several MATLAB releases/versions. In recent releases it has been possible to put the ‘kinetics’ function at the end of a script (avoiding the need to put everything in a function), however the current appropach works in releases previous to that so I continue to use it for similar problems.

[Yasin Islam]

Hello Star Strider, i used a lot of your code and you already helped me a lot with my kinetic modeling!

I still have a question about modeling of chemical kinetics, I'm using different (or better all ODE solvers) and lsqcurvefit to fit my data (starting from 0 for my theta0 values). I also use the GA function in another modeling file also with different ODE solvers. I'm compareing all the models by their RSS or Rsdnrm values to find the best one. Is there another reasonable route for modeling kinetic data? Does it make sense to use the theta values from GA und insert them as theta0 values for lsqcurvefit? So far the lsqcurvefit with theta0=0 give the best RSS values, followed bei lsqcurvefit with theta0 from GA and lastly GA itself.

Thank you very much!

[Star Strider]

Thank you!

‘Is there another reasonable route for modeling kinetic data?’

I believe the Bioinformatics Toolbox (that I don’t have) has specific ODE solvers for them. I’m not aware of any other options.

‘Does it make sense to use the theta values from GA und insert them as theta0 values for lsqcurvefit?’

Using the parameters that ga produces as initial estimates for lsqcurvefit can result in slightly better parpameter estimates. (I believe that there’s an option to do this automatically in ga, however I routinely forget to enable it.)

.

[Yasin Islam]

Thanky you very much for your fast response.

I do have bioinformatics toolbox, i'll take a look for other ode solvers, thanks for the hint.

I just copy the values after GA into lsqcurvefit, but unfortunately the outcome is not as good as starting from theta0=0 instead of GA values.

Once again, thank you very much! I had and propably still have no clue about matlab and was able to write a working code using your code from several posts.

[Star Strider]

My pleasure!

I am somewhat surprised that the parameters estimated by ga as initial paremeter estimates for lsqcurvefit do not result in a better result than other initial estimates.