## 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.

## 1. Introduction

Traditionally deterministic models have formed the bedrock of mathematical epidemiology. The standard approaches such as the susceptible–infected–recovered (SIR) and the susceptible–exposed–infected–recovered (SEIR) models (Anderson & May 1991; Keeling & Rohani 2007) have shaped much of our present day understanding of recurrent epidemics. Stochasticity has long been recognized as an important factor within epidemic modelling, but there is still much debate as to its precise dynamical role and how it can be understood from a theoretical point of view (Bjornstad & Grenfell 2001; Coulson *et al.* 2004).

The first stochastic simulations of an epidemiological system were carried out by Bartlett (1957) on the then new Manchester electronic computer. He was primarily interested in extinction dynamics (Nasell 1999), but also recognized that the simple deterministic models of Hamer (1906) and Soper (1929) would give recurrent dynamics with the addition of stochasticity; in other words, that noise could maintain oscillations by constantly perturbing the system away from its steady state (Aparicio & Solari 2001). The recognition that external forcing is also important to the dynamics of many diseases adds an extra complication (London & York 1973; Bjornstad *et al.* 2002; Altizer *et al.* 2006); trying to understand how this interacts with stochastic effects and thus shapes the dynamics of recurrent epidemics is a challenging problem.

In larger populations, noise has tended to be regarded as a small perturbation to the underlying deterministic dynamics, or as something which can cause switching between different deterministic states (Schwartz 1985; Earn *et al.* 2000). However, the idea that stochasticity plays a more critical role in large populations has been growing over the past few years (Higgins *et al.* 1997; Rohani *et al.* 1999; Alonso *et al.* 2007). For instance, the noise due to demographic stochasticity (noise at the individual level owing to chance events; Nisbet & Gurney 1982) can excite the system's natural frequency, creating a resonance and giving rise to large structured oscillations, even in large populations. The signature of this effect has been detected in real time series data (Bauch & Earn 2003*a*; Grassly *et al.* 2005).

### 1.1. Models

The challenge in epidemic modelling is to form a model which captures the observed dynamics but also elucidates the mechanisms behind them. The most popular approach tends to have two steps: first to create a suitable population-level model (PLM), usually in terms of ordinary differential equations which are deterministic. Next, the corresponding individual-based model (IBM) is formed, and then simulated, to investigate any stochastic properties. Logically, the procedure should be reversed: real populations are finite and the PLM is always an approximation to the underlying IBM. Usually, it is assumed that the PLM will be accurate in the limit of large populations, but what is meant by large must first be defined (Nasell 2002).

Clearly, the IBM, rather than the PLM, should be adopted as the starting point of an investigation. Once an IBM has been defined, it is usually studied using computer simulations. However, simulations are still inferior in at least one respect to the analysis that can be carried out on PLMs: general results valid over a wide range of models and parameters cannot in general be established. In addition, many insights and a deeper understanding can frequently be more readily obtained from analytical studies than from computer simulations. Knowledge of the mathematics required to analyse stochastic models has lagged behind that used to study nonlinear differential equations. Recently, more effort has been put into this area (Nasell 2002; Alonso *et al.* 2007; Kuske *et al.* 2007; Keeling & Ross 2008; Ross *et al.* 2009), although the lack of analytical studies of IBMs has held back the study of stochastic and other effects in models of epidemics.

The modelling of whooping cough provides a good illustration of these points. Whooping cough is a childhood disease (Cherry 2005; Crowcroft & Pebody 2006), which can be described by the SEIR formalism (Hethcote 1997). In England and Wales, before mass vaccination, the time series of case reports shows dynamics which are strongly multiennial, but after vaccination in the 1950s quite regular 3.5–4 year epidemic cycles occur (Fine & Clarkson 1986; Rohani *et al.* 1999). All suitably parametrized deterministic PLMs fail to capture the correct dynamics, but IBMs produce qualitatively correct patterns (Hethcote 1998; Rohani *et al.* 1999). Stochasticity is therefore important, but its precise role has been hard to quantify (Coulson *et al.* 2004). The work which has been carried out has been an amalgam of analytical work on the PLM and simulations of the IBM (Keeling *et al.* 2001; Rohani *et al.* 2002). While this has yielded valuable insights, fundamentally one is left with the problem of interpreting stochastic simulations in terms of deterministic results. More recently, Nguyen & Rohani (2008) have proposed that an SEIR model, with realistically distributed exposed and infectious periods, can capture the qualitative dynamics of whooping cough. The mechanism behind the dynamics is multiple coexisting attractors (Earn *et al.* 2000); noise plays a secondary role of switching the system between these different deterministic states.

