# Modeling experimental data with differential equations

I have the following two sets of experimental data, which show the dependencies of two quantities, namely, $$S$$ and $$B$$, on time ($$0$$ h, $$3$$ h, $$6$$ h, $$9$$ h, $$15$$ h, $$18$$ h, $$21$$ h, and $$24$$ h):

Sdata = {{0, 9.74},{3, 4.92},{6, 8.29},{9, 5.54},{15, 2.08},{18, 1.38},{21, 1.99},{24, 0.893}};

Bdata = {{0, 0.915094},{3, 0.736097},{6, 0.793694},{9, 0.833664},{15, 1},{18, 0.99578},{21, 0.897964},{24, 0.214499}};


I'm trying to model the dynamics of the above data by the following differential equations:

$$\frac{dS(t)}{dt} = 0.31 S(t) \Big( 1 - \frac{S(t)}{6.19} \Big) - \frac{a}{1 + b B(t)} S(t),$$

$$\frac{dB(t)}{dt} = c (1 - B(t)) - d B^2(t) \Big( \frac{1 - B(t)}{B(t)} \Big)^{0.96},$$

where $$a$$, $$b$$, $$c$$, and $$d$$ are some real, preferably positive, constants to be determined from the data.

S'[t] == 0.31 S[t] (1 - S[t]/6.19) - a/(1 + b B[t]) S[t]

B'[t] == c (1 - B[t]) - d B[t]^2 ((1 - B[t])/B[t])^0.96


Now, I'm trying to obtain the constants of my model using Mathematica. I found this answer might be helpful; but, this answer only deals with one set of data.

Any help or hint is appreciated.

EDIT

With the help of @ydd answer, I could modify the model and could do fitting with Mathematica:

Modified model:

$$\frac{dS(t)}{dt} = - \frac{a}{1 + B(t)} S(t),$$

$$\frac{dB(t)}{dt} = \frac{c}{1 + S(t)} B(t) - d B^2(t) \Big( \frac{1 - B(t)}{B(t)} \Big)^n,$$

where $$a$$, $$c$$, $$d$$, and $$n$$ are constants to be determined from the data.

We have:

