The role of stochasticity and its interplay with nonlinearity are central current issues in studies of the complex population patterns observed in nature, including the pronounced oscillations of wildlife and infectious diseases. The dynamics of childhood diseases have provided influential case studies to develop and test mathematical models with practical application to epidemiology, but are also of general relevance to the central question of whether simple nonlinear systems can explain and predict the complex temporal and spatial patterns observed in nature outside laboratory conditions. Here, we present a stochastic theory for the major dynamical transitions in epidemics from regular to irregular cycles, which relies on the discrete nature of disease transmission and low spatial coupling. The full spectrum of stochastic fluctuations is derived analytically to show how the amplification of noise varies across these transitions. The changes in noise amplification and coherence appear robust to seasonal forcing, questioning the role of seasonality and its interplay with deterministic components of epidemiological models. Childhood diseases are shown to fall into regions of parameter space of high noise amplification. This type of ‘endogenous’ stochastic resonance may be relevant to population oscillations in nonlinear ecological systems in general.
A wide range of temporal behaviours, including annual, biennial, multi-annual and irregular oscillations, have been described for the extensive time-series data on childhood diseases, such as measles and whooping cough (Hethcote & Levin 1989; Earn et al. 2000; May 2000; Grenfell et al. 2002). One key test of these models has been their ability to explain the observed transition from regular periodic cycles to irregular aperiodic fluctuations with the advent of mass vaccination and changes in birth rates. The standard paradigm, emerging from a largely deterministic framework, is that seasonal forcing and the endogenous frequency of disease oscillations can dynamically interact to produce large fluctuations and different types of dynamics (Earn et al. 2000). In simpler terms, the dynamics of disease would be similar to those of a seasonally forced pendulum. In this view, away from questions of extinction and associated critical population sizes (e.g. Bolker & Grenfell 1995), noise plays a role, but one that is typically understood through its interaction with the deterministic trajectories of the system. Specifically, the known effects of noise include sustaining damped cycles that would otherwise approach a stable equilibrium (Barlett 1957), allowing the system to remain in a transient regime near an attractor or a repellor (Rand & Wilson, 1991; Sidorowich 1992; Bauch & Earn, 2003a,b), or causing it to shift between different attractors, such as cycles of different periods, leading to irregular behaviour (Earn et al. 2000; Keeling et al. 2001). In these dynamics, the population-level deterministic framework, including transients near asymptotic behaviours, provides the blueprint upon which the effect of stochastic fluctuations can be understood. Here, we propose a more fundamental role for noise in the patterns of disease cycles, with a theory that suggests a simple explanation for the major dynamical transitions in epidemics from regular to irregular fluctuations. This explanation relies on the amplification of demographic stochasticity (McKane & Newman 2005) arising from the discrete nature of individuals (Durrett & Levin, 1994).
Seasonal forcing and nonlinear dynamics do interact in complex ways (Keeling et al. 2001), producing, in some cases (Dushoff et al. 2004), resonant oscillations of significant amplitude. However, the importance of this resonant effect will depend on both the magnitude of seasonal forcing and the endogenous sensitivity of the system to amplify fluctuations. While the former has been intensively investigated (Olsen & Schaffer, 1990; Keeling et al. 2001; Bjørnstad et al. 2002; Grenfell et al. 2002; Bauch & Earn, 2003a,b), the latter has received comparatively less attention. We describe analytically the fluctuations produced by demographic stochasticity in epidemic models, and consider specifically the susceptible–infectious–recovered (SIR) model for infectious diseases with immigration. From the stochastic model, we derive an analytical expression for the power spectral densities (PSD) of the number of infected and susceptible individuals, describing how the temporal variability of these quantities is distributed over different frequencies. Armed with this expression, we show that (i) qualitative changes in disease cycles are consistent with changes in the way disease dynamics amplify stochasticity, (ii) childhood diseases fall in a region of parameter space of strong noise amplification, leading to highly coherent fluctuations and (iii) that the dominant period of the oscillations need not be that of the oscillatory approach to equilibrium in the deterministic system, as typically assumed. The proposed explanation is robust to low levels of spatial coupling, a realistic process in disease systems. It is also robust to a significant degree of seasonal forcing, suggesting a need to re-evaluate the role of this factor in our understanding of childhood diseases. The theory is general and can be applied to study the sensitivity to demographic stochasticity of other ecological systems.
2.1 The open susceptible–infectious–recovered model
Inspired by the previous work (Nisbet & Gurney 1982; McKane & Newman 2005), our results are derived from a detailed analysis of a simple individual-based stochastic SIR model, in which individuals are discrete (Durrett & Levin, 1994) and belong to one of the three possible classes (S, I and R). The number of individuals belonging to each of these classes changes in time according to three different processes—infection, recovery and demography (birth and death). The transitions associated with these events are illustrated graphically in figure 1 (see the electronic supplementary material for details). Susceptible individuals can acquire infection from two routes: internal frequency-dependent transmission (Bjørnstad et al. 2002) at a rate βI/N (where N is the maximum number of individuals that the system can support, or more practically, the city size) and external infection at a small constant rate η. This is the simplest possible way to handle spatial coupling (Xia et al. 2004). Infected individuals recover at a constant rate γ and, therefore, 1/γ is the average time spent in the infectious class. Regardless of their class, individuals die at a constant rate. To consider a constant community size, births and deaths are tightly coupled, so that an ‘empty site’ created by a death is instantly replenished by the birth of a susceptible individual. The death (or birth) rate δ therefore specifies an average turnover time 1/δ or individual lifetime of susceptibles.
2.2 The large N expansion
The stochastic dynamics are fully described by a multivariate master equation (van Kampen 1992). At the leading order, a formal large N expansion of this equation gives rise to the so-called macroscopic law, which corresponds to a deterministic description of the temporal evolution of the macroscopic variables, i.e. the fraction of susceptible and infected individuals in a population of constant size (ϕ(t) and ψ(t), respectively),(2.1)
A description of the stochastic fluctuations of the system requires consideration of higher-order terms in this expansion. In particular, a very good approximation is obtained only if the next-to-leading order is considered. In this way, we obtain a set of Langevin equations for the temporal evolution of the normalized fluctuations of susceptible and infectious individuals around equilibrium values (x and y, respectively),(2.2)where τ is a scaled time and where ηi are the white noises with zero mean and their cross-correlation structure is determined by the expansion as well. The elements aij are the entries of the Jacobian matrix evaluated at the stable fixed point. By Fourier transforming these equations, we are able to analytically calculate the PSD corresponding to the normalized fluctuations, which, therefore, no longer depend on the community size N. The general form of these PSD is given by(2.3)for susceptibles or infectives. The parameter α in the numerator of the PSD is a function of both aij, the entries of the Jacobian matrix, and Bij, those of the noise covariance matrix. Thus, both the disease rates and the noise correlation structure influence this parameter. The parameters in the denominator are the two frequencies Ω and Γ, which depend exclusively on the aij entries, and can therefore be determined from the linearization close to equilibrium. We note that the expression of the PSD differs for susceptible and infectious individuals in the value of α and B. For a comprehensive and detailed description of this difference and the parameters arising in the PSD, see the electronic supplementary material.
As in Nisbet & Gurney (1982), this approach is based on three points: a formal treatment of demographic stochasticity; the presence of stationary stable fixed points; and the power spectral analysis. Here, we apply this three-step approach to infectious diseases for the first time using a powerful technique based on a large N expansion due to van Kampen (1992). It is important to note that this expansion is not limited to small perturbations around equilibrium values, as deterministic perturbations are classically studied. In fact, it is not performed for such small perturbations, but with the inverse of the system size N as our expansion parameter. In this way, by maintaining the essential nonlinear dynamics of disease transmission, we can evaluate the effect of finite community sizes along with the discrete nature of individual interactions on oscillatory patterns. The linearity of resulting equations emerges from the expansion itself for this particular case, and not from linearizing the system close to equilibrium. The stochastic model and the steps of the analysis are described in more detail in the electronic supplementary material.
2.3 Overall amplification and coherence: analytical formulae
In general, the total spectral power of a time-series can be calculated by integrating its PSD over all frequencies. Fortunately, it turns out that the PSD P(ω) given in equation (2.3) can be integrated analytically using standard methods,(2.4)This quantity has been called overall amplification, mapped into the parameter space in figure 4a, and labelled A0 in figure 3 (legend). The properties of the PSD and the way fluctuations have been normalized ensure that this quantity is equal to the mean squared deviation of the original time-series from the equilibrium values(2.5)where N(t) and N* are the time-dependent time-series variable and its equilibrium value, respectively. This measure depends only on the model parameters and quantifies the model's intrinsic sensitivity to fluctuations or, more precisely, the ability of the model to amplify and sustain fluctuations produced by demographic stochasticity.
A measure of how well-structured stochastic oscillations are around the dominant spectral frequency can be given by integrating equation (2.3) between definite limits around the spectral peak frequency. This quantity has been labelled Ap in figure 3,(2.6)where a factor 2 has been included to take into account that the power spectra are normalized differently in the numerical work (the range of ω is assumed to be (0,∞), rather than (−∞,∞)). It turns out that the integral in equation (2.6) can also be evaluated analytically. In fact, the definition we have adopted for coherence, c=Ap/A0, which has been mapped into the parameter space in figure 4c in percentage values, takes the form of the following expression:(2.7)where F(ω) is given by(2.8)and a, and are given by(2.9)
As a check, we can see that if we let and in equation (2.7), then c→1, as expected.
It is well known that the deterministic counterpart of this model shows damped oscillations approaching a stable fixed point (figure 1b), and it is generally assumed that the natural period of the associated stochastic oscillations coincides with the period of the deterministic damped oscillations. To determine analytically the periodicity of the stochastic system, we derived the PSD of the fluctuations around equilibrium using a large N expansion method due to van Kampen (1992) (see the electronic supplementary material). The PSD gives a complete representation of how the variance of the fluctuation is distributed among different frequencies. From the expression for the PSD, we observe that demographic noise can shift the deterministic period of the damped oscillations. Specifically, the stochastic endogenous dominant period—the peak of the PSD—does not necessarily correspond to that of the decaying deterministic oscillations, and, therefore, it cannot generically be obtained from the deterministic system itself as is typically assumed. Conditions are derived for this correspondence to hold (see the electronic supplementary material).
The derived PSD further allows us to examine how noise amplification varies with changes in the transmission rate (β). The series of bifurcations induced by changes in this important parameter have been considered before in deterministic models, in order to address the pronounced change in the character of disease cycles with the advent of mass vaccination and also changing birth rates (Earn et al. 2000). The dynamics of measles have been described by a transition from regular annual or biennial cycles to more irregular, less coherent, fluctuations of smaller amplitude (see fig. 2 in Earn et al. (2000) and fig. 1 in Rohani et al. (1999)). As shown by Earn et al. (2000), changes in vaccination and birth rates can be studied by considering instead changes in an effective transmission rate, . In our model, it can be easily shown that a proportion, p, of vaccination changes both contact rates, and therefore β, and the external infection parameter, η, but for simplicity we first consider changes only in the former and fix the latter to a small constant value (see the electronic supplementary material). Figure 2a shows that for parameters typical of measles, a lower transmission rate leads to a longer period and reduces the overall amplitude of the oscillations. Furthermore, the PSD curve becomes flatter as β decreases, indicating that more frequencies are involved in the stochastic fluctuations, and that the overall variance is more evenly distributed among these frequencies. As a result, the dominant period increases while the cycles become less coherent, more irregular and less pronounced in their amplitude. Numerical calculation of the PSD from an ensemble of simulations confirms this analytical result (figure 2a compared with b). However, this general picture is parameter dependent. For example, for whooping cough, the change in period and the decrease in amplitude are clear but not the broadening of the PSD curve (see fig. 1 in Bauch & Earn (2003b)).
Our analysis so far ignores the seasonality of transmission rates. In childhood diseases, seasonality in transmission is assumed to follow from alternating school–holiday terms. High contact rates characterize school terms, while low contact rates would occur during the holidays. Figure 2c shows numerically that the qualitative changes described with the analytical PSD for the non-seasonal case are comparable with those obtained for the seasonally forced system, even for moderate values of this forcing (β1=0.25). The only difference is the addition of the expected annual peak with seasonal transmission. The pattern is also robust to lower values of the external infection rate and to decreases of both external and internal transmission, as would result from vaccination.
To examine the strength of noise amplification, we introduce two quantitative measures, derived from the predicted PSD, that characterize respectively the total amplification and the coherence of the stochastic fluctuations (figure 3; see also §2). The parameter space of the model depends on three independent dimensionless parameters. For clarity, we have mapped these two measures on a plane to show how sensitivity to stochastic amplification and coherence change depending on the values of and a measure of the degree of acuteness of the infection—the relative duration of the infectious phase in relation to the average lifetime of individuals (figure 4). We observe that typical model parameters for different childhood diseases lie within a region of strong stochastic amplification and, particularly, strong coherence close to an instability boundary. This boundary separates the region where the disease is self-maintained—the endemic or ‘source’ phase—from the region where the disease is maintained at very low average incidence values by external infections—the epidemic or ‘sink’ phase. The role of the third parameter, η/γ, is also qualitatively important, increasingly as the average infectious period of the disease lengthens. Moreover, external infections are responsible for the hyperbolic shape of this instability boundary in the parameter space. We may approach this boundary not only by decreasing the average infectious period, which reduces R0, but also, as a consequence of this hyperbolic shape, by keeping R0 constant and increasing the average lifespan. In fact, the inclusion of external infection defines a boundary in the parameter space conceptually different from the well-established concept of a critical community size. This threshold is defined as the community size above which endemic dynamics is self-maintained. By holding the intrinsic parameters of the disease constant, below the critical community size, epidemic outbreaks and fade-outs are observed due to repeated stochastic disease extinction and re-introduction. Here, regardless of the community size, a boundary in the parameter space separates a sink phase from a source phase. If the system is in the sink phase, the disease has to be necessarily maintained by constant external infections, and, thus, it will be prone to present epidemic outbreaks and fade-outs, while if the system is in the source phase, the dynamics is characterized by endemic strong fluctuations without fade-outs. On our phase diagram, community size would introduce a third dimension, and parameter-dependent critical community sizes could and should be defined for each point of the parameter space.
We have described analytically how noise amplification by the nonlinear dynamics of an infectious disease with permanent immunity changes when a key epidemiological parameter is varied. This pattern is consistent with the qualitative transitions observed for childhood diseases with the advent of vaccination. Noise amplification provides a possible explanation for qualitative changes from regular to irregular oscillations of lower amplitude (see fig. 2 in Earn et al. (2000) or fig. 1 in Keeling et al. (2001)). Another explanation was proposed earlier based on the bifurcation analysis of the deterministic susceptible–exposed–infections–recovered (SEIR) model with seasonal forcing (Earn et al. 2000; Bauch & Earn, 2003b). Specifically, at low transmission rates (or high vaccination), irregular cycles in measles would result from coexisting stable limit cycles with intertwined basins of attraction. We find, however, that this bifurcation diagram is not robust to external infection (and, therefore, spatial coupling between populations); even for very low rates of external transmission, the multiple coexisting cycles are no longer present (see the electronic supplementary material). This observation supports the conclusion that for low transmission rates (β), stochasticity and spatial coupling would play a leading role, while for high transmission rates, stochasticity becomes less important and internal transmission becomes dominant, making a deterministic description and the assumption of a closed system sufficient. In this latter case, our approach reproduces previous results.
The importance of noise to disease dynamics has been recognized since the pioneering work of Bartlett (1957) in phenomena such as the critical population size for disease extinction (e.g. Bolker & Grenfell 1995), the spatio-temporal modelling of coupled cities of different size (Bjørnstad et al. 2002; Xia et al. 2004) and the role of transients in nonlinear systems (Hastings 2001, 2004; Keeling et al. 2001; Rohani et al. 2002; Bauch & Earn, 2003b). Our results extend the fundamental role of noise to endemic cycles away from patterns of extinction, small population sizes and proximity to stable or unstable trajectories of the seasonally forced deterministic system. They raise questions concerning the role of seasonality itself (Hethcote 1998), an area of relevance beyond childhood diseases, as seen in the recent result that a small level of seasonality can resonate with the nonlinear dynamics of influenza (Dushoff et al. 2004).
Our analytical results pertain to disease dynamics in the absence of seasonality. The robustness of the PSD, and hence of the pattern of noise amplification when seasonality is added, extends these results to the seasonal case. The generality of this correspondence in the PSD was not demonstrated analytically; therefore, it raises interesting open questions on the relationship between the dynamics of the forced and unforced system, particularly in the presence of noise (Bauch & Earn, 2003b). The dynamic transitions between the two equilibria of the unforced system for ‘high’ and ‘low’ transmission, respectively (Keeling et al. 2001), have provided a powerful approach to the understanding of school-term forcing, its interplay with noise and extensions to more general forms of seasonality. With this approach, Keeling et al. (2001) have argued that when seasonality is weak and transmission levels are low, stochastic processes should dominate the dynamics. A similar approach may prove fruitful to examine the implications of our results for the seasonal case, perhaps beyond weak seasonality. It may also provide the means to understand how the fully stochastic perspective we have developed here fits with the existing explanations based on the study of transient fluctuations close to oscillatory attractors of seasonally forced systems. This is important for a full understanding of the role of seasonality and noise, and has implications for the development and fitting of nonlinear time-series models to disease data (e.g. Bjørnstad et al. 2002; de Valpine & Hastings 2002; Grenfell et al. 2002).
The patterns of noise amplification we have described here are clearly not restricted to small community sizes; they are also critical for low transmission rates or, equivalently, high vaccination rates in large cities (Jansen et al. 2003). While demographic noise is indeed negligible when considering population fractions for large communities, the variations for actual population numbers can be sizable; fluctuations in the number of cases can exhibit large amplitudes in large cities, and the extremely low values during the troughs of the cycles make externally driven infections potentially critical. Many of the steps of our analysis can also be applied to understand the role of environmental noise (e.g. Nisbet & Gurney 1982), but in this case, the starting point is already a set of deterministic equations with specific assumptions on the noise distribution and associated power spectrum.
By contrast, here, we have started with stochastic processes at the individual level. In this case, the form of demographic noise, whose correlation structure reflects the nonlinear nature of the interaction, is not assumed on an ad hoc basis (see equations (2.2)), but naturally emerges from the expansion we have performed. We have used van Kampen's (1992) method to approximate the master equation of a stochastic system only to next-to-leading order. High-order contributions yield a nonlinear Fokker–Planck equation. Our Langevin equations (see equations (2.2)) are just equivalent to the linear equation we obtain when these higher-order corrections are neglected (see the electronic supplementary material for further details). In this linear regime, the stochastic variables are well described by their averages and respective variances. We have reported and quantified an increase in the variance of the noise when we approach a critical threshold. However, it is known that this description strictly breaks down close enough to critical points. In order to study the dynamics very close to (or crossing) the source–sink critical boundary described here, more and more terms in classical van Kampen's expansion are needed or other methods become more appropriate (Dekker 1980a,b).
By accounting for stochasticity and external infections, we have uncovered the result that some infectious diseases lie in a region of the parameter space close to an instability boundary characterized by well-structured strong oscillations, which raises interesting open questions concerning the coevolutionary forces acting on host–parasite systems. This type of ‘endogenous’ stochastic resonance (McKane & Newman 2005) is quantitatively described here for infectious diseases for the first time. These phenomena may also be relevant to population oscillations in nonlinear ecological systems in general.
D.A. would like to thank I. Pagonabarraga for his initial advice and A. King for his very constructive comments. This work benefited from discussions and was conducted in part, at the Working Group on Seasonality and Population Dynamics of Infectious Diseases supported by the National Center of Ecological Analysis and Synthesis, a centre funded by the NSF, and the UC at Santa Barbara. Further support was provided by the James S. McDonnell Foundation through a Centennial Fellowship in Global and Complex Systems to M.P.