# Solving a delay-differential equation related to epidemiology

For some inexplicable reason, I have recently been interested in epidemiology. One of the classical and simplistic models in epidemiology is the SIR model given by the following system of first-order autonomous nonlinear ordinary differential equations in the real variables $$s,i,r$$ satisfying $$s,i,r\geq 0$$ and $$s+i+r=1$$: \begin{align} s' &= -\beta i s\\ i' &= \beta i s - \gamma i\\ r' &= \gamma i \end{align} \tag{*} (where prime denotes derivative w.r.t. time, $$s,i,r$$ represent the proportion of “susceptible”, “infected” and “recovered” individuals, $$\beta$$ is the kinetic constant of infectiousness and $$\gamma$$ that of recovery; the “replication number” here is $$\kappa := \beta/\gamma$$). It is easy to see that for every initial value $$(s_0,i_0,r_0)$$ at $$t=0$$ the system has a unique $$\mathscr{C}^\infty$$ (even real-analytic) solution (say $$r_0=0$$ to simplify, which we can enforce by dividing by $$s_0+i_0$$); it does not appear to be solvable in exact form in function of time, but since $$s'/r' = -\kappa\cdot s$$ and $$i'/r' = \kappa\cdot s - 1$$ (where $$\kappa := \beta/\gamma$$) it is possible to express $$s$$ and $$i$$ as functions of $$r$$, viz., $$s = s_0\,\exp(-\kappa\cdot r)$$ and of course $$i = 1-s-r$$. This makes it possible to express, e.g., the values of $$s,i,r$$ at peak epidemic (when $$i'=0$$): namely, $$s = 1/\kappa$$ and $$r = \log(\kappa)/\kappa$$); or when $$t\to+\infty$$: namely, $$s = -W(-\kappa\,\exp(-\kappa))/\kappa$$ where $$W$$ is (the appropriate branch of) Lambert's transcendental W function; this is upon taking $$i_0$$ infinitesimal and $$r_0=0$$.

