Reaction kinetics parameter determination

Hi, I am new to MATLAB and I have a series of data points to fit and then calculate the rate constants.

here are the equations:

dcdt(1)=-k(1)*c(1)-k(2)*c(1)-k(3)*c(1)

dcdt(2)=k(2)*c(1)-k(6)*c(2)-k(9)*c(2)

dcdt(3)=k(3)*c(1)-k(5)*c(4)-k(7)*c(6)

dcdt(4)=k(1)*c(1)-k(4)*c(4)-k(5)*c(4)

dcdt(5)=k(6)*c(2)+k(4)*c(4)-k(10)*c(5)

dcdt(6)=k(7)*c(3)-k(8)*c(6)

dcdt(7)=k(8)*c(6)+k(9)*c(2)+k(10)*c(5)-k(11)*c(8)-k(12)*c(9)

dcdt(8)=k(11)*c(7)

dcdt(9)=k(12)*c(8)

and the data:

Time=[0 2 4 10 15 25 40]

c= [100 0 0 0 0 0 0

0 4 37 15 5 4 2

0 10 20 8 4 0.2 0.1

0 5 10 16 4 2 1

0 4 5 11 5 4 3

0 5 6 5 1 0.1 0.04

0 6 12 32 64 44 31

0 0.0071 0.06 0.14 0.3000 3.8500 7.4440

0 0 0 0.0340 0.340 0.5000 0.630]

I would really appreciate the help! Thanks in advance

Answer 1

[Star Strider]

See Parameter Estimation for a System of Differential Equations for a working example of how to do this correctly.

Unfortunately, since you have 7 data rows (governed by the number of elements in the ‘Time’ vector, since that is the independent variable), and are estimating 12 parameters (my version estimates 21 parameters, however the last 9 are the initial conditions of the differential equations, so they don’t count), this will fail.

This runs without error, however because of the mismatch between the data and the parameters being estimated, the results will not be reliable.

Adapting that approach to your code —

function c = kinetics(k,t)
        c0 = k(13:21);
        [t,Cv] = ode15s(@DifEq,t,c0);
        
    function dcdt=DifEq(t,c)
            dcdt = zeros(9,1);
            dcdt(1)=-k(1)*c(1)-k(2)*c(1)-k(3)*c(1);
            dcdt(2)=k(2)*c(1)-k(6)*c(2)-k(9)*c(2);
            dcdt(3)=k(3)*c(1)-k(5)*c(4)-k(7)*c(6);
            dcdt(4)=k(1)*c(1)-k(4)*c(4)-k(5)*c(4);
            dcdt(5)=k(6)*c(2)+k(4)*c(4)-k(10)*c(5);
            dcdt(6)=k(7)*c(3)-k(8)*c(6);
            dcdt(7)=k(8)*c(6)+k(9)*c(2)+k(10)*c(5)-k(11)*c(8)-k(12)*c(9);
            dcdt(8)=k(11)*c(7);
            dcdt(9)=k(12)*c(8);
    end
        c = Cv;
end

Time=[0 2 4 10 15 25 40];
c= [100  0  0  0  0  0  0
      0  4  37 15  5  4  2   
      0  10  20 8  4   0.2  0.1
      0    5   10   16   4   2  1
      0    4    5   11  5  4  3
      0    5    6    5    1  0.1  0.04    
      0    6   12   32   64   44  31
      0    0.0071    0.06   0.14    0.3000    3.8500   7.4440
      0    0    0    0.0340    0.340    0.5000    0.630];
  
  c = c.';
  
  k0 = rand(21,1)*1E-3;

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,Time,c,zeros(size(k0)));

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

tv = linspace(min(Time), max(Time));
Cfit = kinetics(k, tv);

figure(1)
hd = plot(Time, 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')
lgd = legend(hlp,compose('C_{%d}(t)',1:9), 'Location','N');
lgd.NumColumns = 3;

This will work best if you put all of it inside another function, and then call the external function. See the link I referred to for details.

However more data (at least 12 rows of the transposed ‘c’ matrix) are absolutely necessary if you want to estimate all 12 parameters accurately.

EDIT — (3 Mar 2021 at 4:02)

Added plot —

.

[iPanda]

Thank you for sharing this. I am having trouble in running the code in my software, its taking forever. But as I understood, I need more points to estimate my 12 parameters?

[Star Strider]

As always, my pleasure!

Note that I deleted my first Answer and replaced with a second corrected Answer. Change the solver to ode15s and the integration proceeds quickly. (The parameters vary sufficiently in magnitude that the system becomes ‘stiff’.)

‘But as I understood, I need more points to estimate my 12 parameters?’

Yes, definitely. This code will estimate them, however the estimates are unreliable. If you have the Statistics and Machine Learning Toolbox, and you use nlparci to calculate confidence intervals on the fitted parameters, many of the confidence intervals will include zero. (I have not done this, however I am confident of the results.) More data are always better. I suggest that if you run the experiment again that you simply sample it more frequently, for example every 2 time units. That will generate 20 data points for all the compartments, more than enough to accurately estimate the 12 parameters. (I do not suggest attempting to interpolate the existing data, since they are too noisy to provide for reliable interpolation.)

[iPanda]

Will try that, Thank you!

[Star Strider]

As always, my pleasure!