### 1.2. Methodology

A method which begins from a mathematical formulation of the IBM, and, through a simple and straightforward procedure, gives analytical results would be a valuable tool in the investigation of the SIR, SEIR and related models. Such a method exists, and was developed by van Kampen (1992) in statistical physics nearly 50 years ago. It has been applied to understand stochastic amplification and other effects in ecological models (Alonso & McKane 2002; McKane & Newman 2005), while other authors have used similar methods in a more mathematically rigorous context to investigate population dynamics (Pollett 1990; Ross 2006). Within epidemiology, it has only been used for investigating relatively simple models (Nasell 1999; Alonso *et al.* 2007; Black *et al.* 2009; Dangerfield *et al.* 2009).

In this paper, we apply the van Kampen method to the SEIR model with distributed exposed and infectious periods (Lloyd 2001; Keeling & Grenfell 2002; Wearing *et al.* 2005), which has recently been employed as a more realistic model of epidemics (Conlan & Grenfell 2007; Ferrari *et al.* 2008; Nguyen & Rohani 2008). The essence of van Kampen's method is an expansion in , where *N* is the number of individuals in the population. In the limit where *N* is very large, stochastic fluctuations become unimportant and the equations describing the model approach the PLM. The method determines these, rather than postulating them phenomenologically, but it also provides a stochastic correction to the equations of the PLM which takes the form of a set of linear stochastic differential equations. It is the analysis of these equations which allows one to accurately interpret the results of simulations and to trace back the effects seen to the mechanisms involved, something which has previously proved difficult to do.

This method is especially powerful when studying the staged SEIR model. By deriving the theoretical power spectrum, we can distinguish the effects of the more realistically distributed periods from the other factors such as immigration and finite size effects. This allows us to give quantitative results for a broad range of parameters very quickly and easily, which would otherwise need large amounts of simulation time to achieve.

The rest of this paper is split into four sections: in §2, we give an overview of the stochastic SEIR model with distributed exposed and infectious periods and our method of analysis. Section 3 details our results and we discuss the impact of the effect of changing the variances of the distributions on the form of the power spectrum. In §4, we use our methods to reinterpret the dynamics of a seasonally forced model of whooping cough studied by Nguyen & Rohani (2008), and we conclude in §5.

## 2. The staged SEIR model

In the standard SEIR model, each individual within the population belongs to one of four classes: susceptible, exposed (infected but not yet infectious), infectious and recovered (Anderson & May 1991; Keeling & Rohani 2007). The usual assumption is that the exposed and infectious periods are exponentially distributed; thus, the probability of moving from one class to another is independent of the time spent in that class. This is unrealistic and in general the probability of recovery will depend quite strongly on the time since infection (Hope-Simpson 1952; Bailey 1956). This can be remedied by splitting the exposed and infectious classes into a number of subclasses or stages, which are traversed sequentially (Andersson & Britton 2000; Lloyd 2001; Keeling & Grenfell 2002; Nguyen & Rohani 2008). If the number of stages is *M* and *L* in the exposed and infectious classes, respectively, with rates *M**σ* and *L**γ* between them, then we obtain gamma-distributed periods (Cox & Miller 1965; Anderson & Watson 1980; Wearing *et al.* 2005).

