How to estimate parameters for a system of ODE's by using replicates data ?

Hello

I have a conducted a triplicate experiment for oxidation of two compounds over couple of months. I have to estimate the parametrs (reaction orders, rate reactions, and stoichiometric coefficients). I have a system of three ODEs. Can you please advise what matlab tool to use? (using R2019a).

Thanks

Answer 1

[Star Strider]

There are several examples, two being: Parameter Estimation for a System of Differential Equations and of course Monod kinetics and curve fitting.

Note that non-specific Questions result in non-specific Answers.

[gina_et_berg]

thanks for yu reply. I have actually used Monod Kinetics (great and very detailed) one to write my codes. However, it does not produce output. My code is attached for only useing one set of data (not all three) and even with this one I do not get results. The odes are:

dc1/dt = -beta2*k2*c1^ox2*c2^m2 - beta3*k3*c1^ox3*c3^m3

dc2/dt = -k2*c1^ox2*c2^m2

dc2/dt = -k3*c1^ox3*c2^m3

The initial concentration of all compounds are known.

I have two more set of data for all three compounds (running 3 similar experiments at the same time).

Thanks so much for your help.

[Star Strider]

This runs, however you will need to experiment with different initial parameter estimates because the initial parameter estimates you are currently using are apparently inappropriate, and because of that, lsqcurvefit cannot estimate the parameters correctly. All nonlinear parameter estimation routines are very sensitive to the initial estimates, so experiment with them.

The (Revised) Code —

% DO NOT USE GLOBAL VARIABLES!
observations = [ 0	1	2	3	4	5	6	7	11	14	21	35	56 % time points
4135.1	2934.6	785.4	1677.0	2047.1	1608.7	1576.8	1380.9	1353.3	1251.3	1007.4	1024.0	772.9 % F2 concentration
4158.4	2816.8	2026.8	1701.7	1794.2	1736.1	1286.2	1484.8	1072.7	1302.8	1173.1	1166.9	850.5 % F3 concentration
3791.0	2410.1	2062.9	1372.0	2087.4	1708.3	1509.5	1390.0	1162.2	1316.6	1193.5	1017.2	751.0]'; % OX concentration

t_Obs = observations(:,1);
F2_Obs = observations(:,2);
F3_Obs = observations(:,3);
OX_Obs = observations(:,4);
Obs = [OX_Obs F2_Obs F3_Obs];

%parameters
ox_F2 = 4; %exponent of concentration of ox in chemical reaction with F2
m_F2 = 1; %exponent of concentration of F2 in chemical reaction with ox
ox_F3 = 4; %exponent of concentration of ox in chemical reaction with F3
m_F3 = 1; %exponent of concentration of F3 in chemical reaction with ox
beta2 = 1; %number of ox moles consumed for each mole of F2 reacting
beta3 = 1; %number of ox moles consumed for each mole of F3 reacting
k2 = 0.1;
k3 = 0.1;

%Initial concentrations
cin_ox = 40; %g/L initial concentration of added oxidant 
cin_2 = 4.14*1e-3; %g/L detected concentration of F2 in collected groundwater
cin_3 = 4.36*1e-3; %g/L detected concentration of F3 in collected groundwater

%Initial conditions (Fill array with correct values)
cstart = zeros(3,1); %Create array of ALL initial concentrations
cstart(1) = cin_ox; 
cstart(2) = cin_2;
cstart(3) = cin_3; 



X0 = [m_F2; m_F3; ox_F2; ox_F3; k2; k3; beta2; beta3; cstart];

% %durataion of the experiment
% tstart = 0;
% tend = 56; %days
% tspan = [tstart tend];


options=optimoptions(@lsqcurvefit,'Algorithm','levenberg-marquardt','StepTolerance',1e-16)
B = lsqcurvefit(@(b,t)odefit(b,t,cin_ox,cin_2,cin_3,cstart), X0, t_Obs, Obs, zeros(size(X0)), inf(size(X0)), options);
fprintf(1, '\n\tKinetic Orders:\n\t\tm_F2 = %11.3E\n\t\tm_F3 = %11.3E\n\t\tox_F2 = %11.3E\n\t\tox_F3 = %11.3E\n\t\tKinetic Constants:\n\t\tk2 = %11.3E\n\t\tk3 = %11.3E\n\t\tcoefficients:\n\t\tbeta2 = %11.3E\n\t\tbeta3 = %11.3E\n', B);

dcdt = @(t,x,b) [ -b(5)*b(7)*(x(2)^b(1))*(x(1)^b(2)) - b(6)*b(8)*(x(3)^b(3))*(x(1)^b(4)); %oxidant rate reaction
          -b(7)*(x(2)^b(1))*(x(1)^b(2)); % F2 rate reaction
          -b(8)*(x(3)^b(3))*(x(1)^b(4))]; % F3 rate reaction
%[T,X] = ode23s(@(t,x) dcdt(t_Obs,x,B(1:8)), t_Obs, B(9:11));
[T,X] = ode23s(@(t,x) dcdt(t_Obs,x,B(1:8)), t_Obs, cstart);

figure
plot(observations(:,1), observations(:,2:end),'p')
hold on
plot(observations(:,1), odefit(B,observations(:,1),cin_ox,cin_2,cin_3,cstart))
hold off
grid


function X = odefit(b,t,cin_ox,cin_2,cin_3,cstart)
% global cin_ox cin_2 cin_3 tspan cstart
      % b(1:8) = Parameters,  b(9:11) = Initial Conditions
      % MAPPING: x(1) = c1,  x(2) = c2, x(3) = c3,  b(1) = m_F2,  b(2) =  ox_F2
      % b(3) = m_F3, b(4) = ox_F3, b(5) = beta2, b(6) = beta3, b(7) = k2,
      % b(8) = k3, b(9) = cin_ox, b(10) = cin_2, b(11) = cin_3
      
      dcdt = @(t,x,b) [ -b(5)*b(7)*(x(2)^b(1))*(x(1)^b(2)) - b(6)*b(8)*(x(3)^b(3))*(x(1)^b(4)); %oxidant rate reaction
          -b(7)*(x(2)^b(1))*(x(1)^b(2)); % F2 rate reaction
          -b(8)*(x(3)^b(3))*(x(1)^b(4))]; % F3 rate reaction
      [~,X] = ode23s(@(t,x) dcdt(t,x,b), t, cstart);
 end

I eliminated the global variables and passed the necessaary parameters. See the documentation section on Passing Extra Parameters for an extended discussion.

See this Comment for an example of using the ga (genetic algorithm) function to estimate the parameters. It searches the entire parameter space to find an acceptable set of parameters (although sometimes it may take more than one run if the problem is not well-posed).