Recurrent and chaotic outbreaks in SIR model

Abstract

We examine several extensions to the basic Susceptible-Infected-Recovered model, which are able to induce recurrent outbreaks (the basic Susceptible-Infected-Recovered model by itself does not exhibit recurrent outbreaks). We first analyse how slow seasonal variations can destabilise the endemic equilibrium, leading to recurrent outbreaks. In the limit of slow immunity loss, we derive asymptotic thresholds that characterise this transition. In the outbreak regime, we use asymptotic matching to obtain a two-dimensional discrete map which describes outbreak times and strength. We then analyse the resulting map using linear stability and numerics. As the frequency of forcing is increased, the map exhibits a period-doubling route to chaos which alternates with periodic outbreaks of increasing frequency. Other extensions that can lead to recurrent outbreaks include the addition of noise, state-dependent variation and fine-graining of model classes.

1. Introduction

An SIRS model (Susceptible-Infected-Recovered) is the simplest epidemic model describing the transmission of an infectious disease within a population and is the basic building block of mathematical epidemiology. In its most simple form, the basic SIR model is

(1.1) \begin{align} S_{t} & =-\beta SI+\varepsilon R\nonumber \\ I_{t} & =\beta SI-\gamma I\\ R_{t} & =\gamma I-\varepsilon R.\nonumber \end{align}

Here, $\beta$ is the transmission rate, $\gamma$ is the rate of recovery from infection, and $\varepsilon$ is the rate of loss of immunity. Typically, the loss of immunity happens on a much slower timescale (e.g. years) than infection (e.g. days), and for this reason, we are primarily interested in the regime of small $\varepsilon .$

The behaviour of the basic SIR model (1.1) – with constant parameters $\beta,\varepsilon,\gamma$ – is well understood. Despite its simplicity, the basic SIR model and its variants are able to reproduce the progression of epidemic, at least on relatively short timescales [1–5]. On short timescales, the loss of immunity is often discarded (i.e. $\varepsilon =0$ ); the result is typically an epidemic outbreak, where most (but not all) of susceptibles become infected, followed by a recovery period. When loss of immunity is allowed $\left ( \varepsilon \gt 0\right )$ , several outbreaks are possible (especially for very small $\varepsilon$ ), but the dynamics eventually converge to a stable endemic equilibrium, as illustrated in Figure 1(row 1).

Figure 1. Simulations of (2.1) with $\beta (t)=\beta _{0}+a\cos (\omega \varepsilon t)$ for several values of $a$ as shown. Other parameters are $\beta _{0}=1,\omega =2\pi,$ $\gamma =0.5,\ \varepsilon =0.001.$ Initial conditions were taken to be $I(0)=0.01,$ $S(0)=0.99.$ Row 1: constant $\beta$ , damped oscillations towards steady state are observed. Row 2: Convergence towards an adiabatic state. Row 3: Convergence towards adiabatic equilibrium with damped oscillations. Row 4: Onset of small periodic oscillations. Row 5: Large periodic outbreaks. $J_{0}$ is defined in (2.7).

Ultimately, the basic SIR model with constant parameters is unable to reproduce repeated epidemic outbreaks [3]. In this paper, we examine several mechanisms which are able to reproduce multiple outbreaks, including both periodic and chaotic outbreaks. These are:

1. (a) Seasonal variation, such as time-periodic infection rate $\beta =\beta (t)$ . It has been well documented that infectious diseases often have infection rates $\beta \!\left ( t\right )$ which are correlated to the seasons, and adding this variation can lead to very complex outbreak dynamics [2, 6–10].

2. (b) Introduction of sufficient noise into the system or stochastically driven parameters [10–12].

3. (c) State-dependent parameters, such as an infection rate which is suppressed in the presence of infection (e.g. through public health measures), but is relaxed back to a natural rate $\beta _{0}$ after the outbreak passes.

4. (d) Multiple components for the recovered class and reinfection from the recovered class [2, 13].

In this paper, we mostly concentrate on seasonal parameter variations (a). The study of periodic outbreaks has a long history; here we mention but a few references. In [2], London and Yorke estimated a seasonally varying contact rate $\beta (t)$ to simulate annual outbreaks of chickenpox and mumps and biennial outbreaks of measles. Paper [8] considered sinusoidal forcing for a similar model in the high-frequency regime (large $\omega$ ) and found both periodic and chaotic outbreaks. In paper [12], the authors use data-fitting algorithms to ‘discover’ time-periodic parameters which best fit the incidence data for Measles, Varicella and Rubella. A recent paper [9] studied the periodic solution using asymptotic methods. We will be using a similar approach in Section 3.

Let us summarise the results. Consider sinusoidal infection rate of the form

