lsqcurvefit Fitting kinetic model to estimate kinetic parameter.

I found this code in the community and used it for my own data(Thanks to the provider of the code), but no matter how I ran it, the data didn't fit well. How can I find the right initial value? Or what other function should I use to estimate my parameters?

These are the parameters I got from running.

		Theta( 1) =  0.00221
		Theta( 2) =  0.00427
		Theta( 3) =  0.00042
		Theta( 4) =  0.00816
		Theta( 5) =  0.00076
		Theta( 6) =  0.00000
		Theta( 7) =  0.00000
		Theta( 8) =  0.00000
		Theta( 9) =  0.00002
		Theta(10) =  2.94054
		Theta(11) =  0.02476

The code is as follows:

clc;
clear all;

t=[ 0
    10
    20
    30
    40
    50
    60];
c=[ 100 100 0 0 0 0 0
    94.77333333 93.73333333 0.400145927	0.27701074 0.107176434 0.158130584 0.000505722
    84.48 88.87333333 0.604820791 0.42039699 0.148464431 0.359954181 0.001077632
    73.91333333 83.23333333 0.589530533 0.461026019 0.17071451 0.508774341 0.001205872
    68.80666667 76.43333333 0.529979 0.48559649 0.180350987 0.635601375 0.001720902
    63.56666667 73.92 0.502282172 0.52315421 0.198918834 0.700114548 0.002064255
    59.76666667 70.59333333 0.478273783 0.574401193 0.207380132 0.845269187 0.002895749
  ]; 

theta0=rand(11,1)*1E-3;
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));




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


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

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
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, 'c(1)', 'c(2)', 'c(3)', 'c(4)', 'c(5)','c(6)','c(7)', 'Location','best');




function C=kinetics(theta,t)
c0=[100; 100; 0; 0; 0; 0; 0];
[T,Cv]=ode45(@DifEq,t,c0);
    function dC=DifEq(t,c)
        dcdt=zeros(7,1);
        dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
        dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
        dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
        dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
        dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
        dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
        dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
    end
C=Cv;
end

Answer 1

[Star Strider]

That does not appear to me to be the same data.

I get decent results except for and , leading me to believe that there is something wrong with their differential equations. Please check them.

% clc;
% clear all;

t=[ 0
    10
    20
    30
    40
    50
    60];
c=[ 100 100 0 0 0 0 0
    94.77333333 93.73333333 0.400145927	0.27701074 0.107176434 0.158130584 0.000505722
    84.48 88.87333333 0.604820791 0.42039699 0.148464431 0.359954181 0.001077632
    73.91333333 83.23333333 0.589530533 0.461026019 0.17071451 0.508774341 0.001205872
    68.80666667 76.43333333 0.529979 0.48559649 0.180350987 0.635601375 0.001720902
    63.56666667 73.92 0.502282172 0.52315421 0.198918834 0.700114548 0.002064255
    59.76666667 70.59333333 0.478273783 0.574401193 0.207380132 0.845269187 0.002895749
  ]; 

theta0 = [rand(11,1)*1E-3; c(1,:).'];

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

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function 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.00256
		Theta( 2) =  0.00146
		Theta( 3) =  0.00056
		Theta( 4) =  0.00230
		Theta( 5) =  0.00094
		Theta( 6) =  0.00000
		Theta( 7) =  0.00000
		Theta( 8) =  0.00000
		Theta( 9) =  0.00005
		Theta(10) =  2.86843
		Theta(11) =  0.04557
		Theta(12) = 101.54545
		Theta(13) = 100.16133
		Theta(14) =  0.00000
		Theta(15) =  0.00000
		Theta(16) =  0.00036
		Theta(17) =  0.00000
		Theta(18) =  0.00000

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

figure(1)
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, 'c(1)', 'c(2)', 'c(3)', 'c(4)', 'c(5)','c(6)','c(7)', 'Location','best');

figure
for k = 1:size(c,2)
    subplot(4,2,k)
    plot(t,c(:,k),'bp', 'MarkerFaceColor','b')
    hold on
    plot(tv, Cfit(:,k))
    hold off
    grid
    title("C_"+k)
end

