Invariant predictions of epidemic patterns from radically different forms of seasonal forcing

TL;DR

This paper shows that key epidemic model behaviors are unaffected by the exact seasonal forcing pattern, allowing reliable predictions even with approximate seasonal data.

Contribution

It demonstrates the invariance of epidemic bifurcations to the shape of seasonal forcing, simplifying modeling under uncertain seasonal patterns.

Findings

Bifurcations are invariant to forcing shape when amplitude is adjusted.

Qualitative epidemic pattern changes can be predicted without precise seasonal data.

Invariance observed in predator-prey models suggests broader applicability.

Abstract

Seasonal variation in environmental variables, and in rates of contact among individuals, are fundamental drivers of infectious disease dynamics. Unlike most periodically-forced physical systems, for which the precise pattern of forcing is typically known, underlying patterns of seasonal variation in transmission rates can be estimated approximately at best, and only the period of forcing is accurately known. Yet solutions of epidemic models depend strongly on the forcing function, so dynamical predictions---such as changes in epidemic patterns that can be induced by demographic transitions or mass vaccination---are always subject to the objection that the underlying patterns of seasonality are poorly specified. Here, we demonstrate that the key bifurcations of the standard epidemic model are invariant to the shape of seasonal forcing if the amplitude of forcing is appropriately adjusted. Consequently, analyses applicable to real disease dynamics can be conducted with a smooth, idealized sinusoidal forcing function, and qualitative changes in epidemic patterns can be predicted without precise knowledge of the underlying forcing pattern. We find similar invariance in a seasonally forced predator-prey model, and conjecture that this phenomenon---and the associated robustness of predictions---might be a feature of many other periodically forced dynamical systems.

Tables

Table 1. Table 1: Invariance of fold bifurcations [ 22 ] at different ℛ 𝟎 subscript ℛ 0 \bm{\mathcal{R}_{0}} values when the principal period doubling (PD) bifurcation at ℛ 𝟎 = 15.12 subscript ℛ 0 15.12 \bm{\mathcal{R}_{0}=15.12} is matched. Fold ( n 𝑛 n ) refers to a fold bifurcation that gives rise to a period n 𝑛 n attractor (an n 𝑛 n -year epidemic cycle) as ℛ 0 subscript ℛ 0 \mathcal{R}_{0} is increased. The PD does not occur at precisely the same ℛ 0 subscript ℛ 0 \mathcal{R}_{0} for each ( p , α ) 𝑝 𝛼 (p,\alpha) pair due to slight inaccuracies of the numerical continuation software [ 23 , 37 ] . The relative difference refers to max ⁡ [ ( x − x tt ) / x tt ] 𝑥 subscript 𝑥 tt subscript 𝑥 tt \max[(x-x_{\rm tt})/x_{\rm tt}] , where x 𝑥 x is the value of ℛ 0 subscript ℛ 0 \mathcal{R}_{0} at the bifurcation of interest and x tt subscript 𝑥 tt x_{\rm tt} is its value for term-time forcing ( p = − 1 𝑝 1 p=-1 ). All the data in this table are plotted in figure  5 .

Forcing Function		${𝓡}_{\mathrm{𝟎}}$ Bifurcation Point
$p$	$\alpha$	PD (2)	Fold (3)	Fold (4)	Fold (5)	Fold (6)	Fold (7)
$-1.00$	$0.2500$	$15.1199$	$8.8991$	$6.2374$	$4.9755$	$4.2623$	$3.8094$
$-0.50$	$0.1012$	$15.1254$	$8.9014$	$6.2363$	$4.9713$	$4.2554$	$3.7995$
$0.00$	$0.0782$	$15.1146$	$8.8946$	$6.2247$	$4.9554$	$4.2357$	$3.7758$
$0.25$	$0.0839$	$15.1248$	$8.8992$	$6.2292$	$4.9599$	$4.2401$	$3.7800$
$1.00$	$0.1000$	$15.1234$	$8.8976$	$6.2269$	$4.9571$	$4.2370$	$3.7767$
$2.00$	$0.1182$	$15.1240$	$8.8972$	$6.2261$	$4.9558$	$4.2354$	$3.7749$
Relative Difference		$0.0004$	$0.0005$	$0.0020$	$0.0040$	$0.0063$	$0.0090$

