A simple influenza model with complicated dynamics

Abstract

We propose and analyse a model for the dynamics of a single strain of an influenza-like infection. The model incorporates waning acquired immunity to infection and punctuated antigenic drift of the virus, employing a set of differential equations within a season and a discrete map between seasons. We show that the between-season map displays a variety of qualitatively different dynamics: fixed points, periodic solutions, or more complicated behaviour suggestive of chaos. For some example parameters we demonstrate the existence of two distinct basins of attraction, that is the initial conditions determine the long term dynamics. Our results suggest that there is no reason to expect influenza dynamics to be regular, or to expect past epidemics to give a clear indication of future seasons’ behaviour.

1 Introduction

In mathematical epidemiology, infectious disease dynamics are typically studied using continuous time models, such as the Susceptible–Infectious–Recovered (SIR) model (Diekmann et al. 2013). Given a supply of susceptible hosts, gained through loss of immunity, population turnover or both, SIR models support endemic equilibria. Motivated by the complex interacting cycles of seasonal influenza A epidemics in temperate climates, Andreasen (2003) took the novel step of combining a continuous time SIR model for the ‘within-season’ dynamics, with a discrete map to update the population’s susceptibility between seasons. As is well known for iterative maps, the dynamics can be particularly complicated (Strogatz 2001).

While seasonal influenza epidemics are caused by multiple different influenza viruses (see Kucharski et al. 2016 for a recent review of the modelling efforts in studying multi-strain epidemic dynamics], there is limited cross immunity between these groupings. While cross immunity may act to shorten infection, it does not typically provide protection against infection (Brown and Kelso 2009). In contrast, within a subtype (e.g. influenza A(H3N2)), infection elicits a strong antibody response providing a high level of sterilising immunity. In consequence, infection provides long-lasting protection against homologous re-infection (Brown and Kelso 2009). This strong immunity also provides a selection advantage for antigenically novel viruses, leading to the rapid evolution of the virus through a process referred to as antigenic drift. As a result there is a rapid loss of protection in previously infected hosts. While the tempo and detailed characteristics of antigenic drift differ for the four influenza viruses most commonly associated with infection (A(H3N2), A(H1N1), B(Yam) and B(Vic)), all have this feature, allowing for their continued persistence in the population. Antigenic cartography, first used to study influenza A (H3N2) (Smith et al. 2004), demonstrates that as population-level immunity to a particular strain accumulates, the selection advantage for novel viruses increases dramatically, and the dominant circulating strain abruptly changes. This punctuated drift occurs every 2 to 8 years.

Here, motivated by the characteristic behaviour of seasonal diseases such as influenza, we explore a one-strain model of year-on-year infection dynamics. We exploit the richer dynamical behaviour that can emerge from discrete maps (compared to continuous time dynamical systems). Within-season dynamics are studied through a generalisation of the final-size equation, and between-season dynamics are, as in Andreasen (2003), modelled using a discrete map. Our model incorporates waning immunity. Having been infected in one epidemic season an individual remains fully protected for the next, then becomes susceptible but partially protected for another year, before becoming fully susceptible to reinfection. We also model changes in the host–virus interaction due to punctuated evolution of the virus (Adams and McHardy 2011; Koelle et al. 2006; Smith et al. 2004). In Sect. 2 we introduce the philosophy of the method, basing our presentation on the SIR model without waning immunity, and the assumption that changes in the virus reflected in a loss of host protection have a constant effect over years. In Sect. 3 we introduce the model with waning immunity, and derive semi-analytic expressions for the existence and stability of fixed points of the system. We discuss the within-season dynamics in Sect. 3.1 and the between-season dynamics in Sect. 3.2, presenting and discussing a selection of computational results in Sect. 3.3. We assume throughout that there is one influenza season each year. In Sect. 3 we assume that the between-season changes are the same every year. Then, in Sect. 4 we generalise the model to capture the effects of punctuated drift in which the between-season changes are different in different years. We present results where an inter-season change in susceptibles due to antigenic drift occurs every five years, but with population turnover held fixed each year. We compare these results to those where drift is held fixed each year, and briefly to results where the drift occurs every three or seven years.

2 Preliminary model

To explain the basis for the method, consider the well-known SIR model

$$\begin{aligned} \frac{\text{ d }S}{\text{ d }t} = -\beta S I \qquad \frac{\text{ d }I}{\text{ d }t} = \beta S I -\gamma I \end{aligned}$$

where S is the proportion of the population that is susceptible to infection, and I is the proportion infected and infectious. The basic reproduction number is defined by $\mathcal {R}_0=\beta /\gamma$. If $\mathcal {R}_0S(0)>1$ there is an epidemic, during which the proportion of the population susceptible reduces from $S(0)=S_0$ to $\lim _{t\rightarrow \infty } S(t)=S_\infty$, where $S_0e^{-\mathcal {R}_0S_0} =S_\infty e^{-\mathcal {R}_0S_\infty }$. The proportion of the population infected during the epidemic is $\mathcal {P}=S_0-S_\infty$, which solves the equation $\mathcal {P}=S_0\left( 1-e^{-\mathcal {R}_0\mathcal {P}} \right)$ (Diekmann et al. 2013; Ma and Earn 2006).

