# Substituting ParametricNDSolve for NDSolve

This is my first post so let me know if I need to format things differently. I have been using NDSolve for most problems which optimize parameters for ODEs. I believe that ParametricNDSolve would be more robust, but cannot make it work with some of our problems. Here is an example problem:

Et = 0.05;
S0 = 50;
k1 = 60.;
k2 = 600.;
k3 = 10;
k4 = 60.;
k5init = 1000.;
k6init = .05;

dataP = {{0, 0.001, 20.001, 0.837}, {0, 6, 26, 0.810}, {0, 12, 32,
0.825}, {0, 18, 38, 0.825}, {0, 24, 44, 0.880}, {0, 30, 50,
0.822}, {1, 0.001, 20.001, 0.814}, {1, 6, 26, 0.808}, {1, 12, 32,
0.804}, {1, 18, 38, 0.80124}, {1, 24, 44, 0.7694}, {1, 30, 50,
0.776}, {3, 0.001, 20.001, 0.866}, {3, 6, 26, 0.791}, {3, 12, 32,
0.752}, {3, 18, 38, 0.755}, {3, 24, 44, 0.742}, {3, 30, 50,
0.675}, {10, 0.001, 20.001, 0.816}, {10, 6, 26, 0.752}, {10, 12,
32, 0.719}, {10, 18, 38, 0.653}, {10, 24, 44, 0.626}, {10, 30, 50,
0.572}, {30, 0.001, 20.001, 0.778}, {30, 6, 26, 0.661}, {30, 12,
32, 0.612}, {30, 18, 38, 0.560}, {30, 24, 44, 0.513}, {30, 30, 50,
0.444}, {100, 0.001, 20.001, 0.700}, {100, 6, 26, 0.591}, {100,
12, 32, 0.514}, {100, 18, 38, 0.469}, {100, 24, 44, 0.403}, {100,
30, 50, 0.344}};


