Solving coupled ODEs with NDSolve to fit experimental data

Question

I need to fit a multi-exponential data (Intensity x Time) using a model consisting of three differential equations. I want to find the parameters that gives the best fit to the experimental data.

Because the set o differential equations have to be numerically solved, I tried the following approach: I used a initial guess and the _?NubmerQ command as being the numerical inputs for each parameter. The function is defined as a block. The function is then evaluated with the NDsolve command to find A1[t], A2[t] and A3[t] with the supplied values (initial guess) of the parameters.

Although the solution of the differential equations are A1[t], A2[t] and A3[t], I need some normalization constants (l, m and n) multiplying the respective solutions, which should also be evaluated through the fit routine. This make my final fitting function to be (l* A1[t] + m* A2[t] + n*A3[t]).

So, finally I tried to calculate the least squares difference between the evaluated fitting function and the experimental data.

sse[l_?NumberQ, m_?NumberQ, n_?NumberQ, kb_?NumberQ, kab_?NumberQ, 
kbc_?NumberQ, kcb_?NumberQ, kba_?NumberQ] := 
Block[{solTotal, A1, A2, A3}, solTotal = NDSolve[{
  A1'[t] == -(ka + kab)*A1[t] + kba*A2[t], 
  A2'[t] == -(kb + kba + kbc)*A2[t] + kab*A1[t] + kcb*A3[t], 
  A3'[t] == -(kc + kcb)*A3[t] + kbc*A2[t],
  A1[0] == 0.23*10^4, A2[0] == 0.048*10^4, A3[0] == 0},
 {A1[t], A2[t], A3[t]}, {t, 0, 90.0*10^-6}];
Apply[Plus, (fexp - ((l*A1[t] + m*A2[t] + n*A3[t]) /. 
     solTotal[[1]] /. t -> texp))^2]];

fexp and text are the y-values and x-values respectively, of the experimental data. ka and kc are known parameters. Although I may have more parameters to be fitted than functions, some of them are linked to other through known constants and I have specific range of constraints for each one of them. But first I would like to have a more general fitting function with independent parameters.

Finally, I want to use FindMinimum to find the minimum value of sse[l, m, n, kb, kab, kbc, kcb, kba], given the initial guess and constraints.

bf = FindMinimum[{sse[l, m, n, kb, kab, kbc, kcb, kba], 
0 < l < 1*10^-5 && 0 < m < 1*10^-5 && 0 < n < 1*10^-5}, {{l, 
0.1*10^-5}, {m, 0.1*10^-5}, {n, 0.1*10^-5}, {kb, 6*10^6}, {kab, 
2.5*10^7}, {kbc, 1.5*10^7}, {kcb, 1.5*10^6}, {kba, 3*10^6}}]

Although I can make the code to work (kind of...) it only minimizes the linear parameters l, m and n. It doesn't minimizes the parameters kb, kab, kbc, kcb, kba, which would be the goal here. Doesn't matter what is the given initial guess, it maintain the provided values and only minimizes l, m, n.

(*output: {4.80374*10^-6, {l -> 5.82322*10^-6, m -> 2.90073*10^-6, 
n -> 3.63747*10^-6, kb -> 6.*10^6, kab -> 2.5*10^7, kbc -> 1.5*10^7,
kcb -> 1.5*10^6, kba -> 3.*10^6}}*)

Alternatively to the NDSolve and FindMinimum approach, I also tried using ParametricNDSolve to evaluate A1[t], A2[t] and A3[t] separately to latter use the NonlinearModelFit to fit the experimental data.

difeqA1 = 
ParametricNDSolve[{A1'[t] == -(ka + kab)*A1[t] + kba*A2[t], 
A2'[t] == -(kb + kba + kbc)*A2[t] + kab*A1[t] + kcb*A3[t], 
A3'[t] == -(kc + kcb)*A3[t] + kbc*A2[t], A1[0] == 0.23*10^4, 
A2[0] == 0.048*10^4, A3[0] == 0}, 
A1, {t, 0, 90*10^-6}, {kb, kab, kbc, kcb, kba}];