(1.2) \begin{equation} \beta (t)=\beta _{0}+a\cos \!\left ( \omega \varepsilon t\right ) ; \end{equation}

this is meant to model a seasonal variation in the infection rate. The scaling $\omega \varepsilon t$ is chosen so that timescale of seasonal change is comparable with the timescale of immunity loss $\varepsilon t.$ Figure 1 shows what happens when $a$ is used as a control parameter. When $a=0$ , the endemic equilibrium state is stable. In fact as we will see in Section 2, slowly decaying oscillations converge to endemic equilibrium in the limit $\varepsilon \rightarrow 0.$ Increasing $a$ may eventually destabilise the equilibrium. A simple criterion for the threshold value at which this destabilisation occurs is derived in Section 2 (see (2.9)). Instability of endemic state results in a wave of infections as illustrated in Figure 1(bottom). In Section 3, we use asymptotic matching to obtain a two-dimensional discrete map which describes outbreak times and strength. In Section 4, we analyse the resulting discrete map using linear stability and numerics. As the frequency of forcing is increased, the map exhibits a period-doubling route to chaos which alternates with periodic outbreaks of increasing frequency. In Section 6, we use numerics to illustrate how extensions (b-d) can also lead to complex dynamics. Some open questions are discussed in Section 7.

2. Onset of recurrent outbreaks with periodic $\beta$

As shown in Figure 1, large periodic outbreaks appear when $\beta$ has ‘sufficient variation’ (row 4). This is in contrast to the ‘slow variation’ of the adiabatic steady state (time-dependent equilibrium) as illustrated in Figure 1 (row 2). Our first task is to quantify this dichotomy.

We start with a standard simplification of the system (1.1) by noting that the total population $S+I+R$ is conserved in time. For convenience, we rescale the total population to one, so that $S+I+R=1.$ Eliminating $R$ from (1.1), we then obtain an equivalent two-dimensional system

(2.1) \begin{align}S_{t} & =-\beta SI+\varepsilon \!\left ( 1-S-I\right ),\nonumber\\[4pt]I_{t} & =(\beta S-\gamma )I. \end{align}

First consider the case of constant $\beta,\gamma .$ The endemic equilibrium is then given by

(2.2) \begin{equation} S_{eq}=\frac{\gamma }{\beta };\text{ }I_{eq}=\frac{\varepsilon }{\beta }\frac{\left ( \beta -\gamma \right ) }{\gamma +\varepsilon }. \end{equation}

Here and below, we will assume that $\beta \gt \gamma ;$ on the contrary case, there is no (positive) endemic equilibrium and the system converges to the disease-free equilibrium $I=0,\ S=1.$

Linearisation around the endemic equilibrium yields the Jacobian matrix

\begin{equation*} M=\left [ \begin {array}{c@{\quad}c}-\beta I_{eq}-\varepsilon & -\gamma -\varepsilon \\[5pt] \beta I_{eq} & 0 \end {array} \right ]. \end{equation*}

Since $tr(M)=-\varepsilon \frac{\beta +\varepsilon }{\gamma +\varepsilon }\lt 0$ and $\det \!\left ( M\right ) =\varepsilon \!\left ( \beta -\gamma \right ) \gt 0$ , the endemic equilibrium is stable for all parameter values. Nonetheless, a sequence of decaying outbreaks is apparent in Figure 1 (row 1). This is due to the fact that for small $\varepsilon$ , the eigenvalues of $M$ are approximately given by

\begin{equation*} \lambda \sim \sqrt {\varepsilon }\left ( \pm i\sqrt {\beta -\gamma }-\frac {\sqrt {\varepsilon }}{2}\frac {\beta }{\gamma }\right ),\ \ \ \varepsilon \ll 1. \end{equation*}

and are dominated by their purely imaginary part, with a small negative real part. This explains the slowly decaying oscillations towards the endemic equilibrium in this regime.

To explain the transition towards periodic outbreaks, we consider the slow time

(2.3) \begin{equation} y=\varepsilon t. \end{equation}

In addition, the endemic equilibrium has $I_{eq}=O(\varepsilon )$ so we also rescale

\begin{equation*} I=\varepsilon J \end{equation*}

to obtain the system

(2.4) \begin{align} S_{y} & =({-}\beta J-1)S+1-\varepsilon J\nonumber\\[5pt] \varepsilon J_{y} & =\left ( \beta S-\gamma \right ) J \end{align}

Assuming that the system remains close to the endemic equilibrium, we expand in powers of $\varepsilon$ as

(2.5) \begin{equation} J=J_{0}+\varepsilon J_{1}+\ldots ;\ \ \ S=S_{0}+\varepsilon S_{1}+\ldots. \end{equation}

At leading order we then obtain

