9
$\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!

Improve this question
$\endgroup$
10
  • 2
    $\begingroup$ You might want to look into ParametricNDSolveValue[]. $\endgroup$
    – J. M.'s eventual burnout ♦
    May 22, 2015 at 0:21
  • 2
    $\begingroup$ I strongly second @J.M.'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, 2015 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, 2015 at 1:57
  • $\begingroup$ As @MichaelE2 also pointed out, Oleksandr's answer to a past question seems extremely relevant. You obviously are aware of it, since you seem to be using the same code structure he proposed. Could you show us perhaps some annotated fitting results, or otherwise give us an idea of what you are not satisfied with? Additionally, you could try and generate a "fake" data set to feed to your fitting procedure, to see if you can recover the known parameters you used to generate the dataset. $\endgroup$
    – MarcoB
    May 22, 2015 at 3:06
  • 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, 2015 at 16:18

1 Answer 1

Reset to default
5
$\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

Improve this answer
$\endgroup$
4
  • $\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, 2015 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, 2015 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, 2015 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, 2015 at 0:23

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.