Solving for Parameters (k5 and k6), we use data P[t] (I needed to use a different variable for t (te i.e. P[te]) to make it work - I do not know why.

model[k5_?NumericQ, k6_?NumericQ, I0_?NumericQ, tS_, te_] :=
P[te] /. NDSolve[{
Eu'[t] == -k1 Eu[t] S[t] + k2 ES[t] + k3 ES[t] - k4 Eu[t] Iu[t] +
k5 EI[t],
ES'[t] == k1 Eu[t] S[t] - (k2 + k3) ES[t],
EI'[t] == -(k5 + k6) EI[t] + k4 Eu[t] Iu[t],
E2'[t] == k6 EI[t],
S'[t] == -k1 Eu[t] S[t] + k2 ES[t],
P'[t] == k3 ES[t],
Iu'[t] == -k4 Eu[t] Iu[t] + k5 EI[t],
S[0] == 0,
Eu[0] == Et,
ES[0] == 0,
EI[0] == 0,
P[0] == 0,
Iu[0] == I0,
E2[0] == 0,
WhenEvent[
t == tS, {S[t] -> S0, Eu[t] -> Eu[t]/10, ES[t] -> ES[t]/10,
EI[t] -> EI[t]/10, E2[t] -> E2[t]/10, P[t] -> P[t]/10,
Iu[t] -> Iu[t]/10}]},
{Eu, ES, EI, E2, S, P, Iu}, {t, 0, 60}, MaxSteps -> 100000,
PrecisionGoal -> \[Infinity]];

fit = NonlinearModelFit[
dataP, {model[k5, k6, I0, tS, te]}, {{k5, k5init}, {k6,
k6init}}, {I0, tS, te}, Weights -> (1/#4 &),
PrecisionGoal -> \[Infinity], MaxIterations -> 10000];

In[4549]:= fit["BestFitParameters"]

Out[4549]= {k5 -> 609.524, k6 -> 0.0247421}


When I try to use ParametricNDSolve I have errors:

model[k5_?NumericQ, k6_?NumericQ][I0_?NumericQ, tS_, te_] :=
P[te] /. ParametricNDSolve[{
Eu'[t] == -k1 Eu[t] S[t] + k2 ES[t] + k3 ES[t] - k4 Eu[t] Iu[t] +
k5 EI[t],
ES'[t] == k1 Eu[t] S[t] - (k2 + k3) ES[t],
EI'[t] == -(k5 + k6) EI[t] + k4 Eu[t] Iu[t],
E2'[t] == k6 EI[t],
S'[t] == -k1 Eu[t] S[t] + k2 ES[t],
P'[t] == k3 ES[t],
Iu'[t] == -k4 Eu[t] Iu[t] + k5 EI[t],
S[0] == 0,
Eu[0] == Et,
ES[0] == 0,
EI[0] == 0,
P[0] == 0,
Iu[0] == I0,
E2[0] == 0,
WhenEvent[
t == tS, {S[t] -> S0, Eu[t] -> Eu[t]/10, ES[t] -> ES[t]/10,
EI[t] -> EI[t]/10, E2[t] -> E2[t]/10, P[t] -> P[t]/10,
Iu[t] -> Iu[t]/10}]},
{Eu, ES, EI, E2, S, P, Iu}, {t, 0, 60}, {k5, k6},
MaxSteps -> 100000, PrecisionGoal -> \[Infinity]];

fit = NonlinearModelFit[
dataP, {model[k5, k6][I0, tS, te]}, {{k5, k5init}, {k6,
k6init}}, {I0, tS, te}, Weights -> (1/#4 &),
PrecisionGoal -> \[Infinity], MaxIterations -> 10000];


Error:Initial condition I0 is not a number or a rectangular array of numbers. >>

I am obviously missing something. Thanks - Ken

Sometimes you need to be careful with the semantics. I0 and tS are parameters and not variables for the ODEs:

pe = ParametricNDSolve[{
Eu'[t] == -k1 Eu[t] S[t] + k2 ES[t] + k3 ES[t] - k4 Eu[t] Iu[t] + k5 EI[t],
ES'[t] ==  k1 Eu[t] S[t] - (k2 + k3) ES[t],
EI'[t] == -(k5 + k6) EI[t] + k4 Eu[t] Iu[t],
E2'[t] ==  k6 EI[t],
S'[t] == -k1 Eu[t] S[t] + k2 ES[t],
P'[t] ==  k3 ES[t],
Iu'[t] == -k4 Eu[t] Iu[t] + k5 EI[t],
(* initial conditions *)
S[0]== 0, Eu[0]== Et, ES[0]== 0, EI[0]== 0, P[0]== 0, Iu[0]== I0, E2[0]== 0,
WhenEvent[
t == tS, {S[t] -> S0, Eu[t] -> Eu[t]/10, ES[t] -> ES[t]/10,
EI[t] -> EI[t]/10, E2[t] -> E2[t]/10,
P[t]  -> P[t]/10,  Iu[t] -> Iu[t]/10}]},
(* functions *)
{Eu, ES, EI, E2, S, P, Iu},
(* vars & intervals *)
{t, 0, 60},
(* paramters (these don't take intervals since we are not integrating over them) *)
{k5, k6, I0, tS},
MaxSteps -> 100000,  PrecisionGoal -> ∞];

h[k5_?NumericQ, k6_?NumericQ, I0_?NumericQ, tS_?NumericQ,te_?NumericQ] :=
(P /. pe)[k5, k6, I0, tS][te]

fit = NonlinearModelFit[dataP, h[k5, k6, I0, tS, te],
{{k5, k5init}, {k6, k6init}}, {I0, tS, te}];

fit["BestFitParameters"]
(* {k5 -> 610.104, k6 -> 0.0247282} *)

• Fantastic! I still do not understand the inclusion of I0 and tS as parameters. I0 is the initial concentration of I and tS is the time that S is added to the experiment. I would consider these to be variables, with k5 and k6 parameters that define model behavior and need to be optimized. Does including I0 and tS in the parameters in ParametricNDSolve negate any advantage over using NDSolve? It appears that I0 is treated differently than tS since I0_?NumericQ is necessary when using NDSolve but tS_ is sufficient. I0 is in the actual differential equations and tS is part of WhenEvent. Thanks! Jun 16, 2015 at 21:50
• @flyfishken Well, the rule for distinguishing parameters from variables in NDSolve[ ]is very easy. If you're taking derivatives respect to it, it's a variable. Otherwise it's a parameter Jun 17, 2015 at 1:53
• Thank you. I will find out if ParametricNDSolve will make some of our more difficult problems easier. Jun 17, 2015 at 12:54