(2.6) \begin{align} S_{0} & =\gamma/\beta ; \qquad\qquad\qquad\qquad\end{align}

(2.7) \begin{align} J_{0} & =\frac{1-S_{0}-S_{0}^{\prime }}{\beta _{0}S_{0}}=\frac{1}{\gamma }-\frac{1}{\beta }+\frac{\beta ^{\prime }}{\beta ^{2}}. \end{align}

It is here that the effect of the slow adiabatic change becomes apparent: the first two terms of $J_{0}$ are the same as the time-independent equilibrium (2.2) but the term $\beta ^{\prime }$ can make $J_{0}$ negative. When this happens, the adiabatic approximation breaks down and instability ensues, leading to periodic outbreaks. The onset of periodic outbreaks therefore corresponds to the double-root condition:

(2.8) \begin{equation} J_{0}(y)=0=J_{0}^{\prime }(y). \end{equation}

A simple Maple computation, eliminating $y$ from (2.8), yields the following cubic equation in $a^{2}$ for the threshold values of the parameters:

(2.9) \begin{align} 0 & =16\omega ^{2}a^{6}+\left ( \left ( 8\omega ^{4}-20\omega ^{2}-1\right ) \gamma ^{2}+48\beta _{0}\omega ^{2}\gamma -48\beta _{0}^{2}\omega ^{2}\right ) a^{4}+ \left ( \left ( \omega ^{2}+1\right ) ^{3}\gamma ^{4}\right.\nonumber \\ & \quad \left. -2\!\left ( 10\omega ^{2}+1\right ) \left ( \omega ^{2}+1\right ) \beta _{0}\gamma ^{3}+2\beta _{0}^{2}\!\left ( 10\omega ^{4}+35\omega ^{2}+1\right ) \gamma ^{2}-96\beta _{0}^{3}\omega ^{2}\gamma +48\beta _{0}^{4}\omega ^{2}\right ) a^{2}\\ & \quad -\left ( \omega ^{2}+1\right ) ^{2}\beta _{0}^{2}\gamma ^{4}+2\beta _{0}^{3}\!\left ( \omega ^{4}+10\omega ^{2}+1\right ) \gamma ^{3}-\beta _{0}^{4}\!\left ( \omega ^{4}+50\omega ^{2}+1\right ) \gamma ^{2}+48\beta _{0}^{5}\omega ^{2}\gamma -16\beta _{0}^{6}\omega ^{2}.\nonumber \end{align}

For example when $\beta _{0}=1,\,\gamma =0.5,$ and $\omega =2\pi,$ (2.9) yields $a=0.1446.$ This is in agreement with the observed values of the onset of the periodic outbreaks in Figure 1 with $\varepsilon =0.001$

Figure 2(a) shows the comparison between the threshold boundary (2.9) and full numerics. Here, we took $\beta _{0}=1,\gamma =0.5.$ and plotted the implicit curve (2.9) in the $a-\omega$ plane (shown in dashed line). To compare with full numerics, for each choice of $a=0,0.01, 0.02, \ldots, 0.5$ and $\omega =0, 0.4, 0.8, \ldots, 20$ , the solution was computed up to time $t=30,000$ . The solution was deemed to have an outbreak if $\frac{\mbox{max}I(t)}{\mbox{min}I(t)} \gt 1000$ . The choice of ‘1000’ is rather arbitrary and the overall picture barely changes if it is replaced by ‘10,000’ or ‘100’.

Figure 2. Comparison of asymptotic prediction for the outbreak boundary and full numerics. (a) Periodic $\beta (t).$ Dashed curve shows the threshold boundary in the $a-\omega$ parameter space given by (2.9), with $\beta =1+a\cos (\omega \varepsilon t),\gamma =0.5.$ Orange and blue regions correspond to outbreak (unstable) and endemic (stable) regions, respectively, as computed using the full numerical simulations of the full model (1.1) with $\varepsilon =0.001.$ See text for details. (b) Periodic $\gamma (t)$ . Dashed curve shows the threshold boundary as given by (2.11) with $\beta =1,\ \ \gamma =0.5+a\cos (\omega \varepsilon t) .$ Full numerics are as in part (a), computed with $\varepsilon =0.001$ .

Overall, the agreement looks very good, although the asymptotic prediction is degraded for sufficiently large $\omega$ and small $a$ (top left corner of Figure 2). A different asymptotic scaling may be required to more accurately capture the numerics in that regime; we leave it as an open problem for future work.

The method for computing this threshold is rather general. For example, a simpler computation is to take $\beta$ constant and vary $\gamma$ periodically:

(2.10) \begin{equation} \beta (t)=\beta,\ \ \ \gamma (t)=\gamma _{0}+a\cos \!\left ( \omega \varepsilon t\right ). \end{equation}

