Parameter Estimation for a System of Differential Equations

Hello. Ok, so I'm new to matlab and I've got a question regarding parameter estimation for a kinetic model. I have 4 different reactants and their concentrations are c1, c2, c3 and c4. I also have 4 differential equations, each one related to a concentration (c1, c2, c3 and c4, respectively -see below-) and experimental data for all these concentrations on 12 different times plus the initial condition. The k's are the rate coefficients. I want to solve this system of ODE's using ode45 and then use the output to compute the experimental data minus the observed data and use these results to estimate the values of k's using lsqnonlin, but apparently I can't solve these ODE's without numerical values for k -which is what I want to know-. Any help on how to set up the command to solve this?

    function dcdt=batch(t,c,k)
dcdt=zeros(4,1);
dcdt(1)=-k(1)*c(1)-k(2)*c(1);
dcdt(2)= k(1)*c(1)+k(4)*c(3)-k(3)*c(2)-k(5)*c(2);
dcdt(3)= k(2)*c(1)+k(3)*c(2)-k(4)*c(3)+k(6)*c(4);
dcdt(4)= k(5)*c(2)-k(6)*c(4);
end

Data:
t  c1  c2  c3  c4
0  1  0  0  0
0.1  0.902  0.06997  0.02463  0.00218
0.2  0.8072  0.1353  0.0482  0.008192
0.4  0.6757  0.2123  0.0864  0.0289
0.6  0.5569  0.2789  0.1063  0.06233
0.8  0.4297  0.3292  0.1476  0.09756
1  0.3774  0.3457  0.1485  0.1255
1.5  0.2149  0.3486  0.1821  0.2526
2  0.141  0.3254  0.194  0.3401
3  0.04921  0.2445  0.1742  0.5277
4  0.0178  0.1728  0.1732  0.6323
5  0.006431 0.1091  0.1137  0.7702
6  0.002595 0.08301 0.08224 0.835

Thanks in advance!

[Alex Sha]

the results below are much better:

k1 0.756713424357579

k2 0.246983152463589

k3 0.100637714178661

k4 0.24638788199029

k5 0.582535403727263

k6 -0.0269993611093679

[Star Strider]

@Alex Sha — Show the code you used to get those results. The results themselves are otherwise meaningless.

[Torsten]

A negative reaction rate constant (k6) does not make sense.

You will have to constrain the ki to be non-negative.

[Alex Sha]

if want to all parameters are positive， the results will be:

k1 0.757805051886583

k2 0.248942776988651

k3 0.182367452149268

k4 0.462511556999762

k5 0.620286121513934

k6 0

Star Strider has done this kind of job perfectly, I just used another package to obtain the result above, This result can be used to compare with Matlab results。

[Star Strider]

Please do not post new problems to this thread unless they relate directly to it. Post them instead as new Questions.

Votes are always appreciated!

.

[Budda]

Hi Star Srider, How would I change the code if only say the data for first variable ( equation) available in the Igor_Moura.m?

Thank you in advance

[Star Strider]

Budda —

Make appropriate changes to the ‘kinetics’ function (specifically to the ‘C’ variable) to output only the variable you need. In that code, there are 4 columns, so address only the column you want returned in ‘C’.

Answer 1

[Star Strider]

See if the techniques in Monod kinetics and curve fitting will do what you want.

[Igor Moura]

