8
$\begingroup$

So, I have a ODE system, it is a complex biochemical kinetic mechanism with six species changing over time.

S'[t] == -k1 Eu[t] S[t] + k2 ES[t],
Eu'[t] == -k1 Eu[t] S[t] + k6 EP[t] + k2 ES[t],
ES'[t] == k1 Eu[t] S[t] - (k2 + k3) ES[t],
EP'[t] == k3 ES[t] - (k4 + k6) EP[t],
Ec'[t] == k4 EP[t],
P'[t] == k6 EP[t],

with the initial conditions:

S[0] == 100, Eu[0] == 0.5, ES[0] == 0, EP[0] == 0, Ec[0] == 0, 
P[0] == 0

I can solve the ODE system using NDSolve and manipulate it to "manually" fit some experimental data. Now, I have data for two species, and I want to numerically fit my ODE to those. I know three constants k1 (20),k2 (200) and k3 (0.03). I followed the approach described elsewhere, transforming my data in this way:

data = List[dataEP, dataEc];

transformedData = {ConstantArray[Range@Length[data], Length[time]] //Transpose, ConstantArray[time, Length[data]], data}~Flatten~{{2, 3}, {1}};

and then:

Sol = model[k3_?NumericQ, k4_?NumericQ, k6_?NumericQ, i_, te_] := ({EP[te], Ec[te]} /. First[NDSolve[ {
S'[t] == -k1 Eu[t] S[t] + k2 ES[t],
Eu'[t] == -k1 Eu[t] S[t] + k6 EP[t] + k2 ES[t],
ES'[t] == k1 Eu[t] S[t] - (k2 + k3) ES[t],
EP'[t] == k3 ES[t] - (k4 + k5 + k6) EP[t],
Ec'[t] == k4 EP[t],
Ed'[t] == k5 EP[t] ,
P'[t] == k6 EP[t], 
S[0] == 100, Eu[0] == 0.5, ES[0] == 0, EP[0] == 0, Ec[0] == 0, P[0] == 0}, {S, Eu, ES, EP, Ec, P}, {t, 0, 2000}, 
  Method -> Automatic, MaxSteps -> Infinity, 
  PrecisionGoal -> Infinity]])

and then using NonlinearModelFit as following:

fit = NonlinearModelFit[transformedData, {model[k3, k4, k6][i, t]},{k3, k4, k6}, {i, t}]

However, the fitting is really bad. I think the problem is that a) the fitting is not passing through the solver; b) maybe the fitting protocol is not identifying correctly EP and Ec. Another issue is that is not possible to get RSquared and other fitting options. Any help? I tried a lot of different setting and scripts, mostly following this forum. Thanks!!

Here an example of transformed data (i=1 is Ec and i=2 EP):

{{1, 0., 0.00001}, {1, 60.782, 0.01839}, {1, 121.43, 0.0273516}, {1, 
  182.062, 0.05744}, {1, 242.684, 0.066366}, {1, 303.31, 
  0.0834534}, {1, 363.983, 0.0966352}, {1, 424.626, 0.109041}, {1, 
  485.294, 0.124628}, {1, 545.964, 0.129099}, {1, 606.626, 
  0.133582}, {1, 667.293, 0.131262}, {1, 727.959, 0.142481}, {1, 
  788.619, 0.147817}, {1, 849.291, 0.145241}, {1, 909.936, 
  0.14883}, {1, 970.61, 0.154498}, {1, 1031.34, 0.151261}, {1, 
  1092.01, 0.155667}, {1, 1152.71, 0.15563}, {1, 1213.45, 
  0.148236}, {1, 1274.18, 0.15006}, {1, 1334.93, 0.161015}, {1, 
  1395.76, 0.158383}, {1, 1456.59, 0.167894}, {1, 1517.42, 
  0.165273}, {1, 1578.28, 0.170253}, {1, 1639.24, 0.166955}, {1, 
  1700.05, 0.160558}, {1, 1760.98, 0.161363}, {2, 0., 0.00001}, {2, 
  60.782, 0.233408}, {2, 121.43, 0.259436}, {2, 182.062, 
  0.224185}, {2, 242.684, 0.210032}, {2, 303.31, 0.175457}, {2, 
  363.983, 0.169942}, {2, 424.626, 0.163133}, {2, 485.294, 
  0.137899}, {2, 545.964, 0.116932}, {2, 606.626, 0.126436}, {2, 
  667.293, 0.108688}, {2, 727.959, 0.101772}, {2, 788.619, 
  0.0972984}, {2, 849.291, 0.0936195}, {2, 909.936, 0.0893072}, {2, 
  970.61, 0.0889732}, {2, 1031.34, 0.0737908}, {2, 1092.01, 
  0.0348883}, {2, 1152.71, 0.0796826}, {2, 1213.45, 0.0529935}, {2, 
  1274.18, 0.046321}, {2, 1334.93, 0.0341308}, {2, 1395.76, 
  0.0511362}, {2, 1456.59, 0.0326164}, {2, 1517.42, 0.0315381}, {2, 
  1578.28, 0.017776}, {2, 1639.24, 0.0254979}, {2, 1700.05, 
  0.00924619}, {2, 1760.98, 0.0225616}}

I also tried with ParametricNDSolveValue, in this way:

Sol = ParametricNDSolveValue[{
   S'[t] == -k1 Eu[t] S[t] + k2 ES[t],
  Eu'[t] == -k1 Eu[t] S[t] + k6 EP[t] + k2 ES[t],
   ES'[t] == k1 Eu[t] S[t] - (k2 + k3) ES[t],
   EP'[t] == k3 ES[t] - (k4 + k5 + k6) EP[t],
   Ec'[t] == k4 EP[t],
   P'[t] == k6 EP[t], S[0] == 100, Eu[0] == 0.5, 
   ES[0] == 0, EP[0] == 0, Ec[0] == 0, P[0] == 0}, {S, Eu,
    ES, EP, Ec, P}, {t, 0, 2000}, {k3,k4,k6}, MaxSteps -> Infinity, 
  PrecisionGoal -> Infinity]

followed by:

model[k3_,k4_, k6_][i_, t_] := 
  Through[Sol[k3,k4,k6][t], List][[i]] /;
   And @@ NumericQ /@ {k3, k4, k6,i, t};

Fitting again does not make any sense. Constraints also do not help. I tried with just k4>0, I left it overnight but NO fitting at all. I went through other questions, as I mentioned before, Manipulate my model gives reasonable "manual" fitting. Thanks!

$\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, and 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign. $\endgroup$ – bbgodfrey May 21 '15 at 23:56
  • 2
    $\begingroup$ You might want to look into ParametricNDSolveValue[]. $\endgroup$ – J. M. is away May 22 '15 at 0:21
  • 2
    $\begingroup$ I strongly second @Guess's suggestion. I think you need ParametricNDSolve[] or the *Value version for your system. Additionally, would you be able to post some sample data, or at least the reasonable values for the kinetic constants from your manual fit? It's hard to troubleshoot without trying to run some code. $\endgroup$ – MarcoB May 22 '15 at 0:23
  • 3
    $\begingroup$ Multiple potential duplicate: mathematica.stackexchange.com/q/21774; and mathematica.stackexchange.com/q/28461 and its many linked duplicates mathematica.stackexchange.com/q/34807, mathematica.stackexchange.com/q/56318, etc. Are you sure your question is not addressed by any of these? $\endgroup$ – Michael E2 May 22 '15 at 1:57
  • 1
    $\begingroup$ Hi! Brief update: I used ParametricNDSolveValue to solve the ODE system, then I generated some fake data putting some random error. From k4=0.0013, k5=0.0001, k6=0.06, my NonLinearModelFit procedure gave me k4=0.21, k5=0.181 and k6=3.46. Very far from my initial ones!! Hope you can give me a hint! :) $\endgroup$ – BlueOysterCult May 22 '15 at 16:18
4
$\begingroup$

This worked for me. I hope it helps.

I used ParametricNDSolveValue

k1 = 20; k2 = 200; k3 = 0.03;
tmax = 2000;
ode = {S'[t] == -k1 Eu[t] S[t] + k2 ES[t], 
   Eu'[t] == -k1 Eu[t] S[t] + k6 EP[t] + k2 ES[t], 
   ES'[t] == k1 Eu[t] S[t] - (k2 + k3) ES[t], 
   EP'[t] == k3 ES[t] - (k4 + k5 + k6) EP[t], Ec'[t] == k4 EP[t], 
   P'[t] == k6 EP[t],
   S[0] == 100, Eu[0] == 0.5, ES[0] == 0, EP[0] == 0, Ec[0] == 0, 
   P[0] == 0};
paramSOL = ParametricNDSolveValue[ode,
   {Ec, EP, S, Eu, ES, P}, {t, 0, tmax}, {k4, k5, k6}];

Then, define

model[k4_, k5_, k6_][i_, t_] := 
  Through[paramSOL[k4, k5, k6][t], List][[i]] /; And @@ NumericQ /@ {k4, k5, k6, i, t};

And using NonlinearModelFit...

fitted = NonlinearModelFit[data, model[k4, k5, k6][i, t],
        {{k4, 0.1}, {k5, 0.1}, {k6, 0.1}}, {i, t}] // Quiet;

fitted["RSquared"]
fitted["ParameterTable"]

RSquared = 0.990764

enter image description here

Plot of result:

dataEc = Take[data, 30][[All, 2 ;; 3]];
dataEP = Drop[data, 30][[All, 2 ;; 3]];
Show[
 ListPlot[{dataEc, dataEP}, PlotLegends -> {"Ec", "EP"},Frame -> True],
 Plot[ {fitted[1, t], fitted[2, t]}, {t, 0, tmax}] ]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks @Ivan, that's fantastic! I am running it right now, it is taking some time...it seems to me that the error was in model definition, right? Thanks again, hope it will work now! $\endgroup$ – BlueOysterCult May 25 '15 at 22:02
  • $\begingroup$ @BlueOysterCult That's weird. It takes like 2 seconds on my crappy machine. Note I made an edit, I forgot to put the value for tmax. One of the problems I saw was that you weren't specifying any starting points for the NonlinearModelFit to work properly. That might be why you were getting odd results. $\endgroup$ – Ivan May 26 '15 at 0:05
  • 1
    $\begingroup$ Hi @Ivan, it is working fine know, even with more complex model. And no, I tried with initial guesses too and it was the same. I think that the order of the different variable in ParametricNDSolveValue makes a difference. I mean, you specify ParametricNDSolveValue[ode, {EP, Ec, etc...}...] and your data need to be data=List[EP, Ec], in the same order!! I changed the order in your script, and it went odd again! :) Thanks!! – $\endgroup$ – BlueOysterCult May 26 '15 at 0:16
  • $\begingroup$ @BlueOysterCult That's exactly right! The order of the solutions of any solver is given by the order specified in the input. And in the model function the i variable is the index of the solution you are extracting, so 1 is Ec and 2 is EP. $\endgroup$ – Ivan May 26 '15 at 0:23

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.