We assume a proportion b of the population is replaced with susceptibles each year due to births, deaths, immigration and emigration. This population turnover ensures that following an epidemic year a proportion of those protected from infection in the population are replaced with new susceptibles. If the susceptible population is increased sufficiently, an epidemic in the subsequent year is then a possibility. In addition, the strains of the virus that have potential to cause an epidemic change from year to year (see for example Smith et al. 2004). If a proportion $d\leqslant 1$ of the immune population are protected from the second year’s epidemic virus, and a proportion $1-d$ are not, then we can define a so-called dilution constant $c=\left( 1-b\right) d$. The constant c is the probability that an immune individual remains in the population between seasons, and remains immune (assuming these are independent). Hence $1-c$ is the probability that an individual is replaced, or their immunity wanes. The proportion immune at the end of the season is equal to the proportion immune at the beginning, $1-S_0$, plus the proportion infected during the season, $\mathcal {P}$. Hence, if $S_1$ is the proportion susceptible at the beginning of the next season, then

$$\begin{aligned} S_1=S_\infty +\left( 1-c\right) \left( 1-S_0+\mathcal {P}\right) =1-c+c\left( S_0-\mathcal {P}\right) \end{aligned}$$

(1)

We repeat this process over successive years, defining a between-season map $S_n\rightarrow S_{n+1}$. A fixed point of the map is $S^*=1-c\mathcal {P}^*/\left( 1-c\right)$ where $\mathcal {P}^*$ solves

$$\begin{aligned} \mathcal {P}^*=\frac{\left( 1-c\right) \left( 1-e^{-\mathcal {R}_0\mathcal {P}^*}\right) }{1-ce^{-\mathcal {R}_0\mathcal {P}^*}} \end{aligned}$$

It is shown in “Appendix A” that $0>\text{ d }S_1/\text{ d }S_0>-1$. Hence the fixed point is stable whenever it exists, which is when $\mathcal {R}_0>1$.

3 The model

3.1 Within-season dynamics

Let the proportion of the population that is fully susceptible to infection at time t be $S^{\emptyset }(t)$, the proportion partially susceptible be $S^{1}(t)$, the proportion infected (prevalence of infection) be I(t), and the proportion recovered and immune be R(t). The population size is constant, and $S^{\emptyset }(t)+S^{1}(t)+I(t)+R(t)=1$ for all t. We allow $t\in \left( 0,\infty \right)$ to represent an influenza season. The within-season dynamics are fully described by the equations

$$\begin{aligned} \frac{\text{ d }S^{\emptyset }}{\text{ d }t}&= -\beta S^{\emptyset } I \\ \frac{\text{ d }S^1}{\text{ d }t}&= - k\beta S^1 I \nonumber \\ \frac{\text{ d }I}{\text{ d }t}&= \beta \left( S^{\emptyset }+kS^1\right) I-\gamma I \nonumber \end{aligned}$$

(2)

where $k< 1$ is the relative susceptibility of those in the $S^1$ compartment compared to those in the $S^\emptyset$ compartment. The model is summarised in Fig. 1a. The epidemic takes off if $\mathcal {R}>1$, where $\mathcal {R}=\left( S^{\emptyset }(0) +kS^1(0)\right) \mathcal {R}_0$ and $\mathcal {R}_0=\beta /\gamma$. It follows from Eq. 2 that during an epidemic

$$\begin{aligned} \frac{\text{ d }S^{\emptyset }}{\text{ d }t} + \frac{\text{ d }S^1}{\text{ d }t}+ \frac{\text{ d }I}{\text{ d }t}=\frac{1}{\mathcal {R}_0}\frac{\text{ d }}{\text{ d }t}\log S^{\emptyset }=\frac{1}{k\mathcal {R}_0}\frac{\text{ d }}{\text{ d }t}\log S^1 \end{aligned}$$

(3)

Define $\lim _{t\rightarrow \infty }\left( S^{\emptyset }(t),S^{1}(t),I(t),R(t)\right) =\left( S^{\emptyset }_\infty ,S^{1}_\infty ,0,R_\infty \right)$. Then, integrating Eq. 3 over $\left( 0,\infty \right)$, we obtain

$$\begin{aligned} \mathcal {P}= \frac{1}{\mathcal {R}_0}\log \left( S^{\emptyset }(0)/S^{\emptyset }_\infty \right) = \frac{1}{k\mathcal {R}_0}\log \left( /S^{1}(0)/S^{1}_\infty \right) \end{aligned}$$

where $\mathcal {P}=S^\emptyset (0)+S^1(0)-S^{\emptyset }_\infty -S^{1}_\infty$ is the proportion of the population infected in an epidemic, $S^{\emptyset }_\infty =S^{\emptyset }(0)e^{-\mathcal {R}_0\mathcal {P}}$ and $S^{1}_\infty =S^1(0)e^{-k\mathcal {R}_0\mathcal {P}}$. If $\mathcal {R}<1$, $\mathcal {P}=0$. If $\mathcal {R}>1$, $\mathcal {P}$ solves

$$\begin{aligned} \mathcal {P}= S^{\emptyset }(0)\left( 1-e^{-\mathcal {R}_0\mathcal {P}} \right) +S^1(0)\left( 1-e^{-k\mathcal {R}_0\mathcal {P}} \right) \end{aligned}$$

An example of the within-season model dynamics is shown in Fig. 1b.

Fig. 1