Hi, Star Rider, thanks for your help! Both problems are very alike, those techniques might be what I'm looking for, but I just can't get the code to work :( I created a kinetics.m file with the following code:

function C=kinetics(t,theta)
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(:,1);
end

and then I called it using:

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

but I keep getting the following errors:

Error using odearguments (line 83)
The last entry in tspan must be different from the first entry.
%
Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
%
Error in kinetics (line 3)
[T,Cv]=ode45(@DifEq,t,c0);
%
Error in lsqcurvefit (line 202)
          initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
%
Error in cu (line 26)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
%
Caused by:
  Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

[Star Strider]

I apologise for the delay. It’s been four years since I dealt with the Monod kinetics problem and I had to go back to remember what I’d done.

There were two problems:

1. Your ‘kinetics’ function has to have the parameter variable (here ‘theta’) as the first argument, and ‘t’ as the second;
2. You need to return all the values of your ‘Cv’ vector as ‘C’ since you are fitting all of them.

I made those changes and put all the corrected code, data, and plot calls in the attached ‘Igor_Moura.m’ file, so you can run it by typing Igor_Moura in your Command Window. I won’t post it here. Open it in your Editor to see the changes and my added code.

If you want to return ‘theta’, the concentrations, and other variables to your workspace, the easiest way is to add them as outputs to the ‘Igor_Moura’ function.

The Estimated Rate Constants:

  Rate Constants:
    Theta(1) =  0.89197
    Theta(2) =  0.26069
    Theta(3) =  0.24870
    Theta(4) =  0.48183
    Theta(5) =  0.60977
    Theta(6) = -0.01301

The Plot:

Make appropriate changes in the legend code to describe the concentrations.

[Khan Jan]

Hello. My problem is much similar to this example. I want to plot different markers for each experimental data columns c1 to c4 (square, cross, circle and diamond). How to plot them? My second question:Is it possible to constrain the rate constant value between 0 and 0.80 using LSQCURVEFIT function, if yes how to do it?please guide.

[Star Strider]

‘How to plot them’

See the code in the ‘Igor_Moura.m’ file I attached to my previous Comment.

‘Is it possible to constrain the rate constant value between 0 and 0.80 using LSQCURVEFIT function, if yes how to do it?’

Yes. See the documentation for lsqcurvefit.

[Taimur Akhtar]

Thank you Star strider for quick feedback. I got plots and answer to my 2nd question is solved. However, my point relevant to markers is still unclear. I will try to clarify my ambiguous question. Based on code from Igor_Moura.m, i got single marker 'pentagram' with different colors for different concentrations (attachment 2), however my target is to get more than one marker assigned to C1, C2, C3 and C4, respectively (attachment 1). Is it possible to get 4 different markers with 4 different colors assigned to each 4 concentrations? In simple words, is it possible to replace attachment 2 by attachment 1?

[Khan Jan]

I tried to solve them for different markers, but its seems i can get only pentagram with different colors. Star strider, i will appreciate your guidline for solving it.

[Khan Jan]

Both answers Solved, thanks star.

[Taimur Akhtar]

I got plots after solving ODE'S. Now, I need to fit a power law line in the plot. Any idea.

[Wouter Van Hecke]

Hi all,

Great post.

Do you have an idea on how to deal with incomplete data sets? Assume matrix c in the "Igor_Moura.m" file is incomplete and contains some unknown "NaN" values. In that case, lsqcurvefit gives the following error:

"Objective function is returning undefined values at initial point. lsqcurvefit cannot continue. " One option is to remove the entire row(s) with the unknown values, but aren't there more elegant approaches for this problem?

Many thanks for your suggestions!

Best regards, Wouter

[Star Strider]

‘One option is to remove the entire row(s) with the unknown values, but aren't there more elegant approaches for this problem?’

Not that I’m aware of. The only other approach would be to use interpolation, and that has the extreme hazard of creating data that don’t actually exist, especially if the data are noisy or changing rapidly, making the parameter estimates essentially worthless. The Statistics and Machine Learning Toolbox functions ignore the rows with NaN values. I would follow that convention.

[SHUBHAM PATEL]

Hello everyone, I had a very similar problem and used the above code by Star Strider for the fit. One of the problem i am facing is that, i am getting negative values for the parameters which for my specfic problem is not physical.

How do i constrain my parameters to give only postive values.

Thank you.

[Star Strider]

Set the ‘lb’ (lower bound) argument in lsqcurvefit (or other optimization routines that allow bounded parameters) to a zeros vector the size of your parameter vector.

In the context of this code, that would be:

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

[Gustavo Lunardon]

Could someone show me how to do the same minimization, (same data) but with fminsearch? I was trying to adapt the example for a while and tried to pass the variables and parameters in different ways but I keep getting :

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

My attempt is in attachment.

[Star Strider]

Use this:

[theta,fval,exitflag] = fminsearch(@(theta)norm(c - kinetics(theta,t,c0)),theta0);

and:

Cfit = kinetics(theta, tv, c0);

I suggest optimising the initial conidtions as well, rather than using fixed initial conditions.

[Gustavo Lunardon]

Thank you very much Star Strider!

I accompany many of your answers and they are always helpful, but answering this question in so little time exceeded my expectations!

[Star Strider]

My pleasure!

Thank you!

I cannot always respond quickly, however I do my best!

[Daphne Deidre Wong]

Hi Star Strider.

My the set of equations that i need to solve are very much similar to the file attached above 'Igor_Moura.m'. However I was wondering how should i implement the code if i want the reaction to occur over a range of temperature?

or would i need to put in the experimental results data at the certain temperature to get my theta value @ that temperature?

Kind regards, Daphne

[Star Strider]

It depends.

If they are functions of temperature only and not of time, then use temperature as the independent variable and adjust the code accordingly.

If they are functions of time and temperature simultaneously, and both vary, that requires a completely different approach since the system is one of partial differential equations, not ordinary differential equations. The parameters would then not be scalars, but functions of temperature. That would be very difficult. One option would be to solve them with respect to time at each temperature (in a loop) as ordinary differential equations, and then later use appropriate parameter estimation techniques to characterise each parameter as a function of temperature. (That would likely be the approach I would use, at least at the outset, to determine if it was appropriate.)

[Daniel Zhou]

Hi, Star.

I'm using the same code to fit data now. But I want to print this Rsdnrm value. How can I print that value?

[Star Strider]

Daniel Zhou — The ‘Rsdnrm’ output is a scalar double value, and exists in the workspace after lsqcurvefit completes its estimation. Just use fprintf or sprintf to print it.

[Star Strider]

ana ramirez de las heras’s Answer moved here —

Hi Star Strider,

Do I have to know c() values to use this method? I have a similar problem and I have data from the exercise statement but without knowing the values of theta() I don't know how to obtain c() values...

Thank you in advance!

[Star Strider]

ana ramirez de las hera —

In order to estimate the parameters, the ‘c’ data (including the time vector associated with it) have to exist in your workspace. The code in my Answer estimates the kinetic parameters ‘theta’ vector from the data, so knowing them in advance is not necessary, however approximate estimates of their amplitudes are necessary, since the parameter estimation routines require an initial estimate of the parameters.

If you already have the parameters and are simply simulating the system, the ‘c’ values are not necessary.

If this Comment does not address your problems, please post a new Question with more detail as to what you are starting with and the result you want.

[ana ramirez de las heras]

Thank you for answering, Star. I have posted a new Question: https://es.mathworks.com/matlabcentral/answers/685053-ode-system-solution-with-unknown-constant

[Star Strider]

My pleasure!

[Yasin Islam]

Hi Star Strider,

thank you very much for the usefull solution to the problem! I adapted the code to fit for my 6 diff. eq. and 12 k-values. it works.

unfortunately i not only have unknown k-values but also unknown concentrations.

The problem in my case extends up to 18 concentrations with 18 differential equations (i do know time dependend concentration values for 6 of them)

How can I update the code to not only look for my k-values but also for the unknown concentrations in order to represent the experimental concentrations best? Although alot of the unknown concentrations will be equal to 0 and/or neglectable.

[Star Strider]

Yasin Islam — You can obviously fit only the concentrations you have data for. If the other 12 concentrations (actually, differential equations describing them) can be expressed in terms of the parameters (‘k values’) you can estimate, simply simulate them outside the ‘kinetics’ function using the parameters and plot all of them.

If the other 12 differential equations involve parameters that you cannot estimate (or cannot derive from the parameters you can estimate), then you cannot simulate them. There is no way around that unfortunate situation.

[si heon Seong]

@Star Strider

Hi Star Strider,

first, thank you for introducing this page! It was great help to me

One thing I have a additional question about this example

What should I do to give each parameter the constraint I want? For instance, thetas must be bigger than 0 etc..

Thank you

[Star Strider]

si heon Seong — The lsqcurvefit and many other optimisation functions allow parameter constraints, at least setting the upper and lower parameter bounds. Use those arguments.

Here, that would be:

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

with the zeros call setting the lower limit. (The upper limit is not set here, since it is not necessary.)

[si heon Seong]

Dear @Star Strider

I really appreciate your help! Can I ask the last one?

I have a model of 6 coupled differential equations for c1, c2, c3, c4, c5 and c6

but there are only two data set for c1 and c2

Can I still do the fitting job?

Always thank you

[Star Strider]

si heon Seong — You can only fit ‘c1’ and ‘c2’ to your data, so in ‘kinetics’:

C = Cv(:,1:2);

(assuming that ‘c1’ and ‘c2’ are the first 2 differential equations in your system).

There have to be more data (more rows in the data matrix and time vector) than there are parameters to estimate.

My pleasure!

[Khoi Ly]

Hi @Star Strider,

I was directed by your answer in another post regarding a similar problem.

https://www.mathworks.com/matlabcentral/answers/751489-parameter-estimation-for-a-linear-differential-equation-from-multiple-measured-time-series-data?s_tid=srchtitle

In that post, I was trying to estimate the parameters of the differential equations using more than one time series data sets.

In the example above, given.

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

we can estimate the parameters theta. Since the example above uses only one data set, what should I do if I have multiple data sets? Can I just concatenate the time-series data sets together and run the lsqcurvefit as usual?

Also, in my problem, each different time-series data set contains a different initial value. How do I modify the example to allow for different initial values depending on the chosen time-series data set?

Thank you

[sumit kumar]

@Star Strider Thankyou so much. Your replies on various problems have helped me a lot. I was stuck in a parameter estimation problem and your answer is the exact thing i was looking for. You are a saviour man..

[Star Strider]

My pleasure!

A Vote is always appreciated!

[Star Strider]

UPDATE

Slightly revised code with some aethetic tweaks so that the lines and markers colours now match —

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 = 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:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

		Theta( 1) =  0.76483
		Theta( 2) =  0.23492
		Theta( 3) =  0.20878
		Theta( 4) =  0.49179
		Theta( 5) =  0.62211
		Theta( 6) =  0.00000
		Theta( 7) =  0.90288
		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(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

.

[uki71319]

Dear @Star Strider, thanks for your detailed explanations. I've followed your codes, but i still got some problmes. Thus, i posted a new question under this link and hope to find out the solution for my lsqcurvefit problem. Would you mind taking a look and see where the issues are?

https://www.mathworks.com/matlabcentral/answers/1833763-parameter-estimation-for-a-set-of-odes-with-structs-input

Any advice and help would be greatly appreciated. Thanks in advance!

[Star Strider]

I looked at it just now, however it has already been answered.

Note that initial parameter estimates of zeros is never a good approach. Choose something nonzero and preferably not the same values for both parameters.

[Maria Teresa Aguirrezabala Campano]

@Star Strider may I also ask for your help with my problem? https://www.mathworks.com/matlabcentral/answers/1900830-double-monod-maximum-rate-estimation

[Jan]

@Maria Teresa Aguirrezabala Campano: Imagine what happens if all users advertise their questions and ask specific members to answer them. These member will be cluttered by messages and they will find less time to answer questions. Please consider this.

[Continuum]

Hi @Star Strider, thanks for your activity on the platform. I have question;

I have a set of ODEs (attached), I have been able to solve them using ode45, however, my issue now is my experimental results don't match the integrated values of the equations. So, I am looking to fit only the solution for epsilon with it's experimental results to find the best parameters A, B, (A0/alpha), k0, Q, and QG. Is this achievable using the Igor_Moura.m script? If so, which modifications are necessary? Like I said, I have already solved for G, theta, and epsilon, but I only have experimental data for epsilon. Thank you.

[Star Strider]

If you only have one variable (column) to fit, that being epsilon, choose that variable as ‘C’.

So if is ‘dcdt(3)’ then:

C = Cv(:,3);

That should work, providing that you have more (preferably at least twice as many) data than the number of parameters (6) you want to estimate.

[Star Strider]

You have to write the ‘kinetics’ function.

First, write all the differential equations, and then verify that the differential equation system works with parameters of the appropriate magnitudes (they do not have to be exact, since that is the point of the parameter estimation problem). It would likely be best to use ode15s or ode23s to integrate the differential equation system, since it is likely a ‘stiff’ system (most kinetics systems are ‘stiff’).

Writing the rest of the ‘kinetics’ code after that should be straightforward. Include the temperature in your differential equation system as independent variable data. It would likely have to be passed in the same independent variable argument as the time vector:

function C = kinetics(theta, tT)

where:

tT = [t T];

(they are both column vectors), then separate them inside the ‘kinetics’ function. This is your project, so you get to make all the design decisions. Use the ‘kinetics’ function here as a guide.

The parameters you want to estimate have to be presented to lsqcurvefit as a single vector, although you can assign them meaningful names inside the ‘kinetics’ function. That might be something like this:

ko = theta(1);

QG = theta(2);

and so for the rest of them.

I plotted the data, including the extra ‘exp_epsilon’ vector. I do not see any problems with it.

[Continuum]

Sorry, I meant to include my code attached. Thanks.

[Star Strider]

You need to put your code in the same forma as the ‘kinetics’ example (although you can call it whatever you like, providing the naming is consistent throughout your code). As I mentioned previously, put ‘t’ and ‘T’ together in the same independent variable matrix, use that as the independent variable in your call to lsqcurvefit, then separate them in the ‘kinetics’ function into the time vector and the matching temperature vector.

It may be necesary to interpolate the temperature vector to match the time value that the ODE solver generates. See the documentation section on ODE with Time-Dependent Terms (this is in the ode45 documentation, however it works for all the solvers) for an esample on how to do that.

Choose the appropriate vector for ‘epsilon’ that you want to use as the dependent variable, and then return as ‘C’ the appropriate column of the integrated differential equations to fit to it.

[Continuum]

Thanks @Star Strider. In the case of the Igor_Moura.m, what is the essence of t and c values from lines 20 to 44? Are they the solution/output from ode45? If so, why do you need to explicitly copy them over within the function? Is it possible to use the outputs t and c directly from the ode45 in the 'lsqcurvefit'?

[Continuum]

I have tried to write it in your format, but the code just runs continuously which makes me think something is not right. Could it be that the there are too many parameters to fit? Attached is what I tried to run.

P.S. To my earlier question, I figured the t and c were from experiments or what the model is calibrating against?

[Star Strider]

This should actually be a new question. If you post it as such, and post the link here, I will post my edited version of you code.

I finally got it to work (run without error), however the fit is not good, and that will likely require using the genetic algorithm or other Global Optiomization Toolbox functions that can change the parameters (or initial parameter estimates) until it can get a good fit.

[Continuum]

OK @Star Strider, here is a link to my question;

https://www.mathworks.com/matlabcentral/answers/2074946-how-to-find-fitting-optimized-parameters-to-fit-a-system-of-non-linear-odes-to-experiment

Thank you.