function C=kinetics(theta,t)
% c0=[100; 100; 0; 0; 0; 0; 0];
c0 = theta(end-6:end);
[T,Cv]=ode45(@DifEq,t,c0);
    function dcdt=DifEq(t,c)
        dcdt=zeros(7,1);
        dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
        dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
        dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
        dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
        dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
        dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
        dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
    end
C=Cv;
end

I changed the code slightly to add the initial conditions as parameters to be estimated, since in my experience, that generally produces better results, even if the initial conditions themselves do not change much from their original values. (I added the subplot array because otherwise the last 5 compartments are squashed together at the lower end of the first plot.)

.

[Kaihua Liu]

Thank you for your answer! I'll check my equation again.

[Star Strider]

As always, my pleasure!

I will help as much as I can if you have other problems with it (since I wrote the original version of the code you are using, and feel responsible for being certain that it works correctly). I do not know what you are calculating, and since I do not have access to those differential equations, I cannot check it myself.

[Kaihua Liu]

I'm sorry to bother you again, but as a newcomer to MATLAB, how should I change the code to output the model's goodness-of-fit index R2 and the 95% confidence interval (95%CI) of parameter estimates? Thank you again!

[Star Strider]

I was in the process of responding to this when Micro$oft crashed my computer again with more of this azure and autopilot s**t. I am never buying a Windows machine again, and I am going to seriously consider archiving my files, wiping Win 11 from this machine, and installing Linux (that I have to learn about first, however that errort will be more than worth the effort).

Anyway, my intense hatred for Micro$oft aside, the only option fot the confidence intervals is nlpredci. (Use nlparci to get confidence intervals on the parameters.) I can implement that here by tweaking calls to it. I am not certain how to calculate the R² or other goodness-of-fit statistics for a regression such as this. (Occasonally an svd error message appears in the nlpredci call. If that occurs, just run the code again.)

Try this —

% clc;
% clear all;

t=[ 0
    10
    20
    30
    40
    50
    60];
c=[ 100 100 0 0 0 0 0
    94.77333333 93.73333333 0.400145927	0.27701074 0.107176434 0.158130584 0.000505722
    84.48 88.87333333 0.604820791 0.42039699 0.148464431 0.359954181 0.001077632
    73.91333333 83.23333333 0.589530533 0.461026019 0.17071451 0.508774341 0.001205872
    68.80666667 76.43333333 0.529979 0.48559649 0.180350987 0.635601375 0.001720902
    63.56666667 73.92 0.502282172 0.52315421 0.198918834 0.700114548 0.002064255
    59.76666667 70.59333333 0.478273783 0.574401193 0.207380132 0.845269187 0.002895749
  ]; 

theta0 = [rand(11,1)*1E-3; c(1,:).'];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,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.

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.00284
		Theta( 2) =  0.00021
		Theta( 3) =  0.00051
		Theta( 4) =  0.00006
		Theta( 5) =  0.00260
		Theta( 6) =  0.00013
		Theta( 7) =  0.00012
		Theta( 8) =  0.00008
		Theta( 9) =  0.00007
		Theta(10) =  0.17364
		Theta(11) =  0.02283
		Theta(12) = 101.54879
		Theta(13) = 100.39359
		Theta(14) =  0.00000
		Theta(15) =  0.00000
		Theta(16) =  0.00000
		Theta(17) =  0.06798
		Theta(18) =  0.00000

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

for k = 1:size(c,2)
    [Ypred,delta] = nlpredci(@(theta,t)kinetics1(theta,t,k),tv,theta,Rsd,'Jacobian',Jmat);
    ypred(:,k) = Ypred;
    CI(:,k) = delta;
end

figure(1)
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, 'c(1)', 'c(2)', 'c(3)', 'c(4)', 'c(5)','c(6)','c(7)', 'Location','best');

figure
for k = 1:size(c,2)
    subplot(4,2,k)
    hp1 = plot(t,c(:,k),'bp', 'MarkerFaceColor','b', 'DisplayName','Data');
    hold on
    hp2 = plot(tv, ypred(:,k), '-r', 'DisplayName','Regression');
    hp3 = plot(tv, ypred(:,k)+CI(:,k)*[-1 1], ':r', 'DisplayName','±95% Conofidence Intervals');
    hold off
    grid
    title("C_"+k)
