Fitting kinetic model to estimate kinetic parameter

Hello, I want to ask a question. I want to curve fit the data with the kinetic model. I am able to run the code (written and attached below), but the result is still far from my expectation.

function trycurvefitlunaflores 
  
function C=kinetics(theta,t) 
c0=[0.6178; 28.90; 0.4589]; 
[T,Cv]=ode45(@DifEq,t,c0); 

    function dC=DifEq(t,c) 
        dcdt=zeros(3,1); 
        mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2)^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5)); 
        dcdt(1) = mu.*c(1); 
        dcdt(2) = -c(1).*((mu./(theta(6))+theta(7))); 
        dcdt(3) = c(1).*((theta(8).*mu)+theta(9)); 
        dC=dcdt; 
    end 

C=Cv; 

end 

t = [0 
12 
24 
36 
48 
60 
72]; 

c = [0.6178489703 28.90160183 0.7586206897 
0.823798627 28.83295195 4.045977011 
2.402745995 21.62471396 6.827586207 
5.972540046 13.04347826 18.5862069 
8.306636156 6.178489703 34.01149425 
9.885583524 2.402745995 44.88505747 
9.542334096 2.059496568 50.44827586]; 

theta0 = [1;1;1;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)', 'Location','N') 

end 

In the normal circumstance, the plot line (from kinetic model) should follow the data's tren, like this figure below.

Could you help me to fix the plot to make it more similar to the data? Any little advice or help will be much appreciated. Thank you.

[Alex Sha]

Hi, Salman, the results below are obtained by using a package other than Matlab, you may have a check:

Root of Mean Square Error (RMSE): 0.732135213542845
Sum of Squared Residual: 9.64839547636969
Correlation Coef. (R): 0.997606342175452
R-Square: 0.995218413948685

Parameter          	Best Estimate
--------------------	-------------
theta1	-0.0242147478503888
theta2	-12.5872182600746
theta3	-28.8723626819411
theta4	50.506381835912
theta5	0.152709651271307
theta6	0.348919829234199
theta7	0.00280925362885251
theta8	1.98890080928454
theta9	0.0951796615605864

[Nadya Shalsabila salman]

Hi @Alex Sha! thank you for the answer.. the curve looks great! may I know what package are you using and also can I have the full code? Thank you again

[Alex Sha]

1stOpt, an package with great advance features on global optimization, the code is simple as fellow:

ConstStr  mu = ((theta1*c2)/(theta2+c2+(c2^2/(theta3))))*((1-c3/theta4)^theta5);
Parameter theta(9);
Variable t,c(3);
ODEFunction
c1' = mu*c1;
c2' = -c1*((mu/(theta6)+theta7));
c3' = c1*((theta8*mu)+theta9);
Data;
0	0.6178489703	28.90160183	0.7586206897
12	0.823798627	28.83295195	4.045977011
24	2.402745995	21.62471396	6.827586207
36	5.972540046	13.04347826	18.5862069
48	8.306636156	6.178489703	34.01149425
60	9.885583524	2.402745995	44.88505747
72	9.542334096	2.059496568	50.44827586

[Star Strider]

@Nadya Shalsabila salman — For what it’s worth, I gave this my best shot with a number of Global Optimization Toolbox functions over a couple days, and could never get a good fit to with any of them, although I had no problems with and . (I have no idea how ‘1stOpt’ solves these problems. Perhaps that information can be shared.)

[Nadya Shalsabila salman]

Hi @Star Strider! Thanks for trying the problem I had. if you don't mind, may i see the code you have tried?

[Nadya Shalsabila salman]

Is it possible that this difference is due to different units of the two plots? because if I get a suitable plot, the parameter value is very far from what it should be

the substrate, product, and biomass data should be as follows

but in the plotting, I used the initial product (karotenoid) = 0.75, not 75.86 mu g/L

and the parameter should be like this

and ms 0.005 g g-1 h-1

I'm still confused about this unit difference

[Star Strider]

Parameter Estimation for a System of Differential Equations

.

Answer 1

[Alex]