The threshold is then obtained as previously, by setting $J_{0}=0=J_{0}^{\prime }$ and eliminating $y.$ In this case, the threshold is given by a much simpler formula:

(2.11) \begin{equation} a^{2}\!\left ( 1+\omega ^{2}\right ) =\left ( \beta -\gamma _{0}\right ) ^{2} \end{equation}

as shown in Figure 2(b).

3. Return map

Beyond the onset threshold described in Section 2, recurrent outbreaks are observed. These outbreaks are ‘localised’ in time, followed by a long relaxation period. They may be periodic (commensurate with periodic variation in $\beta$ ) or chaotic, depending on parameter values. Such separation of scales allows for matched asymptotic analysis in order to reduce the system to a discrete ‘return map’ describing the timing and size of the subsequent outbreaks.

Figure 3 shows typical ‘outbreak’ dynamics and separation of scales. The dynamics have two timescales: the outbreak (‘inner’ region), characterised by relatively large values of infection $I$ and a sharp drop in $S,$ followed by the ‘outer’ region, characterised by exponentially small values of $I$ and a gradual increase of $S$ due to a slow immunity loss. The inner region occurs at $O(1)$ timescale whereas the outer region occurs on $O(1/\varepsilon )$ timescale. Let $t_{i}$ be the time of outbreak $i,$ and define $y_{i}=\varepsilon t_{i}.$ Let $S_{i}^{+}$ (resp. $S_{i}^{-}$ ) be the susceptible population at the onset (resp. after) of outbreak $i$ . Refer to Figure 3. We now construct both inner and outer regions, then glue them together to obtain the desired discrete return map.

Figure 3. Inner-outer structure of the recurrent outbreaks solution and definition of various quantities that define the discrete return map.

Inner region. Assuming that $\varepsilon$ is very small, on the scale of the outbreak, the immunity loss terms $\varepsilon R$ can be disregarded. Therefore during the time of the outbreak $i$ , the ‘inner region’ becomes just the simplest SI model. Namely, we let

\begin{equation*} t=t_{i}+\tau,\ \ \ \beta _{i}\,:\!=\,\beta (t_{i}); \end{equation*}

Then to leading order, we have the following governing equations for the ‘inner’ region:

\begin{align*} \frac{dS}{d\tau } & =-\beta _{i}SI,\ \ \ \frac{dI}{d\tau }=\beta _{i}SI-\gamma I\\ \lim _{\tau \rightarrow \pm \infty }S(\tau ) & \rightarrow S_{i}^{\mp };\ \ \lim _{\tau \rightarrow \pm \infty }I(\tau )=0 \end{align*}

This is the usual SI model, where the recovered class is removed from the system. We then have $\frac{dI}{dS}=-1+\frac{\gamma }{\beta _{i}}\frac{1}{S}$ which gives the first integral

\begin{equation*} I=S_{i}^{+}-S+\frac {\gamma }{\beta _{i}}\log \frac {S}{S_{i}^{+}}\end{equation*}

and correspondingly, we obtain following relation between $S_{i}^{\pm }:$

(3.1) \begin{equation} 0=S_{i}^{+}-S_{i}^{-}+\frac{\gamma }{\beta _{i}}\log \frac{S_{i}^{-}}{S_{i}^{+}}. \end{equation}

Outer region. Next, we look at the ‘outer’ region $t\in \left ( t_{i},t_{i+1}\right )$ , characterised by exponentially small amount of infection $I\ll 1.$ The extent of this latent region is of $O(1/\varepsilon ).$ It is therefore appropriate to use the scaling

\begin{equation*} y_{i}=\varepsilon t_{i}\end{equation*}

and furthermore, we define a variable $z$ to be an ‘outer’ time variable shifted so that ‘initial time’ $z=0$ corresponds to $t=t_{i}:$

\begin{equation*} t=t_{i}+\varepsilon ^{-1}z. \end{equation*}

Then, $y_{i},z=O(1)$ . Equations (1.1) then become

(3.2) \begin{align} \varepsilon S_{z} & =-\hat{\beta }(z)SI+\varepsilon R \end{align}

(3.3) \begin{align} \varepsilon I_{z} & =\left ( \hat{\beta }(z)S-\gamma \right ) I \end{align}

(3.4) \begin{align} \varepsilon R_{z} & =\gamma I-\varepsilon R \end{align}

where

\begin{equation*} \hat {\beta }\!\left ( z\right ) =\beta (y_{i}+z)=\beta _{0}+a\cos \!\left ( \omega \!\left ( y_{i}+z\right ) \right ). \end{equation*}

In the outer region, we assume that $I$ is exponentially small compared to $S$ and $R$ so that to leading order, (3.2) and (3.4) become

