Kinetic Models (Monod Equations)

Hello,

1. I saw the equations below in litrature. This is a process model for cultivation/fermentation process. Question is, being new to Matlab, how do I code this and fit it to experiemtal data?
2. Can these equations be implemented as is in Simulink and the unknown parameters estimated based on the test/experiment data? like can Simulink solve these equations? I am not sure how to implement this and would any guiance you can provide.
3. Example video from MathWorks is How to Estimate Model Parameters from Test Data with Simulink Video - MATLAB & Simulink (mathworks.com).

Thanks!

[Ghazwan]

Hello! Yes, you can use both Matlab and Simulink. With Matlab, you can use ode45 and fminsearch to fit the equations to the experimental and find the parameters. Basically, the optimization tool functions.

Answer 1

[Star Strider]

See: Monod kinetics and curve fitting

A somewhat improved version of that is in: Parameter Estimation for a System of Differential Equations

.

[Learning]

As always, I appreciate your response! I don’t understand all aspects of the code in the link your provided. In the screenshot I provided, for example how would I code it and solve it using ODE45? Also if I have 4 differential equations, each for a specific concentration, how would I make predictions when I get new data to see if the model is working?

[Star Strider]

As always, my pleasure!

It might be necessary for you to re-code your differential equations if they differ from those in my original code (that simply illustrates a way to solve the problem).

The code in my first reference deals specifically with Monod kinetics. My code in the second reference demonstrates an improved version of the code and displays the fittted parameter responses with the data. If you want to make predictions for times beyond the original tme vector, run the same code with the estimated parameters with a new time vector that extends the original time vector. That would require saving the estimated parameters and using them with the existing code, with a new time vector.

If you output the result of the differential equation integration to a structure instead of the way I have the results, you can use the deval function with a new time vector to plot predictions beyone the existing time vector. That would require a slight re-write of my existing code to display the result with the existing time vector and data, however it is certainly an option.

[Learning]

Thank you again! 1. I see in your code you have Theta0 =[1;1;1;1;1;1]. Is there a reason why all of them are 1s? Aso can the theta or constants be estimated with real data using parameter estimator app in Matlab?

2. Also, how does one ‘validate’ the Monod equation/model with a different experimental data? How does one use this Monod model in real time to estimate concentrations of analytes (assuming each equation is for a certain analyte such as biomass) of interest?

Thank you!

[Learning]

Also, if the Monod model is to simulate the fermentation process and the concentrations of certain materials of interest at specific time point (t), is there a way to tell the algorithm/code to pause or produce a prediction say every 5 minutes (of course depending on the sampling interval of the process) so another measurement from another instrument can be compared with the value estimated by the Monod model instead of the model outputting all the estimated values at ones for all the time interval given to it?

For the current code, all estimated concentrations are produced at once. Say I have an instrument the produces concentrations of the material of interest (represented by one of the differential equations) but the instrument output is every 5 min. Hence can the code be made such that the estimated concentration by the Monod code be output every 5 min so the instrument value can be compared to it?

Thanks!

[Star Strider]

As always, my pleasure!

The ‘[1;1;1;1]’ are simply starting parameter estimates. They can be anything reasonable, although they should have a lower bound of zeros to prevent them from becoming negative (since kinetic parameters are always positive by defintion in most models). My code most recently estimates the initial concentrations as well as the kinetic constants, since that produces a better fit to the data. They should also be bounded as positive.

I am not certain what you mean by ‘validating it with different experimental data’. It is likely best to estimate the parameters for each set of data, although it is certainly possible to use a given set of parameters with different data. To use the same parameters, use the new data with a given set of parameters and see how well they estimate the new data. I would need a specific example (and data) to expand on this in order to understand what you want to do.

It is certainly possible to use a different time vector with a given set of estimated parameters, and the code does that in order to produce a smooth curve for the plots at the end. That time vector does not have to end when the original time vector ends, and can easily be longer. This extrapolates beyond the region-of-fit, and while that is not recommended (because that amounts to ‘guessing’ what the model will do beyond the actual data), it is certianly possible.

The original parameter estimation uses only the data at the time instants specified, and those can be produced at the end of the parameter estimation process. Interpolating data at new times within the original interval is permitted (as opposed to extrapolating beyond the original time interval that can be done but is not recommended), and the code does that to draw the plots at the end.

