# Fitting experimental data by using ParametricNDSolveValue and NonlinearModelFit

I am a newcomer to Mathematica. Basically I just want to fit the data (enzyme kinetic data) shown below to a system of odes' by using NonlinearModelFit. The data, an 101x2 array, contains time in the first column and substrate concentration in the second. The model contains 4 species (e[t], es[t], s[t] and p[t]) and species s[t] is the one I want to fit the data to. I get the following error:

"The function value {<<1>>} is not a list of real numbers with dimensions {101} at {k1,k2,k3} = {4.,3.,1.}". >>

Clicking on ">>" sends me to the basic NonlinearModelFit help page and I'm stuck there. The entire code is shown below. Thanks if anybody can help.

Francesco

(*some data*)
data={{0, 4.9112}, {20., 4.75011}, {40., 4.43818}, {60., 4.28744}, {80.,
3.97296}, {100., 3.86888}, {120., 3.69122}, {140., 3.59596}, {160.,
3.22247}, {180., 2.85438}, {200., 2.81939}, {220., 2.88236}, {240.,
2.49125}, {260., 2.55379}, {280., 2.33662}, {300., 2.34136}, {320.,
1.88169}, {340., 1.9444}, {360., 1.73578}, {380., 2.04545}, {400.,
1.74068}, {420., 1.70471}, {440., 1.37455}, {460., 1.35169}, {480.,
1.29391}, {500., 1.35778}, {520., 1.1509}, {540., 1.18335}, {560.,
0.846087}, {580., 0.957338}, {600., 0.855021}, {620.,
0.727364}, {640., 0.886429}, {660., 0.817111}, {680.,
0.748117}, {700., 0.569694}, {720., 0.77641}, {740.,
0.661459}, {760., 0.561378}, {780., 0.56037}, {800.,
0.500522}, {820., 0.322087}, {840., 0.44058}, {860.,
0.359604}, {880., 0.31989}, {900., 0.278633}, {920.,
0.318697}, {940., 0.150813}, {960., 0.427698}, {980.,
0.364589}, {1000., 0.292937}, {1020., 0.27481}, {1040.,
0.182754}, {1060., 0.349605}, {1080., 0.220416}, {1100.,
0.149073}, {1120., 0.343196}, {1140., 0.173815}, {1160.,
0.126286}, {1180., 0.145337}, {1200., 0.0800335}, {1220.,
0.043485}, {1240., 0.399296}, {1260., 0.303941}, {1280.,
0.161308}, {1300., -0.00255049}, {1320., 0.0296389}, {1340.,
0.0919508}, {1360.,
0.182537}, {1380., -0.0356638}, {1400., -0.140977}, {1420.,
-0.0581143}, {1440., 0.115227}, {1460., 0.116371}, {1480.,
0.118025}, {1500., 0.0556984}, {1520., 0.0831993}, {1540.,
0.0135393}, {1560., 0.143889}, {1580., -0.0817538}, {1600.,
0.0968327}, {1620., -0.0364522}, {1640., 0.0121839}, {1660.,
0.0983604}, {1680., 0.144547}, {1700., -0.0734307}, {1720.,
0.162225}, {1740., 0.100122}, {1760., 0.0253859}, {1780.,
0.0108251}, {1800., 0.00686486}, {1820., -0.00330938}, {1840.,
0.0277739}, {1860., 0.0291533}, {1880., 0.105267}, {1900.,
0.174073}, {1920., 0.0668537}, {1940., -0.00195318}, {1960.,
0.080458}, {1980., 0.0352437}, {2000., -0.0870161}};
(* Dimensions[data] {101, 2} *)

tmax = Max[data[[All, 1]]];

(* the model, k1, k2, k3 are the parameters *)

model = ParametricNDSolveValue[{e'[t] == (k2 + k3) es[t]
- k1 e[t] s[t], es'[t] == -e'[t], s'[t] == k2 es[t] - k1 e[t] s[t],
p'[t] == k3 es[t], e[0] == 0.001, s[0] == 5, es[0] == 0,
p[0] == 0}, s[t], {t, 0, tmax}, {k1, k2, k3}]