\begin{equation*} R_{z}\sim -R,\ \ \ S_{z}\sim R, \end{equation*}

having solutions

(3.5) \begin{equation} R\sim (1-S_{i}^{-})e^{-z};\ \ \ S=1-R. \end{equation}

Substituting (3.5) into (3.3) yields the equation for $I$ in the outer region:

(3.6) \begin{equation} I\sim C\exp \!\left \{ \frac{1}{\varepsilon }\int\limits _{0}^{z}\left \{ \hat{\beta }(z)\left ( 1-(1-S_{i}^{-})e^{-z}\right ) -\gamma \right \} dz\right \} \end{equation}

Matching. The outer region ends (and the next outbreak begins) when $I$ changes from exponentially small to O(1). This yields $y_{i+1}$ as a function of $y_{i}$ :

\begin{equation*} y_{i+1}=y_{i}+d_{i}, \end{equation*}

where

(3.7) \begin{equation} \int\limits_{0}^{d_{i}}\left \{ \beta _{0}+a\cos \!\left ( \omega _{0}y_{i}+\omega _{0}z\right ) \left ( 1-(1-S_{i}^{-})e^{-z}\right ) -\gamma \right \} dz=0. \end{equation}

Finally, from (3.5) we obtain

(3.8) \begin{equation} S_{i+1}^{+}=1-(1-S_{i}^{-})e^{-d_{i}} \end{equation}

Together, equations (3.1) (3.7) and (3.8) define a ‘return map’.

Figure 4 shows the comparison between the full numerical solution of the ODEs and the asymptotics (3.1) (3.7) (3.8). Overall, a good agreement is observed. Note that $S_{i}^{\pm }$ , as defined via the asymptotic discrete return map, is independent of $\varepsilon$ , whereas the ODEs do depend on $\varepsilon .$

Figure 4. Comparison between the full ODE (2.1) simulation (left) and the discrete return map (right). The ode (2.1) was solved with $\beta (t)=1+a\cos (4\pi \varepsilon ),$ $\gamma =0.5,\ \varepsilon =0.001.$ The parameter $a$ was gradually ramped up using the formula $a=0.1+10^{-7}t$ with $t=0\ldots 10^{6}.$ Forward Euler method was used with timestep $dt=0.1$ Maxima of $S(t)$ as plotted versus $a$ . Right: Iterations of the return map (3.1) (3.7) and (3.8). Here, 4000 iterations are plotted, while gradually increasing $a$ according to the formula $a_{i}=0.1+0.1\frac{i}{i_{\max }}$ ; here, $i$ is the iteration number and $i_{\max }=4000.$ .

Figure 5. (a) Transitions between $k=1$ periodic solution ( $\omega =3$ ) and $k=2$ solution $(\omega =4)$ . Parameters are $\beta _{0}=1,\ a=0.5,\ \gamma =0.7,\ \varepsilon =0.001$ , and $\omega$ as indicated $.$ (b) Plot of the shift $s$ versus $\omega .$ Other parameters are $\beta _{0}=1,\ a=0.5,\ \gamma =0.7$ . Red dots denote the locations of the fold points. (c) Bifurcation diagram of the discrete return map, as $\omega$ is slowly ramped up. Dashed curves overlay the location of the fixed points corresponding to $k=1\ldots 4$ from (b) The stable branch is between the red point $\omega _{f}$ and the cyan point $\omega _{c}.$ The solutions in (a) correspond to the values of $\omega$ denoted by vertical lines in (c).

4. Periodic outbreaks and stability

A periodic train of outbreaks such as shown in Figure 5(a) (top or bottom) corresponds to a fixed point of the discrete return map derived in Section 3. Since $\beta (y)=\beta _{0}+a\cos \!\left ( \omega y\right )$ is periodic with frequency $\omega$ , any such periodic solution must have its period commensurate with this frequency. Such solution has the form

\begin{equation*} S_{i}^{\pm }=S^{\pm };\ y_{i}=s+2\pi ki/\omega, \,k\in \mathbb {N} \end{equation*}

where the integer $k$ counts the number of periods of $\beta (y)$ per outbreak (i.e. $k$ ‘years’ per single outbreak), and $s$ is the ‘shift’ between the location of the outbreak and the maximum of $\beta .$ Substituting these conditions into the discrete map yields the following conditions for the fixed point of period $2\pi k/\omega :$

(4.1) \begin{align} 0 & =\int\limits _{0}^{2\pi k/\omega }\left \{ \beta (s+y)\left ( 1-(1-S^{-})e^{-y}\right ) -\gamma \right \} dy; \end{align}

(4.2) \begin{align} S^{+} & =1-(1-S^{-})e^{-2\pi k/\omega } \end{align}