Now one of the many ways in which this SIR model is simplistic is that it assumes that recovery follows an exponential process (hence the $$-\gamma i$$ term in $$i'$$) with characteristic time $$1/\gamma$$. A slightly more realistic hypothesis is recovery in constant time $$T$$. This leads to the following delay-difference equation: \begin{align} s' &= -\beta i s\\ i' &= \beta i s - \beta (i s)_T\\ r' &= \beta (i s)_T \end{align} \tag{\dagger} where $$(i s)_T = i_T\cdot s_T$$ and $$f_T(t) = f(t-T)$$. When I speak of a “$$\mathscr{C}^\infty$$ solution” of ($$\dagger$$) on $$[0;+\infty[$$ I mean one where the functions $$s,i,r$$ are $$\mathscr{C}^\infty$$ and satisfy the equations whenever they make sense (so the second and third are only imposed for $$t\geq T$$), although one could also look for solutions on $$\mathbb{R}$$.

I am interested in how ($$\dagger$$) behaves with respect to ($$*$$). Specifically,

• Does ($$\dagger$$) admit a $$\mathscr{C}^\infty$$ (or better, real-analytic) solution on $$[0;+\infty[$$ for every initial value $$(s_0,i_0,r_0)$$ at $$t=0$$? Or maybe even on $$\mathbb{R}$$? Is this solution unique? (Note that we could try to specify a solution by giving initial values on $$[0;T[$$ and working in pieces, but this does not answer my question as it would not, in general, glue at multiples of $$T$$ to give a $$\mathscr{C}^\infty$$ solution.)

• Can we express this solution in closed form or, at least, express $$s$$ and $$i$$ in function of $$r$$?

• Qualitatively, how do the (at least continuous!) solutions of ($$\dagger$$) differ from those of ($$*$$), and specifically:

• How do the behaviors compare when for $$t\to 0$$?

• How do the values compare around peak epidemic $$i'=0$$?

• How do the limits when $$t\to+\infty$$ compare?

• I stumbled upon the SEIR model which differentiates between infected and exposed people - there is only a German Wikipedia page (de.m.wikipedia.org/wiki/SEIR-Modell). The model does not have delay, though.
– Dirk
Commented Apr 1, 2020 at 18:42
• Nice question. Regarding Equation $(*)$, how does the parameters, particularly $\beta$ and $\gamma$ control the "speed" (such as half life) of $I$ approaching the asymptotic state for time infinity? Do you have a reference on this issue?
– Hans
Commented Apr 10, 2020 at 5:27
• +1. There is an exact closed form solution to the SIR model arxiv.org/pdf/1403.2160.pdf.
– Hans
Commented Apr 10, 2020 at 5:53
• I was unable to make your model work (numerically, I haven’t tried your solutions yet). It seems that a regular SIR with a delay in the terms that have gamma works best with actual data ... that is, instead of gamma i, you would have gamma i_T in both the second and third equation of (*). I don’t know how to solve it analytically though. It would be nice. I will keep trying to make your model work ...
– EGME
Commented Apr 19, 2020 at 19:51
• By the way, it might be worth pointing out that there a variety of models called SIR models and that the equations which you describe are specifically known as the Kermack-McKendrick equations. Commented Aug 20, 2021 at 0:07

After trying many Ansätze in the form of series, I stumbled upon the fact (which I find rather remarkable) that ($$\dagger$$) admits an exact solution in closed form. To express this, let me first introduce the following notations:

• $$\kappa := \beta T$$ (reproduction number), which I assume $$>1$$;

• $$\Gamma := -W(-\kappa \exp(-\kappa))/\kappa$$ the solution in $$]0,1[$$ to $$\Gamma = \exp(-\kappa(1-\Gamma))$$, which will be the limit of $$s$$ when $$t\to+\infty$$ (both in ($$*$$) and in ($$\dagger$$), starting with $$i$$ and $$r$$ infinitesimal);

• $$X := \exp(\beta(1-\Gamma) t)$$, a change of variable on the time parameter, in which $$s,i,r$$ will be expressed.

The solution is then given by: \begin{align} s &= \frac{(1-\Gamma)^2+\Gamma c X}{(1-\Gamma)^2+c X}\\ i &= \frac{(1-\Gamma)^4 c X}{((1-\Gamma)^2+c X)((1-\Gamma)^2+\Gamma c X)}\\ r &= \frac{\Gamma (1-\Gamma) c X}{(1-\Gamma)^2+\Gamma c X} \end{align} where c is an arbitrary positive real parameter which merely serves to translate the solution.

From this we can deduce the following about the behavior at start, peak and end of the epidemic, and how it compares to ($$*$$) (for which I let $$\kappa := \beta/\gamma$$ and $$\Gamma$$ defined by the same formula):

• Initially, $$i$$ grows like $$c\,\exp(\beta(1-\Gamma) t)$$ (and $$r$$ like $$\frac{\Gamma}{1-\Gamma}$$ times this). This is in contrast to ($$*$$) where $$i$$ initially grows like $$c\,\exp((\beta-\gamma) t)$$ (and $$r$$ like $$\frac{\gamma}{\beta-\gamma}$$ times this), i.e., for a given reproduction number $$\kappa$$, solutions of ($$\dagger$$) grow faster than those of ($$*$$).

• Peak epidemic happens for $$X = \frac{(1-\Gamma)^2}{c\sqrt{\Gamma}}$$, at which point we have $$s = \sqrt{\Gamma}$$ and $$i = (1-\sqrt{\Gamma})^2$$ and $$r = \sqrt{\Gamma}(1-\sqrt{\Gamma})$$. This is in contrast to ($$*$$) where we have $$s = \frac{1}{\kappa}$$ and $$i = \frac{\kappa-\log\kappa-1}{\kappa}$$ and $$r = \frac{\log\kappa}{\kappa}$$; so, for a given reproduction number $$\kappa$$, solutions of ($$\dagger$$) have a peak epidemic with fewer uninfected ($$s$$) than those of ($$*$$).

• The limit when $$t\to+\infty$$ is the same in ($$\dagger$$) as in ($$*$$), namely $$s \to \Gamma$$.