Monod kinetics and curve fitting

Hi to everyone!

I think i should start immediately with a short description of what i' m dealing with. I have a system of two (2) equations. Generally, these equations are expressed as following:

1. dS/dt = (mmax * X * S)/(Y * (Ks + S)
2. dX/dt = (mmax * X * S)/ (Ks + S) - (b * X)

S = So and X = Xo

In the laboratory, i' ve measured concentration values versus time, e.g. (t, S). Moreover, i know the initial value Xo. As a result, i' m trying to estimate all the other parameters.

So the problem is how do i solve the differential equation system and then use this solution for the curve fitting procedure? I tried simbiology and the absence of a constraint search for the best solutions leads to estimates with no physical meaning.

What should i try next? Perhaps the solution is associated with non linear fitting tools, but i would like a more specific answer by an experienced user, before i start the effort.

Thank you in advance!

[Star Strider]

PLEASE do not post new Questions here as Comments or Answers related to this Question.

Post new Questions, and copy the URL here as a reference.

Answer 1

[Star Strider]

The differential equations may have long since been integrated and published:

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC291754/?tool=pubmed

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC240226/?tool=pubmed

so all you need to do is fit them to your data. (They are both free PDFs.) I didn't look through the lists in the articles to be sure they have specifically integrated the equations you're interested in fitting.

I can't claim to be an experienced user with respect to fitting chemical kinetics data, since I haven't done anything with chemical kinetics in a while, but I'm familiar with MATLAB's nonlinear curve fitting routines. Either ‘nlinfit’, ‘lsqcurvefit’, or others should be able to estimate the parameters you want. They can also provide confidence intervals on the parameter estimates with ‘nlparci’ and the fitted estimates with ‘nlpredci’.

[kknd]

At first, thank you for your question.

Yes, you are correct. This type of kinetics has an integrated form when expressed simply as in my question. However, if we decide to describe more complicated microbial mechanisms (toxicity, inhibition etc), the two differential equations become quite complex. Therefore, they have been solved numerically only. Not to mention when we have to describe the growth of more than one bacteria.

This is why i insisted on how we could use the numerical solution of the ode in a curve fitting process of our data.

In any case, thanks again!

[Star Strider]

It was easy for me to code your differential equations and get ‘ode45’ to integrate them, and I found some example data to test. Passing all the parameters through ‘lsqcurvefit’ so the objective function will integrate the equations and give ‘lsqcurvefit’ the output it needs to fit the data is proving to be more difficult. The MATLAB documentation is unfortunately not helpful. The problem seems to be passing the time vector and the parameter vector (and other parameters such as the ODE initial conditions) from ‘lsqcurvefit’ to the objective function that also does the integration, so that the objective function will return a vector of estimates of ‘S(t)’. It’s likely not a difficult problem, but one I haven’t been able to solve in the few hours I’ve been working on it.

[Star Strider]

The solution was as obvious as I thought it would be. (I wasn't passing the arguments in the correct order.) The function ‘MonodKinetics1’ goes into its own file. The call to ‘lsqcurvefit’ goes in your main routine. ‘Time’ is obvious and ‘Sdata’ are the ‘S(t)’ data you collected. Note that I wrote my code for ‘Time’ and ‘Sdata’ being column vectors.

I don’t have your data to test, but this code works with some data I found online, and ‘lsqcurvefit’ works well with it.

The calling statements:

B0 = rand(4,1) * 100;

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@MonodKinetics1,B0,Time,Sdata);

The objective function ‘MonodKinetics1.m’:

function S = MonodKinetics1(B, t)
% MONODKINETICS1 codes the system of differential equations describing one
%   version of Monod kinetics:
%       dS/dt = (mmax * X * S)/(Y * (Ks + S));
%       dX/dt = (mmax * X * S)/ (Ks + S) - (b * X);
%   with:
%       Variables:  x(1) = S,  x(2) = X
%       Parameters: mmax = B(1),  Y = B(2),  Ks = B(3),  b = B(4)

x0 = rand(2,1);

[T,Sv] = ode45(@DifEq, t, x0);

      function dS = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = (B(1) .* x(2) .* x(1)) ./ (B(2) .* (B(3) + x(1)));
          xdot(2) = (B(1) .* x(2) .* x(1)) ./ (B(3) + x(1)) - (B(4) .* x(2));
          dS = xdot;
      end

S = Sv(:,1);

end

The ‘Answers’ site wrapped some lines, so be sure to check for that and make the appropriate changes if you copy-paste the code I provided above.

[kknd]

Star Strider you were a real help here! This was exactly what i needed. Nicely done!!!

[Star Strider]

My pleasure! I suggest three changes I thought of overnight. First, add a third argument to create:

function S = MonodKinetics1(B, t, ~)

Second, change B from a 4x1 vector to a 6x1 vector in your calling program. Then change the x0 assignment to:

x0 = B(5:6);

Third, replace the output assignment for S to:

if nargin == 2
    S = SX(:,1);
elseif nargin == 3
    S = SX;
end

(I changed ‘Sv’ to ‘SX’ because it's the [n x 2] matrix consisting of [S X].) Then if you want to see what both S(t) and X(t) are, add a third argument to the call in your main program to something like:

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@MonodKinetics1, B0, Time, Sdat, [], [], MonodKin1Opts);

FitData = MonodKinetics1(B,Time,1);

so that:

S = FitData(:,1);

X = FitData(:,2);

These changes don't affect the function of ‘MonodKinetics1’, but make it easier to use and more robust in its parameters estimation.

[hna]

Hello,

I have just read this old question dated from 6 years ago. How would you plot the data to have the Sdata against Time and the curve fitted line generated from the lsqcurvefit() funtion ?

Thanks in adance,

[Star Strider]

With the estimated parameters ‘B’:

tv = linspace(min(t), max(t), 250);
Mfit = MonodKinetics(B,t);

figure
plot(t, Mfit)
grid

[hna]

Thank you very much for this quick answer

++

[Star Strider]

My pleasure!

A vote for my Answer would be appreciated!

[Michele Notarantonio]

Hello,

I have another type of problem similar to the one explained here. My problem is " I have a set of equations (biochemical reactions) and, initially, I want to simulate the mathematical model with parameters choosen by literature" How the code should change? Pratically I need to solve the combined system of algebraic and differentially equation but I'm finding a lot of trouble. I have written all the equations that I need to simulate the model (11 differential equation and 3/6 algebraic equations to define kinetics rate). I tried to simulate them with ode15s and ode 45 but the results haven't physycally meaning (i'm finding some negatives values). What should I do?

Thanks for the help!

[Star Strider]

Note that lsqcurvefit (and several other optimisation functions) can use parameter constraints, so setting the lower bound to a vector of zeros will prevent negative parameters.

[TheGymnopedist]

Hi, thanks for the useful solution; however I am getting an error.

"Error using lsqcurvefit. Function value and YDATA sizes are not equal."

This is because although the 2 matricies have the same number of columns (23), but different numbers of rows between ydata (8) and the function value (number varies due to ode solver).

Do you have any suggestions as to how to solve this issue?

[Star Strider]

@Ben Mayoh —

The result of the integration has to match the size of the dependent variable matrix. The ODE solver should not have any effect on this, however if the solver detects a singularity, it will stop and the resulting matrix sizes (specifically the row sizes) will not be equal. The solution for that is to keep trying different initial parameter estimates to lsqcurvefit until you get a set that works. (Functions in the Global Optimization Toolbox can make this easier.) Fitting an (Nx23) data matrix is likely going to be a problem.

Later versions of this approach (for example Parameter Estimation for a System of Differential Equations) estimate the initial conditions as well, so consider that option.

Without knowing more about your particular data and model, I cannot offer specific guidance. For that, it would be best for you to provide that in a new question post rather than pursue it here, posting he code for your differential equations and attaching your data. Also consider using one of the stiff solvers (for example ode15s) in the event that your paraameters have widely-varying magnitudes (several orders-of-magnitude difference). And of course, check your differential equations code for errors, since that can be a source of problems.