Fitting and parameter estimation of kinetic model ODE system

Hello, I am trying to fit experimental data with simultaneous ODE like shown below:

clc
global t c 
t = [0
24 
48 
96 
144 
192
240 
288]; 
c = [0.08 22.55530474 0.003
0.508650519 22.84424379 0.1349480969 
1.089965398 23.13318284 0.9653979239 
2.584775087 20.78555305 2.169550173 
4.951557093 16.1986456 6.197231834 
5.906574394 13.30925508 10.01730104 
6.321799308 12.62302483 13.38062284 
6.737024221 13.23702032 14.79238754]; 
theta0 = [0.05;1.66;110;76;3.66;53.26;0.5;0.0128;2.35;0.02;0.0017;63.95;19.1;72.5]; 
options = optimset('display','iter');
options.MaxFunEvals = 10e8;
options.TolX = 10e-8;
options.TolFun = 10e-7;
options.MaxIter = 10e6;
[theta] = fminsearch(@minfunction, theta0, options);
tv = linspace(min(t), max(t)); 
Cfit = kinetics(theta, tv, c(1,:)); 

figure
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit(:,[1 2 3]));
hold off
grid
xlabel('Waktu(jam)')
ylabel('Biomassa(g/L),Substrat(g/L),Astaksantin(mg/L)')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'Location','N')

for k1 = 1:length(theta) 

    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1)) 

end 

function C=kinetics(theta, t, c0)
[T,Cv]=ode45(@DifEq,t, c0); 
    function dC=DifEq(t,c) 
        dcdt=zeros(3,1); 
        I=80;
        mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2).^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5)).*(I/(theta(6)+I)); 
        dcdt(1) = mu.*c(1); 
        dcdt(2) = -c(1).*((mu./(theta(7))+theta(8))); 
        dcdt(3) = c(1).*((theta(9).*mu)+theta(10)).*(c(2)/(theta(11)+c(2)+(c(2)^2/theta(12)))).*(c(3)/(theta(13)+c(3))).*(I/(theta(14)+I));
        dC=dcdt; 
    end 
C=Cv;
end 

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

However, in that code, I have stated that "I=80", where actually I want to make the value of "I" is 35 (from t 0 to 96 hours) and 80 (t more than 96 hours). "I" is the light intensity.

Could you help me what should I do to change the code according to my request? Thank you for the help!

Answer 1

[Star Strider]

Optimising 14 parameters is asking a lot of fminsearch! I switched to lsqcurvefit here because fminsearch is not likely to be useful beyond about 7 parameters. Also, don’t clear at the beginning, and definitely do not use global variables! Pass the other arguments as extra parameters, if necessary.

This works, although it will be necessary to let it run longer than I have here so that lsqcurvefit can converge on a better sollution. I added the initial conditions as the last 3 elements of ‘theta’ so they are now parameters-to-be-estimated. In my experience, this works better than fixing them at specific values in the code.

The equation system can reasonably be called ‘stiff’ so I switched from ode45 toode15s with a significant speed imprivement.

The numeric ODE solvers have problems with abrupt, non-differentialbe ‘step’ transitions, and there are two ways to deal with that. One is to interrupt the integration with the first value of ‘I’ then save the last values and time and re-start it with those as the intial conditions and re-start the integration. The other is to use the tanh function to privide a differentiable transition at the appropriate time. That’s what I did here.

I also updated the plot calls so that the colours match and the plot makes more sense.

% clc
% global t c                                      % Don't Use 'global' Variables!  
t = [0
    24 
    48 
    96 
    144 
    192
    240 
    288]; 
c = [0.08 22.55530474 0.003
    0.508650519 22.84424379 0.1349480969 
    1.089965398 23.13318284 0.9653979239 
    2.584775087 20.78555305 2.169550173 
    4.951557093 16.1986456 6.197231834 
    5.906574394 13.30925508 10.01730104 
    6.321799308 12.62302483 13.38062284 
    6.737024221 13.23702032 14.79238754]; 
