how to model kinetic reaction system

Question

Anastasya Silitonga 2020-5-24

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/532118-how-to-model-kinetic-reaction-system

评论： Sharath Ankathi 2024-5-2

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;

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sharath Ankathi 2024-5-2

@Anastasya Silitonga Can you share your complete solution for estimating kinetic parameters. Thanks in advance

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-5-24

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/532118-how-to-model-kinetic-reaction-system#answer_437648

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.

13 个评论
显示 11更早的评论隐藏 11更早的评论

Anastasya Silitonga 2020-5-25

在 MATLAB Online 中打开

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 2020-5-25

在 MATLAB Online 中打开

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.

Yasin Islam 2021-5-1

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 2021-5-1

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.

请先登录，再进行评论。

how to model kinetic reaction system

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

13 个评论
显示 11更早的评论隐藏 11更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

how to model kinetic reaction system

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

13 个评论 显示 11更早的评论隐藏 11更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论