10
$\begingroup$

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]
Improve this question
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Aug 14, 2016 at 17:16
  • $\begingroup$ 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 $\endgroup$
    – Irrrr
    Jun 20, 2019 at 14:15

2 Answers 2

Reset to default
12
$\begingroup$

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]

enter image description here

Addendum

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.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Jun 20, 2019 at 15:34
  • $\begingroup$ Thanks, @JimB. This is a very useful observation. $\endgroup$
    – bbgodfrey
    Jun 20, 2019 at 15:59
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Jun 20, 2019 at 16:09
5
$\begingroup$

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

enter image description here

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.

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.