# ParametricNDSolve with a delay differential equation

I have a set of delay differential equations that I solve numerically from 0 < t < T. y[T] is then used as the initial condition for z to start the next round. I would like to keep track of the evolution of y[T] for different values of ta. I planned to do this by generating the list ytotal for different values of ta using Table.

x0 = 10^8; a = 10^-8;  dv = 2; d = 0.5; T = 4; b = 200;
sol = ParametricNDSolve[{
x'[t] == -d x[t] - a *x[t]*z[t],
y'[t] == a b Exp[-d ta]*x[t - ta]*z[t - ta] - dv y[t],
z'[t] == -dv z[t],
x[t /; t <= 0] == x0, y[t /; t <= 0] == 0,
z[t /; t <= 0] == z0},
{x, y, z}, {t, 0, T}, {z0, ta}, MaxStepSize -> 10^1000,
Method -> {"StiffnessSwitching"}, MaxSteps -> 10^6,
WorkingPrecision -> MachinePrecision];
{z0 = 1; ytotal[ta] = {z0};
For[i = 0, i < 100, i++,
{z0 = Evaluate[y[z0, ta][T] /. sol];
ytotal[ta] = Append[ytotal[ta], z0];
zeq[ta] = Last[ytotal[ta]]}]}


My current approach is not working. I thought it might be because of an issue arising from using ParametricNDSolve to vary the delay ta. Also, I am using a For loop which I know is advised against in Mathematica but I couldn't figure out how to iterate sol without one.

• To be sure, you want to solve until t == T for a given ta, then change ta and do it again? Commented Jun 10, 2020 at 15:20
• I want to keep the same ta for all 100 iterates. The goal is to generate lists of y[T] for different values of ta over 100 iterates. Commented Jun 10, 2020 at 15:25
• It seems that restarting the numerical integration with only t == 0 initial conditions loses the history between iterates then. Do you even need to stop & restart NDSolve? Commented Jun 10, 2020 at 15:32
• Yes, I think I do need to restart numerical integration because I need to reset y[0] == 0 and z[0] == z0. (z0 = y[T]) Commented Jun 10, 2020 at 15:36
• Hmm, it's more complex than I thought! Could you add some explanation of what you're modeling? When z0 becomes y[T], should it also inherit y[t]'s history? Commented Jun 10, 2020 at 15:39

I think the easiest way to get what you want is with Table:

z0 = 1; ta = 0.01;
ytotal[ta] = Table[z0 = y[z0, ta][T] /. sol, {i, 100}];


For those parameter values, it looks like the population quickly goes extinct:

ListPlot[ytotal[ta], PlotRange -> All]


If you change parameters they can persist.

I'm not an expert in delay differential equation models, but I'm a bit concerned about how the initial history of z[t] might introduce artifacts in your results. A priori, I imagined that increasing the developmental delay ta would hurt the parasite but instead it seemed to help. But then I figured that for t < ta, y[t] would be increasing due to the a b Exp[-d ta]*x[t - ta]*z[t - ta] term, but were there really any first generation parasites z[t] running around when t < 0? Seems delicate. I tried a different initial history z[t /; t <= 0] == If[t < 0, 0, z0] but that often led to numerical problems (ParametricNDSolve::ndsz -- step size is effectively zero; singularity or stiff system suspected). What do you think?

• I think I should be able to get this work. Thanks for your help! And yes, longer delays help the parasite in this case. Commented Jun 10, 2020 at 21:30
• Glad it helps. I suppose d < dv means the parasite would rather hang out in the host? Commented Jun 10, 2020 at 22:50
• Yes, that's exactly it Commented Jun 10, 2020 at 22:54