Within-season and between-season dynamics of the model. a Flow between compartments within season. b Numerical example showing proportions in each compartment during an epidemic. Parameter values $\mathcal {R}_0=3.5$ and $k=0.5$. c Summary of the dynamics. Within season those initially susceptible, $S^\emptyset _0$ and $S^1_0$, may be infected during an epidemic, moving to $R_\infty$ (red arrows); or remain uninfected forming $S^\emptyset _\infty$ and $S^1_\infty$ (dashed red arrows). Those initially immune to infection, $R_0$, take no part in an epidemic, hence $R_\infty =\mathcal {P}+R_0$. Between seasons a proportion $1-c$ of those in $R_\infty$ go to $S^\emptyset$ (dashed blue arrows). A proportion c of those in $R_\infty$ remain immune, going to $R_1$ or to $S^1_1$ depending on whether they have just been infected ($\mathcal {P}$) or not ($R_0$) (blue arrows). Those in $S^1_\infty$ go to $S^\emptyset _0$ (black arrow). Note that $R_0$ in b and c denotes the value of R(t) at $t=0$, not the basic reproduction number $\mathcal {R}_0$ (color figure online)

3.2 Between-season dynamics

As we noted in Sect. 2, the combined effect of antigenic drift and population turnover can be mimicked by replacing a proportion $1-c$ of the population with susceptibles between influenza seasons. We now have two susceptible classes, and the initial conditions for the epidemic season are constructed as follows. Those in the R compartment at the beginning of the previous season move to $S^1$, and those infected during the previous season remain in R. Both of these transitions are diluted due to drift and population turnover. Those in the $S^1$ compartment at the beginning of the previous season that did not get infected during that season move to $S^\emptyset$, and those in $S^\emptyset$ remain in $S^\emptyset$. If we use a subscript 0 to denote the initial condition for a season (for example $S^\emptyset _0=S^\emptyset (0)$), the subscript $\infty$ to denote the value of that state variable at the end of the season (as before), and a subscript 1 to denote the initial condition for the subsequent season, then we map the initial conditions from one season to the next by $S^\emptyset _1=1-cR_\infty$, $S^1_1=cR_0$ and $R_1=c\left( R_\infty -R_0\right)$. In other words, if at the beginning of season zero $S^\emptyset (0)=S^\emptyset _0$ and $S^1(0)=S^1_0$, and during season zero a proportion $\mathcal {P}$ of the population is infected, then

$$\begin{aligned} S^\emptyset _1=1-c+c\left( S^\emptyset _0e^{-\mathcal {R}_0\mathcal {P}}+S^1_0e^{-k\mathcal {R}_0\mathcal {P}} \right) =1-c+c\left( S^\emptyset _0+S^1_0-\mathcal {P}\right) \end{aligned}$$

with $S^1_1=c\left( 1-S^\emptyset _0-S^1_0\right)$ and $R_1=1-S^\emptyset _1-S^1_1$. The within-season and between-season dynamics of the model are shown in Fig. 1c.

Using the subscript n to define the initial condition of a state variable at the beginning of season n, or the proportion of the population infected during season n, we have the inter-season map

$$\begin{aligned} \mathbf {s}_{n+1}= \mathbf {C}\mathbf {s}_n+\mathbf {q}(\mathcal {P}_n) \end{aligned}$$

(4)

where

$$\begin{aligned} \mathbf {s}_n=\left( \begin{array}{c} S^{\emptyset }_n \\ S^1_n \end{array}\right) \qquad \mathbf {C}=\left( \begin{array}{cc} c &{} c\\ -c &{} -c\end{array}\right) \qquad \mathbf {q}(\mathcal {P}_n)=\left( \begin{array}{c} 1-c-c\mathcal {P}_n \\ c \end{array}\right) \end{aligned}$$

and

$$\begin{aligned} \mathcal {P}_n= S^{\emptyset }_n \left( 1-e^{-\mathcal {R}_0\mathcal {P}_n}\right) +S^{1}_n\left( 1-e^{-k\mathcal {R}_0\mathcal {P}_n}\right) \end{aligned}$$

(5)

if $\left( S^{\emptyset }_n+k S^{1}_n \right) \mathcal {R}_0>1$, and $\mathcal {P}_n=0$ otherwise. Fixed points of the map (4) solve

$$\begin{aligned} \mathbf {s}^*=\left( \begin{array}{c} 1-c\mathcal {P}^*-c^2\mathcal {P}^* \\ c^2\mathcal {P}^* \end{array}\right) \end{aligned}$$

where

$$\begin{aligned} \mathcal {P}^*= S^{\emptyset *} \left( 1-e^{-\mathcal {R}_0\mathcal {P}^*}\right) +S^{1*}\left( 1-e^{-k\mathcal {R}_0\mathcal {P}^*}\right) \end{aligned}$$

Combining these results, $\mathcal {P}=\mathcal {P}^*$ solves

$$\begin{aligned} \mathcal {P}=\dfrac{1-e^{-\mathcal {R}_0\mathcal {P}} }{\left( 1+c \right) \left( 1-ce^{-\mathcal {R}_0\mathcal {P}} \right) +c^2e^{-k\mathcal {R}_0\mathcal {P}}} \end{aligned}$$

The Jacobian matrix of the map (4) is

$$\begin{aligned} \mathbf {J}\left( S^{\emptyset }_n,S^1_n\right) =\left( \begin{array}{cc} c- c\dfrac{\partial \mathcal {P}_n}{\partial S^{\emptyset }_n} &{} c- c\dfrac{\partial \mathcal {P}_n}{\partial S^1_n} \\ -c &{} -c\end{array}\right) \end{aligned}$$

Differentiating $\mathcal {P}_n$ we obtain