Equations

$$\frac{\mathrm{d}S}{\mathrm{d}t}$$

$$\frac{\mathrm{d}I}{\mathrm{d}t}$$

$$\frac{\mathrm{d}R}{\mathrm{d}t}$$

$${\beta }_{\cos }(t)=⟨\beta ⟩[1+\alpha \cos (2\pi t)],$$

$${\beta }_{\mathrm{tt}}(t)=\{\begin{array}{ll}{\beta }_{\mathrm{high}}& \text{mbox}schooldays,\\ {\beta }_{\mathrm{low}}& \text{mbox}otherdays,\end{array}$$

$$⟨\beta ⟩=p_{\mathrm{s}}{\beta }_{\mathrm{high}}+(1-p_{\mathrm{s}}){\beta }_{\mathrm{low}},$$

$$\alpha =\frac{1}{2}\left(\frac{{\beta }_{\mathrm{high}}-{\beta }_{\mathrm{low}}}{⟨\beta ⟩}\right),$$

$${\mathrm{osc}}_{\mathrm{tt}}(t)=\{\begin{array}{ll}2(1-p_{\mathrm{s}})& \text{mbox}schooldays,\\ -2p_{\mathrm{s}}& \text{mbox}otherdays,\end{array}$$

$${\beta }_{\mathrm{tt}}(t)=⟨\beta ⟩[1+\alpha {\mathrm{osc}}_{\mathrm{tt}}(t)].$$

$$0\le \alpha \le \frac{1}{2p_{\mathrm{s}}},$$

$${\mathcal{R}}_0=\frac{⟨\beta ⟩}{\gamma +\mu }.$$

Full text

*Corresponding author (email: ip98@cornell.edu)

Keywords  Infectious disease, seasonal forcing, SIR epidemic model, predator-prey model, nonlinear oscillators, dynamical systems

Cite this article: Papst I, Earn DJD. 2019 Invariant predictions of epidemic patterns from radically different forms of seasonal forcing. J. R. Soc. Interface 16: 20190202. \urlhttp://dx.doi.org/10.1098/rsif.2019.0202

1 Introduction

Periodic forcing of infectious disease transmission arises from a number of sources, including seasonal weather changes[1, 2], annual cycles in birth rates[3, 4], and school terms[5, 6, 7, 8]. While the period of forcing is always one year, the pattern of forcing depends strongly on the source (e.g., nearly sinusoidal for weather or births, but sharp and asymmetric for school terms[9, 6]) and the forcing function is never known exactly.

The detailed structure of the seasonal forcing function strongly affects the resulting pattern of disease incidence. For example, figure 1 compares solutions of the standard seasonally forced susceptible-infectious-recovered (SIR) epidemic model [10, 11] (see §2) with different forcing patterns (school terms or a sinusoid). In each case, the basic reproduction number ($\mathcal{R}_{0}$) and mean infectious period ($1/\gamma$) are the same, and the amplitude of seasonality ($\alpha$) has been chosen to ensure that the system displays a strictly biennial cycle. The differences between the time series in figure 1 emphasize that the solution of a model can match the details of an observed temporal pattern of prevalence only if the pattern of forcing is known (or estimated accurately).

However, in situations in which we are primarily interested in the length of recurrent cycles of epidemics—or whether or not the recurrent epidemic pattern exhibits dynamical chaos—it is not clear that knowing the detailed forcing pattern is critical. Indeed, it was noted in ref. 6 (endnote 13) that in a parameter regime relevant to measles dynamics, qualitatively similar bifurcation structure can be obtained with different seasonal forcing functions. In particular, the structure of the bifurcation diagram that facilitates predictions of observed dynamical transitions (bifurcation parameter $\mathcal{R}_{0}$ with fixed seasonal forcing amplitude $\alpha$) is qualitatively similar for school term forcing with amplitude estimated from data ($\alpha=0.25$) or sinusoidal forcing with much lower amplitude ($\alpha=0.08$).

Here, for a wide range of forcing functions, we examine the quantitative and qualitative bifurcation structure of the SIR model (below) and a seasonally forced predator-prey model (in \hrefrun:rsif20190202supp1.pdf\hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material, §LABEL:sec:supppredprey). Revealing invariance of dynamical transitions in these systems deepens our understanding of epidemic and predator-prey dynamics, and suggests a more general type of invariance in the dynamical structure of forced nonlinear oscillators.

2 Methods

2.1 Seasonally forced SIR model

The standard seasonally forced SIR model for a closed population can be expressed as a system of three ordinary differential equations,

[TABLE]

where $S$, $I$ and $R$ denote the proportions of the population that are susceptible, infectious and recovered from the infection, $\mu$ is the per capita birth and death rate, $\gamma$ is the recovery rate, and $\beta(t)$ is the time-dependent transmission rate. The first two equations above do not depend on the variable $R$, and $R=1-S-I$, so the third equation is not required to specify the dynamics.

In most previous work, the transmission rate $\beta(t)$ has been taken to vary sinusoidally for simplicity[12, 13, 14, 15],

[TABLE]

where $\langle\beta\rangle$ is the mean transmission rate, $\alpha$ is the amplitude of seasonality, and time is measured in units of the forcing period (one year). For childhood diseases such as measles, mumps, rubella or whooping cough, the primary source of seasonal forcing is typically school terms[5]. Consequently, a substantial body of work has employed a forcing function that is high in “term-time” and low otherwise[9, 16, 6, 17, 18, 19],

[TABLE]

where $\beta_{\rm high}$ and $\beta_{\rm low}$ are constant transmission rates, with $\beta_{\rm high}>\beta_{\rm low}$. If $p_{\rm{s}}$ is the proportion of the year that children spend in school then

[TABLE]

so if we define the amplitude of seasonality to be

[TABLE]

and denote the term-time “oscillation” function to be

[TABLE]

we can rewrite equation 2.3 in terms of the two parameters $\langle\beta\rangle$ and $\alpha$,

[TABLE]

The amplitude of forcing, $\alpha$, is constrained to lie in the range[18]

[TABLE]

to ensure that the forcing function is always non-negative. Regardless of the form of periodic forcing, the basic reproduction number (the average number of secondary cases per primary case in a wholly susceptible population) is[20]

[TABLE]

2.2 Family of forcing functions

In order to better understand how the dynamical structure of the SIR model depends on the pattern of seasonality in transmission, we constructed a continuous family of forcing functions that connects term-time and sinusoidal forcing. The family is parameterized with a shape parameter $p$, such that the forcing follows the school term schedule for $p=-1$, is square and symmetric for $p=0$, and is sinusoidal for $p=1$. The forcing function, which we denote $\operatorname{osc}_{p}(t)$, is symmetric for $p\geq 0$ and asymmetric for $p<0$ (more time above than below the mean). For $p>1$, the sine wave is squashed and in the limit $p\to\infty$ approaches two delta function impulses (up at integral times of $t=0,1,2,...$ years, and down at half-integral times of $t=0.5,1.5,2.5,...$ years). Several members of this family of forcing functions are plotted in figure 2; details of the construction can be found in \hrefrun:rsif20190202supp1.pdf\hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material, §LABEL:sec:defin-family-forc.

2.3 Stroboscopic map

Perturbative analysis of periodic orbits is usually simplified by considering a Poincaré map, which discretely samples continuous trajectories each time they pass through a particular hypersurface in the state space [21, §8.7]. For periodically forced systems, the natural Poincaré map is the stroboscopic map, which samples trajectories once per forcing period (cf. Example 8.7.2 in [21]). Our bifurcation analysis is based on the one-year stroboscopic map for the SIR system (2.1), so $n$-year periodic attractors are represented by $n$ discrete points.

2.4 Continuation of bifurcations

It is rarely possible to find explicit analytical expressions for bifurcation points as functions of model parameters. Instead, the conditions that define a given bifurcation are typically solved numerically and then “continued” in parameter space (i.e., starting from one bifurcation point, a curve in parameter space is computed on which the bifurcation conditions are satisfied to some level of precision[22]). We used standard open-source software (XPPAUT[23]) to continue bifurcation points in one and two model parameters.

The one-parameter bifurcation diagrams shown in figure 3 were constructed by finding equilibria and cycles of the annual stroboscopic map of the seasonally forced SIR model, and then continuing those points as functions of $\mathcal{R}_{0}$ (details in \hrefrun:rsif20190202supp1.pdf\hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material, §LABEL:sec:comp-bif-diags).

To find the sinusoidal forcing amplitude that yields a bifurcation diagram (figure 3, bottom panel) that is similar to that of the term-time forced model (figure 3, top panel), we proceeded as follows. Starting from the stable period doubling (PD) bifurcation point highlighted with a square in the term-time bifurcation diagram (figure 3, top panel, $\mathcal{R}_{0}=15.12$), we continued the bifurcation in the two-dimensional $(p,\alpha)$ parameter space. At each point on the resulting curve, shown in figure 4, the value of $\mathcal{R}_{0}$ for the PD remains the same. Consequently, this curve can be thought of as a function, $\alpha(p)$, which specifies the forcing amplitude ($\alpha$) that yields a PD at the same $\mathcal{R}_{0}$ value for any given shape of forcing pattern ($p$). In particular, the second point marked with a square in figure 4 shows that the forcing amplitude that yields a PD at $\mathcal{R}_{0}=15.12$ with sinusoidal forcing ($p=1$) is $\alpha=0.1$, the amplitude used to make the bifurcation diagram shown in the bottom panel of figure 3.

3 Results and Discussion

3.1 Quantitative invariance

Figure 3 displays unexpected similarity of bifurcation diagrams for SIR models with radically different forcing functions. The close correspondence between the positions of bifurcations in the top and bottom panels of figure 3 is surprising given how the diagrams were constructed. We first created the top panel (term-time forced SIR model bifurcation diagram as a function of $\mathcal{R}_{0}$). We then made the middle panel, using all the same parameter values but with sinusoidal rather than term-time forcing. To make the bottom panel, we used sinusoidal forcing again but adjusted the forcing amplitude so that the position ($\mathcal{R}_{0}$ value) of the principal period doubling bifurcation matched the position of the same bifurcation in the top panel ($\mathcal{R}_{0}=15.12$; see §2). Remarkably, with no further adjustments, the five fold bifurcations[22] listed in table 1 match to similar precision. These fold bifurcations occur in different parts of the parameter space on disconnected branches; yet their positions are conserved for a continuous family of forcing functions between term-time and sinusoidal (see figure 5, table 1, and \hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material), suggesting that the global dynamical structure of the model is determined by a property of the forcing function that is independent of its shape.

3.2 Qualitative invariance

Some aspects of the invariance are qualitative rather than quantitative. While the fold bifurcations at the left of the five “fold branches” in the $\mathcal{R}_{0}$ bifurcation diagrams (figure 3) appear to be invariant provided the position of the principal period doubling is conserved, the folds at the right edges of these branches do not line up precisely. Thus, the observed quantitative invariance is restricted to the “births” of these branches and not their “deaths”. Nevertheless, the main qualitative structure is preserved.

3.3 Unknown source of invariance

What the bifurcation-conserving property might be is not clear. Preliminary work indicated that fixing the average spectral power in the forcing function might be sufficient to conserve the bifurcation structure[17], but detailed analysis shows this to be false (see \hrefrun:rsif20190202supp1.pdf\hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material, figure LABEL:fig:compalphap-full). With the hope of identifying another spectral feature that might be associated with invariance of bifurcations, we also considered the shape of the full power spectrum of the forcing function. Unfortunately, as the pattern of forcing is smoothly changed from term-time to sinusoidal while keeping the bifurcation structure constant (see \hrefrun:rsif20190202supp1.pdf\hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material, figure LABEL:fig:powspecs), the full spectrum varies dramatically and no specific spectral feature underlying the observed bifurcation invariance was apparent to us.

3.4 Summary of results

Figure 5 and table 1 summarize our main results. Our continuation of the stable period doubling (PD) bifurcation (figure 4) determines the forcing amplitude that yields the same $\mathcal{R}_{0}$ value for the PD for each member of the family of forcing functions, including term-time ($p=-1$, $\alpha=0.25$) and sinusoidal ($p=1$, $\alpha=0.1$). This fact is indicated in figure 5 by the vertical alignment of the black squares (each square occurs at $\mathcal{R}_{0}=15.12$). The filled circles of each colour are also aligned vertically in figure 5, showing that as $\mathcal{R}_{0}$ increases the positions of all the folds that create attractors with periods longer than two years are invariant with respect to the forcing pattern. table 1 shows that these “birth folds” are invariant to similar precision as the PD as it is continued. The open circles and squares in figure 5 show that the positions of other bifurcations—the corresponding “death folds” and some intermediate PDs that appear in the longest period fold branches as $p$ is increased—are not invariant to the shape of the forcing function.

4 Conclusion

The bifurcation invariance that we have observed in the seasonally forced SIR model has significant practical importance, since the form of forcing in real infectious disease systems is never known precisely.

From a theoretical perspective, valuable rigorous analyses (e.g., proving the co-existence of multiple attractors[24] or the existence of chaos[25]) can be performed more easily if sinusoidal forcing is assumed; given the bifurcation invariance, such analyses are likely less dependent on the pattern of forcing than might previously have been thought.

From an applied perspective, predicting transitions in epidemic patterns [6, 17, 26, 19] depends on estimates of the pattern of forcing [27, 28, 29, 30], which can change in both amplitude and shape over decades or centuries [29, 19]; the bifurcation invariance helps explain why transition analyses have worked so well [6, 17, 29, 19], in spite of crude estimates of transmission forcing.

Why such bifurcation invariance occurs, and the extent to which it can be generalized to other periodically forced dynamical systems[21, 31, 32], remains to be seen. The SIR model is one of the very simplest nonlinear systems, having only a single quadratic nonlinearity. As a first step towards generalization, we have found that the type of invariance exhibited by the SIR model is also displayed by a periodically forced predator-prey model (\hrefrun:rsif20190202supp1.pdf\hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material, §LABEL:sec:supppredprey). It seems likely that similar invariances occur in many other forced dynamical systems. If a property can be identified that produces the observed invariance, this may lead to important general insights about periodically forced dynamical systems.

In the meantime, dynamical phenomena observed in idealized population models employing sinuoidal forcing[24, 25, 15, 33, 34, 35, 36] can reasonably be assumed to be relevant to epidemics and ecosystems in the real world, and epidemiological transition analyses [6, 17, 29, 26, 19] can be conducted with confidence using imperfectly estimated transmission forcing.

Data accessibility

Code required to replicate the results presented in this paper is available in \hrefrun:rsif20190202supp1.pdf\hrefhttps://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2019.0202&file=rsif20190202supp1.pdfelectronic supplementary material, §LABEL:sec:supp-code.

Authors’ contributions

IP and DJDE designed and performed research, analysed data, and wrote the paper; IP wrote computer programs and performed simulations. Both authors gave final approval for publication and agree to be held accountable for the work performed therein.

Competing interests

We declare we have no competing interests.

Funding

IP was supported by NSERC postgraduate scholarships (CGS-M, PGS-D) and an Ontario Graduate Scholarship. DJDE was supported by NSERC.

Acknowledgements

We are grateful to the Earn and Ellner research groups, Sigal Balshine, Ben Bolker, Stephen Ellner and Steven Strogatz, for discussions and comments.