Using fminsearch, defining a custom cost-function and using the initial values of Alex Sha finally did the trick.

I seperated the program into seperate files. The main file:

clear

close all

clc

global t c

t = [0

12

24

36

48

60

72];

c = [0.6178489703 28.90160183 0.7586206897

0.823798627 28.83295195 4.045977011

2.402745995 21.62471396 6.827586207

5.972540046 13.04347826 18.5862069

8.306636156 6.178489703 34.01149425

9.885583524 2.402745995 44.88505747

9.542334096 2.059496568 50.44827586];

theta0 = [-0.024;-12.55;-28.80;50.63;0.168;0.3501;0.003146;1.9802;0.0954];

options = optimset('display','iter');

options.MaxFunEvals = 10e8;

options.TolX = 10e-5;

options.TolFun = 10e-5;

options.MaxIter = 10e3;

[theta] = fminsearch(@minfunction, theta0, options);

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv, c(1,:));

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)', 'Location','N')

the kinetics function:

function C=kinetics(theta, t, c0)

[T,Cv]=ode45(@DifEq,t, c0);

function dC=DifEq(t,c)

dcdt=zeros(3,1);

mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2)^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5));

dcdt(1) = mu.*c(1);

dcdt(2) = -c(1).*((mu./(theta(6))+theta(7)));

dcdt(3) = c(1).*((theta(8).*mu)+theta(9));

dC=dcdt;

end

C=Cv;

end

and the cost function:

function error = minfunction(theta)

global t c

c0 = c(1,:);

c_estim = kinetics(theta,t,c0);

w1 = 1;

w2 = 1;

w3 = 1;

error_temp = w1*abs(c(:,1)-c_estim(:,1)) + ...

w2*abs(c(:,2)-c_estim(:,2)) + ...

w3*abs(c(:,3)-c_estim(:,3)) ;

error = sum(error_temp);

end

In the cost function one could use different weights for the error of c1, c2 and c3. This could be of interest if one would worry more about the error for example of c1 than c2.

fminsearch now tries to find parameters for the differential equation which minimizes:

Hope this helps.

[Nadya Shalsabila salman]

Hi @Alex, thank you very much for your answer! I wish you always health and happiness :)

[Alex]

You are very welcome. Would you mind marking the answer as "accepted" if it solved your problem?

[Nadya Shalsabila salman]

Sure! @Alex

[Nadya Shalsabila salman]

Hi @Alex, I'm sorry to bother you again.. Could you explain to me what the value 1, 12, and 20 in error code? The result of error like this picture below. I don't understand why the R2 is not close to 1

[Alex]

The idea behind the weights was to give you the option to say: "i worry more about the error of c3, than c2 and therefore want to give the error a higher weight".

Looking again at the cost function makes me think about an error I made there. Maybe it is better to define the cost function like:

function error = minfunction(theta)

global t c

c0 = c(1,:);

c_estim = kinetics(theta,t,c0);

w1 = 1;

w2 = 1;

w3 = 1;

error_temp = w1*abs(c(:,1)-c_estim(:,1)) + ...

w2*abs(c(:,2)-c_estim(:,2)) + ...

w3*abs(c(:,3)-c_estim(:,3)) ;

error = sum(error_temp);

end

The error is now defined the following way:

(Before it would have been possible that errors cancel out each other, which could have been the reason why I had to add the weights. This is of course not the prefered way.). As you can see I now have set the weights to 1, since it looks as if they are not really necessary anymore. But if you want to play around with the weights, you are of course free to go if it improves the results.

I will update the original answer.

[Nadya Shalsabila salman]

Okay.. thank you very much @Alex. Hereby I attach the latest code I have tried, this code is using the latest error code you provided and the result is pretty decent.

Answer 2

[Alex]

Your basic structure seems to be ok. Here are some hints you could consider:

• Try playing around with the solver options, like using different algorithms, different initial points and tolerances:

theta0 = [1;1;1;1;1;1;1;1;1];

options = optimoptions('lsqcurvefit', 'Algorithm','trust-region-reflective');