$$\begin{aligned} \frac{\partial \mathcal {P}_n}{\partial S^{\emptyset }_n}&=\frac{1- e^{-\mathcal {R}_0\mathcal {P}_n}}{1-\mathcal {R}_0S^{\emptyset }_n e^{-\mathcal {R}_0\mathcal {P}_n}-k\mathcal {R}_0S^1_ne^{-k\mathcal {R}_0\mathcal {P}_n}} \\ \frac{\partial \mathcal {P}}{\partial S^1_n}&=\frac{1- e^{-k\mathcal {R}_0\mathcal {P}_n}}{1-\mathcal {R}_0S^{\emptyset }_n e^{-\mathcal {R}_0\mathcal {P}_n}-k\mathcal {R}_0S^1_ne^{-k\mathcal {R}_0\mathcal {P}_n}} \end{aligned}$$

The fixed point is stable if the spectral radius of $\mathbf {J}$ is less than one.

The fixed point of the between-season map corresponds to a periodic solution with period $p=1$. We now investigate period-p solutions of the model, defined by $\mathcal {P}_{n+p}=\mathcal {P}_n$ for all n and integer $p>1$. For example, a period-2 solution is a fixed point of the map

$$\begin{aligned} \mathbf {s}_{n+2}=\mathbf {C}\left( \mathbf {C}\mathbf {s}_{n}+\mathbf {q}(\mathcal {P}_{n})\right) +\mathbf {q}(\mathcal {P}_{n+1}) \end{aligned}$$

but $\mathbf {C}^2=\mathbf {0}$ so

$$\begin{aligned} \mathbf {s}_{n+2}=\mathbf {C}\mathbf {q}(\mathcal {P}_{n})+\mathbf {q}(\mathcal {P}_{n+1})=\left( \begin{array}{c} 1 - c\mathcal {P}_{n+1} - c^2\mathcal {P}_n \\ c^2\mathcal {P}_n \end{array}\right) \end{aligned}$$

(6)

There are two possibilities, we either have an epidemic in alternate years, or an epidemic in both years.

• 1. A special case of a period-2 solution is when there are epidemics in alternate years, hence $\mathcal {P}_{n+1}=0$. We then have $S^{\emptyset }_{n} = 1-c^2\mathcal {P}_n$ and $S^1_{n} = c^2\mathcal {P}_n$. Substituting in Eq. 5 we obtain

  $$\begin{aligned} \mathcal {P}_n =\frac{1-e^{-\mathcal {R}_0\mathcal {P}_n}}{1+c^2\left( e^{-k\mathcal {R}_0\mathcal {P}_n}-e^{-\mathcal {R}_0\mathcal {P}_n}\right) } \end{aligned}$$

• 2. More generally, a period-2 solution would have epidemics in every year, with $\mathcal {P}_{n+2}=\mathcal {P}_{n}$ but $\mathcal {P}_{n+1}\ne \mathcal {P}_{n}$ and $\mathcal {P}_{n+1}\ne 0$. From Eq. 6 we have $S^1_{n} = c^2\mathcal {P}_n$ and substitution into Eq. 5, followed by some minor rearranging yields

  $$\begin{aligned} S^{\emptyset }_{n}=\frac{1-c^2\left( 1-e^{-k\mathcal {R}_0\mathcal {P}_n}\right) }{1-e^{-\mathcal {R}_0\mathcal {P}_n}}\mathcal {P}_n \end{aligned}$$

  Also from Eq. 6, $S^{\emptyset }_{n} = 1-c\mathcal {P}_{n+1}-c^2\mathcal {P}_n$, so for consistency

  $$\begin{aligned} \mathcal {P}_{n+1}=\frac{1}{c}\left( 1- \frac{1+c^2\left( e^{-k\mathcal {R}_0\mathcal {P}_n}-e^{-\mathcal {R}_0\mathcal {P}_n}\right) }{1-e^{-\mathcal {R}_0\mathcal {P}_n}}\mathcal {P}_n \right) \end{aligned}$$

If it exists, the period-2 solution is stable when the spectral radius of the product $\mathbf {J}\left( S^{\emptyset }_{n+1},S^1_{n+1}\right) \mathbf {J}\left( S^{\emptyset }_n,S^1_n\right)$ is less than one.

Equation 6 provides a relationship between the current size of the susceptible populations and the final size of the two previous epidemics. This relationship holds even if c is not a constant, but changes from year to year. The dynamics that follow from a variable c, used to model punctuated antigenic drift, are explored in Sect. 4.

3.3 Results

For fixed values of the dilution constant c, we calculated the fixed points of the map, $\mathbf {s}^*$, the associated final size of the epidemic, $\mathcal {P}^*$, and their stability over a range of values of $\mathcal {R}_0$ and k. Examples are shown in Figs. 2, 3 and 4. In Fig. 2a, b the blue shaded regions correspond to parameter values where the fixed point is locally stable, and the red regions to values where it is unstable. Also shown on the figure are horizontal black lines corresponding to the parameter ranges used for the bifurcation diagrams in Fig. 2c, d and the orbit diagrams in Fig. 3.

Figure 2a shows that for $c=0.7$ the fixed point is stable apart from a region of parameter space with high values of k. This is illustrated in Fig. 2c, which shows that when $\mathcal {R}_0=2$ a stable steady state bifurcates to a period-2 solution at $k=0.694$. The periodic solution has no epidemic $\left( \mathcal {P}=0\right)$ in alternate years. Figure 2b shows that for $c=0.9$ there are two regions of parameter space where the fixed point is unstable. Figure 2d shows that when $\mathcal {R}_0=8$ a stable steady state bifurcates to a period-2 solution at $k=0.344$. The periodic solution has epidemics in every year, alternating in size, until $k=0.431$ where it becomes a solution with epidemics in alternate years. At $k=0.703$ the fixed point regains stability, with an unstable period-2 solution appearing through a pitchfork bifurcation. Numerical results (not shown) have confirmed the existence of the two attractors and the bistable behaviour in this parameter region.