theta0 = [0.05;1.66;110;76;3.66;53.26;0.5;0.0128;2.35;0.02;0.0017;63.95;19.1;72.5]; 
theta0 = [theta0.*rand(size(theta0)); c(1,:).'];
% theta0 = rand(17,1);
% % options = optimset('display','iter');
% % options.MaxFunEvals = 10e8;
% % options.TolX = 10e-8;
% % options.TolFun = 10e-7;
% % options.MaxIter = 10e6;
[theta] = lsqcurvefit(@objfunction, theta0, t, c, zeros(size(theta0)))

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.700000e+03.

theta = 17×1

0.0329 0.6924 14.6130 42.2253 3.4208 0.6521 0.2545 0.0003 0.0337 0.3789

tv = linspace(min(t), 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')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

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

		Theta( 1) =   0.032858
		Theta( 2) =   0.692390
		Theta( 3) =  14.612964
		Theta( 4) =  42.225345
		Theta( 5) =   3.420782
		Theta( 6) =   0.652061
		Theta( 7) =   0.254480
		Theta( 8) =   0.000271
		Theta( 9) =   0.033661
		Theta(10) =   0.378941
		Theta(11) =  49.187783
		Theta(12) =  11.247759
		Theta(13) =   0.845725
		Theta(14) =  13.764339
		Theta(15) =   0.216072
		Theta(16) =  22.779756
		Theta(17) =   0.189724

function C=kinetics(theta, t)
c0 = theta(end-2:end);                              % Last Three 'theta' Elements Are Initial Conditions (To Be Estimated As Parameters)
I = @(t) 35+22.5*(1+tanh(t-96));                         % Differentiable Function For 'I' That 'Steps' At The Appropriate Time!
[T,Cv]=ode15s(@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)).*(I(t)/(theta(6)+I(t))); 
        dcdt(1) = mu.*c(1); 
        dcdt(2) = -c(1).*((mu./(theta(7))+theta(8))); 
        dcdt(3) = c(1).*((theta(9).*mu)+theta(10)).*(c(2)/(theta(11)+c(2)+(c(2)^2/theta(12)))).*(c(3)/(theta(13)+c(3))).*(I(t)/(theta(14)+I(t)));
        dC=dcdt; 
    end 
C=Cv;
end 

% function error = minfunction(theta,t,c)
% % global t c                                      % Don't Use 'global' Variables!  
% 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);

function c_estim = objfunction(theta,t)
% global t c                                      % Don't Use 'global' Variables!  
% c0 = c(1,:);
c_estim = kinetics(theta,t);
% 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

Let it run a bit longer to get a better fit by increasing the funciton evaluation limit and other options as necessary. Another option is to sue the genetic algorithm to estimate the parameters, however it is likely that lsqcurvefit can estimate them if given a longer time to do so.

.

[Nadya Shalsabila salman]

Hi! Thank you very much @Star Strider! However, could you explain why the I function is written like this "I = @(t) 35+22.5*(1+tanh(t-96))"? I do not understand about this tanh thing. And I also wonder why I get different result everytime I run the code?

[Star Strider]

My pleasure!

As I mentioned in my original Answer:

‘The numeric ODE solvers have problems with abrupt, non-differentialbe ‘step’ transitions, and there are two ways to deal with that. One is to interrupt the integration with the first value of ‘I’ then save the last values and time and re-start it with those as the intial conditions and re-start the integration. The other is to use the tanh function to provide a differentiable transition at the appropriate time. That’s what I did here.’

So I created the ‘I’ function to provide this solution —

t = linspace(0, 288, 500);
I = @(t) 35+22.5*(1+tanh(t-96));

figure
plot(t,I(t))
grid
xlabel('t')
ylabel('I(t)')
xlim([0 288])
ylim([30 90])

So ‘I’ is 35 until the function gets to 93 then shifts at 96 to become 80 after that. The time at which it shifts can be adjusted by changing the 96 in the tanh argument to something else. Everything else about the function remains the same.

For example, changing the function slightly to add an additional argument —

t = linspace(0, 288, 500);
I = @(t,tau) 35+22.5*(1+tanh(t-tau));

figure
plot(t,I(t,90),'-b')
hold on
plot(t, I(t,96),'-g')
plot(t, I(t,102),'-r')
hold off
grid
xlabel('t')
ylabel('I(t)')
legend('\tau =   90','\tau =   96','\tau = 102', 'Location','best')
xlim([0 288])
ylim([30 90])

Choose whatever value of τ works best (that is, wherever the transition should be, either before (), at (), or after () 96 hours). The differentiable function provides a smooth transition, so it should work regardless of the value of τ, providing that the value is within reason and that it meets the criteria set out in the assumptions of the differential equation system.

.

[Nadya Shalsabila salman]

You've been very helpfu! Thank you @Star Strider

[Star Strider]

As always, my pleasure!

.

[Nadya Shalsabila salman]

sorry to bother you again @Star Strider, what should I change if I want the "I" value becomes more than 80? and how to decide the number that is going to be changed?

[Star Strider]

No worries.

Use this version of ‘I’ —

I = @(t,tau,lo,hi) lo+(hi-lo)/2*(1+tanh(t-tau));    % Arguments: 't'=Time Vector, 'tau'=Delay, 'lo'=Low Value (Before 'tau'), 'hi'=High Value (After 'tau')

Copy and paste the entire line (including the comment) so that the comment-documentation stays with it. It now permits varying all three parameters, the delay (switching time), low value, and high value. I tested it and it appears to work for every real value of every parameter, regardless of sign.

.

[Nadya Shalsabila salman]

