Stochasticity in staged models of epidemics: quantifying the dynamics of whooping cough

Abstract

Although many stochastic models can accurately capture the qualitative epidemic patterns of many childhood diseases, there is still considerable discussion concerning the basic mechanisms generating these patterns; much of this stems from the use of deterministic models to try to understand stochastic simulations. We argue that a systematic method of analysing models of the spread of childhood diseases is required in order to consistently separate out the effects of demographic stochasticity, external forcing and modelling choices. Such a technique is provided by formulating the models as master equations and using the van Kampen system-size expansion to provide analytical expressions for quantities of interest. We apply this method to the susceptible-exposed-infected-recovered (SEIR) model with distributed exposed and infectious periods and calculate the form that stochastic oscillations take on in terms of the model parameters. With the use of a suitable approximation, we apply the formalism to analyse a model of whooping cough which includes seasonal forcing. This allows us to more accurately interpret the results of simulations and to make a more quantitative assessment of the predictions of the model. We show that the observed dynamics are a result of a macroscopic limit cycle induced by the external forcing and resonant stochastic oscillations about this cycle.

Figure 1.

A typical stochastic realization of the staged SEIR model. The dashed curve shows the deterministic result; see the electronic supplementary material for these equations. Parameters are β = 3.4 d⁻¹, σ = 1/8 d, M = 8, γ = 1/5 d, L = 5, μ = 5.5 × 10⁻⁵ d⁻¹, η = 10⁻⁶ d⁻¹ and N = 10⁶.

Figure 2.

Theoretical power spectrum (solid curve) for the staged SEIR model parametrized for measles, with results from simulations (open circles). The dashed curve is the standard SEIR result, assuming exponentially distributed exposed and infectious periods. The shift in period is from 1.9 to 2.1 years. Transient periods are discarded before evaluating the numerical spectrum; parameters are as in figure 1.

Figure 3.

The (a) amplitude and (b) peak frequency of the analytical power spectrum as a function of M, the exposed period variance parameter. Increasing M (decreasing the variance of the exposed period distribution) increases the amplitude of the spectrum, while having negligible effect on the peak frequency. The intrinsic parameters are typical of measles as given in figure 1. In this example, the exposed period is longer than the infectious period so the spectra with larger L have smaller amplitudes. Filled diamond, L = 10; filled square, L = 3; filled circle, L = 1.

Figure 4.

The (a) amplitude and (b) peak frequency of the analytical power spectrum as a function of L, where the infectious period is longer than the exposed period. Increasing L leads to an increase in the peak frequency and amplitude of the spectrum. The change in frequency is independent of M. Parameters are typical of whooping cough: β = 1.2 d⁻¹, σ = 1/8 d, γ = 1/14 d, μ = 5.5 × 10⁻⁵ d⁻¹ and η = 10⁻⁶ d⁻¹. Filled diamond, M = 10; filled square, M = 3; filled circle, M = 1.

Figure 5.

Comparison of the analytical power spectra for whooping cough using the SEIR (dashed curves) and staged SEIR (solid curves) models. Parameters are given in the main text. Using the staged version shifts the peak periods from 2.7 to 2.4 years pre-vaccination and from 4.5 to 4 years post-vaccination. The vertical lines show the peak frequencies from an analysis of just the deterministic model.

Figure 6.

Theoretical and numerical power spectra for the whooping cough model: (a) pre-vaccination and (b) post-vaccination. The filled circles are the average of 2000 realizations of the stochastic simulation, including seasonal forcing. The dark curves are the theoretical predictions. The forcing slightly shifts the frequency and amplitude of the stochastic peaks from the analytical prediction, but is otherwise good. The value of N is 5 × 10⁶, with other parameters being given in the main text. The spectra are normalized to N^1/2 as we are predominantly interested in the stochastic parts and comparing with analytical calculations. The annual peaks in the numerical spectra are cropped for clarity in the comparison.