0
$\begingroup$

The following code attempts to fit a system of differential equations (enzyme kinetics) to some data, which for now happens to be fake data extracted from a solution with a particular set of parameters. The approach is almost identical to the one presented here by Oleksandr R.

ClearAll["Global`*"];

rates = {
R1[t] == 
Subscript[r, COD]*X[t]*
 Subscript[S, S][t]/(Subscript[k, S] + Subscript[S, S][t])*
 Subscript[S, Mox][t]/(Subscript[k, Mox] + Subscript[S, Mox][t]), 
R2[t] == 
Subscript[r, NO3]*X[t]*
 Subscript[S, NO3][t]/(Subscript[k, NO3] + Subscript[S, NO3][t])*
 Subscript[S, Mred][
   t]/(Subscript[k, Mred1] + Subscript[S, Mred][t]), 
R3[t] == 
Subscript[r, NO2]*X[t]*
 Subscript[S, NO2][t]/(Subscript[k, NO2] + Subscript[S, NO2][t])*
 Subscript[S, Mred][
   t]/(Subscript[k, Mred2] + Subscript[S, Mred][t]), 
R4[t] == 
Subscript[r, NO]*X[t]*
 Subscript[S, NO][t]/(Subscript[k, NO] + Subscript[S, NO][t])*
 Subscript[S, Mred][
   t]/(Subscript[k, Mred3] + Subscript[S, Mred][t]), 
R5[t] == 
Subscript[r, N2O]*X[t]*
 Subscript[S, N2O][t]/(Subscript[k, N2O] + Subscript[S, NO2][t])*
 Subscript[S, Mred][
   t]/(Subscript[k, Mred4] + Subscript[S, Mred][t])};

conc = {Subscript[S, NO3]'[t] == -R2[t], 
Subscript[S, NO2]'[t] == R2[t] - R3[t], 
Subscript[S, NO]'[t] == R3[t] - R4[t], 
Subscript[S, N2O]'[t] == 0.5 R4[t] - R5[t], 
Subscript[S, N2]'[t] == R5[t], Subscript[S, S]'[t] == -R1[t], 
Subscript[S, Mox]'[t] == -(1 - Subscript[Y, H]) R1[t] + R2[t] + 
 0.5 R3[t] + 0.5 R4[t] + R5[t], X'[t] == Subscript[Y, H]*R1[t], 
Subscript[S, Mox][t] + Subscript[S, Mred][t] == Subscript[c, tot]};

Subscript[r, NO3] = 0.045; Subscript[r, NO2] = 0.059; 
Subscript[r, NO] = 0.56; Subscript[r, N2O] = 0.23;

Subscript[k, S] = 0.1; Subscript[k, NO3] = 0.018; 
Subscript[k, NO2] = 0.0041;
Subscript[k, NO] = 0.000011;
Subscript[k, N2O] = 0.0025;

Subscript[k, Mox] = 0.0001;
Subscript[Y, H] = 0.6;
Subscript[c, tot] = 0.01;

initc = {WhenEvent[t == 0.5, Subscript[S, NO2][t] -> 5/14], 
Subscript[S, NO3][0] == 1, Subscript[S, NO2][0] == 0, 
Subscript[S, NO][0] == 0, Subscript[S, N2O][0] == 0, 
Subscript[S, N2][0] == 0, Subscript[S, S][0] == 170/14, 
Subscript[S, Mox][0] == 0.005, X[0] == 70};
Subscript[t, 0] = 1.5;

params = {Subscript[r, COD], Subscript[k, Mred1], 
Subscript[k, Mred2],
Subscript[k, Mred3], Subscript[k, Mred4]};

sol = ParametricNDSolveValue[
Join[rates, conc, initc], {Subscript[S, NO3], Subscript[S, NO2], 
    Subscript[S, NO], Subscript[S, N2O], Subscript[S, N2], 
Subscript[S, S], 
    Subscript[S, Mox], Subscript[S, Mred], X}, {t, 0., 
Subscript[t, 0]}, params, 
  Method -> {"TimeIntegration" -> "StateSpace"}]


abscissae = Range[0., 1.5, 0.1];
ordinates = With[{r = 0.064, kM1 = 0.0015, kM2 = 0.00058,
kM3 = 0.000010, kM4 = 0.00024}, 
Through[sol[r, kM1, kM2, kM3, kM4][abscissae], List]];

data = ordinates + 
RandomVariate[NormalDistribution[0, 0.1^2], 
Dimensions[ordinates]];
ListLinePlot[data, DataRange -> {0, 1.5}, PlotRange -> All, 
AxesOrigin -> {0, 0}];

transformedData = {ConstantArray[Range@Length[ordinates], 
  Length[abscissae]] // Transpose, 
ConstantArray[abscissae, Length[ordinates]], ordinates}~
Flatten~{{2, 3}, {1}};

model[r_, kM1_, kM2_, kM3_, kM4_][i_, t_] := 
Through[sol[r, kM1, kM2, kM3, kM4][t], List][[i]] /; 
And @@ NumericQ /@ {r, kM1, kM2, kM3, kM4, i, t};

fit = NonlinearModelFit[
transformedData, {model[r, kM1, kM2, kM3, kM4][i, t], {0. < r < 1.,
  0. < kM1 < 0.01, 0. < kM2 < 0.05, 0. < kM3 < 0.001, 
 0. < kM4 < 0.01}}, {{r, 0.1}, {kM1, 0.001}, {kM2, 0.01}, {kM3, 
 0.001}, {kM4, 0.001}}, {i, t}, MaxIterations -> 10, 
Method -> {"NMinimize", 
 Method -> {"SimulatedAnnealing", "PerturbationScale" -> 3}}, 
EvaluationMonitor :> 
Print["r=", r, " kM1=", kM1, " kM2=", kM2, " kM3=", kM3, " kM4=", 
 kM4]];
fit["BestFitParameters"]

However, it appears that after a couple of dozens of evaluations, NonlinearModelFit stalls on an iteration of ParametricNDSolveValue called with a certain "bad" set of parameters, and is unable to get out of it. Is there a work-around to abort ParametricNDSolveValue in order for NonlinearModelFit to move onto the next set of parameters?

Note: Method -> "NMinimize in NonlinearModelFit is the only main method that allows for fitting to occur, and Method -> {"TimeIntegration" -> "StateSpace"} in ParametricNDSolveValue is the only one that seems to give the desired numerical solutions.

$\endgroup$
  • 1
    $\begingroup$ Does a plot of the transformedData suggest something either wrong with the data or that the model is not capable of a good fit to that data? Using ListPointPlot3D[transformedData] shows that when i=6 there is a very different response which is orders of magnitude larger than all of the other values of i. $\endgroup$ – JimB Oct 26 '17 at 0:11
  • $\begingroup$ transformedData is actually fake data from a particular solution that doesn't cause any problems. With so many dependent variables in the model, some of them are expected to be of different order of magnitude than others. It indeed should be that i=6,9 are larger than the rest. The initial parameters in the fit are also somewhat close to the original data. $\endgroup$ – HappyTreeFriend Oct 26 '17 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.