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
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 [Reference Alonso, McKane and Pascual1–Reference Greenman, Kamo and Boots5]. 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).
Ultimately, the basic SIR model with constant parameters is unable to reproduce repeated epidemic outbreaks [Reference Dietz3]. In this paper, we examine several mechanisms which are able to reproduce multiple outbreaks, including both periodic and chaotic outbreaks. These are:
-
(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 [Reference Bertozzi, Franco, Mohler, Short and Sledge2, Reference Hethcote6–Reference Kolokolnikov and Iron10].
-
(b) Introduction of sufficient noise into the system or stochastically driven parameters [Reference Kolokolnikov and Iron10–Reference Schwartz12].
-
(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.
-
(d) Multiple components for the recovered class and reinfection from the recovered class [Reference Bertozzi, Franco, Mohler, Short and Sledge2, Reference Stollenwerk, Spaziani, Mar, Arrizabalaga, Knopoff, Cusimano, Anam, Shrivastava and Aguiar13].
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 [Reference Bertozzi, Franco, Mohler, Short and Sledge2], 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 [Reference Horrocks and Bauch8] considered sinusoidal forcing for a similar model in the high-frequency regime (large $\omega$ ) and found both periodic and chaotic outbreaks. In paper [Reference Schwartz12], 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 [Reference Kermack and McKendrick9] 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
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
First consider the case of constant $\beta,\gamma .$ The endemic equilibrium is then given by
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
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
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
In addition, the endemic equilibrium has $I_{eq}=O(\varepsilon )$ so we also rescale
to obtain the system
Assuming that the system remains close to the endemic equilibrium, we expand in powers of $\varepsilon$ as
At leading order we then obtain
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:
A simple Maple computation, eliminating $y$ from (2.8), yields the following cubic equation in $a^{2}$ for the threshold values of the parameters:
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’.
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:
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:
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.
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
Then to leading order, we have the following governing equations for the ‘inner’ region:
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
and correspondingly, we obtain following relation between $S_{i}^{\pm }:$
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
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}:$
Then, $y_{i},z=O(1)$ . Equations (1.1) then become
where
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
having solutions
Substituting (3.5) into (3.3) yields the equation for $I$ in the outer region:
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}$ :
where
Finally, from (3.5) we obtain
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 .$
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
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 :$
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:
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
where $S^{\pm }$ are as given by (4.3). Linearising the return map we obtain, after some algebra,
where
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:
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.
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
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
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.
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:
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.
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:
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 [Reference Stollenwerk, Spaziani, Mar, Arrizabalaga, Knopoff, Cusimano, Anam, Shrivastava and Aguiar13]. It would be interesting to apply some of the ideas here to the resulting continuum limit to predict destabilisation and outbreaks in this model.
Competing interests
None.