Using ParametricNDSolveValue and MultiNonlinearModelFit to fit an ODE system to datasets

I asked a question here about fitting an ODE system to given datasets. The great answer of @ydd solved the problem nicely. In the mentioned answer, the initial values are taken as the initial points in the datasets. Now, I want to make the initial values also parameters to be determined by fitting, so I'm using @Sjoerd Smit's answer for this, but it doesn't produce results.

My datasets are:

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}};


The datasets with standard deviations:

SdatawithSD = {{0, Around[9.74, 2.89]}, {3, Around[4.92, 1.65]}, {6,
Around[8.29, 4.04]}, {9, Around[5.54, 2.45]}, {15,
Around[2.08, 1.91]}, {18, Around[1.38, 0.962]}, {21,
Around[1.99, 2.41]}, {24, Around[0.893, 0.359]}};

BdatawithSD = {{0, Around[0.915094, 0.1]}, {3,
Around[0.736097, 0.091668]}, {6, Around[0.793694, 0.082575]}, {9,
Around[0.833664, 0.070242]}, {15, Around[1, 0.002851]}, {18,
Around[0.99578, 0.015591]}, {21, Around[0.897964, 0.04783]}, {24,
Around[0.214499, 0.01231]}};


My ODE system is:

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

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

where $$a$$, $$c$$, $$d$$, and $$n$$ are constants.

The Mathematica code:

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] == -a/(1 + B[t]) S[t],
B'[t] == c (1 - B[t])/(1 + S[t]) - d B[t]^2 ((1 - B[t])/B[t])^n,
S[0] == s0, B[0] == b0}

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

eval[{s0_?NumericQ, b0_, a_, c_, d_, n_}, t_?NumericQ] :=
Through[Once[sol[s0, b0, a, c, d, n]][t]];

ResourceFunction["MultiNonlinearModelFit"][{Sdata, Bdata},
eval[{s0, b0, a, c, d, n},
t], {{s0, 9}, {b0, 0.9}, {a, 0.18}, {c, 0.58}, {d, 0.25}, {n,
0.75}}, t]


For starting values of the parameters $$a$$, $$c$$, $$d$$, and $$n$$, I used the final results in the @ydd's answer which give the values of these parameters. But still, using these values, this approach doesn't work.

Any help is appreciated!

I split the data values and uncertanties into two separate lists, and also grab the time values:

sVals = SdatawithSD[[All, 2, 1]];
bVals = BdatawithSD[[All, 2, 1]];

sErr = SdatawithSD[[All, 2, 2]];
bErr = BdatawithSD[[All, 2, 2]];

tPoints = SdatawithSD[[All, 1]];

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



I tried using MultiNonlinearModelFit, but the solutions were exploding off to massive numbers when I did that. I tried adding constraints as shown in the documentation but this just produced an error. So in the meantime, I will use my previous method of using a linear interpolation of the data, plugging the interpolation into the system of equations, and finding the parameters that minimize the squared difference between the lhs and rhs of each equation in the system (The first equation determines a, then the value of everything else is determined from the second equation). I randomly sample the data each time assuming a normal distribution with the same mean and standard deviation as the data:

sDists = MapThread[NormalDistribution[#1, #2] &, {sVals, sErr}];
bDists = MapThread[NormalDistribution[#1, #2] &, {bVals, bErr}];
nps = 100;

SeedRandom[1234];

AbsoluteTiming[
randParamsFitted =
Table[
{randS, randB} = {Thread[{tPoints, RandomVariate /@ sDists}],
initRule = {s0 -> randS[[1, 2]], b0 -> randB[[1, 2]]};
ClearAll[fitB0, fitS0, fitS, fitB];

(*{fitS0[t_],fitB0[t_]}={Fit[randS,powers,t],Fit[randB,powers,
t]};*)

{fitS0, fitB0} = {Interpolation[randS, InterpolationOrder -> 1],
Interpolation[randB, InterpolationOrder -> 1]};
fitB[t_] := Clip[fitB0[t], {0, 1}];
fitS[t_] := Clip[fitS0[t], {0, Infinity}];

eqsWithFit = eqs /. {S -> fitS, B -> fitB};

residatT = MapApply[(#1 - #2)^2 &, eqsWithFit];
aEq = residatT[[1]];
restEq = residatT[[2]];
params = {a, c, d, n};
restParams = Rest@params;

objA = Sum[aEq, {t, tPoints}];
aMin = NArgMin[{objA, a > 0}, a];

restEq = restEq /. {a -> aMin};
restObj = Sum[Simplify[residatT[[2]]], {t, tPoints}];
restMin =
NArgMin[{restObj, restParams \[VectorGreaterEqual] 0},
restParams];

paramRule = Thread[params -> Join[{aMin}, restMin]];
Join[initRule, paramRule]
, nps];
]


Then we solve the differential system with each parameter set. Some of the solutions produce non-real solutions, or solutions with singularities which I elminate by finding the solutions that have the full domain from t=0 to 24. I then calculate the std dev. of the real solution using Around and plot. The solutions for S are fine, but the fits to B have very large uncertainty relative to the data magnitude:

allSolns =
Table[NDSolveValue[eqs /. i, {S[t], B[t]}, {t, 0, 24}], {i,
randParamsFitted}];

domainList = Table[(#[[0]]["Domain"] & /@ i), {i, allSolns}];
cands = Position[domainList, {{{0., 24.}}, {{0., 24.}}}] // Flatten;
{realS, realB} = Transpose@allSolns[[cands]];

{aroundS, aroundB} = {Table[{t, Around[realS // Re]}, {t, 0, 24}],
Table[{t, Around[realB // Re]}, {t, 0, 24}]};

ListPlot[{aroundS, SdatawithSD}, PlotRange -> All,
Joined -> {True, False},
PlotStyle -> {Blue, Black}, PlotLegends -> {"fit", "data"},
PlotLabel -> "S(t)"]

ListPlot[{aroundB, BdatawithSD}, PlotRange -> All,
Joined -> {True, False}, PlotStyle -> {Orange, Black},
PlotLegends -> {"fit", "data"}, PlotLabel -> "B(t)"]


• The warnings are from NDSolve not being able to have small enough step size near singularities which occur from bad parameters. These are thrown out. It took 5 minutes total to run all the code on my device (Apple M1 Pro Chip on macOS Sonoma 14.1.1, 16 GB RAM)
– ydd
Nov 26, 2023 at 20:48
• Ah nevermind, you can't run it because I forgot to include the definitions of sDists and bDists, it should work now
– ydd
Nov 26, 2023 at 21:09