Solving ODE System and Estimating Parameters with Experimental Data

Question

I am new to Mathematica. I am trying to find out the best-fit values of the parameters appearing in an ODE system of three equations. The three ODE equations are:

{a'[t] == -k[1]*a[t], b'[t] == -k[2]*b[t] + k[1]*a[t], c'[t] == k[2]*b[t]}

The experimental data set for the three equations with respect to time (in minutes) is:

adata = {{0, 1}, {2, 0.88}, {6, .69}, {10, .53}, {20, .28}, {30, .15}, {50, .043}, {70, .012}, {90, 0}, {120, 0}, {150, 0}, {200, 0}}
bdata = {{0, 0}, {2, 0.12}, {6, .29}, {10, .42}, {20, .56}, {30, .57}, {50, .46}, {70, .33}, {90, 0.22}, {120, 0.12}, {150, .06}, {200, 0.02}}
cdata = {{0, 0}, {2, 0.003}, {6, .030}, {10, .050}, {20, .16}, {30, .28}, {50, .50}, {70, .66}, {90, 0.78}, {120, 0.88}, {150, .94}, {200, 0.98}}

(First entry in each pair is time)
The goal here is to solve ODE, and then optimize the parameters (k1, k2) such that the equations can best describe the experimental data. I would be very thankful to you, if you can spare some of your valuable time and help me with this query. the objective fuction to be minimzed is as follows

But it is giving me the error like 'not a real number' etc.

Answer 1

ClearAll["Global`*"]

adata = {{0, 1}, {2, 0.88}, {6, .69}, {10, .53}, {20, .28}, {30, .15}, {50, .043}, {70, .012}, {90, 0}, {120, 0}, {150, 0}, {200, 0}};
bdata = {{0, 0}, {2, 0.12}, {6, .29}, {10, .42}, {20, .56}, {30, .57}, {50, .46}, {70, .33}, {90, 0.22}, {120, 0.12}, {150, .06}, {200, 0.02}};
cdata = {{0, 0}, {2, 0.003}, {6, .030}, {10, .050}, {20, .16}, {30, .28}, {50, .50}, {70, .66}, {90, 0.78}, {120, 0.88}, {150, .94}, {200, 0.98}};

eqns = {a'[t] == -k1 a[t], b'[t] == -k2 b[t] + k1 a[t],  c'[t] == k2 b[t]};
Thread[{aa[t_], bb[t_], cc[t_]} = DSolveValue[{eqns, a[0] == 1, b[0] == 0, c[0] == 0}, {a[t], b[t],  c[t]}, t]];



model[k1_, k2_] := Sum[(aa[adata[[i, 1]]] - adata[[i, 2]])^2 + (bb[bdata[[i, 1]]] - 
    bdata[[i, 2]])^2 + (cc[cdata[[i, 1]]] - cdata[[i, 2]])^2, {i, Length@adata}]        


 fit = Last@ NMinimize[model[k1, k2], {k1, k2}]

{k1 -> 0.0632016, k2 -> 0.0211064}

Thread[{k1, k2} = Values@fit];

Show[Plot[{aa[t], bb[t], cc[t]}, {t, 0, 200}, Frame -> True], ListPlot[{adata, bdata, cdata}]]

Answer 2

First thing first, when you have data points, always plot and have a look at it.

I replaced your k[1] and k[2] with $ \delta $ and $ \alpha $.

eqns = {a'[t] == -k[1]*a[t], b'[t] == -k[2]*b[t] + k[1]*a[t], c'[t] == k[2]*b[t]} /. {k[1] -> \[Delta], k[2] -> \[Alpha]} 

ListLinePlot[{adata, bdata, cdata}, ImageSize -> Large]

Looks like a bunch of $a e^{x}$ functions to me...And your system of equations is small, they may have an analytical solution.

DSolve[{eqns, a[0] == 1, b[0] == 0, c[0] == 0}, {a[t], c[t], b[t]}, t]

$$ \left\{\left\{a(t)\to e^{\delta (-t)},b(t)\to -\frac{\delta e^{\alpha (-t)-\delta t} \left(e^{\delta t}-e^{\alpha t}\right)}{\alpha -\delta },c(t)\to \frac{e^{\alpha (-t)-\delta t} \left(\alpha e^{\alpha t+\delta t}-\delta e^{\alpha t+\delta t}-\alpha e^{\alpha t}+\delta e^{\delta t}\right)}{\alpha -\delta }\right\}\right\}. $$

Oh, they do. We're in luck! We can write a simple manipulate for fitting.

Manipulate[Plot[{E^(-t \[Delta]), -((E^(-t \[Alpha] - t \[Delta]) (-E^(t \[Alpha]) + E^(t \[Delta])) \[Delta])/(\[Alpha] - \[Delta])), (E^(-t \[Alpha] - t \[Delta]) (-E^(t \[Alpha]) \[Alpha] + E^(t \[Alpha] + t \[Delta]) \[Alpha] + E^(t \[Delta]) \[Delta] - E^(t \[Alpha] + t \[Delta]) \[Delta]))/(\[Alpha] - \[Delta])}, {t, 0.1, 199}, Epilog -> {{Red, Point[adata]}, {Blue, Point[bdata]}, {Green, Point[cdata]}}, ImageSize -> Large], {{\[Delta], 0.0613}, 0.0001, 0.08}, {{\[Alpha], 0.02085}, 0.0001, 0.03} ]

We've found by hand some good initial values for $\delta$ and $\alpha$

We can use this as initial values for a nonlinear model fit.

m1 = NonlinearModelFit[
  adata, {a[t]} /. sol, {{\[Delta], 0.0613}, {\[Alpha], 0.02085}}, t]
m2 = NonlinearModelFit[
  bdata, {b[t]} /. sol, {{\[Delta], 0.0613}, {\[Alpha], 0.02085}}, t]
m3 = NonlinearModelFit[
  cdata, {c[t]} /. sol, {{\[Delta], 0.0613}, {\[Alpha], 0.02085}}, t]

Plot the result

Plot[{m1[t], m2[t], m3[t]}, {t, 0, 200}, PlotRange -> All, 
 Epilog -> {{Red, Point[adata]}, {Blue, Point[bdata]}, {Green, 
    Point[cdata]}}]

and then get some infos about our parameters that NonlinearModelfit found.

m1["ParameterTable"]
m2["ParameterTable"]
m3["ParameterTable"]

$$ \begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline \delta & 0.063221 & 0.000225209 & 280.722 & \text{8.092528367214482$\grave{ }$*${}^{\wedge}$-21} \\ \alpha & 0.02085 & 0. & \infty & \text{0$\grave{ }\grave{ }$323.6072453387798} \\ \end{array} $$

$$ \begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline \delta & 0.0629746 & 0.000365423 & 172.333 & \text{1.0634929126057997$\grave{ }$*${}^{\wedge}$-18} \\ \alpha & 0.0210366 & 0.0000895516 & 234.91 & \text{4.805391875551799$\grave{ }$*${}^{\wedge}$-20} \\ \end{array} $$

$$ \begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline \delta & 0.0640477 & 0.00221151 & 28.961 & \text{5.61521578981255$\grave{ }$*${}^{\wedge}$-11} \\ \alpha & 0.0210622 & 0.000314458 & 66.9795 & \text{1.3404461637525215$\grave{ }$*${}^{\wedge}$-14} \\ \end{array} $$

Histogram[m1["FitResiduals"]]
Histogram[m2["FitResiduals"]]
Histogram[m3["FitResiduals"]]

Answer 3

Using the language resources

adata = {{0, 1}, {2, 0.88}, {6, .69}, {10, .53}, {20, .28}, {30, .15}, {50, .043}, {70, .012}, {90, 0}, {120, 0}, {150, 0}, {200, 0}};
bdata = {{0, 0}, {2, 0.12}, {6, .29}, {10, .42}, {20, .56}, {30, .57}, {50, .46}, {70, .33}, {90, 0.22}, {120, 0.12}, {150, .06}, {200, 0.02}};
cdata = {{0, 0}, {2, 0.003}, {6, .030}, {10, .050}, {20, .16}, {30, .28}, {50, .50}, {70, .66}, {90, 0.78}, {120, 0.88}, {150, .94}, {200, 0.98}};
adata0 = Table[Flatten[{1, adata[[k]]}], {k, 1, Length[adata]}];
bdata0 = Table[Flatten[{2, bdata[[k]]}], {k, 1, Length[adata]}];
cdata0 = Table[Flatten[{3, cdata[[k]]}], {k, 1, Length[adata]}];
transformedData = Join[Join[adata0, bdata0], cdata0];

tmax = 200;
eqns = {a'[t] == -k1*a[t], a[0] == a0, b'[t] == -k2*b[t] + k1*a[t], b[0] == b0, c'[t] == k2*b[t], c[0] == c0};
sol = ParametricNDSolveValue[eqns, {a, b, c}, {t, 0, tmax}, {k1, k2, a0, b0, c0}];

model[k1_, k2_, a0_, b0_, c0_][i_, t_] := Through[sol[k1, k2, a0, b0, c0][t], List][[i]] /; And @@ NumericQ /@ {k1, k2, a0, b0, c0, i, t};
fit = NonlinearModelFit[transformedData, model[k1, k2, a0, b0, c0][i, t], {k1, k2, a0, b0, c0}, {i, t}] // Quiet;

fit["ParameterTable"]

obtaining the results

$$ \left[ \begin{array}{ccccc} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \text{k1} & 0.0630922 & 0.000256218 & 246.244 & \text{1.358897$\grave{ }$*${}^{\wedge}$-52} \\ \text{k2} & 0.0210129 & 0.0000735043 & 285.874 & \text{1.3335103$\grave{ }$*${}^{\wedge}$-54} \\ \text{a0} & 0.99959 & 0.00160667 & 622.151 & \text{4.545291$\grave{ }$*${}^{\wedge}$-65} \\ \text{b0} & 0.000233603 & 0.00172373 & 0.135522 & 0.893075 \\ \text{c0} & 0.00247035 & 0.00120983 & 2.04189 & 0.0497497 \\ \end{array} \right] $$