(4.3) \begin{align} 0 & =S^{+}-S^{-}+\frac{\gamma }{\beta (s)}\log \frac{S^{-}}{S^{+}} \end{align}

An example of such a fixed point with $k=1$ and $2$ is shown in Figure 5(a).

We eliminate $S^{\pm }$ from the system (4.1), (4.2) as follows:

(4.4) \begin{equation} S^{-}=1-\frac{\int _{0}^{2\pi k/\omega }\left \{ \beta (s+y)-\gamma \right \} dy}{\int _{0}^{2\pi k/\omega }\beta (s+y)e^{-y}dy};\ \ \ \ S^{+}=1-\frac{\int _{0}^{2\pi k/\omega }\left \{ \beta (s+y)-\gamma \right \} dy}{\int _{0}^{2\pi k/\omega }\beta (s+y)e^{-y}dy}e^{-2\pi k/\omega }. \end{equation}

This results in a single algebraic equation for $s$ (4.3). Solving the resulting equation numerically, we found at most two fixed points for each $k=1,2,\ldots .$ Figure 5(b) shows the bifurcation diagram $\omega$ vs. $s$ for $k=1,\ldots,7.$ For a given $k$ , there is a fold point $\omega _{k}$ such that the two branches exist if and only if $\omega \gt \omega _{k}.$

Note how multiple possible fixed points can coexist for a single value of $\omega ;$ for example there are four when $\omega =5$ (two with $k=1$ and two with $k=2$ ). However, only one solution appears to be stable (the upper branch of $k=2$ when $\omega =5$ ). Figure 5(c) shows the dynamics of the return map as well as the fixed points, as $\omega$ is gradually increased. As $\omega$ is increased, successive steady states $k=1,2,3\ldots$ lose stability, and the solution transitions to the next value of $k.$

5. Stability

Each branch (indexed by $k$ ) in Figure 5(b) starts at the fold point $\omega =\omega _{f}$ where the two solutions in the $s-\omega$ plane meet. This is indicated by red circles. As $\omega$ is increased, eventually each branch loses its stability at some $\omega =\omega _{c}$ , indicated by cyan circle in the figure. To compute $\omega _{c}$ , we linearise around periodic states (indexed by $k$ ) of the return map. Let

\begin{equation*} y_{i}=s+2\pi ki/\omega +\xi _{i},\ \ S_{i}^{\pm }=S^{\pm }+\phi _{i}^{\pm }, \end{equation*}

where $S^{\pm }$ are as given by (4.3). Linearising the return map we obtain, after some algebra,

\begin{equation*} \left ( \begin {array}{c}\xi _{i+1}\\[3pt] \phi _{i+1}^{+}\end {array} \right ) =\left ( \begin{array}{c@{\quad}c}A & B\\[3pt] C & D \end {array} \right ) \left ( \begin{array}{c}\xi _{i}\\[3pt] \phi _{i}^{+}\end {array} \right ), \end{equation*}

where

(5.1) \begin{align} A & =1-\frac{\int _{0}^{\frac{2\pi k}{w}}\beta ^{\prime }(s+z)\left ( 1-(1-S^{-})e^{-z}\right ) dz}{\beta (s)\left ( 1-(1-S^{-})e^{-\frac{2\pi k}{w}}\right ) -\gamma }+\frac{\int _{0}^{\frac{2\pi k}{w}}\beta (s+z)e^{-z}dz}{\beta (s)\left ( 1-(1-S^{-})e^{-\frac{2\pi k}{w}}\right ) -\gamma }\frac{\gamma \beta ^{\prime }(s)\log \!\left ( \frac{S^{-}}{S^{+}}\right ) }{\beta ^{2}(s)-\frac{\gamma \beta (s)}{S^{-}}}, \end{align}

(5.2) \begin{align} B & =-\frac{\int _{0}^{\frac{2\pi k}{w}}\beta (s+z)e^{-z}dz}{\beta (s)\left ( 1-(1-S^{-})e^{-\frac{2\pi k}{w}}\right ) -\gamma }\frac{\beta (s)-\frac{\gamma }{S^{+}}}{\beta (s)-\frac{\gamma }{S^{-}}},\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \end{align}

(5.3) \begin{align} C & =e^{-\frac{2\pi k}{w}}\left ( S^{-}-1-\frac{\gamma \beta ^{\prime }(s)\log \!\left ( \frac{S^{-}}{S^{+}}\right ) }{\beta ^{2}(s)-\frac{\gamma \beta (s)}{S^{-}}}\right ) -e^{-\frac{2\pi k}{w}}(S^{-}-1)A,\qquad\qquad\qquad\qquad\qquad\qquad\qquad \end{align}

