ODE fitting to dataset

Question

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!

Answer 1

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

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