end
hl = legend([hp1 hp2 hp3(1)]);
subplot(4,2,8)
Ax = gca;
Ax.Visible = 0;
hl.Position = Ax.Position;

function C=kinetics(theta,t)
% c0=[100; 100; 0; 0; 0; 0; 0];
c0 = theta(end-6:end);
[T,Cv]=ode45(@DifEq,t,c0);
    function dcdt=DifEq(t,c)
        dcdt=zeros(7,1);
        dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
        dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
        dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
        dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
        dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
        dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
        dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
    end
C=Cv;
end

function C=kinetics1(theta,t,k)
% c0=[100; 100; 0; 0; 0; 0; 0];
c0 = theta(end-6:end);
[T,Cv]=ode45(@DifEq,t,c0);
    function dcdt=DifEq(t,c)
        dcdt=zeros(7,1);
        dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
        dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
        dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
        dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
        dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
        dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
        dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
    end
C=Cv(:,k);
end

NOTE — I added the ‘kinetics1’ function since that is necessary to run the nlpredci function with individual vectors, since it does not accept matrix arguments. That also explains the for loop that calls it.

EDIT — Corrected typographical errors.

.

[Kaihua Liu]

I agree with your point of view, WIN11 is indeed not the most perfect version at present. First of all thanks for your answer, I ran the code and it seems to be the 95%CI of the whole curve, how can I get the 95%CI of the estimated parameter 'theta'?

[Star Strider]

For the parameter estimates, use the nlparci function.

That call would be:

ci = nlparci(theta,Rsd,'jacobian',Jmat)

The two-sided confidence intervals will be in the order of the parameter vector ‘theta’.

The parameters are significant if both confidence intervals have the same sign. If they are opposite signs, that means the 95% confidence intervals include zero, and that parameter is not necessary for the model.

.

[Kaihua Liu]

I have set the upper limit ub and lower limit lb of the parameter estimation, why there are still negative numbers? And the interval range is far beyond the calculated value. Something about parameter estimation is really bothering me.

[Star Strider]

The way I have them set (‘lb’ is zero, ‘ub’ not specified, so essentially infinite), all the parameter values are positive. Check to be certain that all the arguments are in the correct locations in the argument list. Otherwise, copy my lsqcurvefit call and it should work.

EDIT — (5 Aug 2023 at 16:59)

Thinking about this further, I believe that you are referring to the 95% confidence limits of the parameters, not the parameters themselves, none of which should be negative in this particular model.

As I mentioned in my previous Comment:

• The parameters are significant if both confidence intervals have the same sign. If they are opposite signs, that means the 95% confidence intervals include zero, and that parameter is not necessary for the model.

So wide confidence intervals are occasionally to be expected, and the parameter value should be midway between them.

If the confidence intervals are of the same sign (either positive or negative), the parameter is significantly different from zero with 95% probability, and so contributes to the regression. If the confidence intervals are of opposite signs, that means that they include zero (that parameter is not statistically significantly different from zero with 95% probability), so that specific parameter is actually not necessary in the regression.

.

[Kaihua Liu]

I'm sorry it took so long to reply to your message. I've been busy with the experiment recently. Recently, I have been perplexed by this problem. There must be some errors in parameter estimation. How can I get the standard errors of these parameters theta? The previous program seems to have only got the confidence interval of the concentration c.

[Star Strider]

Getting the standard errors for the parameters is likely going to require writing your own function (to either calculate them directly or back-calculate from the confidence intervals), and probably is not necessary anyway. The nlparci function (discussed earlier here and here) will give the 95% confidence intervals for the parameters, and that is actually the most significant information.

I agree that there is a problem with the code, particularly with and . Look at them specifically.

[Kaihua Liu]

Thank you very much for your answer! I'll check the code carefully.

[Star Strider]

As always, my pleasure!

[dinesh kumar s]

@Star Strider Hi, I was working on similar problem like this and ofcourse I was using the same code written by you for the kinetic modeling. I hope the parameters are good. I just want to thank you @Star Strider for the code. It works well and good. Thanks.

[Star Strider]

@dinesh kumar s — My pleasure!

A Vote would be appreciated!