How to speed up NonlinearModelFit applied on a NDSolve solution?

Question

I'm trying to perform fittings of a model defined through a system of ODEs to data consisting of time courses of a measured quantity. As described in a previous question, the system of ODEs is as follows:

X'[t]:= m[t].X[t]
X[t_] := {h[t], r[t], rh1[t], rh2[t], rh3[t]}
m[t_] := {{-k1*r[t], 0, ki1, 0, 0}, {0, -k1*h[t], ki1, 0, 0}, {0,   k1*h[t], -(ki1 + k2), ki2, 0}, {0, 0, k2, -(ki2 + k3), ki3}

and the actual quantity I need to evaluate is:

F:= aF.X[t]
aF:={0,aR,aRH1,aRH2,aRH3}

After the incorporation of some suggestions given in this forum, I end up with the following procedure to fit the model:

Model to be fitted:

modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt] =
Module[{X, m, aF, sol},
X[t_] := {h[t], r[t], rh1[t], rh2[t], rh3[t]};
m[t_] := {{-k1*r[t], 0, ki1, 0, 0}, {0, -k1*h[t], ki1, 0, 0}, {0,     k1*h[t], -(ki1 + k2), ki2, 0}, {0, 0, k2, -(ki2 + k3), ki3}, {0, 0, 0, k3, -ki3}};
aF = {0, aR, aRH1, aRH2, aRH3};
sol = NDSolve[{X'[t] == m[t].X[t], X[0] == {Ht, Rt, 0, 0, 0}}, X[t], {t, 0, 500}];
Function[{tu}, Evaluate[aF.Flatten[X[t] /. sol] - aR*Rt] /. t :> tu]])

Data:

data = Table[{t, .22 (1 - E^(-7.2 t)) + 0.10 (1 - E^(-0.084 t)) + 0.15 (1 - E^(-0.027 t))}, {t, 0, 500, 0.1}]

Fitting routine (leaving some parameters fixed):

k1 = 0.0012;
k2 = 0.43;
ki2 = 0.0055;
k3 = 0.01;
ki3 = 0.0096;
aR = 0.0034;
Ht = 800;
Rt = 62; (*The last three quantities are not fitting parameters but fixed ones, whose values are known. *)
fit = NonlinearModelFit[data, {modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt][t], {2 <= ki1 <= 20, 0.001 <= aRH1 <= 0.1, 0.001 <= aRH2 <= 0.1, 0.001 <= aRH3 <= 0.1}}, {{ki1, 5}, {aRH1, .032}, {aRH2, .01}, {aRH3, .012}}, t, Method -> {NMinimize, Method -> "NelderMead"}]]

Here is the problem. The model works and the fitting code too, but it is too slow. For instance, it takes about 110 s (in Linux Mint MATE 1.8.1 on a dual core AMD A4 machine, with 3 GiB RAM) for a 2 iterations run and a data set of 5000 points. Given that I pretend to use the program to fit series of ~ 8 curves of ~ 14000 points, this performance is obviously too slow (it is much faster in other programs such us COPASI). Please, let me know if you have any idea on how to speed up the code.

Edit 1:

To asses the quality of the fitting I add here a new data set obtained from the model with a given set of parameter values:

newdata = Table[{t, Evaluate[With[{k1 = 0.0012, ki1 = 5, k2 = 0.43, ki2 = 0.0055, k3 = 0.01, ki3 = 0.0096, aR = 0.0034, aRH1 = .032, aRH2 = .01,    aRH3 = .012},modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt]]][t]}, {t, 0, 500}]

Answer 1

After trying many ways to fit the -numerically solved- model described above, with lots of help from fellows, it turned out that -for this kind of model at least- fitting is much faster if ParametricNDSolve is used, instead of NDSolve, to solve the ODEs of the model. This is just what 'Guess who it is' slightly suggested in a comment to my original question.

I just give an example based on the original model comparing the two procedures.

Model to be fitted:

Basically, the model is:

 X'[t]=m[t].X[t],
 X[t_]={h[t],r[t],rh1[t],rh2[t],rh3[t]} (* vector of concentrations of species of the modeled chemical reaction *)
 m[t_] := {{-k1*r[t], 0, ki1, 0, 0}, {0, -k1*h[t], ki1, 0, 0}, {0, 
k1*h[t], -(ki1 + k2), ki2, 0}, {0, 0, k2, -(ki2 + k3), ki3}, {0, 
0, 0, k3, -ki3}} (* the so called 'kinetic matrix' of the system *)

