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?
-
Best Answer
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$.