ParametricNDSolve with a delay differential equation

Question

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.

Answer 1

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?