I can provide more specific comments with specific data. I will then have a better idea of what the specific issues are.

.

[Learning]

Thank you again for the detailed response!

Final question: what if I already know those paratmeters in your code? i.e theta(1), theta(2) etc.

so far I can get it to work in this way,

c0 = [2];

tspan = 0:0.5:6

[t,c] = ode45('batchfun',tspan,c0) %batchfun is saved as an .m file

plot (t,c)

function dcdt = batchfun(t,c)

theta(1) = 0.76483; %from your code

theta(2) = 0.23492;

dcdt(1) = -theta(1).*c(1)-theta(2).*c(1);

I however get errors when I try to include dcdt(2), dcdt(3) etc. anything I'm ddoing wrong? I also gets an error when c0 = [1;0;0;0] with the current code I have shown here. Basically if I know all the constants, how would your code looks like.

Thanks!

[Star Strider]

If you already know the parameters, there is no need to do the estimation, so comment out the lsqcurvefit call. Just load them as the ‘theta’ vector (by readint a file containing them or simply pasting the vector into the code) to pass them to the function, similar to what I did after the lsqcurvefit call.

That might be something like this (using your own ‘kinetics’ code) —

function C=kinetics(theta,t)
    c0=theta(end-3:end);
    [T,Cv]=ode45(@DifEq,t,c0);
    %
    function dC=DifEq(t,c)
        dcdt=zeros(4,1);
        dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
        dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
        dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
        dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
        dC=dcdt;
    end
    C=Cv;
end

% LOAD THE 'theta' VECTOR AND THE 'time' VECTOR (OR AT LEAST A TWO-ELEMENT VERSION OF THE 'time' VECTOR) HERE

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
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, compose('C_%d',1:size(c,2)), 'Location','best')
end

I have not specifically tested it, however this should work.

.

[Learning]

Thank you!!

[Star Strider]

As always, my pleasure!

[Learning]

So far been successful using your code. I have tried other equations and having some issues...in your

• Igor_Moura.m file, t and c values are provided in the body of the script/code. t and c are then subsequently provided in the lsqcurvefit call. What if I want to call it from the workspace and not provide t and c in the body of the code? How do i do that? example for the lsqcurvefit function, [theta,Rsdnrm, Rsd, ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c). how do I call t and c from the workspace?

Thank you!

[Star Strider]

Those variables must exist in your workspace. If you have them in a file, read the file to import the data and then use that matrix and time vector in the code. (Beyond that, I am not certain what you are referring to.)

[Learning]

So I have c and t in the workspace and have deleted them from the body of the code m (in your Igor_Moura.m) t and c data is provided in the body of the code. I have deleted them from the body of the code and have created t and c in the workspace. However it give error when I run the code. Not sure why it doesn’t recognize t and c in the workspace. I’m sure I’m doing something wrong or not calling them well into the code but don’t know what. Example of how to call them? Or read them from a file?

[Star Strider]

I need to know more.

If you can copy the code and paste it into a comment here I can then understand the problem. If it involves reading files, I would need to have them provided as well, to be certain they are being imported and referenced correctly.

[Learning]

Sure. Thanks for your patience for this. your code is below:

function Igor_Moura

% 2016 12 03

% NOTES:

%

% 1. The ‘theta’ (parameter) argument has to be first in your

% ‘kinetics’ funciton,

% 2. You need to return ALL the values from ‘DifEq’ since you are fitting

% all the values

function C=kinetics(theta,t)

c0=[1;0;0;0];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);

dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);

dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);

dcdt(4)= theta(5).*c(2)-theta(6).*c(4);

dC=dcdt;

end

C=Cv;

end

t=[0.1

0.2

0.4

0.6

0.8

1

1.5

2

3

4

5

6];

c=[0.902 0.06997 0.02463 0.00218

0.8072 0.1353 0.0482 0.008192

0.6757 0.2123 0.0864 0.0289

0.5569 0.2789 0.1063 0.06233

0.4297 0.3292 0.1476 0.09756

0.3774 0.3457 0.1485 0.1255

0.2149 0.3486 0.1821 0.2526

0.141 0.3254 0.194 0.3401

0.04921 0.2445 0.1742 0.5277

0.0178 0.1728 0.1732 0.6323

0.006431 0.1091 0.1137 0.7702

0.002595 0.08301 0.08224 0.835];