The gamma distribution is convenient as, with this particular formulation, the mean exposed and infectious periods remain fixed at 1/*σ* and 1/*γ*, while the variances can be changed with the parameters *M* and *L*. If we set *M* = *L* = 1, we regain the exponential distributions of the standard SEIR model. As the number of stages is increased, the gamma distribution becomes more central (smaller variance) and an individual remains on average exposed/infectious for a more constant amount of time. For large values of *M* or *L*, the gamma distribution is approximately normal, and the periods approach fixed duration as *M*, *L* → ∞ (Grossman 1980).

Birth and death rates are set equal to *μ*, and these events are linked, so that the total population *N* remains constant. This allows for the elimination of the variable relating to recovered individuals. Frequently, birth and death processes are assumed to happen at the same rate, but remain as distinct events. This still results in fluctuations in the total population size for finite systems. By linking these events at the stochastic level, the population size remains constant at any system size, so that we can still eliminate the variable relating to recovered individuals. This only has an effect on the dynamics for very small populations.

In the analytical approach, seasonal forcing is not explicitly included (Fine & Clarkson 1982; Aron & Schwartz 1984; Keeling *et al.* 2001) and *β*, the contact rate between individuals, is taken as a constant. In §4, we show how a model with forcing can be analysed using a suitable approximation scheme. We also include a small immigration term *η*, allowing for infectious imports (Hagenaars *et al.* 2004; Alonso *et al.* 2007). The formulation we use corresponds to a commuter model, where a susceptible has the possibility of catching the disease while away from the main population, and then returning (Engbert & Drepper 1994). Vaccination is simple to incorporate, as it can be shown in the deterministic version of the model that vaccinating a proportion *p* of the population at birth scales *β* by (1 − *p*), and this remains true for the stochastic version (Earn *et al.* 2000). The processes which define the stochastic system are given in table 1.

The model can now be investigated by two methods: simulation with Gillespie's (1976) algorithm, and analytically by constructing the master equation corresponding to the processes given in table 1 and performing van Kampen's system-size expansion (van Kampen 1992; Black *et al.* 2009). Both methods allow the power spectrum of the fluctuations of the infective time series to be found. Full technical details of the model and analysis are given in the electronic supplementary material.

## 3. Results

### 3.1. The power spectrum

For orientation, we first show in figure 1 a typical simulation of the stochastic model with parameters corresponding to measles as given by Keeling & Grenfell (2002). We see strong stochastic oscillations caused by the amplification of the demographic stochasticity. The theoretical power spectrum for the same parameters is shown in figure 2. The effect of using the more realistic gamma distributions, as compared with the exponentially distributed model, is to shift the peak frequency and increase the amplitude of the spectrum. These effects were first described for the SIR model with realistic recovery profiles by Simoes *et al.* (2008) and Black *et al.* (2009).

The deterministic version of this model (derived in the limit *N* → ∞) shows weakly damped oscillations tending to a stable fixed point (Hethcote & Tudor 1980; Lloyd 2001; see figure 1). A linear stability analysis yields equations for the eigenvalues which allows for the approximate calculation of the damping time and frequency of these oscillations (Grossman 1980; Lloyd 2001). By contrast, the van Kampen method gives exact analytical expressions (in terms of the frequency, *ω*) for the power spectrum, allowing precise calculation of the peak frequency and amplification. It is possible to give these expressions explicitly in terms of the original parameters of the model, but they would be so unwieldy when given in this way that they would be of little use. They can be constructed from the information given in the electronic supplementary material, but we prefer to represent them graphically where the nature of the oscillations is clearer to see.

Although the eigenvalues of the deterministic system, to a large extent, determine the form (peak frequency and amplitude) of the power spectrum of the fluctuations (see the electronic supplementary material), this is not exact, as the demographic stochasticity can shift the peak frequency away from that predicted from the deterministic equations (Alonso *et al.* 2007; Keeling & Ross 2008). It does, however, form a good approximation for epidemic models of this kind.

The spectrum is normalized so that square root of the amplitude is equal to the root mean square of the fluctuations. Essentially, the spectrum is seen to be very similar to that of the SIR model; the parameter range for which large peaks in the power spectrum occur is also where the deterministic model shows large transients (Alonso *et al.* 2007).

### 3.2. Effect of the variance of the distributions on the power spectrum

In this section, we examine the effect of the parameters *M* and *L* on the form of the power spectrum, which control the variances of the gamma distributions. The behaviour is obviously parameter dependent, but, if we hold the intrinsic parameters {*β*, *σ*, *γ*, *δ*, *η*} fixed, then we can make a number of observations.

Independent of the other parameters, changing *M* has a simple effect on the power spectrum. Increasing *M* (decreasing the variance of the exposed period distribution) increases the amplitude of the spectrum with a negligible change to its frequency. This is illustrated in figure 3. The longer the exposed period, as compared with the infectious period, and the larger *L*, the larger the increase in amplitude. The largest changes happen for the smaller values of *M*, typically less than 10, where the exposed period distribution changes the most.

The effect of *L* on the power spectrum is more complicated and can change both the peak frequency and amplitude. Firstly, independent of *M*, increasing *L* shifts the peak frequency of the spectrum to higher values. The effect on the amplitude depends on whether the infectious period is longer than, approximately equal to or shorter than the exposed period. When the infectious period is longer than the exposed period (*σ* > *γ*), then the results are similar to that of the SIR model with gamma-distributed infectious periods (Lloyd 2001; Black *et al.* 2009). Thus, increasing *L* increases both the amplitude and the peak frequency of the spectrum (figure 4). Again, the largest changes are seen for smaller values of *L*.

If the exposed period is longer than the infectious period (*γ* > *σ*), as with measles (Anderson & May 1991), then there is still an increase in frequency with increasing *L*, but a decrease in amplitude (figure 3). When the two average periods are approximately equivalent (*γ* ≈ *σ*), the behaviour is not as easy to quantify and the amplitude is not necessarily an increasing function of *L* (Lloyd 2001); however, any shifts tend to be small.

## 4. Whooping cough model

In this section, we apply our analysis to the seasonally forced, staged SEIR model, proposed by Nguyen & Rohani (2008) as an explanation for the long-term dynamics of whooping cough before and after mass vaccination. Their model includes seasonal forcing by assuming that the contact rate, *β*(*t*), follows a term-time pattern (Schenzle 1985; Keeling *et al.* 2001),
4.1
where *β*_{0} is the baseline contact rate, *β*_{1} the magnitude of forcing and term(*t*) is a periodic function that switches between 1 during school terms and −1 during holidays. We use the England and Wales term dates as set down by Keeling *et al.* (2001). The problem now explicitly includes time dependence, which is simple to include in simulations with appropriate changes to the Gillespie algorithm (Anderson 2007). It is also possible to include time dependence in the master equation, but this introduces considerable complications in the derivation of the power spectrum, which we will consider in a forthcoming publication. Instead, here we replace *β*(*t*) in the master equation with the average effective *β*, defined as
4.2
where *p*_{s} is the fraction of the time spent in school, as opposed to on holiday; for the term dates, we use *p*_{s} = 0.75. As we show later, simulations confirm this is a good approximation.

All other parameters are taken from Nguyen & Rohani (2008). Before vaccination *R*_{0} = 17, from which we find 〈*β*〉 ≈ *γ**R*_{0} = 1.21. After vaccination, *p*=0.6 so 〈*β*〉^{vac} ≡(1 − *p*)〈*β*〉 = 0.48 (Earn *et al.* 2000). Assuming *β*_{1} = 0.25 and using equation (4.2), we find that *β*_{0} = 1.075 pre-vaccination and . These are lower than 〈*β*〉 because children spend more time in school when *β*(*t*) is higher. The intrinsic parameters, estimated from household incubation data, are 1/*σ* = 8 days, *M* = 1, 1/*γ* = 14 days, *L* = 4 and *μ* = 5.5 × 10^{−5} per day. Note that there is a typographical error in the caption of fig. 6 of Nguyen & Rohani (2008), which should read *n* = 4, *m* = 1. These authors assume a small rate of infectious imports of *δ* = 10 per million per year. We can convert this to our commuter model formulation of immigration by noticing that to a good approximation *η* ≈ (*δ*/*N*)*R*_{0} (Keeling & Rohani 2007, p. 210). This gives *η* = 5 × 10^{−7} per day.