(5.4) \begin{align} D & =e^{-\frac{2\pi k}{w}}\frac{\beta (s)-\frac{\gamma }{S^{+}}}{\beta (s)-\frac{\gamma }{S^{-}}}\left \{ 1-\frac{\int _{0}^{2\pi k/\omega }\left \{ \beta (s+z)-z\right \} dz}{\beta (s)\left ( 1-(1-S^{-})e^{-\frac{2\pi k}{w}}\right ) -\gamma }\right \} .\qquad\qquad\qquad\qquad\qquad\qquad\quad \end{align}

The eigenvalues $\lambda _{\pm }$ of the matrix $\left ( \begin{array}{c@{\quad}c}A & B\\[4pt] C & D \end{array} \right )$ determine the stability of the system. A steady state is stable if $\left \vert \lambda _{\pm }\right \vert \lt 1.$ At the threshold $\omega =\omega _{c},$ the Hopf bifurcation is observed corresponding to the point in parameter space where $\left \vert \lambda _{\pm }\right \vert =1.$

Figure 5(c) shows that $\omega _{c}$ for branch $k$ occurs before $\omega _{f}$ for branch $k+1.$ So while multiple equilibria coexist (corresponding to different values of $k$ ), it appears that at most one equilibrium is stable.

6. Other models with periodic and/or chaotic outbreaks

6.1. Stochastically induced outbreaks

Instead of periodically varying parameters, recurrent outbreaks can be triggered simply when there is enough noise in the system. The precise nature of noise is not as important as the magnitude of the noise relative to the immunity loss rate $\varepsilon$ . Every biological system has inherent noise; the question is how much noise is enough to cause the outbreak. To be concrete, consider the parameter $\beta (t)$ which evolves according to an Ornstein–Uhlenbeck process:

(6.1) \begin{equation} \frac{d\beta }{dt}=a\!\left ( \beta _{0}-\beta \right ) +b\frac{dW}{dt}. \end{equation}

Here, $dW=\xi \sqrt{dt}$ where $\xi$ is the iid normal distribution; $a$ is the strength of attraction of $\beta$ towards the ‘base’ state $\beta _{0}$ and $b$ is the standard deviation of the Brownian drift. Figure 6 shows simulations of (2.1), (6.1) for several values of $\varepsilon,$ with other parameters as indicated. For each run, we used the same trajectory for $\beta$ (which here is independent of $\varepsilon$ ). For larger values of $\varepsilon,$ the stochastic fluctuations are observed, but overall the endemic state ‘follows’ $\beta (t)$ without producing big outbreaks: the blue curve follows the black. However for smaller $\varepsilon,$ the equilibrium state has difficulty ‘catching up’ to the fluctuations in $\beta,$ and this tends to produce large recurrent outbreaks. Overall, it is the ratio of noise level $b$ compared to immunity loss rate $\varepsilon$ that determines whether recurrent outbreaks appear. The analysis of this phenomenon (how much noise is needed to trigger outbreaks) is left for future work.

Figure 6. Effect of stochasticity on recurrent outbreaks. Simulation of (2.1) and (6.1) with $\beta _{0}=1,\gamma =0.5,$ $a=0.001,\ b=0.01$ and with $\varepsilon$ as indicated. Same stochastic path for $\beta$ was used in all three simulations. Top: the endemic state $S$ ‘follows’ $\gamma/\beta$ without major outbreaks. Middle: decreasing $\varepsilon$ induces some outbreaks. Bottom: even smaller $\varepsilon$ results in recurrent outbreaks.

6.2. Public health measures

Consider a model where the infection rate $\beta$ is temporarily adjusted in response to an outbreak. This can be due to, e.g., stay-at-home measures, social distancing etc. These restrictions are eventually lifted when the outbreak passes. For concreteness, we consider $\beta$ to depend on $S$ via an ODE

(6.2) \begin{equation} \frac{d\beta }{dt}=a\!\left ( F(S)-\beta \right ) \end{equation}

where $F(S)$ is a sharp transition function such that $F(S)\sim \beta ^{+}$ for $S\gt S_{0}$ and $F(S)\sim \beta ^{-}$ for $S\lt S_{0}.$ More concretely, take

(6.3) \begin{equation} \gamma =0.5,\ \varepsilon =0.01,\ \ \ F(S)=\frac{1}{2}+\frac{1}{2}\tanh (10\!\left ( S-0.7\right ) )),\ \ \ a=1. \end{equation}

Then, $\beta ^{+}=1,\beta ^{-}=0.$ This control means that $\beta$ quickly approaches $0$ when 30% of the population becomes infected, but goes back to 1 as epidemic passes. Figure 7 shows the resulting dynamics with periodic outbreaks.

Figure 7. Periodic outbreaks for model (6.2) using parameters (6.3).