Fig. 2

a and b Maps of the $\left( k,\mathcal {R}_0\right)$ parameter space with a$c=0.7$ and b$c=0.9$. In blue areas the fixed point is stable, in red areas it is unstable. c and d Bifurcation diagrams, the final size $\mathcal {P}$ against k with c: $\mathcal {R}_0=2$ and $c=0.7$ corresponding to the horizontal black line in a; and d: $\mathcal {R}_0=8$ and $c=0.9$ corresponding to the upper horizontal black line in b. Stable fixed points and period-2 solutions are shown with unbroken lines, unstable fixed points and period-2 solutions with broken lines (color figure online)

Fig. 3

Orbit diagrams with $\mathcal {R}_0=2$ and $c=0.9$. a and b the final size $\mathcal {P}$ against k. The unstable fixed point is shown by a broken black line. c and d the proportions in the susceptible compartments $S^\emptyset$ and $S^1$ against k. Initial conditions are $I(0)=10^{-5}$ with: a and c, $S^\emptyset _0=1.0-I(0)$, $S^1_0=0$; b and d, $S^\emptyset _0=0.5-I(0)$, $S^1_0=0.5$

Fig. 4

The $\left( S^\emptyset ,S^1\right)$ plane when $\mathcal {R}_0=2$ and $c=0.9$ showing a the attractors when $k=0.2$; b the basin of attraction for the attractors when $k=0.2$; c: the attractor when $k=0.28$; and d: the attractors when $k=0.18$. In a and d a period-3 attractor is shown with blue circles. In a, c and d the regions in red are a high-period attractor, the thin black lines show a sample of transitions between subsequent points on the attractor, and the yellow dots are arranged at every fifth point. In a, c and d the broken black line is at $\mathcal {R}_0\left( S^\emptyset +kS^1\right) =1$. In b, regions in blue are initial conditions for the period-3 attractor (blue circles in a), and regions in red are initial conditions for the high-period attractor (red structure in a) (color figure online)

The dynamics along the lower of the horizontal lines in Fig. 2b, at $\mathcal {R}_0=2$ and $c=0.9$ for $0.1<k<0.6$, are shown in Fig. 3. Orbit diagrams were constructed by fixing k, iterating the map one thousand times, then plotting the results of the next one hundred iterations. The value of k was then incremented and the process repeated. For Fig. 3a, c the initial conditions were: $S^\emptyset =1-10^{-5}$, $S^1=0$; for Fig. 3b, d: $S^\emptyset =0.5-10^{-5}$, $S^1=0.5$. The orbit diagrams show a stable fixed point for $0.281<k<0.482$, and a period-2 solution with epidemics in alternating years for $k>0.482$. The attractor for small values of k is a period-3 solution: two years without an epidemic $\left( \mathcal {P}=0\right)$ followed by a single year with $\mathcal {P}\simeq 0.8$, most clearly seen by inspecting results for $S^\emptyset$ in Fig. 3c. Between these values and the value where the fixed point becomes stable is a region with more complicated dynamics. Comparing Fig. 3a, c with Fig. 3b, d, there is a range of k where some initial conditions generate dynamics tending to the period-3 attractor, and other initial conditions that lead to an attractor with complicated dynamics, either period-p with p large or chaotic. We refer to this as a high-period attractor. These dynamics are explored further in Fig. 4, where results are presented on the feasible region of the $\left( S^\emptyset ,S^1\right)$ plane. Figure 4a shows the two attractors present when $k=0.2$. The period-3 attractor is shown with blue circles, and the high-period attractor is shown in red. Although the pentagon-like shape of the attractor is suggestive of a period-5-like structure, the yellow points plotted for every fifth iteration are on every ‘side’ of the pentagon, indicating that this is not the case. Figure 4b shows the basins of attraction for the two attractors in Fig. 4a. The regions in blue are initial conditions for the period-3 attractor, and the region in red corresponds to initial conditions for the high-period attractor. For comparison with the results shown in Fig. 4a, c shows a single high-period attractor for $k=0.28$, and Fig. 4d shows a period-3 attractor and a high-period attractor for $k=0.18$. The dynamics when $k=0.28$ (Fig. 4c) are particularly intriguing; the attractor has five isolated points, presumably a remnant of the pentagonal shape seen for other values of k, and every fifth iterate (yellow dots) lies on a distinct segment of the attractor’s boundary. In Fig. 4a, c, d we show a broken black line defined by $\mathcal {R}_0\left( S^\emptyset +kS^1\right) =1$. The regions of the $\left( S^\emptyset ,S^1\right)$ plane to the right of that line correspond to values of $S^\emptyset$ and $S^1$ for which there is an epidemic, $\mathcal {P}\ne 0$.

The results presented in Fig. 4 are strikingly similar to those obtained for the example analysed in detail by Andreasen (2003), and summarised in Fig. 6 of that reference. The model analysed in Andreasen (2003) was one where previous infection reduced subsequent infectiousness rather than reduced susceptibility. However, for some parameter values a stable limit cycle existed, whereas other values gave rise to either a stable limit cycle or chaotic dynamics, with different basins of attraction. The dynamics of both models are typical of systems with border-collision bifurcations (Simpson 2016 and Simpson pers. comm.), and should be the subject of further investigation.

