## Fitting Monod Equation with ODE45 to data using lsqcurvefit function

### StarSign1997 (view profile)

on 22 Apr 2019
Latest activity Answered by Alex Sha

### Alex Sha (view profile)

on 21 Oct 2019 at 8:42

### David Wilson (view profile)

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

Star Strider

on 22 Apr 2019

### David Wilson (view profile)

on 22 Apr 2019
Edited by David Wilson

### David Wilson (view profile)

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: Vinoj Liyanaarachchi

### Vinoj Liyanaarachchi (view profile)

on 19 Oct 2019 at 4:30
First of all, I would like to thank StarSign1997 and David Wilson
StarSign1997 Can you kindly provide the final code.... It would be a great help.
I also have to model the same set of equations. Experimental data are as follows:
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]';
I used the code provided by David Wilson (I don't know how to disperse the intial values, So I used same code provided)
Curves were not fitted to the experimental data well. StarSign1997

### StarSign1997 (view profile)

on 19 Oct 2019 at 13:27
Ahh yes. I tried your data and we have the same results. Here is the script I used in the command window. The function parafun1 is the same as the one posted.
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]';
plot(t,[x,s,p],'o-')
umax = 0.5;
ks = 2.5;
Yxs = 5;
Yps = 1.2;
xt0 = [0.06;9.54;0.00]; % read off from data plot [x0 ,S0, p0]';
p0 = [umax; ks; Yxs; Yps; xt0]
ypred = paramfun1(p0,t); % measurement predictions
options = optimoptions ('lsqcurvefit','FiniteDifferenceStepSize',1e-4,'FiniteDifferenceType','central');
[pfit, presnorm, presidual] = lsqcurvefit(@paramfun1,p0,t,[x,s,p],[],[],options);
pfit = lsqcurvefit(@paramfun1,p0,t,[x,s,p])
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 16],xt0);
plot(tspan,a(:,1),'-',tspan,a(:,2),'--',tspan,a(:,3),'k-.', ...
t,x,'o',t,s,'+',t,p,'s')
title('Monod Equation Simulation')
xlabel('Time')
ylabel('Concentration')
legend('Simulated Biomass','Simulated Substrate','Simulated Ethanol','Biomass (Actual)','Substrate (Actual)','Ethanol (Actual)')
Vinoj Liyanaarachchi

### Vinoj Liyanaarachchi (view profile)

on 19 Oct 2019 at 16:17
Thank you StarSign1997 !!
Are there any techniques to improve the parameter estimation? Above shapes of the curves are not acceptable in my case (microalgae cultivation).

### Alex Sha (view profile)

on 21 Oct 2019 at 7:33

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

### Vinoj Liyanaarachchi (view profile)

on 21 Oct 2019 at 7:50
Thank you Alex Sha. Actually I'm new to matlab.
Can you share the code..