Figure 5 shows the predicted analytical spectra before and after vaccination. The results of using a gamma-distributed infectious period, as compared with a standard exponential, is to shift the peak frequencies, while only slightly increasing the amplitudes of the power spectra. Further increases in *L* have much smaller effects on the power spectrum (figure 4). Figure 5 also shows the peak frequencies from simply carrying out a deterministic analysis. As can be seen for this choice of parameters, demographic stochasticity causes only very small shifts in frequency away from those predicted by the deterministic theory (Alonso *et al.* 2007).

Figure 6 shows power spectra from simulations of the staged SEIR model which include the time-dependent *β*, along with the analytical predictions. The simulations were initialized near to the fixed point of the unforced system, and transients were allowed to damp down before the data for the power spectra were collected. Other initial conditions were tried, but it was found that the dynamics always settled into the annual attractor. As the transients have essentially the same frequency as the stochastic oscillations, it is important to make sure these have damped down, otherwise the stochastic peak in the power spectrum will appear enhanced. The length of each time series was about 700 years; this is certainly very large, but was chosen to facilitate a good resolution on the power spectrum without aliasing (Priestley 1982). The simulations show an annual peak at 1 year, which comes from the macroscopic limit cycle induced by the external seasonal forcing, and a stochastic peak at lower frequencies. In general, there are more peaks at 1/2, 1/3, etc. years, but for clarity these are not shown.