Thank you very much @Star Strider. May I ask you again? So I have tried this code to plot in two y axis.. but the program is error. Could you help me to solve my problem?

The code from this below (1 y axis)

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')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

to (2 y axis) I have tried

figure
yyaxis left
hd = plot(t,c(:,[1 2]),'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
ylabel('Biomassa dan Substrat(g/L)')
yyaxis right
plot(t,c(:,3),'p')
hold on
hlpp=plot(tv,Cfit(:,3));
hold off
ylabel('Karotenoid (\mug/L)')
xlabel('waktu (jam)')
grid
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')
xlabel('Waktu (jam)')

[Star Strider]

I got it to plot correctly, however I had to abandon the custom colours and the legend, partially because that tends to confuse the interpretation of the yyaxis plot. (The key to making it work with yyaxis is completely copying the plotting loops under both yyaxis calls.)

t = [0
    24 
    48 
    96 
    144 
    192
    240 
    288]; 
c = [0.08 22.55530474 0.003
    0.508650519 22.84424379 0.1349480969 
    1.089965398 23.13318284 0.9653979239 
    2.584775087 20.78555305 2.169550173 
    4.951557093 16.1986456 6.197231834 
    5.906574394 13.30925508 10.01730104 
    6.321799308 12.62302483 13.38062284 
    6.737024221 13.23702032 14.79238754]; 
theta0 = [0.05;1.66;110;76;3.66;53.26;0.5;0.0128;2.35;0.02;0.0017;63.95;19.1;72.5]; 
theta0 = [theta0.*rand(size(theta0)); c(1,:).'];
% theta0 = rand(17,1);
% % options = optimset('display','iter');
% % options.MaxFunEvals = 10e8;
% % options.TolX = 10e-8;
% % options.TolFun = 10e-7;
% % options.MaxIter = 10e6;
[theta] = lsqcurvefit(@objfunction, theta0, t, c, zeros(size(theta0)));

Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.

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

figure
yyaxis left
hdl = plot(t, c(:,[1 2]), 'p');
% for k1 = 1:size(c,2)
%     CV(k1,:) = hd(k1).Color;
%     hd(k1).MarkerFaceColor = CV(k1,:);
% end
hold on
hlpl = plot(tv, Cfit(:,[1 2]));
% for k1 = 1:size(c,2)
%     hlp(k1).Color = CV(k1,:);
% end
hold off
ylabel('Concentration C_1, C_2')

yyaxis right
hdr = plot(t, c(:,3), 'p');
% for k1 = 1:size(c,2)
%     CV(k1,:) = hd(k1).Color;
%     hd(k1).MarkerFaceColor = CV(k1,:);
% end
hold on
hlpr = plot(tv, Cfit(:,3));
% for k1 = 1:size(c,2)
%     hlp(k1).Color = CV(k1,:);
% end
hold off
grid
xlabel('Time')
ylabel('Concentration C_3')

% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')


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

		Theta( 1) =   0.035400
		Theta( 2) =   0.105053
		Theta( 3) =  82.419693
		Theta( 4) =  18.813437
		Theta( 5) =   5.176098
		Theta( 6) =  48.379487
		Theta( 7) =   0.446478
		Theta( 8) =   0.006028
		Theta( 9) =   0.758765
		Theta(10) =   0.147805
		Theta(11) =   0.001606
		Theta(12) =  85.192877
		Theta(13) =   7.038255
		Theta(14) =  32.922152
		Theta(15) =   0.205821
		Theta(16) =  22.386365
		Theta(17) =   0.102030

function C=kinetics(theta, t)
c0 = theta(end-2:end);                              % Last Three 'theta' Elements Are Initial Conditions (To Be Estimated As Parameters)
I = @(t) 35+22.5*(1+tanh(t-96));                         % Differentiable Function For 'I' That 'Steps' At The Appropriate Time!
[T,Cv]=ode15s(@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)).*(I(t)/(theta(6)+I(t))); 
        dcdt(1) = mu.*c(1); 
        dcdt(2) = -c(1).*((mu./(theta(7))+theta(8))); 
        dcdt(3) = c(1).*((theta(9).*mu)+theta(10)).*(c(2)/(theta(11)+c(2)+(c(2)^2/theta(12)))).*(c(3)/(theta(13)+c(3))).*(I(t)/(theta(14)+I(t)));
        dC=dcdt; 
    end 
C=Cv;
end 

% function error = minfunction(theta,t,c)
% % global t c                                      % Don't Use 'global' Variables!  
% 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);

function c_estim = objfunction(theta,t)
% global t c                                      % Don't Use 'global' Variables!  
% c0 = c(1,:);
c_estim = kinetics(theta,t);
% 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

.

[Nadya Shalsabila salman]

Thank you very much @Star Strider!

[Star Strider]

As always, my pleasure!

.