More complex models can be considered where $\beta$ is controlled by $I$ and/or $R,$ but generally similar phenomenon results. While this model is able to generate periodic outbreaks, we were unable to get chaotic outbreaks, even after playing with more complex models of this type.

6.3. Multiple recovered classes

Recurrent infections can occur if we refine either infection and/or recovery class. Consider the following model with $n$ recovery classes:

(6.4) \begin{equation} \left \{ \begin{array}{l}\dfrac{dS}{dt}=-\beta _{s}IS+\mu _{n}R_{n}\\[12pt] \dfrac{dI}{dt}=\beta _{s}IS+\sum \beta _{i}IR_{i}-\gamma I\\[12pt] \dfrac{dR_{1}}{dt}=-\beta _{1}IR_{1}+\gamma I-\mu _{1}R_{1}\\[12pt] \dfrac{dR_{i}}{dt}=-\beta _{i}IR_{i}+\mu _{i-1}R_{i-1}-\mu _{i}R_{i},\ \ \ i=2\ldots n \end{array} \right. \end{equation}

In this model, there are multiple recovered classes representing a more granular view of loss of immunity with time. Reinfection can also occur among the ‘recovered’ classes with different reinfection rates $\beta _{i}.$ When $n=1$ , this model reduces to a 3-component SIR model with reinfection rate $\beta _{1}$ ; such model is too simple to exhibit oscillatory behaviour (similar to analysis in Section 2, it is easy to show that when $n=1,$ the endemic state is always stable when the disease-free equilibrium is unstable). However even when $n=2$ , recurrent outbreaks become possible. As an example, take $n=2,$ $(\beta _{1},\beta _{2},\beta _{s})=\left ( 1,0,1\right ),\ \gamma =0.7,$ $\mu _{i}=0.01,\ i=1,2;$ this can be thought of as a delay in the onset of immunity, followed by eventual loss of immunity. Periodic recurrent outbreaks are observed for these parameters, see Figure 8.

Figure 8. (a) Model (6.4) illustrating sustained persistent outbreaks with $n=2, (\beta _{s},\beta _{1},\beta _{2})=$ $(1,1,0),\ \gamma =0.7,$ $\mu _{1}=\mu _{2}=0.01.$ (b) Model (6.4) with $\gamma$ as indicated and other parameters as in (6.5). Note the period-doubling route to chaos.

Period-doubling and chaotic outbreaks can also occur in this model for larger $n$ and non-monotone $\beta _{i}.$ One example is given in Figure 8, using the following parameter values:

(6.5) \begin{equation} n=100;\ \ \beta _{s}=1,\ \beta _{i}=\left \{ \begin{array}{c}1,\ \ 40\leq i\leq 80\\[6pt] 0,\ \ \ \text{otherwise}\end{array} \right. ;\ \ \mu _{i}=0.5. \end{equation}

The non-monotonicity of $\beta _{i}$ in (6.5) is hard to justify biologically; it would be more realistic to assume that $\beta _{i}$ is increasing with $i$ . While it is easy to find sustained oscillations with such monotone $\beta _{i},$ it is an open question to find examples of chaotic behaviour with this monotone $\beta _{i}.$ The full analysis of system (6.4) is also left for future work.

7. Discussion

In this paper, we studied several extensions of SIR model which can lead to recurrent/chaotic outbreaks. We analysed how seasonal variations can destabilise an endemic equilibrium leading to outbreaks. We then used asymptotic matching to reduce complex dynamics to a two-dimensional discrete map for outbreak timing and strength.

Many open questions remain. We assumed slow immunity loss comparable with timescale of seasonal change ( $\omega =O(1)$ in notation of (1.2)). A different regime would be to assume an even slower immunity loss; this would correspond to $\omega \gg 1$ . As seen in Figure 2, asymptotics start to fail in this regime (i.e. fast but slight variations). It would be interesting to develop a theory that is able to capture this regime.

Simulations in §6.1 illustrate that stochastic parameter variation leads to a similar destabilisation of the endemic equilibrium. It is an open question to ‘quantify’ this phenomenon: how much stochasticity is required to induce big outbreaks (as opposed to small variations around endemic equilibrium)?

The model of public health in Section 6.2 leads to periodic outbreaks. However, we were not able to produce chaotic outbreaks using a similar paradigm. More effort would be welcome in this direction.

In Section 6.3, we provided an example of a model with multiple classes that leads to period-doubling and chaotic outbreaks. The example is hard to justify biologically due to non-uniformity in $\beta _{i}.$ Finding a more realistic example of this is an open problem. Finally, it would be interesting to consider the continuum limit of many-class model. This naturally leads to integro-differential system [13]. It would be interesting to apply some of the ideas here to the resulting continuum limit to predict destabilisation and outbreaks in this model.