Sdata = {{0, 9.74}, {3, 4.92}, {6, 8.29}, {9, 5.54}, {15, 2.08}, {18,
1.38}, {21, 1.99}, {24, 0.893}};
Bdata = {{0, 0.915094}, {3, 0.736097}, {6, 0.793694}, {9,
0.833664}, {15, 1}, {18, 0.99578}, {21, 0.897964}, {24, 0.214499}};
order = 1;
interpolatedData = {intS,
intB} = {Interpolation[Sdata, InterpolationOrder -> order],
Interpolation[Bdata, InterpolationOrder -> order]};
sys = {S'[t] == -a/( 1 + B[t]) S[t],
B'[t] == c B[t]/(1 + S[t]) - d B[t]^2 ((1 - B[t])/B[t])^n};
squareDiffs = MapApply[(#1 - #2)^2 &, sys];
withInt[t_] = squareDiffs /. {S -> intS, B -> intB};
totalSquaredError = Total@Flatten[withInt /@ (Range[0, 24, 3])];
forMin = Join[{totalSquaredError}, restrictions];
restrictions = Thread[{a, c, d, n} > 0];
{resid, bestFitParams} =
NMinimize[forMin, {a, c, d, n}, Method -> "RandomSearch"]
init = {S[0] == Sdata[[1, 2]], B[0] == Bdata[[1, 2]]};
new = Join[sys /. bestFitParams, init];
{sSol[t_], bSol[t_]} = NDSolveValue[new, {S[t], B[t]}, {t, 0, 24}];

lp = ListPlot[Bdata];
p = Plot[bSol[t], {t, 0, 24}, PlotRange -> All];
Show[lp, p]
lpp = ListPlot[Sdata];
pp = Plot[sSol[t], {t, 0, 24}, PlotRange -> All];
Show[lpp, pp]


Now, I have a question which is mostly mathematical: Can one slightly modify these two differential equations in order the fitted curve to the red data, i.e., $$B(t)$$, also captures the bump in the data around $$t = 15$$?

• Quick little experiment. B and c do not depend on S. Find B and c from your data using perhaps that answer. Then start again with known B,c and unknown S,a,b that depends only on your data to find S,a,b. Report if that works. There might be better a,b,c,B,S if you did it all at once, but this should let you try something to get started. If you plot your solution on top of your data does it look good enough?
– Bill
Commented Nov 22, 2023 at 12:27
• That answer has another link to a multiple data set approxiamtion technique with ParametricNDSolveValue mathematica.stackexchange.com/questions/28461/…. By the disproportionate values of 10^7, rescale the system for {S/10^7, B}. A fast simualation with all constants =1 and the given start value at t=0 yields complex values by the irrational power ^0.96 of negative values. Commented Nov 22, 2023 at 12:36
• @Tim Your second ode solves for B[t] with parameters c,d. You should first check your model for this single ode. My simulations don't evaluate in this case (NMinimize, NonlinearModelFit) Commented Nov 22, 2023 at 17:01

## Update

Sdata = {{0, 9.74}, {3, 4.92}, {6, 8.29}, {9, 5.54}, {15, 2.08}, {18,
1.38}, {21, 1.99}, {24, 0.893}};
Bdata = {{0, 0.915094}, {3, 0.736097}, {6, 0.793694}, {9,
0.833664}, {15, 1}, {18, 0.99578}, {21, 0.897964}, {24, 0.214499}};
sys = {S'[t] == -a/(1 + B[t]) S[t],
B'[t] == c B[t]/(1 + S[t]) - d B[t]^2 ((1 - B[t])/B[t])^n};


Using the new system of differential equations, and assuming a value for n we can get a series expansion around $$t=0$$ of the solutions for $$S(t)$$ and $$B(t)$$ using AsymptoticDSolve. AsymptoticDSolve won't work for symbolic n, so I separately found an (approximately) optimal value for $$0 (more on how I got this value of $$n$$ further down). I found n=0.5 is pretty good. I get a 5th order expansion of the solution to our system of diff. equations:

sys0 = sys /. {n -> 0.5};
pars = {a, d, c};
init = {S[0] == Sdata[[1, 2]], B[0] == Bdata[[1, 2]]};
withInit = Join[sys0, init];
{sTrue[t_], bTrue[t_]} =
AsymptoticDSolveValue[withInit, {S[t], B[t]}, {t, 0, 5}];



I realized after writing my last answer that $$S$$ is about an order of magnitude larger than $$B$$ in terms of function value, so summing the total error of the two would bias good $$S$$ solutions, and not so much good $$B$$ solutions. So instead I sort of normalize the error by dividing the norm of the error by the norm of the data:


sVals = Sdata[[All, 2]];
bVals = Bdata[[All, 2]];

normS = Norm@sVals;
normB = Norm@bVals;

tVals = Sdata[[All, 1]];

normErrS = Norm[sTrue /@ tVals - sVals];
normErrB = Norm[bTrue /@ tVals - bVals];

relErrS = normErrS/normS;
relErrB = normErrB/normB;
obj = relErrS + relErrB;

{err, bestParams} = NMinimize[obj, pars, Method -> "RandomSearch"]
(*{0.334796, {a -> 0.185199, d -> 0.259166, c -> 0.583751}}*)


And the $B$ fit is now improved:

{sSol[t_], bSol[t_]} = {sTrue[t], bTrue[t]} /. bestParams;
lp = ListPlot[Bdata, PlotStyle -> Red];
p = Plot[bSol[t], {t, 0, 24}, PlotRange -> All,
PlotStyle -> Directive[Dashed, Orange]];
Show[lp, p, PlotRange -> All, PlotLabel -> "B(t)"]
lpp = ListPlot[Sdata];
pp = Plot[sSol[t], {t, 0, 24}, PlotRange -> All];
Show[lpp, pp, PlotLabel -> "S(t)"]


The fitted 5th-order expansions of $$S(t)$$ and $$B(t)$$ are:

sSol[t]
bSol[t]
(*"
9.74 - 0.941906 t + 0.0415182 t^2 - 0.00119407 t^3 +
0.0000434503 t^4 - 9.82512*10^-7 t^5

0.915094 - 0.0163685 t - 0.00056303 t^2 + 0.00013121 t^3 +
9.24318*10^-6 t^4 - 6.08513*10^-7 t^5
*)


Note that the fits are sensitive to small changes in parameter values. Rationalizing the coefficients to within 10^-6 completely changes the quality of the fit:

sSol[t_] = Rationalize[sSol[t], 10^-6];
bSol[t_] = Rationalize[bSol[t], 10^-6];
lp = ListPlot[Bdata, PlotStyle -> Red];
p = Plot[bSol[t], {t, 0, 24}, PlotRange -> All,
PlotStyle -> Directive[Dashed, Orange]];
Show[lp, p, PlotRange -> All, PlotLabel -> "B(t)"]
lpp = ListPlot[Sdata];
pp = Plot[sSol[t], {t, 0, 24}, PlotRange -> All];
Show[lpp, pp, PlotLabel -> "S(t)"]



The value n=0.5 was determined just by eye by finding the $$n$$ value that minimizes the objective function (n has to be optimized separately because AsymptoticDSolveValue does not produce a solution with symbolic n). I only searched the points 0.1,0.2,,...1. but I figured it was good enough for a quick n fit:

errGivenN[nIn_] := With[{n0 = nIn},
sys0 = sys /. n -> n0;
withInit = Join[sys0, init];
{sTrue[t_], bTrue[t_]} =
AsymptoticDSolveValue[withInit, {S[t], B[t]}, {t, 0, 5}];

normErrS = Norm[sTrue /@ tVals - sVals];
normErrB = Norm[bTrue /@ tVals - bVals];

relErrS = normErrS/normS;
relErrB = normErrB/normB;
obj = relErrS + relErrB;

{err, bestParams} = NMinimize[obj, pars, Method -> "RandomSearch"];
err
]

testVals = {#, errGivenN[#]} & /@ Range[0, 1, 0.1];
ListLinePlot[testVals]


First get a linear interpolation between data points:

order = 1;
interpolatedData = {intS,
intB} = {Interpolation[Sdata, InterpolationOrder -> order],
Interpolation[Bdata, InterpolationOrder -> order]};


Then input the system of differential equations:

sys = {S'[t] == 0.31 S[t] (1 - S[t]/6.19) - a/(1 + b B[t]) S[t],

B'[t] == c (1 - B[t]) - d B[t]^2 ((1 - B[t])/B[t])^0.96};


Let's turn our system into a minimization problem by replacing lhs == rhs with (lhs-rhs)^2:

squareDiffs = MapApply[(#1 - #2)^2 &, sys]

(*{((a S[t])/(1 + b B[t]) - 0.31 (1 - 0.161551 S[t]) S[t] +
Derivative[1][S][t])^2, (-c (1 - B[t]) +
d ((1 - B[t])/B[t])^0.96 B[t]^2 + Derivative[1][B][t])^2}*)


And now let's replace S[t] and B[t] in our system with the linear data interpolations (this also replaces the derivatives as well):

withInt[t_] = squareDiffs /. {S -> intS, B -> intB};


And let's sum all the squared errors at the datapoints t = 0,3,...,24:

totalSquaredError = Total@Flatten[withInt /@ (Range[0, 24, 3])];


And now we can find the positive real {a,b,c,d} that minimizes totalSquaredError. I found RandomSearch works well in this case::

restrictions = Thread[{a, b, c, d} > 0];
forMin = Join[{totalSquaredError}, restrictions];
{resid, bestFitParams} =
NMinimize[forMin, {a, b, c, d}, Method -> "RandomSearch"]

(*{11.5479, {a -> 294399., b -> 1.29843*10^7, c -> 0., d -> 0.345276}}*)


We can now plug these parameter values back into our differential system. We can then solve the system (using Bdata[[1,2]] and Sdata[[1,2]] as initial conditions) and compare the solutions to the numerical data:

init = {S[0] == Sdata[[1, 2]], B[0] == Bdata[[1, 2]]};
new = Join[sys /. bestFitParams, init];
{sSol[t_], bSol[t_]} = NDSolveValue[new, {S[t], B[t]}, {t, 0, 24}];

lp = ListPlot[{Sdata, Bdata}, PlotLegends -> {S, B}];
p = Plot[{sSol[t], bSol[t]}, {t, 0, 24}, PlotRange -> All];
Show[lp, p]


• I will probably be busy rest of day but I will try to think about this. I guess the first thing to ask is where your differential system of equations comes from? I think it would be wrong to modify the system to fit the data. Instead, the system should be based on some theory that can describe the observed phenomena .
– ydd
Commented Nov 22, 2023 at 19:18
• Another thing to ask is can you put any tighter restrictions on the parameter values {a,b,c,d}? This could also potentially help.
– ydd
Commented Nov 22, 2023 at 19:19
• @Tim: Sjoerd Smit does this in their answer (see the parameter values s0 and b0 in ParametricNDSolve. Also, I just updated my answer using your new model.
– ydd
Commented Nov 23, 2023 at 22:23
• My worry about the fit being very sensitive to parameter values is just that it might be hard to find the actual best parameter values, as you'd have make very small step sizes in the solution space to match how rapidly the objective function changes with parameter value.  The model with $c(1-B(t))$ is fine near $n=1$. The fits are ok. I really think you should stick to one model for fitting that can capture the obeservations. Maybe ask someone who knows about how the skin properties represented by $B$ affect the number of bacteria colonies.
– ydd
Commented Nov 24, 2023 at 17:29
• The use of Rationalize[] to gauge the sensitivity of the model to parameters seems to be misleading. It seems to have rationalized the coefficients of powers of time rather than the model coefficients $a,c,d$ and $n$. With this time scale and the associated very-low-magnitude coefficient of thet^5 term, Rationalize[,10^-6] cut out the t^5 term completely. That now-omitted term with its negative coefficient is what keeps the Updated model from exploding at late times.
– EdM
Commented Nov 28, 2023 at 16:06

This should help you underway, but I couldn't find good starting values for the parameters:

Sdata = {{0, 9.74}, {3, 4.92}, {6, 8.29}, {9, 5.54}, {15, 2.08}, {18, 1.38}, {21, 1.99}, {24, 0.893}};
Bdata = {{0, 0.915094}, {3, 0.736097}, {6, 0.793694}, {9,  0.833664}, {15, 1}, {18, 0.99578}, {21, 0.897964}, {24,  0.214499}};

eqs = {
S'[t] == 0.31 S[t] (1 - S[t]/6.19) - a/(1 + b B[t]) S[t],
B'[t] == c (1 - B[t]) - d B[t]^2 ((1 - B[t])/B[t])^0.96,
S[0] == s0,
B[0] == b0
}

sol = ParametricNDSolveValue[
eqs,
{S, B},
{t, 0, 24},
{s0, b0, a, b, c, d}
];

(*Helper function*)
eval[{s0_?NumericQ, b0_, a_, b_, c_, d_}, t_?NumericQ] := Through[Once[sol[s0, b0, a, b, c, d]][t]];

ResourceFunction["MultiNonlinearModelFit"][
{Sdata, Bdata},
eval[{s0, b0, a, b, c, d}, t],
{{s0, 10}, {b0, 1}, a, b, c, d},
t
]


If you find good starting values for a, b, c and d, I think it should work.

• @Tim Honestly, it's different to say without good domain knowledge of what the data represent. Sounds more like a question for another SE site. But at least you should now have the tools for tweaking the model and comparing the results. Commented Nov 23, 2023 at 9:04