# Solving coupled ODEs with NDSolve to fit experimental data

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.23*10^4, A2 == 0.048*10^4, A3 == 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[] /. 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.23*10^4,
A2 == 0.048*10^4, A3 == 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.23*10^4,
A2 == 0.048*10^4, A3 == 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.23*10^4,
A2 == 0.048*10^4, A3 == 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.

• Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Dr. belisarius Oct 1 '15 at 5:55

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.23*10^4, A2 == 0.048*10^4, A3 == 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 :)

• If you provide reasonable starting values for the parameters perhaps it will end faster – Dr. belisarius Oct 1 '15 at 17:31
• thank you so much for your help! I went through your simplified model to better understand it, and it is very "clean" code. I liked it very much. Than I finalized first the fit of your proposed 'data' and plotted it. i.stack.imgur.com/jUmnT.png It fitted very well! What I gonna do next is to adapt your suggestions to my model. I hope it works! :) Thanks once again! – renatons Oct 1 '15 at 18:27