# Fitting Monod Equation with ODE45 to data using lsqcurvefit function

97 views (last 30 days)

Show older comments

StarSign1997
on 22 Apr 2019

Commented: Antara Mazumder
on 6 Jun 2020

Hello!

I am fitting Monod equation to a data containing substrate (s), biomass (x), and ethanol (p) concentration against time. The objective is to get the parameters: 1) umax, 2) ks, 3) Yxs, and 4)Yps that will best represent the data. The differential equations are:

Here is my initial code using assumed values of the four parameters:

umax = 0.5;

ks = 6.5;

Yxs = 0.2;

Yps = 1.2;

%a(1) = x

%a(2) = s

%a(3) = p

f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];

xt0 = [0.0904,9.0115,0.0151];

[tspan,a] = ode45(f,[0 25],xt0);

figure

plot(tspan,a(:,1),tspan,a(:,2),tspan,a(:,3))

Here is the code for trying to fit it into the actual data (script file):

function pos = paramfun1(x,tspan)

umax = x(1);

ks = x(2);

Yxs = x(3);

Yps = x(4);

xt0 = x(5:7);

f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];

[~,pos] = ode45(f,tspan,xt0);

Here is my call function (in the command window):

xt0 = zeros(1,7);

xt0(1) = umax;

xt0(2) = ks;

xt0(3) = Yxs;

xt0(4) = Yps;

data =[0 3 5 8 9.5 11.5 14 16 18 20 25 27, 0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093; 0 3 5 8 9.5 11.5 14 16 18 20 25 27, 9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0; 0 3 5 8 9.5 11.5 14 16 18 20 25 27, 0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869];

%time =[0 3 5 8 9.5 11.5 14 16 18 20 25 27];

[pbest,exitflag,output] = lsqcurvefit(@paramfun,xt0,tspan,data);

fprintf('New parameters: %f, %f, %f, %f',pbest(1:4));

The error is function value not equal to YDATA. Btw, this code was based from an example in MATLAB. (https://www.mathworks.com/help/optim/ug/fit-differential-equation-ode.html)

My data is:

time = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]

x = 0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093

s = 9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0

p = 0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869

Please help! I do not know how to input my data into the lsqcurvefit function.

Thanks in advance!

### Accepted Answer

David Wilson
on 22 Apr 2019

Edited: David Wilson
on 22 Apr 2019

Here's my solution.

First set up the measured data, plot it, and look to see that it is all OK.

t = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]'

x = [0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093]';

s = [9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0]';

p = [0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869]';

plot(t,[x,s,p],'o-')

No problems there. Note that the measured data is in columns.

Now, we need a function to export the predicted measurements as a function of the unknown parameters and initial conditions. My function is:

function ypred = paramfun1(p,t)

umax = p(1); ks = p(2); Yxs = p(3); Yps = p(4); % use disperse here ...

xt0 = p(5:7); % initial conditions

f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); ...

-(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); ...

(1/Yps)*umax*a(1)*a(2)/(ks + a(2))];

[~,ypred] = ode45(f,t,xt0);

end

May I suggest that you think about naming your functions sensible names. Also it is more elegant to use disperse (from the user's group) here.

Now we are ready to test it, and make sure our output is the same shape as the actual data.

xt0 = [0.1;9;0.1]; % read off from data plot [x0 ,S0, p0]';

p0 = [umax; ks; Yxs; Yps; xt0]

Ypred = paramfun1(p0,t); % measurement predictions

Now we are ready for the fit. I'm using your initial conditions, and we can get reasonable start values from the data.

pfit = lsqcurvefit(@paramfun1,p0,t,[x,s,p])

Now plot the fitted curve and the actual data. They should line up.

umax = pfit(1); ks = pfit(2);Yxs = pfit(3); Yps = pfit(4);

f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];

xt0 = pfit(5:7);

[tspan,a] = ode45(f,[0 25],xt0);

plot(tspan,a(:,1),tspan,a(:,2),tspan,a(:,3), '-', ...

t,x,'o',t,s,'+',t,p,'s')

Plot shows the regressed fit is not too bad:

##### 7 Comments

Antara Mazumder
on 6 Jun 2020

### More Answers (2)

Alex Sha
on 21 Oct 2019

The follow results obtained from 1stOpt may be used for comparison:

Code:

Variable t,x,s,p;

ODEFunction x'=umax*s*x/(ks + s);

s'=-(1/Yxs)*umax*s*x/(ks + s);

p'=(1/Yps)*umax*s*x/(ks + s);

Data;

t = [0 2 4 6 8 10 12 14 16];

x = [0.06 0.11 0.46 0.78 1.42 2.36 2.49 2.49 2.54];

s = [9.54 9.33 9.52 9.06 8.05 7.27 6.03 5.80 5.51];

p = [0.00 0.00 0.00 0.01 0.18 0.35 0.49 0.53 0.55];

Results:

Root of Mean Square Error (RMSE): 0.250520910244571

Sum of Squared Residual: 1.50625743527444

Correlation Coef. (R): 0.981509121545543

R-Square: 0.963360155677103

Adjusted R-Square: 0.914727231437843

Determination Coef. (DC): 0.942399473562309

F-Statistic: 14.7361249635671

Parameter Best Estimate

-------------------- -------------

umax 0.117013263727751

ks -5.83263810513223

yxs 0.699511367324087

yps 5.01252941213249

##### 1 Comment

Vinoj Liyanaarachchi
on 21 Oct 2019

Alex Sha
on 21 Oct 2019

##### 0 Comments

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!