Mathematica unable to interpret equations in DSolve

Question

I'm trying to solve a system of seven ODEs with seven conditions.

Since the following approach using DSolve is yielding the same output as the input, Mathematica can't find a solution.

DSolve[
 {EA'[t] == k1 e[t] a[t] + k22 EAB[t] - k2 b[t] EA[t] - k11 EA[t],
  EB'[t] == k4 e[t] b[t] - k44  EB[t] + k33 EAB[t] - k3 a[t] eb[t],
  EAB'[t] == 
   k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] + k33 EAB[t],
  p'[t] == k5 EAB[t],
  e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t] ,
  a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t],
  b'[t] == k44 EB [t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t] ,

  EA[0] == 0,
  EB[0] == 0,
  EAB[0] == 0,
  a[0] == 1,
  e[0] == 1/10,
  b[0] == 100,
  p[0] == 0
  }, {EAB[t], EB[t], EA[t], a[t], e[t], b[t], p[t]}, t]

I then attempted this problem using ParametricNDSolve, and then subsequently fitting to the parameters using NonlinearModelFit and FindFit. With FindFit, I get the error:

"The function value {<<1>>} is not a list of real numbers with \
dimensions "

With NonlinearModelFit, I get

NDSolve Parametric Function is not a list of real \ numbers with dimensions

What I'm doing is here:

pSoln = ParametricNDSolve[
  {EA'[t] == 
    k1* e[t]* a[t] - k11* EA[t] + k22 *EAB[t] - k2* b[t] *EA[t] ,
   EB'[t] == k4 e[t] b[t] - k44  EB[t] + k33 EAB[t] - k3 a[t] EB[t],
   EAB'[t] == 
    k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] + 
     k33 EAB[t],
   p'[t] == k5 EAB[t],
   e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t] ,
   a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t],
   b'[t] == k44 EB [t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t] ,

   EA[0] == 0,
   EB[0] == 0,
   EAB[0] == 0,
   a[0] == 1,
   e[0] == .1,
   b[0] == 100,
   p[0] == 0}, 
  p, {t, 0, 1000}, {k1, k11, k2, k22, k3, k33, k4, k44, k5}]

fit = FindFit[data, 
   pSoln[k1, k11, k2, k22, k3, k33, k4, k44, k5][
    t], {{k1, .15}, {k11, 1}, {k2}, k22, k3, 
    k33, {k4, .75}, {k44, 0}, {k5, .75}}, t];

For the NonlinearModelFit error, I just changed FindFit to NonlinearModelFit.

UPDATE: I think that ParametricNDSolve is the way to go based on the other stack posts I've viewed; however, I'm unable to plot the output of ParametricNDSolve/ParametricNDSolveValue. I had made the following function where I could tune the values of the constants:

    pFunction[t_?NumericQ, k1_?NumericQ, k11_?NumericQ, k2_?NumericQ, 
   k22_?NumericQ, k3_?NumericQ, k33_?NumericQ, k4_?NumericQ, 
   k44_?NumericQ, k5_?NumericQ]  := 
  Evaluate[pSoln][k1, k11, k2, k22, k3, k33, k4, k44, k5][t];

With that, I then wanted to do NMinimize using the values from the manual sliders to find the actual fits:

{rmse, params} = NMinimize[{Sqrt[
    Mean[(pFunction[data[[1]][[All, 1]], k1, k11, k2, k22, k3, 
         k33, k4, k44, k5] - data[[1]][[All, 2]])^2]],
   0 < k1 < 0.2,
   0 < k11 < 2,
   0 < k2 < 0.05,
   0 < k22 < 1,
   0 < k3 < 0.2,
   0 < k33 < 0.2,
   0 < k4 < 2,
   0 < k44 < 3,
   0 < k5 < 1},
  {k1, k11, k2, k22, k3, k33, k4, k44, k5}
  ]

Whether or not I used ?NumericQ or _Real for the inputs in pFunction, NMinimize crashes with

"The function value .... is \
not a number at {k1,k11,k2,k22,k3,k33,k4,k44,k5}"

///

For what it's worth, a subset of the data that I'm working with is here:

data={
    {0.,0.00649550468169123},
    {4.09905486238222,-0.021118728609071737},
    {8.22806574924588,0.002077434829573818},
    {12.3587388435858,0.14235961440477035},
    {16.4890594300343,-0.017660745493849748},
    {20.6200130450588,0.11140389865988036},
    {24.7526023395476,0.06847436631515905},
    {28.8848119119313,0.1442917327119322},
    {33.0176280094345,0.1861025713545064},
    {37.1520539765966,0.1443058843144346},
    {41.2860737678823,0.2289055285168752},
    {45.420674325329,0.19174285727745513},
    {49.5568596712869,0.22341139281802608},
    {53.6926131564865,0.25854671096480975},
    {57.8289224014593,0.13868939793044502},
    {61.96679209018,0.24859665212276605},
    {66.1042050078925,0.35602799421670167},
    {70.24214943697,0.17457668699584183},
    {74.3816307065195,0.20476112596323676},
    {78.5206310725897,0.24270049895368628},
    {82.6601394624533,0.27548916188262884},
    {86.8011618335922,0.32913163018521907},
    {90.9416799472357,0.22361905520264616},
    {95.082683358369,0.2679259771614342},
    {99.2251786358719,0.3272748985310891},
    {103.367147078656,0.3395447955662898},
    {107.509578852036,0.3195026272328241},
    {111.653481119137,0.30730789723336444},
    {115.796834747233,0.5040178412127853},
    {119.940630494433,0.4274042258121084},
    {124.08587610083,0.29313682148998943},
    {128.230552030975,0.3253599695876959},
    {132.375649618136,0.32145305280770675},
    {136.522177162011,0.49417265921299086},
    {140.668114751638,0.4589792754509556},
    {144.814454277751,0.5139683506801757},
    {148.962204582262,0.45910962815361306},
    {153.109345404262,0.5113903786348323},
    {157.256869174255,0.5571819491842587},
    {161.405785258981,0.5745371402536695},
    {165.554073071574,0.4602900338438135},
    {169.702725564653,0.5993781357941224},
    {173.85275261245,0.5814410492884232},
    {178.00213332463,0.5483466054674879},
    {182.151861158352,0.5050097702205707},
    {186.302946478108,0.5611301955777416},
    {190.453368111136,0.5534989203499485},
    {194.604120001686,0.6772346363574281},
    {198.75621298741,0.5516487233158627},
    {202.907625632799,0.467445102907578},
    {207.059352351908,0.6021862725531534},
    {211.212404438627,0.50787909253444},
    {215.36476021308,0.6741722001846713},
    {219.517414542041,0.6841842864704483}
};

PS. The question was edited.

Answer 1

Update

I was a little too clever when I extracted the dependent variables from the equations and made an error to the variable to fit. The 'p' variable to fit occurred in equation 4, but was extracted as the 7th variable. I have corrected the error so that the fit corresponds to the "p" variable.

As I commented in my answer, parameter fitting is frequent topic on MSE. I have been experimenting with StringTemplate as a possible way generalize the fitting to address this common pattern. Parameter fitting is a difficult subject and will depend on your data quality, your model, and your intial guesses. In your particular case, we are trying to fit 9 parameters to noisy data that could be well characterized by a straight line.

Approach

• Use ParametricNDSolveValue to create the model.
• Use StringTemplates to handle lists of parameters and variables.
• Generate a Manipulate slider model to debug model and understand the effects of parameter changes.
• Transfer initial guesses from manipulate to perform a fit.

Implementation

I commented the code so I hope it self explanatory. First assign the constants and prep the data.

(* Data for fitting Column 2 is "p" *)
(* data not reproduced here *)
(* Max Time *)
tmax = 1000;
(* Equations used to extract parameters and variables *)
test = {EA'[t] == 
    k1*e[t]*a[t] - k11*EA[t] + k22*EAB[t] - k2*b[t]*EA[t], 
   EB'[t] == k4 e[t] b[t] - k44 EB[t] + k33 EAB[t] - k3 a[t] EB[t], 
   EAB'[t] == 
    k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] + 
     k33 EAB[t], p'[t] == k5 EAB[t], 
   e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t], 
   a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t], 
   b'[t] == k44 EB[t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t]};

Here is the code to solve the Diffeqs and create the Manipulate slider to improve initial guesses.

(* Generate System of Differential Equations *)
e1 = EA'[t] == 
  k1*e[t]*a[t] - k11*EA[t] + k22*EAB[t] - k2*b[t]*EA[t]; e2 = 
 EB'[t] == k4 e[t] b[t] - k44 EB[t] + k33 EAB[t] - k3 a[t] EB[t]; e3 =
  EAB'[t] == 
  k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] + 
   k33 EAB[t]; e4 = p'[t] == k5 EAB[t];
e5 = e'[t] == 
  k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t]; e6 = 
 a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t]; e7 = 
 b'[t] == k44 EB[t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t];
ics = {EA[0] == 0, EB[0] == 0, EAB[0] == 0, a[0] == 1, e[0] == 1/10, 
   b[0] == 100, p[0] == 0};
