curve fitting using lsqcurvefit on kinetic data for parameter estimation

Hello,

I am fitting some experimental data (protein digestion kinetics) to the following model y = ymax+(ymax-y0)*exp(-k*t) using lsqcurvefit, were t is time (independent variable), y is concentration (dependent variable), and k, ymax and y0 are coefficient representing, respectively, the rate of the reaction, the maximum final concentration and the initial concentration.

The fitting seems to work well but the issue I cannot understand is why I get values of ymax higher than y0 when it should be the opposite.

Below you can find the code I'm using, do you have any idea/suggestion on where the issue could be?

Thanks a lot in advance for the hepl!!

xdata=[0; 30; 60; 90; 120; 180; 240];
ydata=[1.607; 2.346; 2.621; 2.967; 3.238; 3.479; 3.566];
coef = ["ymax","y0","k","R2","R2adj","RMSE"];

fun = @(x,xdata) x(1)+(x(1)-x(2))*exp(-x(3)*xdata)

fun = function_handle with value:
    @(x,xdata)x(1)+(x(1)-x(2))*exp(-x(3)*xdata)

x0 = [1,1,0.01];
lb = [0,0,0];
ub = [10,10,0.5];

[x,resnorm,residual,exitflag,output] = lsqcurvefit(fun,x0,xdata,ydata,lb,ub);

Local minimum found.

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

figure(1);
plot(xdata,ydata,'o',xdata,fun(x,xdata),'-');

SSresid = sum(residual.^2);
SStotal = (numel(ydata)-1) * var(ydata);
R = 1 - SSresid/SStotal;
Radj = 1 - (SSresid/SStotal) * ((numel(ydata)-1)/(numel(ydata)-1-1));

RMSE = rmse(fun(x,xdata),ydata);
r = [x R Radj RMSE];
coef = [coef;r];

figure(2);
scatter(xdata,residual);

coef

coef = 2x6 string array
    "ymax"      "y0"        "k"           "R2"        "R2adj"     "RMSE"    
    "3.6947"    "5.7563"    "0.012023"    "0.9945"    "0.9934"    "0.047899"

Answer 1

[Star Strider]

The model itself is a bit misleading.

A better option might be:

That produces an equivalent fit with parameters that make sense.

xdata=[0; 30; 60; 90; 120; 180; 240];
ydata=[1.607; 2.346; 2.621; 2.967; 3.238; 3.479; 3.566];
coef = ["ymax","y0","k","R2","R2adj","RMSE"];

% fun = @(x,xdata) x(1)+(x(1)-x(2))*exp(-x(3)*xdata)
fun = @(x,xdata) x(2)+(x(1)-x(2))*(1-exp(-x(3)*xdata));
x0 = [1,1,0.01];
lb = [0,0,0];
ub = [10,10,0.5];

[x,resnorm,residual,exitflag,output] = lsqcurvefit(fun,x0,xdata,ydata,lb,ub);

Local minimum found.

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

x

x = 1×3

3.6947 1.6331 0.0120

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure(1);
plot(xdata,ydata,'o',xdata,fun(x,xdata),'-');

SSresid = sum(residual.^2);
SStotal = (numel(ydata)-1) * var(ydata);
R = 1 - SSresid/SStotal;
Radj = 1 - (SSresid/SStotal) * ((numel(ydata)-1)/(numel(ydata)-1-1));

RMSE = rmse(fun(x,xdata),ydata);
r = [x R Radj RMSE];
coef = [coef;r];

figure(2);
scatter(xdata,residual);

coef

coef = 2x6 string array
    "ymax"      "y0"        "k"           "R2"        "R2adj"     "RMSE"    
    "3.6947"    "1.6331"    "0.012023"    "0.9945"    "0.9934"    "0.047899"

.

[federico drudi]

Thanks a lot guys, all the comments were very helpfull. I didn't realise that I've simply typed the wrong model, the correct one is the following y = ymax+(y0-ymax)*exp(-k*t), and gives identical results to the one you suggested.

fun = @(x,xdata) x(1)+(x(2)-x(1))*exp(-x(3)*xdata);

Thanks again for the support!!

Cheers.

[Star Strider]

As always, my pleasure!

Answer 2

[Umang Pandey]

Hi Federico,

Looking at your xdata and ydata, your ydata is increasing with increase in xdata, but the delta for each consecutive increase decreases, implying it is fitting "y = a - bexp(-kx)" where a,b,k are > 0. Since you are expecting "ymax > y0", your ydata should have been decreasing with increase in xdata, with the delta also decreasing for each consecutive decrease.

I have attached the following curves for your reference:

Case 1 : y = 4 - 3*exp(-2x) ; Assumptions : ymax = 4, y0 = 7, k = 2

Case 2 : y = 4 + 3*exp(-2x) ; Assumptions : ymax = 4, y0 = 1, k = 2

As you can see from the curve, your data/curve you obtained fits the first case.

Hope this helps!

Best,

Umang