fit = NonlinearModelFit[data,model[k1, k2, k3][t], {{k1, 4.}, {k2, 3}, {k3, 1}}, t];
plotfit = Plot[model[k1, k2, k3] /. fit, {t, 0, tmax}];
plotdata = ListPlot[data, PlotStyle -> PointSize[0.02]];
Show[plotfit, plotdata]

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Commented Aug 14, 2016 at 17:16
• I am a newcomer to Mathematica also (use ver 9.0). I have a lot of mistakes when I try to reproduce this exemple. data = {{0, 4.9112}, {20., 4.75011}, {40., 4.43818}, {60., 4.28744}, {80., 3.97296}, {100., 3.86888}, {120., 3.69122}, {140., 3.59596}, {160., 3.22247}, {180., 2.85438}, {200., 2.81939}, {220., 2.88236}, {240., 2.49125}, {260., 2.55379}, {280., 2.33662}, {300., 2.34136}, {320., 1.88169}, {340., 1.9444}, {360., 1.73578}, {380., 2.04545}, {400., 1.74068}, {420., 1.70471}, {440., 1.37455}, {460., 1.35169}, {480., 1.29391}, {500., 1.35778}, {520., 1.1509}, {540., 1.18335}, {560., 0.846 Commented Jun 20, 2019 at 14:15

A minor change to model eliminates the error cited in the question. (Replace s[t] by s in the final line of the code immediately below.)

model = ParametricNDSolveValue[{e'[t] == (k2 + k3) es[t] - k1 e[t] s[t],
es'[t] == -e'[t], s'[t] == k2 es[t] - k1 e[t] s[t],
p'[t] == k3 es[t], e[0] == 0.001, s[0] == 5, es[0] == 0,
p[0] == 0}, s, {t, 0, tmax}, {k1, k2, k3}];


Then, using NonlinearModelFit with better initial guesses yields a solution in seconds.

fit = NonlinearModelFit[data, model[k1, k2, k3][t], {{k1, 6}, {k2, 3}, {k3, 3}}, t];


Finally, a corrected use of fit provides the desired result.

plotfit = Plot[fit[t], {t, 0, tmax}];
plotdata = ListPlot[data, PlotStyle -> PointSize[0.01]];
Show[plotfit, plotdata]


For completeness, the fitted values of {k1, k2, k3} are

fit["BestFitParameters"]
{k1 -> 637.856, k2 -> 16540.4, k3 -> 80.546}


So, the initial guesses used above were not so good after all. Evidently, they simply helped NonlinearModelFit to obtain a better start. Other initial guesses, such as {{k1, 12}, {k2, 6}, {k3, 6}}, lead to different fitted values, here {k1 -> 434.784, k2 -> 11248.8, k3 -> 80.5462}, but with a final result indistinguishable to the eye from the plot above, even when displayed on a log scale.

• The variability in parameter estimates depending on starting values is likely due to the model being overparameterized as fit["CorrelationMatrix"] // MatrixForm results in $\left( \begin{array}{ccc} 1. & 1. & -0.765387 \\ 1. & 1. & -0.765385 \\ -0.765387 & -0.765385 & 1. \\ \end{array} \right)$. The estimators of k1 and k2 appear to be perfectly correlated with each other.
– JimB
Commented Jun 20, 2019 at 15:34
• Thanks, @JimB. This is a very useful observation. Commented Jun 20, 2019 at 15:59
• Actually, setting k1=700 and k3=80 and fitting for k2 results in a slightly smaller mean square error (fit["EstimatedVariance"]) and smaller AIC value. That's not to say that the model is "bad" but simply that there isn't the right kind of data to estimate all three parameters. Another indication of trouble is that the fit["ParameterTable"] shows huge standard errors for the full model.
– JimB
Commented Jun 20, 2019 at 16:09

In case no better answer presents itself, this is a manual fit:

Manipulate[f[ti_] := Evaluate[model[k1, k2, k3] /. t -> ti];
Show[lp,
Plot[Evaluate[model[k1, k2, k3]], {t, 0, tmax},
PlotStyle -> {Thick, Red}],
Graphics[Text[
Sqrt[Total[
Table[(f[data[[i, 1]]] - data[[i, 2]])^2, {i,
Length[data]}]]], {1500, 4}]]], {k1, 0.01, 40, 0.01}, {k2,
0.01, 80, 0.01}, {k3, 0.01, 40, 0.01}]
`

The goal is to minimise the numerical value displayed in the top right corner, this is the square root of the total of the square of the distances between the fit and the data points.