options.MaxFunctionEvaluations = 10e5;

options.StepTolerance = 10e-7;

options.FunctionTolerance = 10e-7;

options.MaxIterations = 10e5;

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

• Use the accurate intitial condition for your differential equation:

c0=[0.6178489703 28.90160183 0.7586206897];

• Watch out for local minimum warnings. Matlab provides further thoughts on how to avoid local minimum convergence when clicking on the link in the matlab console (or directly via this link):

>> trycurvefitlunaflores

Local minimum possible.

After playing around a bit I was able to produce the following result:

I am sure that with a little tweaking you will be able to also fit C_1 correctly. Let me know if you found a solution.

[Nadya Shalsabila salman]

hi @Alex Sha and @aleper, thank you for the advice! Yes, after I try and found a solution, I will attached the code in this page. Once again thank's so much and have a nice day!

[Nadya Shalsabila salman]

Hi again! So I have plot again the kinetic model.. In this code below I added code to create two y-axis, on the right and left. but the result (plot and kinetic parameter) is still very far from the desired. Could you help me to solve my problem, please?

function terbaru
  
function C=kinetics(theta,t) 
c0=[0.6178; 28.90; 0.7586]; 
[T,Cv]=ode45(@DifEq,t,c0); 

    function dC=DifEq(t,c) 
        dcdt=zeros(3,1); 
        mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2)^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5)); 
        dcdt(1) = mu.*c(1); 
        dcdt(2) = -c(1).*((mu./(theta(6))+theta(7))); 
        dcdt(3) = c(1).*((theta(8).*mu)+theta(9)); 
        dC=dcdt; 
    end 

C=Cv; 

end 

t = [0 
12 
24 
36 
48 
60 
72]; 

c = [0.6178489703 28.90160183 0.7586206897 
0.823798627 28.83295195 4.045977011 
2.402745995 21.62471396 6.827586207 
5.972540046 13.04347826 18.5862069 
8.306636156 6.178489703 34.01149425 
9.885583524 2.402745995 44.88505747 
9.542334096 2.059496568 50.44827586]; 

theta0 = [1;1;1;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
yyaxis left
plot(t,c(:,[1 2]),'p')
hold on
hlp=plot(tv,Cfit(:,[1 2]));
hold off
ylabel('Concentration (g/L)')
ylim([min(c(:,1)) max(ylim)])
yyaxis right
plot(t,c(:,3),'p')
hold on
hlpp=plot(tv,Cfit(:,3));
hold off
ylabel('Concentration (\mug/L)')
ylim([min(c(:,3)) max(ylim)])
grid
legend([hlp([1 2]); hlpp],'C_1','C_2','C_3', 'Location','')


end 

[Alex Sha]

Hi, you may take the results I provided above as initial start values in your Matlab code, and see what's happen

[Nadya Shalsabila salman]

Hi @Alex Sha, I have tried the initials and the algorithm you suggested. So far the results are better than before, thank you very much. but I don't understand why C1's graph is not plotted well. So this my new code and the result

function terbaru
  
function C=kinetics(theta,t) 
c0=[0.6178489703; 28.90160183; 0.7586206897]; 
[T,Cv]=ode45(@DifEq,t,c0); 

    function dC=DifEq(t,c) 
        dcdt=zeros(3,1); 
        mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2)^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5)); 
        dcdt(1) = mu.*c(1); 
        dcdt(2) = -c(1).*((mu./(theta(6))+theta(7))); 
        dcdt(3) = c(1).*((theta(8).*mu)+theta(9)); 
        dC=dcdt; 
    end 

C=Cv; 

end 

t = [0 
12 
24 
36 
48 
60 
72]; 

c = [0.6178489703 28.90160183 0.7586206897 
0.823798627 28.83295195 4.045977011 
2.402745995 21.62471396 6.827586207 
5.972540046 13.04347826 18.5862069 
8.306636156 6.178489703 34.01149425 
9.885583524 2.402745995 44.88505747 
9.542334096 2.059496568 50.44827586]; 

