Find kinetic constants from differential equations

Hey! So a escription on my problem:

I have a compartimental model that contains 4 different elements that have a kinetic behaviour between thm. The differential equations of this model can be described as:

dydt(1) = k(1)*y(4)+k(2)*y(2)+k(3)*y(1)-(k(4)*y(1)+k(5)*y(1)+k(6)*y(1)+1/25.028)*y(1);

dydt(2) = k(2)*y(1)-(k(5)*y(2)+1/25.028)*y(2);

dydt(3) = k(3)*y(1)-(k(4)*y(3)+1/25.028)*y(3);

dydt(4) = k(1)*y(1)-(k(4)*y(4)+1/25.028)*y(4);

I have some experimental data that I want to use to predict th different k, but I have been struggeling with solving them with the ode45 function or similars. Could someone help me?

Answer 1

[Star Strider]

Some examples: Monod kinetics and curve fitting , Parameter Estimation for a System of Differential Equations , and there are several others.

[Marina Batlló RIus]

Hello! Thanks for you response. I was looking at the Monod kinetics and curve fitting example, and I think I am missing something. When I see the function that you provided as response, there are two variables B0 and x0 (one in the funtion itself and one used when calling) that I do not fully understand. Could you help me?

[Star Strider]

My pleasure!

The ‘B0’ vector are the initial parameter estimates for lsqcurvefit.

The ‘x0’ vector are the initial conditions to differential equations in the ode45 call

The code in the second link are likely more appropriate for your problem. I have also updated it to:

t=[0.1
    0.2
    0.4
    0.6
    0.8
    1
    1.5
    2
    3
    4
    5
    6];
c=[0.902	0.06997	0.02463	0.00218
    0.8072	0.1353	0.0482	0.008192
    0.6757	0.2123	0.0864	0.0289
    0.5569	0.2789	0.1063	0.06233
    0.4297	0.3292	0.1476	0.09756
    0.3774	0.3457	0.1485	0.1255
    0.2149	0.3486	0.1821	0.2526
    0.141	0.3254	0.194	0.3401
    0.04921	0.2445	0.1742	0.5277
    0.0178	0.1728	0.1732	0.6323
    0.006431	0.1091	0.1137	0.7702
    0.002595	0.08301	0.08224	0.835]; 
theta0 = rand(10,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

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

		Theta( 1) =  0.76485
		Theta( 2) =  0.23489
		Theta( 3) =  0.20888
		Theta( 4) =  0.49193
		Theta( 5) =  0.62211
		Theta( 6) =  0.00000
		Theta( 7) =  0.90288
		Theta( 8) =  0.07145
		Theta( 9) =  0.02841
		Theta(10) =  0.00000

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')

function C=kinetics(theta,t)
c0=theta(end-3:end);
[T,Cv]=ode45(@DifEq,t,c0);

    function dC=DifEq(t,c)
        dcdt=zeros(4,1);
        dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
        dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
        dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
        dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
        dC=dcdt;
    end
C=Cv;
end

NOTE — In this version, the initial conditions for the differential equationos are parameters to be estimated. The plot colours also match (data and estimates).

.

[Marina Batlló RIus]

Thank you so much!

[Star Strider]

As always, my pleasure!