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.