Chemical Kinetics: fitting multiple sets of time series data using multiple differential reaction equations

Question

Let me start by saying that this is my first post, so if there are suggestions on proper formatting that is annoying you, please let fly and I will do my best to modify.

I was doing well solving the topic in the title until I attempted to model second order reactions, which seemed to "break" Mathematica. I'll use an example of analogous code to get my point across.

So first I define my time series data as follows, including initial concentrations:

data = {{0., 182.115, 0.}, {1., 166.486, 1.87153}, {5., 136.178, 
   9.76618}, {10., 112.277, 19.1688}, {15., 111.448, 27.2978}, {25., 
   81.5446, 43.3609}, {35., 73.6962, 56.2542}, {45., 62.9892, 
   62.6874}, {55., 50.7161, 68.5274}, {65., 45.2139, 79.0054}};

{time, Xylose, Furfural} = Transpose[data];

Fin = Furfural[[1]];

Xin = Xylose[[1]];

Then I set up the various rate equations for the chemical species, and solve for the rate constants; so far, this is all well and good, and the output can be seen below.

sol = ParametricNDSolveValue[{X'[t] == -a X[t] - b X[t], X[0] == Xin, 
    F'[t] == a X[t] - c F[t], F[0] == Fin}, {X, F}, {t, 0, 70}, {a, b,
     c}];

transformeddata = {ConstantArray[Range@Length[{Xylose, Furfural}], 
      Length[time]] // Transpose, 
    ConstantArray[time, Length[{Xylose, Furfural}]], {Xylose, 
     Furfural}}~Flatten~{{2, 3}, {1}};

model[a_, b_, c_][i_, t_] := 
  Through[sol[a, b, c][t], List][[i]] /; 
   And @@ NumericQ /@ {a, b, c, i, t};

fit = NonlinearModelFit[transformeddata, 
   model[a, b, c][i, t], {a, b, c}, {i, t}];

fit["RSquared"]

fit["ParameterTable"]

Show[ListPlot[Table[Take[data, All, {1, l, l - 1}], {l, 2, 3}], 
  PlotStyle -> PointSize[0.02]], 
 Plot[{fit[1, t], fit[2, t]}, {t, 0, 70}]]

(This works great for modeling first order reaction kinetics)

The "Show stopper" comes when I try to represent a reaction with second order kinetics; all the above are first order kinetics. So for instance when I change the following:

sol = ParametricNDSolveValue[{X'[t] == **-a X[t]** - b X[t], X[0] == Xin, 
    F'[t] == **a X[t]** - c F[t], F[0] == Fin}, {X, F}, {t, 0, 70}, {a, b,
     c}];

into

sol = ParametricNDSolveValue[{X'[t] == **-a X[t] X[t]** - b X[t], X[0] == Xin, 
    F'[t] == **a X[t] X[t]** - c F[t], F[0] == Fin}, {X, F}, {t, 0, 70}, {a, b,
     c}];

then the output I get has the following errors:

In short, I'm not sure why this solve function is having trouble with representing second order reactions in this way. Any advice or information would be greatly appreciated.

Answer 1

You need to provide starting values for NonlinearModelFit.

First, I took your model and played with its parameters untill it looked good:

Show[Plot[model[0., 0.03, 0][1, t], {t, 0, 70}, PlotRange -> All], 
 ListPlot@Transpose@{time, Xylose}]

(You could wrap it in Manipulate.) It's very sensitive to the values of a, moderately sensitive to b and unnoticeably sensitive to c.

Next, I used those parameters as starting values (no errors at all):

fit = NonlinearModelFit[transformeddata, 
   model[a, b, c][i, t], {{a, 0.}, {b, 0.03}, {c, 0}}, {i, t}];

fit["RSquared"]

0.991887

fit["ParameterTable"]

Show[ListPlot[Table[Take[data, All, {1, l, l - 1}], {l, 2, 3}], 
  PlotStyle -> PointSize[0.02]], 
 Plot[{fit[1, t], fit[2, t]}, {t, 0, 70}]]