difeqA2 = 
ParametricNDSolve[{A1'[t] == -(ka + kab)*A1[t] + kba*A2[t], 
A2'[t] == -(kb + kba + kbc)*A2[t] + kab*A1[t] + kcb*A3[t], 
A3'[t] == -(kc + kcb)*A3[t] + kbc*A2[t], A1[0] == 0.23*10^4, 
A2[0] == 0.048*10^4, A3[0] == 0}, 
A2, {t, 0, 90*10^-6}, {kb, kab, kbc, kcb, kba}];

difeqA3 = 
ParametricNDSolve[{A1'[t] == -(ka + kab)*A1[t] + kba*A2[t], 
A2'[t] == -(kb + kba + kbc)*A2[t] + kab*A1[t] + kcb*A3[t], 
A3'[t] == -(kc + kcb)*A3[t] + kbc*A2[t], A1[0] == 0.23*10^4, 
A2[0] == 0.048*10^4, A3[0] == 0}, 
A3, {t, 0, 90*10^-6}, {kb, kab, kbc, kcb, kba}];

nlfit = NonlinearModelFit[data, 
{l*(A1[kb, kab, kbc, kcb,kba] /. difeqA1) + m*(A2[kb, kab, kbc, kcb, kba] /.difeqA2) + n*(A3[kb, kab, kbc, kcb, kba] /. difeqA3), 
4*10^6 < kb < 8*10^6, 1*10^7 < kab < 5*10^7, 1*10^7 < kbc < 5*10^7, 1*10^6 < kcb < 5*10^6, 1*10^6 < kba < 5*10^6}, 
{l, m, n, {kb, 6*10^6}, {kab, 2.5*10^7}, {kbc, 1.5*10^7}, {kcb, 1.5*10^6}, {kba, 3*10^6}}, t]

This latter attempt of using ParametricNDSolve wasn't very helpful since I could not even make the NonlinearModelFit to work (or the FindFit).

So, my best shot so far is the NDSolve with FindMinimum. However, I have the issue that it only minimizes the linear constants l, m, and n.

I'm kind of lost right now, and I would really appreciate any help or hint from you.

Thanks a lot, in advance, for your attention!

Answer 1

I simplified your model just a bit.

sol = ParametricNDSolveValue[{
   A1'[t] == -k1*A1[t] + k2*A2[t],
   A2'[t] == -k3*A2[t] + k4*A1[t] + k5*A3[t],
   A3'[t] == -k6*A3[t] + k7*A2[t],
   A1[0] == 0.23*10^4, A2[0] == 0.048*10^4, A3[0] == 0},
  {A1, A2, A3}, {t, 0, 1}, {k1, k2, k3, k4, k5, k6, k7}]

model[k1_, k2_, k3_, k4_, k5_, k6_, k7_, l_, m_, n_][t_] := 
      {l, m, n}.Through[sol[k1, k2, k3, k4, k5, k6, k7][t], List] /; 
   And @@ NumericQ /@ {k1, k2, k3, k4, k5, k6, k7, l, m, n, t};

Now we generate a dataset:

data = Table[{t, model[2, 3, 1.5, 1, 2, 3, 1, 1, 5, 1][t]}, {t, 0, 1, 10^-2}];
ListLinePlot@data

And now we try Non Linear Model Fit:

Monitor[NonlinearModelFit[
   data, {model[k1, k2, k3, k4, k5, k6, k7, l, m, n][t], 
    And @@ Thread[6 > {k1, k2, k3, k4, k5, k6, k7, l, m, n} > 0]}, {k1, k2, k3, 
    k4, k5, k6, k7, l, m, n}, {t}, 
   Method -> {NMinimize, Method -> "SimulatedAnnealing"}],
  {k1, k2, k3, k4, k5, k6, k7, l, m, n}];

Minimizing ten variables is a lot of work, and I haven't time to finish the calcs to see the result, but at least it's running :)