4 Punctuated antigenic drift

4.1 Between-season dynamics

In Sect. 3.2 we assumed that population turnover and antigenic drift occur by the same amount each year, as measured by their affect on the proportion of the population regarded as susceptible to infection. That is, the values of b and d, and hence c are independent of the year. This is clearly not the case. In particular, antigenic drift has been shown to proceed in almost a step-like manner, with minimal change over several years followed by a more pronounced change over one year (Smith et al. 2004). Here we relax the requirement that the dilution constant c is independent of the year, replacing it with a value $c_n$. The map in Eq. 4 becomes

$$\begin{aligned} \mathbf {s}_{n+1}=\mathbf {C}_n\mathbf {s}_{n} +\mathbf {q}_n(\mathcal {P}_n) \end{aligned}$$

Regardless of the values $c_n$ we still have $\mathbf {C}_n\mathbf {C}_{n+1}=\mathbf {0}$, so the mapping only depends on the current and immediately preceding years:

$$\begin{aligned} \mathbf {s}_{n+2} =\mathbf {C}_{n+1}\mathbf {s}_{n+1} +\mathbf {q}_{n+1}(\mathcal {P}_{n+1})&=\mathbf {C}_{n+1}\mathbf {q}_{n}(\mathcal {P}_n) +\mathbf {q}_{n+1}(\mathcal {P}_{n+1}) \\&=\left( \begin{array}{c} 1- c_{n} c_{n+1}\mathcal {P}_n - c_{n+1}\mathcal {P}_{n+1}\\ c_{n} c_{n+1}\mathcal {P}_n \end{array}\right) \end{aligned}$$

In order to make a comparison with our previous results, we assume that the $c_n$ vary in a cycle with $c_{n+m}=c_n$ for some integer m and all n. We then look for cycles with $\mathcal {P}_{n+m}=\mathcal {P}_n$, hence $\mathbf {q}_{n+m}(\mathcal {P}_{n+m})=\mathbf {q}_n(\mathcal {P}_n)$ and $\mathbf {s}_{n+m}=\mathbf {s}_{n}$. We focus on dynamics with $m=5$. We assume population turnover to be constant ($2\%$ per year), and no antigenic drift for four years followed by some drift following the fifth year. Then the same pattern repeats. To have some basis for comparison with the previous examples we set $0.98^4c_{n+1}=c^5$ when n is an integer divisible by five, and $c_n=0.98$ otherwise, with $c=0.7$ or $c=0.9$. Hence $\prod _{j=1}^5 c_{n+j}=c^5$ for all n. For comparison we also computed some results with $m=3$ and with $m=7$, setting $\prod _{j=1}^m c_{n+j}=c^m$.

Fig. 5

a Orbit diagram for the model with punctuated antigenic drift $\left( m=5\right)$, the final size $\mathcal {P}$ against k with $\mathcal {R}_0=2$, $c_{5m+1}=0.1822$, m integer and $c_n=0.9800$ otherwise, for $0.4<k<1.0$. b Orbit diagram for the model with punctuated equilibrium, the final size $\mathcal {P}$ against k with $\mathcal {R}_0=8$, $c_{5m+1}=0.6402$, m integer and $c_n=0.9800$ otherwise, for $0.1<k<0.8$. c Thirty consecutive values of $\mathcal {P}$ with $k=0.5375$ and other parameters as in a. d Thirty consecutive values of $\mathcal {P}$ with $k=0.625$ and other parameters as in b. The sequence of successive values is colour coded red, magenta, cyan, blue, green (color figure online)

Fig. 6

a Orbit diagram for the model with punctuated antigenic drift $\left( m=5\right)$, the final size $\mathcal {P}$ against k with $\mathcal {R}_0=2$, $c_{5m+1}=0.6402$, m integer and $c_n=0.9800$ otherwise, for $0.1\leqslant k\leqslant 0.6$. The initial conditions are $S^\emptyset _0=1-10^{-5}$ and $S_0^1=0$. The sequence of successive values is colour coded red, magenta, cyan, blue, green. b The $\left( S^\emptyset ,S^1\right)$ plane showing initial conditions for the period-20 attractor (green) and the high-period attractor (red) with $k=0.18$ and other parameter values as in a. c The period-20 attractor (green) and the high-period attractor (red) with all parameter values as in b. d The period-20 attractor (green) and the high-period attractor (red) with all parameter values as in b. The thin black lines show a sample of transitions between subsequent points on the attractor, and the yellow dots are arranged at every fifth point. In c and d the broken black line is at $\mathcal {R}_0\left( S^\emptyset +kS^1\right) =1$ (color figure online)

Fig. 7

a Orbit diagram for the model with punctuated antigenic drift $\left( m=3\right)$, the final size $\mathcal {P}$ against k with $\mathcal {R}_0=2$, $c_{3m+1}=0.7591$, m integer and $c_n=0.9800$ otherwise, for $0<k<1$. The sequence of successive values is colour coded red, magenta, cyan. b Orbit diagram for the model with punctuated equilibrium $\left( m=7\right)$, the final size $\mathcal {P}$ against k with $\mathcal {R}_0=2$, $c_{7m+1}=0.5399$, m integer and $c_n=0.9800$ otherwise, for $0<k<1$. The sequence of successive values is colour coded red, magenta, cyan, blue, green, black, yellow. In a and b the initial conditions are $S^\emptyset _0=1-10^{-5}$ and $S_0^1=0$. c The period-3 attractor (blue) and the high-period attractor (red) with parameter values as in a and $k=0.24$. d The period-28 attractor (green) and the high-period attractor (red) with parameter values as in c and $k=0.3$. In c and d the broken black line is at $\mathcal {R}_0\left( S^\emptyset +kS^1\right) =1$ (color figure online)