Pre-vaccination, the analytical prediction is good at lower frequencies, but there is enhancement, giving a higher amplitude and very slight shift in frequency. This is in contrast to the post-vaccination model where the analytical curve provides a very good fit, but still somewhat enhanced. The amount of power under each peak depends on the population size. The annual macroscopic peak scales with *N*, while the stochastic part scales with *N*^{1/2}. In the post-vaccination case, the stochastic peak is larger and the annual peak is smaller by a factor of a third than in the pre-vaccination simulations.

In the time series, the dominant period can appear to change. This is due to the superposition of the macroscopic limit cycle and stochastic fluctuations. The power of the stochastic oscillations is spread out over a range of frequencies. Although there is a dominant frequency, at certain times the stochastic oscillations will be of a much longer/shorter frequency, with a much reduced amplitude. This allows the macroscopic signal to be seen more strongly, thus generating the effect of a changing period.

Simulating the system at smaller population sizes results in small shifts to the peak frequencies (Simoes *et al.* 2008), but does not alter the qualitative picture. For the whooping cough model, significant shifts occur at *N* < 10^{6}. Taking a smaller forcing magnitude (e.g. *β*_{1} = 0.15) results in significantly better agreement with the analytical predictions at lower population sizes.

## 5. Discussion

We have applied the van Kampen method to the more realistic staged SEIR epidemic model. The ability to derive the exact power spectrum gives considerable insight, especially with regards to the effects of the parameters *M* and *L* on the form of the spectrum. Both Lloyd (2001) and Keeling & Grenfell (2002) show that gamma-distributed models are much more sensitive to seasonal forcing, in the sense that the bifurcation to biennial dynamics happens at lower forcing magnitudes, than the basic exponential model. Lloyd (2001) attributes this increased sensitivity to a destabilization of the model, but our analysis goes further to show that a distributed infectious period also tends to increase the natural frequency of the system. We conjecture that it is this effect rather than the decreased stability that makes these models more sensitive to seasonal forcing. This should be easy to test within this framework as, with the appropriate choice of parameters, we can tune the size and frequency of the power spectrum. It should be possible to see whether a model with only a distributed exposed period is more sensitive to forcing, as this would be destabilized but with the same natural frequency.

It should also be feasible to include other heterogeneities such as spatial structure (Hagenaars *et al.* 2004; Lloyd & Jansen 2004) or age structure (Schenzle 1985; Hethcote 1997), or further complications such as waning immunity (Grenfell & Anderson 1989). The method does start to break down for smaller values of *N*, where deviations from the predicted spectrum are seen. This is due to the proximity of the fade-out boundary at *I* = 0. This could be overcome by going to the next order in the expansion, which would give corrections to the spectrum which depend on the population size (van Kampen 1992).

The review article by Coulson *et al*. (2004) defines two categories for the role of noise: active and passive. These have been characterized as follows (Nguyen & Rohani 2008). Active noise is where ‘stochasticity interacts with the nonlinearity in the deterministic clockwork producing patterns that cannot result from either factor alone’. Passive noise is where ‘stochasticity influences the transition among different deterministic states’. Previous work on whooping cough models formulated with the basic SEIR model (Rohani *et al.* 1999, 2002; Keeling *et al.* 2001) concludes that the noise is of the ‘active’ type and plays a major role in the observed dynamics. This is largely due to the failure of the deterministic model to predict anything other than annual dynamics. In contrast, Nguyen & Rohani (2008) conclude that a realistic deterministic staged SEIR model can account for the qualitative patterns, due to the existence of multiple coexisting attractors, which are not present in the basic SEIR model. This explanation reverts to a passive interpretation of noise, which kicks the system between deterministic states. All these models have been shown to be consistent with the observed whooping cough dynamics, but a persuasive explanation for the mechanisms has been lacking.