And the quantity measured (to be fitted to experimental values) is:

F = {0,aR,aRH1,aRH2,aRH3}.X[t]

Implementation of the model and fitting:

The model is formulated below using NDSolve or ParametricNDSolve.

Clear[F, Xsol, G, Ysol];
X[t_] := {h[t], r[t], rh1[t], rh2[t], rh3[t]};
m[t_] := {{-k1*r[t], 0, ki1, 0, 0}, {0, -k1*h[t], ki1, 0, 0}, {0, k1*h[t], -(ki1 + k2), ki2, 0}, {0, 0, k2, -(ki2 + k3), ki3}, {0, 0, 0, k3, -ki3}};
 (* Fixed parameters for this example *)
aR = 0; k2 = 0.347; ki2 = 0.0235; k3 = 0.0507; ki3 = 0.00645; ht = 800; rt = 20;
(* Solving of ODEs *)
Xsol = ParametricNDSolve[{X'[t] == m[t].X[t], X[0] == {ht, rt, 0, 0, 0}}, X[t], {t, 0, 500, 0.01}, {{k1, 0.0001, 0.1}, {ki1, 0.1, 100}, {k2, 0.0001, 10}, {ki2, 0.00001, 10}, {k3, 0.00001, 10}, {ki3, 0.00001, 10}, {ht, 1,     10000}, {rt, 1, 1000}}]
Ysol := NDSolve[{X'[t] == m[t].X[t], X[0] == {ht, rt, 0, 0, 0}}, X[t], {t, 0, 500, 0.01}]

(* Quantity to be fitted *)
F[k1_?NumberQ, ki1_?NumberQ, k2_?NumberQ, ki2_?NumberQ, k3_?NumberQ,  ki3_?NumberQ, aR_?NumberQ, aRH1_?NumberQ, aRH2_?NumberQ, aRH3_?NumberQ, t_?NumberQ, rt_?NumberQ] := (F[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, ht, rt] = 
Function[{tu}, {0, aR, aRH1, aRH2, aRH3}.Through[(X[t] /. Xsol)[k1, ki1, k2, ki2, k3, ki3, ht, rt]] /. t :> tu])

G[k1_?NumberQ, ki1_?NumberQ, k2_?NumberQ, ki2_?NumberQ, k3_?NumberQ,   ki3_?NumberQ, aR_?NumberQ, aRH1_?NumberQ, aRH2_?NumberQ, aRH3_?NumberQ, ht_?NumberQ, rt_?NumberQ] := (G[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, ht, rt] = 
Function[{tu}, Evaluate[{0, aR, aRH1, aRH2, aRH3}.Flatten[X[t] /. Ysol]] /. t :> tu])

The data for the example was:

data = Table[{t, 0.22 (1 - E^(-7.25 t)) + 0.10 (1 - E^(-0.084 t)) + 0.15 (1 - E^-0.027)}, {t, 0, 500}]

The fitting was run as follows:

 Clear[FitF, FitG]
 Timing[FitF = NonlinearModelFit[data500, F[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, ht, rt][t], {{k1, 0.001}, {ki1, 7}, , {aRH1, .032}, {aRH2, .01}, {aRH3, .012}}, t]]
 Timing[FitG = NonlinearModelFit[data500, G[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, ht, rt][t], {{k1, 0.001}, {ki1, 7}, , {aRH1, .032}, {aRH2, .01}, {aRH3, .012}}, t]]

The output gave the same set of parameter values (along with some NDSolve's warnings):

 {k1 -> 0.000991454, ki1 -> 7.00047, Null -> 1., aRH1 -> 0.117238, aRH2 -> 0.0278981, aRH3 -> 0.0189699}

But very different timings: ~ 23 s for the 'ParametricNDSolve'-model (F), and ~ 1128 s for the 'NDSolve'-model. That is, a 50 fold difference. (I must mention that such difference in timing was also observed in other formulations I tried in which no warnings or errors were prompted)

Using the ParametricNDSolve strategy, the simultaneous fitting of 9 parameters of the model to 9 curves of about 14000 points each, converged very well and it took about 3600 s (1 hour). That's really fine. Now we are talking.

Maybe the convenience of using ParametricNDSolve is obvious for many -it is not for a beginner as me-, but I think it deserves some discussion, at least for future reference for a new beginner.