theta0 = [1;1;1;1;1;1;1;1;1]; 

options = optimoptions('lsqcurvefit', 'Algorithm','trust-region-reflective');
options.MaxFunctionEvaluations = 10e5;
options.StepTolerance = 10e-7;
options.FunctionTolerance = 10e-7;
options.MaxIterations = 10e5;
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c, [],[], options)

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
yyaxis left
plot(t,c(:,[1 2]),'p')
hold on
hlp=plot(tv,Cfit(:,[1 2]));
hold off
ylabel('Concentration (g/L)')
ylim([min(c(:,1)) max(ylim)])
yyaxis right
plot(t,c(:,3),'p')
hold on
hlpp=plot(tv,Cfit(:,3));
hold off
ylabel('Concentration (\mug/L)')
ylim([min(c(:,3)) max(ylim)])
grid
legend([hlp([1 2]); hlpp],'C_1','C_2','C_3', 'Location','')


end 

[Alex Sha]

Hi, taking:

theta0 = [1;1;1;1;1;1;1;1;1];

as:

theta0 = [-0.024;-12.55;-28.80;50.63;0.168;0.3501;0.003146;1.9802;0.0954];

run and see the result

[Nadya Shalsabila salman]

Wow, the results are very surprising.. the graphics look very good.. but what value is put in this tetha? after running the code, I also still get a warning like this:

Local minimum possible.

lsqcurvefit stopped because the size of the current step is less than

the value of the step size tolerance.

Warning: Imaginary parts of complex X and/or Y arguments ignored.

> In latest (line 61)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

> In recent (line 68)

and the parameter (theta) still very far from the value it should be

Value should be:

Miumax=0.145;

Ks=1.6607;

Ki=50.3151;

Kp=8798.15;

gamma=3.6658;

Yxs=0.3941;

alfa=371.45;

betha=1.5502;

ms=0.005;

and value i got from this simulation:

Theta(1) = -0.02189

Theta(2) = -12.84832

Theta(3) = -30.56369

Theta(4) = 48.31902

Theta(5) = 0.06129

Theta(6) = 0.36215

Theta(7) = 0.00873

Theta(8) = 2.14285

Theta(9) = 0.08252

[Alex Sha]

if paameters of [Miumax,Ks,Ki,Kp,gamma,Yxs,alfa,betha,ms] correspond to [theta1,theta2,...,theta9], then you could take:

Miumax=0.145;

Ks=1.6607;

Ki=50.3151;

Kp=8798.15;

gamma=3.6658;

Yxs=0.3941;

alfa=371.45;

betha=1.5502;

ms=0.005;

into your ODE function to test and verify the performance, as my test, the results are much bad.

[Nadya Shalsabila salman]

Hi @Alex Sha, thank you very much for your answer! Anyway, I want to try 1stOpt Software to estimate the parameter. But, when I request the link, the software informer tell me "The download link for requested software cannot be found. We will update our database soon and let you know when 1stOpt is available for download"

Can you tell me how to download the 1st Opt?

[Alex Sha]

Hi, having a look at: 1stOpt

Answer 3

[Manish]

function C=kinetics(theta,t)

c0=[0.6178; 28.90; 0.4589];

[T,Cv]=ode45(@DifEq,t,c0);

function dC=DifEq(t,c)

dcdt=zeros(3,1);

mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2)^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5));

dcdt(1) = mu.*c(1);

dcdt(2) = -c(1).*((mu./(theta(6))+theta(7)));

dcdt(3) = c(1).*((theta(8).*mu)+theta(9));

dC=dcdt;

end

C=Cv;

end

t = [0

12

24

36

48

60

72];

c = [0.6178489703 28.90160183 0.7586206897

0.823798627 28.83295195 4.045977011

2.402745995 21.62471396 6.827586207

5.972540046 13.04347826 18.5862069

8.306636156 6.178489703 34.01149425

9.885583524 2.402745995 44.88505747

9.542334096 2.059496568 50.44827586];

theta0 = [1;1;1;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)', 'Location','N')

end