theta0=[1;1;1;1;1;1];

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

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

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))

end

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

figure(1)

plot(t, c, 'p')

hold on

hlp = plot(tv, Cfit);

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

end

now if I remove write it again an remove t and c data from the body of the code, it looks like this below. I then write t and c to the workspace and I try calling it in the code, it doesn't work. One thing also is that, c0 = [1;0;0;0], is this the initial conitions for C(1), C(2), C(3) and c(4)?

c=[0.902 0.06997 0.02463 0.00218

0.8072 0.1353 0.0482 0.008192

0.6757 0.2123 0.0864 0.0289

0.5569 0.2789 0.1063 0.06233

0.4297 0.3292 0.1476 0.09756

0.3774 0.3457 0.1485 0.1255

0.2149 0.3486 0.1821 0.2526

0.141 0.3254 0.194 0.3401

0.04921 0.2445 0.1742 0.5277

0.0178 0.1728 0.1732 0.6323

0.006431 0.1091 0.1137 0.7702

0.002595 0.08301 0.08224 0.835];

t=[0.1

0.2

0.4

0.6

0.8

1

1.5

2

3

4

5

6];

function Igor_Moura

% 2016 12 03

% NOTES:

%

% 1. The ‘theta’ (parameter) argument has to be first in your

% ‘kinetics’ funciton,

% 2. You need to return ALL the values from ‘DifEq’ since you are fitting

% all the values

function C=kinetics(theta,t)

c0=[1;0;0;0];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);

dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);

dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);

dcdt(4)= theta(5).*c(2)-theta(6).*c(4);

dC=dcdt;

end

C=Cv;

end

theta0=[1;1;1;1;1;1];

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

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

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))

end

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

figure(1)

plot(t, c, 'p')

hold on

hlp = plot(tv, Cfit);

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

end

[Star Strider]

I changed the code slightly to reflect updates to my original code, and added a lower bound of zeros to constrain all the parameters to be positive.

This example seems to use the original data, so pasting your data into it where I put ‘c’ and ‘t’ should work.

Update ‘theta0’ to reflect your differential equation system so the length is (number of parameters + number of differential equations).

This version of the code also estimates the initial conditions of the differrential equations as parameters, as well as the kinetic constants, since that generally provides a better fit to the data. The initial conditions are the last N parameters of theta0’ and ‘theta’ where N is the number of differential equations in the system. The rest of the code adapts to the data, with the plot legend matching the columns of ‘c’.

In the original version of the code, the functions and the calling code all had to be within a separate function (or different function files, however that was inconvenient to post here so I put them all in one function and then called that function to execute the code). Several versions/releases ago MATLAB allowed the functions to be in the script, providing they are all at the end of the script, so this version of my original code follows that later convention.

Try this —

c=[0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835];

t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];


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

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality 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.76482
		Theta( 2) =  0.23492
		Theta( 3) =  0.20877
		Theta( 4) =  0.49179
		Theta( 5) =  0.62211
		Theta( 6) =  0.00000
		Theta( 7) =  0.90287
		Theta( 8) =  0.07146
		Theta( 9) =  0.02840
		Theta(10) =  0.00000

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

figure
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, compose('C_%d',1:size(c,2)), 'Location','best')

function C=kinetics(theta,t)
c0=theta(7:10);
[T,Cv]=ode45(@DifEq,t,c0);

function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dC=dcdt;
end
C=Cv;
end

Put your ‘c’ and ‘t’ in place of these, and provide your own version of ‘kinetics’ and it should work.

.

[Learning]

This is very helpful. You’ve done great and I am only grateful. I cannot get Cfit, theta vector and Tv to be output into the workspace. Logic is after finding the rate constants, I want to save these Cfit, tv and theta vector to the workspace so I can plot new data on it without using the Isqcurvefit (will comment out this call in the code).

But thank you brother. I appreciate it.

[Star Strider]

My pleasure!

I do not understand the problems you are having.