8
$\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]
$\endgroup$
  • $\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 '16 at 17:16
10
$\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.

$\endgroup$
4
$\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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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