Our approach, which allows us to be much more quantitative, can encompass both older SEIR results and the newer staged SEIR model of Nguyen & Rohani (2008). For the unforced system, there are strong stochastic oscillations about a macroscopic fixed point. External forcing creates a macroscopic annual limit cycle with stochastic oscillations about it. For the case of whooping cough, the oscillations about the limit cycle are very similar to those of the unforced model and the two components of the power spectrum, the macroscopic (from the limit cycle) and the stochastic, can be viewed as approximately independent of each other. Using the staged SEIR model changes this picture only in that it shifts the frequency of the stochastic oscillations; it does not change the basic nature of the dynamics.

Our results argue against multiple coexisting attractors in the stochastic whooping cough system. Preliminary results show that the bifurcation diagrams of Nguyen & Rohani (2008) are not robust to the inclusion of immigration. Even small amounts as used in the simulations depress the onset of period 3 attractors which only appear at significantly higher magnitudes of forcing. This effect, whereby immigration can destroy longer period attractors, has been reported before (Engbert & Drepper 1994; Ferguson *et al.* 1996; Alonso *et al.* 2007). One can speculate that the lack of immigration in the deterministic analysis of Nguyen & Rohani (2008) might explain their findings; alternatively there may be other stable solutions of the deterministic equations, but that the noise is never strong enough to move the system towards them.

Our interpretation, in many ways, echoes the findings of Bauch & Earn (2003*b*). They find resonant and non-resonant peaks in their power spectra of real time series which correspond to what we term the stochastic and macroscopic peaks. They also show that the natural damping period of the transients near the annual attractor is well approximated by the damping period from the unforced deterministic model (Bauch & Earn 2003*a*). Our analysis goes further to show that the power spectra are also well approximated by their unforced counterparts.

At the heart of the van Kampen approach is a splitting of the dynamics into two components: a macroscopic part and a stochastic correction to this. In the case of the seasonally forced model, the macroscopic part is a limit cycle and the stochastic corrections are oscillations about this limit cycle. Firstly, we must ask whether the macroscopic part is really a limit cycle, that is, whether the stable solution of the deterministic equations is a limit cycle or whether the multiple attractors, unstable limit cycles, invasion orbits, etc. found or postulated by previous authors have a role or appear in the fully stochastic model. One has to begin from a finite-event-driven model to study these questions, but the problem of disentangling these effects in the time series and power spectra remains. This is where the systematic approach of the system-size expansion is invaluable, because it allows us to separate out the different structures which are seen in the power spectrum. So, for instance, one would expect to see a 3 year peak in the power spectrum if there was a stable 3 year attractor in the model. We also note that the van Kampan expansion frees us from the need to introduce concepts such as invasion orbits (Rand & Wilson 1991; Rohani *et al.* 2002), which reflect the local geometry around an attractor. Such structure is automatically encoded into the van Kampen expansion and so there is no need to explicitly consider it.

It is surprising that the predictions from the unforced model remain so good in the presence of forcing—there is no *a priori* reason to expect this to be true. There are some deviations when the forcing is large, but we expect these would be explained by a theoretical treatment in which the time dependence of *β*(*t*) given by equation (4.1) was not replaced by the effective value (equation (4.2)). This more complete theory would tell us whether these deviations are from stochastic oscillations about the limit cycle or due to some macroscopic effect such as the presence of an unstable multiennial limit cycle (Boland *et al.* 2009). The fact that we get such a good fit using approximation (4.2) suggests the former explanation. We do not expect the current approach to provide a good fit to forced measles spectra, as in that case there is almost certainly a more complex macroscopic dynamics, namely a biennial limit cycle, to take into account. A completely systematic approach to stochastic perturbations of forced systems, which is currently being developed, should be able to elucidate these points. However, we believe that the good agreement which we have found in this paper between the simulations of the fully time-dependent model and the van Kampen expansion using the effective contact rate (4.2) shows the essential correctness of our approach. It also clarifies several aspects of the dynamics of epidemic models of the SEIR type.

## Acknowledgements

A.B. would like to thank the EPSRC for the award of a postgraduate grant.

## Footnotes

- Received November 21, 2009.
- Accepted January 25, 2010.

- © 2010 The Royal Society