4.2 Results

Figures 5, 6 and 7 show results from the model with punctuated antigenic drift. Figure 5a, b present orbit diagrams showing the final size of epidemics for a range of values of k and $m=5$, and Fig. 5c, d show samples of these dynamics for particular values of k.

Figure 5a shows results for $\mathcal {R}_0=2$, $\prod _{j=1}^5 c_{n+j}=0.7^5$ (the product of the $c_n$’s chosen to coincide with those used earlier, see Fig. 2c) and $0.4<k<0.95$. The figure shows a five-cycle for all values of k. For $0.4<k<0.477$ there is an epidemic in three out of five years, for $0.477<k<0.554$ an epidemic in four out of five years, then for $0.554<k<0.95$ an epidemic in three out of five years. The pattern is illustrated for $k=0.5375$ in the bar graph, Fig. 5c. Figure 5b shows results for $\mathcal {R}_0=8$, $\prod _{m=1}^5 c_{n+m}=0.9^5$ (again the product of the $c_n$’s chosen to coincide with those used earlier, see Fig. 2d) and $0.1\leqslant k\leqslant 0.8$. The figure shows a bifurcation from a period-5 solution to a period-10 solution at $k=0.340$. This is very close to the critical value of k for the bifurcation from a fixed point to a period-2 solution seen in Fig. 2d. This pattern is illustrated for $k=0.625$ in the bar graph, Fig. 5d, clearly showing a ten-cycle.

Figure 6a shows an orbit diagram for $\mathcal {R}_0=2$, $\prod _{j=1}^5 c_{n+j}=0.9^5$ and $0.1<k<0.6$. For $0.1<k<0.133$ there is a fifteen-cycle, with epidemics in seven out of fifteen years. This pattern is not visible in the figure as many components of the orbit are superimposed, but has been confirmed by numerical results (not shown). Figure 6a shows complicated behaviour for $0.133<k<0.342$ and $0.499<k<0.6$. For $0.342<k<0.499$ the solution is a five-cycle, with epidemics in two successive years followed by a year with no epidemic, then a year with an epidemic followed by a year without. The pattern then repeats. The pattern is the same (but with different values) either side of the abrupt transition at $k=0.371$. As an example of the complicated dynamics that eventuate, the basins of attraction for the two attractors present when $k=0.275$ is shown in Fig. 6b. Some initial conditions for $S^\emptyset$ and $S^1$ result in a period-20 solution, and others lead to a high-period attractor. The two attractors are shown in the $\left( S^\emptyset ,S^1\right)$ phase plane in Fig. 6c, d. As in Fig. 4d, the black lines in Fig. 6d connect subsequent points on the attractor, with every fifth point marked in yellow. It can be seen that for this example the five disjoint regions of the attractor are generated by the periodicity of the $c_n$’s.

Figure 7 shows results from the model with punctuated antigenic drift and $m\ne 5$. Figure 7a shows an orbit diagram for $m=3$ with $\mathcal {R}_0=2$, $\prod _{j=1}^3 c_{n+j}=0.9^3$ and $0<k<1.0$. The results show a three-cycle for $k<0.224$, with an epidemic every three years, each epidemic followed by two years with no epidemic. There is also a three-cycle for $0.295<k<0.479$, but with an epidemic every year. For other values of $k<0.565$ the results show complicated dynamics, but there is a twelve-cycle for $0.565<k<0.640$, and a period-6 solution for $k>0.640$. Figure 7c shows an example of the complicated dynamics for these parameter values. When $k=0.24$ there is a period-3 attractor and a high-period attractor. Figure 7b, d show similar results for $m=7$. For small values of k there is a period-21 solution, with many years with $\mathcal {P}=0$. For larger values of k (close to one) there is a period-14 solution. Intermediate values show a variety of dynamic behaviour. For example, Fig. 7d shows that when $k=0.3$ there is both a period-28 attractor, and a high-period attractor.

5 Discussion

We have presented and analysed a model for the dynamics of an infectious disease with annual epidemics and waning protection between epidemics, motivated by the epidemiology of a single strain or subtype of influenza. The model incorporates replenishment of susceptibles due to population turnover, the waning of host immunity and antigenic drift. We have examined the situations where antigenic drift is constant from year to year, as measured by the mean retention of immunity in the population, and where drift is punctuated (as it is for influenza A H3N2), essentially zero for some years followed by a more pronounced change (Smith et al. 2004). Our model is equivalent to a special case of the models proposed by Andreasen (2003), but not one that was analysed in detail. The philosophy of combining a differential equation model for within season dynamics with a discrete mapping between seasons was also used for a simple model for parasite populations (Roberts and Heesterbeek 1998). That model was also found to exhibit complicated dynamics for some regions in parameter space.