eqns = {e1, e2, e3, e4, e5, e6, e7}~Join~ics;
(* Extract Variables from test *)
vbles = Cases[test, f_'[t_] -> f, All] // Union;
(*Parameters with target and desired ranges*)
mat = ConstantArray[0, {9, 4}];
(* Extract parameters from test *)
mat[[All, 1]] = Cases[test, _Symbol, Infinity] // Union // Most;
mat = {
   {k1, 0.15, 0.01, 0.4},
   {k11, 1, 0.01, 2},
   {k2, 0.05, 0, 0.1},
   {k22, 0.05, 0, 0.1},
   {k3, 0.05, 0, 0.1},
   {k33, 0.05, 0, 0.5},
   {k4, 0.75, 0.01, 1},
   {k44, 0.5, 0, 1},
   {k5, 0.75, 0, 2}
   };
(* Reduced parameter matrix if we want to fix variables *)
rmat = mat;
(* Build Manipulate sliders *)
sfun =  StringRiffle[(StringTemplate[
         "{{`1`,`2`},`3`,`4`,Appearance\[Rule]\"Labeled\"}"] @@ #) & \
/@ #, ","] &;
sliders = sfun[rmat];
(* Extract Parameters from mat *)
parms = mat[[All, 1]];
(* Create String Representations of parms *)
sparms = StringRiffle[ToString[#] & /@ parms, ","];
(* Create patterns and string reps of parameters *)
pats = Pattern @@@ (#*_ & /@ parms);
spats = StringRiffle[ToString[#] & /@ pats, ","];
(* List Plot of the data *)
lp = ListPlot[data, PlotLegends -> {"data-p"}];
(* ParametricNDSolveValue *)
pfun = ParametricNDSolveValue[eqns, vbles, {t, 0, tmax}, parms];
(*Create an appropriate model function to fit*)
modelstring = "{#[[1]],#[[2]],#[[3]],#[[4]],#[[5]],#[[6]],#[[7]]}&";
(* Create some PlotLegends *)
pl = ",PlotLegends\[Rule]" <> StringRiffle[ToString[#], ","] &@vbles;
(* Build the model expression *)
ToExpression[
 StringTemplate[
   "model[`pats`][t_]:=`ms`@Through[pfun[`params`][t],List]\
/;And@@NumericQ/@{`params`};"][<|"pats" -> spats, "params" -> sparms, 
   "ms" -> modelstring|>]]
(* Create slider model to help with initial guess *)
globalstring = 
  StringTemplate["global={`params`};"][<|"params" -> sparms|>];
mantemp = 
  "Manipulate[`g`\[IndentingNewLine]Show[lp,Plot[Evaluate@(model[`\
params`][t]),{t,0,tmax},PlotRange\[Rule]{0.,.8}`pl`],ImageSize->Large]\
,`sliders`]";
ToExpression@
 StringTemplate[mantemp][<|"sliders" -> sliders, "params" -> sparms, 
   "g" -> globalstring, "pl" -> pl|>]
(*Display global variable of current slider parameters *)
Dynamic@global

We can adjust the sliders to get a better initial guess for the parameters by trying to manually the sliders from the initial values. I intentionally did not manually fit as well as I could have to test how well FindFit change parameters

Here is the code to do the fit and create a slider model initialized at the starting parameters. The "Gradient" method seemed to work reasonably well, but the methods will need to be experimented with to achieve the best results.

(*Initial parameter guesses from global*)
initguess = MapThread[List, {parms, First@Dynamic@global}];
fit = Quiet@
   FindFit[data, 
    model[k1, k11, k2, k22, k3, k33, k4, k44, k5][t][[7]], initguess, 
    t, Method -> "Gradient"];
(*Create new Manipulate Based on Fit*)
mat2 = mat;
mat2[[All, 2]] = fit[[All, 2]];
ToExpression@
 StringTemplate[mantemp][<|"sliders" -> sfun[mat2], 
   "params" -> sparms, "pl" -> pl|>]

I perturbed the initial guesses and re-ran the fit. The fit looks similar, but the fitted parameters are clearly different indicating the non-uniqueness depending on initial guesses. Good initial and long term asymptotic data would aid in the parameter estimation.

Answer 2

You can get a series expansion for small t with AsymptoticDSolveValue; maybe that's good enough if you go to high-enough order in t:

AsymptoticDSolveValue[{
  EA'[t] == k1 e[t] a[t] + k22 EAB[t] - k2 b[t] EA[t] - k11 EA[t], 
  EB'[t] == k4 e[t] b[t] - k44 EB[t] + k33 EAB[t] - k3 a[t] eb[t], 
  EAB'[t] == k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] + k33 EAB[t],
  p'[t] == k5 EAB[t], 
  e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t], 
  a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t], 
  b'[t] == k44 EB[t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t], 
  EA[0] == 0, EB[0] == 0, EAB[0] == 0, a[0] == 1, e[0] == 1/10, b[0] == 100, p[0] == 0},
  {EAB[t], EB[t], EA[t], a[t], e[t], b[t], p[t]},
  {t, 0, 3}]

(output is third-order series expansion in t)