Need help with solve a system of delay differential equations

Question

I am trying to solve and plot delayed SIRD model. My code is:

S4 = 
  NDSolve[
    {x'[t] == - (0.3/(80*10^6)) x[t] y[t], 
     y'[t] == (0.3/(80*10^6)) x[t] y[t] - 
                 (0.3/(80*10^6)) x[t - 14] y[t - 14] - 0.01 y[t], 
     z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - 14], a'[t] == 0.01 y[t],  
     x[0] == 80*10^6, y[t /; t <= 0] == E^t, y[0] == 150, z[0] == 0, 
     a[0] == 0 }, {x, y, z, a}, {t, 0, 200}]`

In the solution I am getting the error message:

NDSolve::ndsz: At t == 164.8072392646089, step size is effectively zero; singularity or stiff system suspected.

And this results in values of graph that shouldnt be possible, active cases can't be negative.

What am I doing wrong?

Answer 1

I think there's a subtle inconsistency in your equations. You have to discount the recovery rate by the survival probability of infecteds, otherwise you double-count removals from y (both due to death and recovery), hence the negative y[t]. Luckily that is easy to do because survivorship is an exponential of the death rate and the infection time:

S4 = NDSolve[{
   x'[t] == -(0.3/(80*10^6)) x[t] y[t], 
   y'[t] == (0.3/(80*10^6)) x[t] y[t] - (0.3/(80*10^6)) x[t - 14] y[t - 14] E^(-0.01*14) - 0.01 y[t], 
   z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - 14] E^(-0.01*14), 
   a'[t] == 0.01 y[t], x[0] == 80*10^6, y[t /; t <= 0] == E^t, 
   y[0] == 150, z[0] == 0, a[0] == 0}, {x, y, z, a}, {t, 0, 200}];

Plot[Evaluate[{x[t], y[t], z[t], a[t]} /. S4], {t, 0, 200}]

Answer 2

Remove y[0] == 150and set $\tau=2$ and replace $y(t \cdots)$ by $y(t-\tau \cdots)$ can remove some warm message.

τ = 2; 
S4 = 
 NDSolve[{x'[t] == -(0.3/(80*10^6)) x[t] y[t - τ], 
   y'[t] == (0.3/(80*10^6)) x[t] y[t - τ] - (0.3/(80*10^6)) x[
       t - 14] y[t - τ - 14] - 0.01 y[t - τ], 
   z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - τ - 14], 
   a'[t] == 0.01 y[t - τ], x[0] == 80*10^6, 
   y[t /; t <= 0] == E^t, z[0] == 0, a[0] == 0}, {x, y, z, a}, {t, 0, 
   200}]
Plot[{x[t], y[t], z[t], a[t]} /. S4 // Evaluate, {t, 0, 200}, 
 AxesOrigin -> {0, 0}, PlotRange -> All, AspectRatio -> 1]