From a mathematical perspective, we were able to find analytical results for the final size of the within-season epidemic, and for the fixed point and period-2 solution of the between-season map when the antigenic drift was constant from year to year. We also found some analytic results for the model with punctuated drift, and were able to compute the dynamics for examples with constant and punctuated drift. We found parameter regions where the between-season map had a fixed point, implying identical epidemics every year. We also found regions with period-2, period-3 and other periodic solutions, and dynamics referred to as high-period which may correspond to the existence of a chaotic attractor. We found examples where the complex attractor coexisted with a stable period-3 solution, the two attractors having distinct basins of attraction. Where antigenic drift was punctuated, with drift being confined to once every m years, the resulting attractor had a period that was an integer multiple of m. A curious feature of the results was the pentagon-like structure of the attractor for examples with $m=1$ (Fig. 4a) and $m=5$ (Fig. 6c), and triple pentagons when $m=3$ (Fig. 7c).

From an epidemiological or public health perspective, our results are of interest because they help to establish that the observed variation in the year-to-year behaviour of influenza epidemics may arise directly from the underlying dynamical structure of the system. Put another way, variation is to be expected, even in the absence of changes in parameters due to the ecological (strain) structure of the virus (Hadfield et al. 2018; Nextstrain Real-time tracking 2018) or, say, incorporation of stochasticity into a model. Given this, there is no reason to expect influenza epidemics to be regular, or for recent past behaviour to give a clear indication of a future season’s behaviour. Of course, this is not to say that short term predictions based on genetic data cannot be made for the purposes of vaccine strain selection (Bedford and Neher 2018), and indeed, combining insights from those approaches with models such as that presented here may provide new insight into the problem of vaccine strain selection.

Returning to the dynamics of the system, if one were willing to consider our model to be a sufficiently realistic model of influenza epidemic activity, and that activity to be well described by the model with parameters that lead to chaotic dynamics, then not only do we expect the variability, but we cannot predict any of the longer term trends in epidemic behaviour. This has potential implications for how influenza epidemics are evaluated and, perhaps, responded to over the longer year-to-year timescale.

Although we have presented a simple model in terms of how past infection induces a heterogeneity of susceptibility in a population, our analysis has revealed a rich behaviour for which different parameter values lead to different dynamics. The between-season map demonstrates dynamical behaviour reminiscent of actual influenza dynamics for specific parameter ranges, in so far as the dynamics are difficult to predict. Other authors have proposed influenza models with different mechanisms for seasonal forcing, and observed a variety of complicated dynamics. It has been observed that with the inclusion of forcing terms (Bacaer and Ouifki 2007; He and Earn 2007) or other non-linear feedback mechanisms (Earn et al. 2000; Viboud et al. 2006), these systems are capable of supporting undamped oscillatory dynamics, and a host of more complex dynamical behaviour. Similarly Mathews et al. (2007), in a study of pandemic influenza in 1918–19, used a pair of iterated discrete maps to study re-infection patterns across multiple waves and troughs, with a focus on transient dynamics.

While we have focused on the seasonal dynamics of a single virus, influenza presents annually with multiple strains or subtypes. Where there is an interaction between these virus variants, the relative timing of their introduction has an influence on their dynamics, see for example Roberts (2012). While we and others have shown that the dynamics of a single virus variant leads to complicated behaviour, the dynamics of interacting variants will be much harder to describe, model and predict. Finding methods for analysing these dynamics, and applying them in the context of public health, is the subject of ongoing research.

Appendix A: Stability of inter-season map for the preliminary model

Here we establish that when a fixed point of the map in Eq. 1 exists, it is stable. The relationship between the proportions susceptible at the beginning and end of the epidemic can be expressed $S_0e^{-\mathcal {R}_0S_0} =S_\infty e^{-\mathcal {R}_0S_\infty }$. Writing $y=\mathcal {R}_0S_\infty$ and $x=\mathcal {R}_0S_0$, we see that $ye^{-y}=xe^{-x}$, and for an epidemic to occur $x>1$ and $y<1$. The between-season map is $S_1=1-c+cS_\infty$, so

$$\begin{aligned} \frac{\text{ d }S_1}{\text{ d }S_0}=c \frac{\text{ d }S_\infty }{\text{ d }S_0} \end{aligned}$$

As $c<1$, for stability we need to prove that $\left| y^{\prime }(x)\right| <1$ for $x>1$. Differentiating

$$\begin{aligned} y^{\prime }(x)=\frac{\left( 1-x\right) e^{-x}}{\left( 1-y\right) e^{-y}}=\frac{\left( 1-x\right) y}{\left( 1-y\right) x}<0 \end{aligned}$$

(A1)

At $x=y=1$ we use l’Hôpital’s rule to get

$$\begin{aligned} y^{\prime }(1)=\lim _{x\rightarrow 1}\frac{\left( 1-x\right) y}{\left( 1-y\right) x}=\lim _{x\rightarrow 1}\frac{\left( 1-x\right) y^{\prime }(x)-y}{1-y-xy^{\prime }(x)} \end{aligned}$$

and hence $y^{\prime }(1)=-1$. Re-arranging Eq. A1

$$\begin{aligned} \left( \frac{1}{y}-1\right) y^{\prime }(x)=\frac{1}{x}-1 \end{aligned}$$

and differentiating

$$\begin{aligned} \left( \frac{1}{y}-1\right) y^{\prime \prime }(x)=\frac{y^{\prime }(x)^2}{y^2}-\frac{1}{x^2} \end{aligned}$$

If we had $y^{\prime }(x)=-1$, then recalling $x>1$ and $y<1$ the right hand side is positive, so $y^{\prime \prime }(x)>0$. Hence for $x>1$ we have $0>y^{\prime }(x)>-1$, and the result is established.