## Abstract

Rubella is a completely immunizing and mild infection in children. Understanding its behaviour is of considerable public health importance because of congenital rubella syndrome, which results from infection with rubella during early pregnancy and may entail a variety of birth defects. The recurrent dynamics of rubella are relatively poorly resolved, and appear to show considerable diversity globally. Here, we investigate the behaviour of a stochastic seasonally forced susceptible–infected–recovered model to characterize the determinants of these dynamics and illustrate patterns by comparison with measles. We perform a systematic analysis of spectra of stochastic fluctuations around stable attractors of the corresponding deterministic model and compare them with spectra from full stochastic simulations in large populations. This approach allows us to quantify the effects of demographic stochasticity and to give a coherent picture of measles and rubella dynamics, explaining essential differences in the recurrent patterns exhibited by these diseases. We discuss the implications of our findings in the context of vaccination and changing birth rates as well as the persistence of these two childhood infections.

## 1. Introduction

Rubella is a completely immunizing, directly transmitted viral infection, generally presenting as a mild and potentially even asymptomatic childhood disease [1]. As a result, rubella tends to be underreported, and its recurrent dynamics are fairly poorly characterized. Nevertheless, because infection during early pregnancy may cause spontaneous abortion or congenital rubella syndrome (CRS), which may entail a variety of birth defects [2], understanding the dynamics of rubella is of considerable public health importance. Dynamical features of rubella may alter the CRS burden via their effects on the average age of infection. Episodic dynamics may increase the average age of infection, as the intervals between larger outbreaks provide the opportunity for individuals to age into later age classes [3,4]. Likewise, local extinction dynamics can allow individuals to remain susceptible as they age into their childbearing years [5,6], resulting in the potential for a considerable CRS burden once rubella is reintroduced.

Empirically, rubella seems to be linked to either (i) annual dynamics, as in Mexico [7], Peru [5] or parts of Africa [8,9]; (ii) spiky dynamics, as in Canada [10]; and (iii) some hint at multi-annual regularity, as in Japan [11], England and Wales [12], and various European countries [13]. In figure 1, we show three time series that represent the range of observed rubella dynamics. Spectral analyses of time series are particularly useful for understanding temporal patterns exhibited by different data [14,15]. The characteristic feature of rubella spectra is an annual peak at 1 year and a multi-annual peak at 5–6 years exhibited by all data in figure 1. Rubella also seems to experience regular fade-outs [7], which is of great epidemiological importance, particularly in the context of increasing global control efforts. The propensity for stochastic extinction is characterized by the critical community size (CCS), below which the infection tends to go extinct in epidemic troughs. Analyses of dynamics in Mexico and Peru suggest a CCS of over 10^{6} for rubella [5,7].

Measles provides a natural comparison for rubella, as it is another viral childhood disease with a very similar life cycle (particularly, direct transmission). In addition, it is perhaps the most extensively studied of the childhood infections, and its dynamics are very well understood [3,16–27]. Before the start of vaccination in England and Wales, both biennial dynamics (e.g. in London) and annual dynamics (e.g. in Liverpool) were observed [18,20,24,26]. The underlying driver of this variability has been identified as differences in birth rate, combined with annual seasonality in transmission driven by school term times [18,19,26,28]. In sub-Saharan Africa, chaotic dynamics have been shown to result from a very high birth rate, combined with extreme seasonal forcing [21]. Both highly irregular dynamics [16,20] (e.g. following vaccination in England and Wales) and triennial dynamics [29] (e.g. in Baltimore between 1928 and 1935) have also been reported. The spectral analyses of measles data exhibiting the described dynamics can be found now in the standard textbook [12]. The CCS of measles is rather smaller than that of rubella, estimated by Bjørnstad *et al*. [18] to be between 3 × 10^{5} and 5 × 10^{5} for England and Wales.

The two key ingredients underlying models of childhood diseases such as rubella and measles are (i) seasonality in transmission owing to schooling patterns and (ii) demographic stochasticity arising from the discrete nature of population [26,28,30,31]. Although various approaches have been used to understand the dynamics of rubella [10,32–34], most of the analyses have been essentially deterministic. Keeling *et al*. [32] considered a term-time forced susceptible–infected–recovered (SIR) model and compared its dynamics with rubella data in Copenhagen. From this, they concluded that the dynamics of rubella may result from switching between two cyclic attractors (annual and multi-annual limit cycles) of the deterministic model. Although the deterministic analysis they present is comprehensive, there is only a limited amount of evidence to suggest that the switching will occur in contexts that include demographic stochasticity. In particular, in this study [32], stochasticity was introduced into the model as multiplicative noise of arbitrary amplitude instead of using, for instance, standard stochastic simulations based on the Gillespie algorithm (for unforced models) [35] and its extensions (for seasonally forced models) [36]. Such simulations produce exact realizations of the stochastic process, whose full dynamics is given by the solution of the master equation as described in §2.2. For large populations, the master equation is approximated by the deterministic model with additive noise [37,38].

Bauch & Earn [10] studied a term-time forced susceptible–exposed–infected–recovered model and showed that frequencies obtained from the linear perturbation analysis of the deterministic model are in good agreement with positions of the peaks in spectra of data records of various childhood infections. The application of this approach to rubella data for Canada predicted two distinct peaks at 1 and 5.1 years, close to what we see in figure 1*f*. With the exception of [10], where stochastic simulations for Canada parameter values were also performed, there has been no work on rubella using a fully stochastic approach dealing with demographic stochasticity.

Here, we use this approach to characterize different rubella dynamics and illustrate patterns by comparison with measles. To this end, we perform the theoretical analysis of spectra of stochastic fluctuations around stable attractors of a seasonally forced deterministic SIR model and compare them with spectra obtained from full stochastic simulations based on a modification [36] of the algorithm by Gillespie [35]. The mathematical techniques used in this study have been developed for ecological and epidemiological models [37,39,40] and applied to model temporal patterns of measles and pertussis [38,41,42]. The picture that emerges to explain rubella dynamics is close to that proposed in reference [10] but goes beyond it, because our analysis allows us to obtain the full structure of a spectrum (as opposed to the deterministic analysis of Bauch & Earn [10] where only frequencies of the spectral peaks could be predicted). By introducing key spectral statistics (described below), we systematically investigate how the dominant period, amplitude and coherence of stochastic fluctuations change across a broad range of epidemiological parameters. We then discuss the implications of our analysis in the context of changing birth rates and vaccination levels, as well as their implications for the persistence of measles and rubella.

## 2. Methods

### 2.1. Model

The individual-based stochastic model, we explore in this paper, follows a standard seasonally forced SIR structure [12,43]. At any time *t*, it consists of a discrete population of constant size *N* divided into compartments of susceptible, *S*(*t*), infected, *I*(*t*) and recovered, *R*(*t*), individuals. Susceptible individuals become infected (and infectious) at a frequency-dependent rate *β*(*t*)*I*(*t*)/*N*, where *β*(*t*) is a seasonally varying transmission rate. For childhood diseases, *β*(*t*) captures an increase in the number of contacts between school children during terms with respect to holidays [28]. Both term-time and sinusoidal forcing have been used to model these changes [10,20,32,37,38,44,45]. Previous studies [10,20] have shown that the type of forcing is not crucial for the essential dynamic structure (the bifurcation diagram) of the sinusoidally and term-time forced models if the seasonal forcing amplitude is adjusted appropriately. Specifically, for rubella, the same dynamics is obtained if the term-time forcing amplitude is 2.7 times larger than the sinusoidal forcing amplitude (see §3.2.1). We therefore focus on a sinusoidally forced *β*(*t*) = *β*_{0}(1 + *ε* cos 2*π**t*), where *β*_{0} is the average transmission rate and *ε* is the amplitude of seasonal forcing, and confirm later that the dynamic temporal patterns observed in the simulations of the term-time forced model are similar. Infected individuals recover at constant rate *ν* (1/*ν* is the average infectious period). As is common in the mathematical epidemiology literature [12,43], we restrict our attention to the case when birth and death rates *μ* (1/*μ* is the average lifetime) are equal, and thus the total population size *N* is constant. This allows us to reduce the number of independent variables to two and define the state of the system as *σ* = {*S*(*t*), *I*(*t*)}. From *β*_{0}, *ν* and *μ*, we can express one of the most important epidemiological parameters [12,43], the basic reproductive ratio *R*_{0} = *β*_{0}/(*ν* + *μ*). *R*_{0} is the average number of secondary cases caused by one infectious individual introduced into a fully susceptible population; *R*_{0} will be used throughout the text.

### 2.2. Theoretical analysis

Two main approaches can be used to investigate the dynamics of the stochastic model formulated above. An analytical approach starts from the formulation of the model as a master equation for the probability of finding the system in state *σ* with *S*(*t*) susceptibles and *I*(*t*) infectives at time *t* [46–49]. Much understanding about the stochastic dynamics relevant for recurrent epidemics can be gained if this equation is expanded in the powers of [46]. An extensive discussion of this approach has already been given at length in the literature in the context of epidemic models, and we refer the reader to [37,38,50] for formal details. Here, we describe only the aspects that are important for this paper. In essence, the method involves the substitutions and in the master equation that can then be expanded to obtain two systems of equations [46]. At the leading order, the expansion gives rise to a set of ordinary differential equations for the *mean* variables, i.e. the fractions (densities) of susceptible and infected individuals, and . These equations are the same as the standard deterministic SIR model with sinusoidal forcing [26]. At next-to-leading order, it gives rise to a set of stochastic differential equations for *fluctuations* of susceptible, *x*_{S}(*t*), and infected, *x*_{I}(*t*), individuals about the mean behaviour given by the deterministic model [46]. From these equations, we are able to analytically calculate power spectra of fluctuations for susceptibles, *P*_{S}(*f*), and infectives, *P*_{I}(*f*), as functions of frequency *f*. We are interested in the endemic behaviour of the model, so the spectra correspond to the fluctuations about *stable attractors* of the deterministic model which for *ε* = 0 and *ε* > 0 are the endemic *fixed point* [50] and stable *limit cycles* with a period that is an integer multiple of a year [37,38], respectively. Further technical details relating to analytical calculations are given in the electronic supplementary material. Throughout the text, we will use the words *theoretical* and *analytical* interchangeably to refer to spectra computed as explained in this section.

In our analysis, we will focus on a spectrum *P*_{I}(*f*), which, for the sake of simplicity, will be denoted as *P*(*f*), and its three characteristics, namely *dominant period*, *amplification* and *coherence* [50,51] (figure 2). We define the dominant period of stochastic fluctuations as the inverse frequency of the maximum of the highest stochastic peak. We also compute the total spectral power that equals the area under a power spectrum curve. This quantity defines the ability of the system to sustain oscillations of all frequencies and shall be referred to as the amplification of stochastic fluctuations. Finally, the coherence is defined as the ratio of spectral power lying within 10% from the dominant period and the total spectral power. It serves to measure how well-structured oscillations about the dominant period are.

As a rule, a theoretical spectrum of unforced epidemic models (*ε* = 0) has one peak [41,50], whereas it can have several peaks of different amplitudes for *ε* > 0 [37,38]. Away from bifurcation points of the deterministic model, one of them is usually much higher than the others. We are not aware of any work assessing the relevance of secondary peaks to recurrent epidemic behaviour seen in the real data. The highest peak, however, has been shown to be important in understanding the interepidemic periods observed in time series of pertussis and measles [38,41,42], and is therefore used in the definition of a spectrum's characteristics in this paper.

### 2.3. Simulations

We simulate the model using an extension of Gillespie's algorithm [35,36] which produces stochastic trajectories for {*S*(*t*), *I*(*t*)} in continuous time. These are processed further to compute numerical spectra and test them against the theoretical prediction for *P*(*f*). The simulation length is 500 years, and the first 50 years are discarded. In numerical work, a time series for fluctuations *x*_{I}(*t*) is obtained as where is the fraction of infectives averaged over many realizations of the model. From *x*_{I}(*t*), we compute a spectrum *P*(*f*) using the discrete Fourier transform. For *ε* > 0, we also present a spectrum of the entire ‘signal’ (scaled by population size *N*), *I*(*t*), which will be referred to as a *full* spectrum. By definition, *P*(*f*) includes only stochastic peaks, whereas the full spectrum includes both deterministic peaks corresponding to a limit cycle and stochastic peaks corresponding to fluctuations about it. For either of these spectra, we will use the words *simulated* and *numerical* interchangeably to emphasize that they were computed using the method described in this section. For each set of parameters, 250 simulations are recorded, and all final spectra are averaged over those where no extinctions occurred during 500 years. The initial conditions for susceptibles and infectives are chosen from *S*(0) = round(*s*_{c}*N*) + *U*(0,30) and *I*(0) = round(*i*_{c}*N*) + *U*(0,30), where round(·) is rounding to the nearest integer, *U* is the uniform distribution and *s*_{c} and *i*_{c} is the fixed point (for *ε* = 0) or a random point on the limit cycle (for *ε* > 0). The random number generators used in the Gillespie algorithm are initialized with unrepeated seeds which guarantees that the simulated stochastic trajectories are all different. We have also checked that with this choice of initial conditions all simulations converge to a stationary state within 50 years (transient period) or die out and so are not taken into consideration.

Both the theoretical analysis (described in §2.2) and the numerical analysis based on simulations (described in §2.3) are suitable for the investigation of temporal patterns in large populations, such as those corresponding to the time series in figure 1, and both have advantages and limitations. The simulations can be quite easily implemented but progressively become computationally intensive as the population size, *N*, increases. In the type of systematic study performed here, the numerical analysis would become very time-consuming for populations larger than one million individuals. However, it is exactly for such large populations that the approximate analytical spectrum, computed from the expansion in the inverse population size, is expected to predict the dynamics very well [38]. In the case when the seasonality is absent, the analytical spectrum is given by a simple formula (see the electronic supplementary material) which can be readily used to compute spectral characteristics. For the seasonally forced model, the spectrum can be written as an analytical formula too [37,38]; however, the calculation is more evolved and has to be done numerically, because no closed-form expression for a limit cycle can be found. The theoretical analysis also helps to understand the mechanisms behind the dynamics such as the change in temporal patterns when approaching bifurcation points of the deterministic model, that are not always clear from the spectral analysis of simulated time series. The limitations of the theory are that (i) it does not allow us to compute the spectrum of fluctuations when the deterministic model has several stable coexisting limit cycles for a given set of parameters, and (ii) in the vicinity of bifurcations the perfect agreement with simulations is achieved for populations larger than one million. Both cases will be discussed in detail in the next section.

The time-series data presented in figure 1 correspond to large populations (Hidalgo's population size was about 2.1 million [7], and the other time series are likely to correspond to even larger populations). In the Results section, we perform the systematic comparison of recurrent patterns for a large range of realistic values of the infectious period and the basic reproductive ratios. For the sake of computational speed, we focus on the population of one million of individuals and we demonstrate how the spectral analysis can be used to predict the dynamics in even larger populations.

## 3. Results

We compare the numerical and theoretical predictions for different spectra, and the three measures we have conventionally chosen to characterize them. In the beginning, we explore a large region of parameter space and later discuss the main findings for rubella and measles.

### 3.1. No seasonal forcing: *ε* = 0

#### 3.1.1. Theoretical and simulation results

We first restrict our attention to the case when there is no seasonality, for which an explicit expression for the analytical spectrum can be found (see the electronic supplementary material). The deterministic SIR model has only one endemic fixed point provided *R*_{0} > 1 [12,43]. Figure 3 shows analytical and numerical results for the dominant period of stochastic fluctuations about it as well as their amplification and coherence for a range of basic reproductive ratios, *R*_{0}, and infectious periods, 1/*ν*. Analytical spectra can be obtained for any *R*_{0} > 1. In practice, long numerical simulations may not be feasible for all parameter combinations where *R*_{0} > 1, because the system experiences frequent extinctions when the infectious period is short. The results are presented in figure 3, and the domain where this happens is shown in grey. In what follows, the line separating the grey region from the rest of the parameter space will be called the extinction boundary. We would like to stress that this concept is used for our convenience and is different from the concept of the CCS. In particular, its location depends on the population size, the length of time series and the number of runs used in simulations. For populations smaller/larger than one million individuals and simulation length larger/smaller than 500 years, the grey region would be extended/abridged, and the extinction boundary would be shifted. In the region of parameter space amenable to the exploration of the long-term dynamics of the model, we observe that the structure of the spectra uncovered by the theory is clearly visible in the simulations too. Small systematic deviations between the two predictions are expected and occur close to the extinction boundary. These are due to the non-Gaussian character of the fluctuations which cannot be explained within the theory used in this study (see the electronic supplementary material). To achieve the perfect agreement with the simulations within this region, the theoretical analysis would require considering the next order corrections of the order to the macroscopic equations. The deviations are mainly reflected in the broadening of a spectrum and appearance of secondary peaks. As a consequence, amplification (figure 3*e*) is slightly increased (coherence, figure 3*f*, is correspondingly decreased) in the simulations in the area adjacent to the grey region. In addition to these changes, an increase of the dominant period of stochastic fluctuations in the simulated spectra may be observed [37,51]. This effect will be discussed in more detail in the next section where seasonality is included as, without seasonality, it is only barely apparent (compare figure 3*a* and *d*).

#### 3.1.2. Implications for rubella and measles dynamics

The structure discovered in figure 3 allows us to derive an initial picture of the dynamics of rubella relative to measles. Based on the estimates of parameters typical of these diseases for the pre-vaccination era [7,10,12,18,26,32], we superimpose their approximate locations in all the panels of figure 3. The rubella estimates are for Mexico and Canada, and the measles estimates are for England and Wales. For rubella, the infectious period, 1/*ν*, is about 18 days and *R*_{0} ranges from 3.4 to 9.5 in Mexico (figure 1*a*) [7]. *R*_{0} in the Canadian province Ontario ranges from 4.6 to 6.5 where the lower bound is the estimate for the years shown in figure 1*c* [10]. For measles, we have taken the most frequently used values for large cities (e.g. London) in England and Wales before vaccination, 1/*ν* about 2 weeks and *R*_{0} around 18.

The results so far ignore the seasonality of transmission rate and so are insufficient to explain the patterns of measles in which it plays a pivotal role [18,26,52]. However, they have important implications for understanding the dynamics of rubella. As we shall confirm shortly, for *ε* > 0, the spectrum of stochastic fluctuations for this disease is close in form to that obtained for the unforced model that correctly predicts a dominant period associated with the stochastic peak of about 5–6 years (figure 3*a*,*d*) as seen from the comparison with the left peaks in the data spectra (figure 1). This period is similar to the natural period of small amplitude perturbations from the endemic fixed point recovered in the purely deterministic setting [32].

Our analysis of the stochastic model allows us to quantify other features of fluctuations using amplification and coherence (figure 3*b,c*,*e*,*f*). For rubella, the amplification is large indicating that the epidemic patterns of the unforced model represent high amplitude oscillations. High coherence suggests that only a few of the frequencies involved in the stochastic fluctuations account for most of the variance of time series. This peculiar type of dynamics sets rubella close to the extinction boundary. Large coherent multi-annual epidemics with troughs deeper than in the region with higher *R*_{0} cause regular extinctions. In §3.2, we discuss how these descriptions of the stochastic dynamics of rubella are changed in the presence of seasonality and compare it with the dynamics of measles.

### 3.2. Seasonal forcing: *ε* > 0

#### 3.2.1. Theoretical and simulation results

For *ε* > 0, the spectra are associated with stochastic fluctuations about stable attractors of the deterministic model, i.e. stable limit cycles of a period in multiples of a year [37,38]. The seasonally forced deterministic SIR model has a complex bifurcation diagram with regimes where multiple limit cycles may coexist [44,45]. For high birth rates and high seasonality, regions corresponding to chaotic dynamics are found [21]. Across most of the range of parameter space we explore here, either annual or biennial limit cycles are present. We performed the theoretical analysis for these attractors for different parameters and found that the agreement between the theory and simulations is, in general, excellent. Nevertheless, small discrepancies are again expected if a limit cycle is not sufficiently stable and/or a population is small. In particular, this happens near the extinction boundary, and is therefore relevant for rubella.

For parameters reflecting rubella, an annual limit cycle is found in the deterministic model. To illustrate the effects of population size on simulated spectra, we show in figure 4 an analytical spectrum about this attractor (dashed red line) and full numerical spectra for *N* = 10^{7} (solid green line) and *N* = 10^{6} (solid black line). For the larger population size (figure 4*a*), the simulated spectrum exhibits two types of peaks (solid green line). There is a dominant annual peak (corresponding to the deterministic annual cycle) and a subdominant broad multi-annual peak (corresponding to stochastic fluctuations about it). The latter is indistinguishable from the theoretical spectrum (dashed red and solid green lines coincide, deviating only for values corresponding to frequencies of around 1 year). The derivation of an approximate theoretical spectrum (dashed red curve) from the expansion in powers of suggests that the stochastic fluctuations, *x*_{I}(*t*), are Gaussian (), see the electronic supplementary material. This allows us to represent the full spectrum of as the sum of two parts: a deterministic part that scales as *N*^{2} and a stochastic part that scales as . The full spectra in our analysis are normalized (divided) by *N* as we mention in §2.3. Therefore, for large populations where there is a perfect agreement between the analytical and simulated spectra as in figure 4*a*, the amplitude of the deterministic peak is proportional to *N*, and the amplitude of the stochastic peak is independent of *N*.

This scaling is captured even for smaller populations where the fluctuations become non-Gaussian. For *N* = 10^{6}, the peak at 1 year becomes subdominant (see the solid black line in figure 4*b*). This indicates that the contribution of an annual component in the time series decreases with decreasing *N*. As for fluctuations beyond the annual component, at least two stochastic peaks at 5.8 and 2.9 years can be clearly seen (solid black line). Although the theory does not capture them in full, the agreement is still good and, more importantly, the systematic qualitative changes can be predicted. For small populations, the dominant period of fluctuations in simulations is slightly increased and their variance is distributed over a larger range of frequencies with respect to theoretical predictions [37]. This is compatible with a general observation of the increased stochasticity and therefore much more irregular dynamics in small populations [24,53]. We would like to point out that the discussed discrepancies are not attributed to sample size effects (the number of runs used to compute the spectrum). The latter may lead to discrepancies only for parameters at the very border with the black line (see figure S1 in the electronic supplementary material). The example, we presented here, was for *R*_{0} = 4, which is two points away from the extinction boundary. Simulations for larger *R*_{0} show smaller deviations from analytical calculations even for populations as small as *N* = 10^{6} (see figure S2 in the electronic supplementary material). As mentioned before, we expect these results to be robust to the form of seasonal forcing. This is confirmed in figure S3 of the electronic supplementary material which shows that the simulated spectra from figure 4 are reproduced by the term-time forced model at a 2.7 larger forcing amplitude.

Another situation which the analytical theory cannot fully account for is stochastic switching between different attractors of the deterministic model [20,32]. The computation of an analytical spectrum about a limit cycle requires the knowledge of its geometric orbit (see electronic supplementary material). In the sinusoidally forced SIR model, several stable attractors coexist in the regions of small *R*_{0} and 1/*ν* [44], and spectra of stochastic fluctuations about each of them can be obtained separately [37]. The theoretical analysis, however, does not allow us to predict which of the attractors will be observed in simulations and what their relative contribution to the total stochastic dynamics is. Previous analysis of the stochastic dynamics of measles and pertussis showed that the only attractors seen in simulations of the seasonally forced SIR model (and other related models of infectious diseases) are annual and biennial cycles [37,38,41]. The stochastic switching was observed to happen exclusively between these attractors and only for measles. This result is, however, of limited value to us because it is restricted to particular parameter choices, and so we cannot assume that the switching does not happen in the broader span of the parameter space.

Spectra from simulations contain complete information about the frequency distribution of oscillations and are thus helpful to identify switching between attractors through the presence of unexpected peaks. Figure 5 shows simulation results for the dominant period of stochastic fluctuations, amplification and coherence for two values of seasonality *ε* and *N* = 10^{6}. In addition to these quantities, we compute the dominant period in the full spectrum that includes both stochastic and deterministic peaks.

To examine the effect of seasonality on stochastic fluctuations, figure 5*a–f* can be directly compared with figure 3*a*–*f* for the unforced model. As *ε* increases, the (grey) domain with frequent extinctions is extended and approaches the measles parameters. For most values of *R*_{0} and 1/*ν*, we have explored (the coloured region) the dynamics of the stochastic model are associated with fluctuations about only few attractors. First, the biennial cycle is found inside the region of increased coherence and amplification in figure 5*b*,*c,e*,*f* which includes measles and is absent in figure 3*b*,*c,e*,*f*. The dominant period in the full spectrum of time series demonstrating such a behaviour is at 2 years (figure 5*h*). Second, stochastic switching between annual and triennial cycles is detected in a small region relevant for measles with low values of *R*_{0}, see an unexpected increase of amplification in figure 5*e* around *R*_{0} = 10. The spectra here have a dominant annual peak (figure 5*h*) and a subdominant triennial peak (figure 5*d*). The amplification of oscillations associated with the latter is, however, much higher than what we would expect to see for fluctuations around an annual cycle (figure 5*b*). In figure 6, we show the full spectrum and a typical time series corresponding to the switching between the annual and triennial cycles. Third, in the rest of the (coloured) region that includes rubella, the spectra are similar to those of the unforced model. The dynamics of the stochastic model here corresponds to fluctuations about an annual cycle and we discuss it first.

#### 3.2.2. Implications for rubella and measles dynamics

From figure 5*a,d*, we see that for relatively small basic reproductive ratios, typical of rubella, the seasonality does not affect the dominant period of stochastic fluctuations which continues to be centred at about 5–6 years. The amplification (coherence) is only slightly increased (decreased) as *ε* increases (figure 5*b*,*c,e*,*f*). The full spectra of rubella resemble that of figure 3 with a sharp peak at 1 year and a broad multi-annual peak.

For future discussion of the implications of vaccination and decline or increase in birth rates, it is useful to investigate how the spectra of rubella change with *R*_{0}. Keeping the amplitude of seasonal forcing and infectious period of rubella fixed and increasing *R*_{0}, we expect the period of stochastic fluctuations as well as their amplification to decrease. This is seen from figure 5 and also illustrated in figure 7*a* where the full spectra for parameters close to rubella estimates are shown. The relative contribution of multi-annual and annual frequency components in model time series can be read from the same figure. For small *R*_{0}, the fluctuations are large, and the multi-annual peak is dominant but for larger *R*_{0,} it becomes subdominant, and the annual peak is enhanced. The increase of *ε* (as well as the increase of population size as discussed before in the text accompanying figure 4) results in the enhancement of the deterministic peak too (compare figure 5*g* and *h*), but does not change the dominant period of fluctuations significantly.

For measles, there are major changes in the behaviour as both coherence and amplification increase drastically for *ε* > 0, see the newly appeared oval-shaped regions near measles parameters in figure 5*b*,*c*,*e*,*f*. To demonstrate that this phenomenon indicates the appearance of a biennial cycle in simulations we show in figure 7*b* the full spectra for fixed infectious period, 1/*ν* = 16 days, and a range of *R*_{0} ∈ [22,29]. These values of 1/*ν* and *R*_{0} are slightly higher than the commonly used estimates for measles in England and Wales, for example, 1/*ν* = 13 days and *R*_{0} = 18 for London before vaccination. The same qualitative dynamics is observed for lower values if the simulation length is shorter than 500 years (results not shown). For *R*_{0} = 22, the spectra are typical of fluctuations about an annual limit cycle with a deterministic peak at 1 year and a broad stochastic peak near 2 years. If *R*_{0} is increased further, then the fluctuations around an annual cycle become macroscopic (the stochastic peak at 2 years becomes much higher) smoothly turning into a biennial limit cycle. This transition corresponds to a period doubling bifurcation in the deterministic model. A strong biennial behaviour with a dominant peak at 2 years and a secondary harmonic at 1 year is observed for example, for *R*_{0} = 26. Finally, for even larger *R*_{0}, we see a transition from biennial to an annual cycle again. The set of transitions seen in figure 7*b* is typical of measles and have been observed in related models of infections dynamics via analytical and numerical studies in other research [20,37,38,45,54,55].

The seasonality would act to change the picture in figure 7*b* in the following way. From the comparison of figure 5*g,h*, the region of parameter space where such a behaviour is seen is expected to get larger with increasing *ε*; in particular, for *ε* > 0.1, the period doubling transition is induced for values of *R*_{0} much smaller than in figure 7*b*.

Previous analysis of measles data from England and Wales and the USA has shown that transitions in the dynamics due to an increase or decline of birth rates as well as the introduction of vaccination are associated with transition between annual and biennial limit cycles [20]. Using a simple mapping from changes in vaccination and birth rates to effective changes in *R*_{0} introduced in reference [20], our results are consistent with this view. For large communities with very high birth rates (high *R*_{0}) such as Liverpool before vaccination, US cities in the period after the Great Depression or cities in developing countries, we would expect to be in the regime with an annual cycle [19,20,23,56]. Other large cities with smaller birth rates such as London are in the regime with a biennial cycle [26]. The corresponding spectrum with narrow and sharp peaks at 2 and 1 years has been the main reason of more regular and thus more predictable patterns of measles epidemics in large cities. The vaccination introduced in UK in 1968 lowered *R*_{0} and induced a transition from the biennial to the annual cycle with large stochastic fluctuations. Our analysis thus offers an insight into the factors responsible for the shift from regular epidemics of measles before vaccination to less irregular in the vaccine era [16].

The last finding deserving a further comment concerns the switching between an annual and triennial cycles found for moderate values of *R*_{0} (figures 5*e* and 6). This behaviour may be responsible for the triennial cycles observed in Baltimore and other US cities during the Great Depression [29] but more thorough analysis is required to confirm this.

## 4. Discussion

In this paper, we have investigated the behaviour of the stochastic seasonally forced SIR model based on the spectra of long time series for a large range of basic reproductive ratios and infectious periods. For relatively low values of *R*_{0} relevant for rubella, the model predicts spectra with a stochastic multi-annual peak at about 5–6 years and a deterministic annual peak. Both peaks are observed in the spectra of rubella data (figure 1). The multi-annual peak stays largely unchanged under the introduction of seasonality (figures 3 and 5*a*,*d*) or population size (figure 4) which explains its presence in time series from different locations.

Using the complementary measures, coherence and amplification, we further studied how the structure of stochastic fluctuations in the model changes with *R*_{0}. By definition of these measures (§2), it is not possible to estimate them from a single (data) time series such as shown in figure 1. We therefore cannot compare our model's predictions for these measures to the data directly. However, it is still possible to match the data and model's full spectra (i.e. not only the positions of the stochastic and deterministic peaks but also their heights) given the information about the rate of reporting and the population size. We did this for the Mexico time series (see figure S4 in the electronic supplementary material) and obtained a very good agreement between the model and the data. For the Canada and Japan time series we lack the aforementioned information, so these analyses yield the full spectra with correct positions of the peaks (as discussed and shown in §3.2.2) but their heights can vary depending on the reporting rate and population size used in the model.

A visual inspection of simulated time series demonstrates intriguing behaviour emerging from the interaction between stochasticity and a deterministic annual cycle. Figure 8 shows typical time series for the unforced (red dashed line) and seasonally forced (blue solid line) cases. If *ε* = 0, then the epidemic patterns represent multi-annual coherent oscillations. As *ε* is increased, we find qualitatively different dynamics all of which correspond to spectra with a multi-annual and an annual peak. Figure 8*a* shows an example of annual epidemics of alternating amplitudes modulated by an oscillation of a long period corresponding to the period of stochastic multi-annual fluctuations. This dynamics is qualitatively similar to the multi-annual regularity observed, for example, for rubella in Japan (figure 1). We also find very large outbreaks followed by outbreaks of much lower amplitude as in figure 8*b*. Such a behaviour may be responsible for the spiky dynamics observed in, for example, Canada (figure 1). Note that the spikes in the data could also arise from spatial effects such as local extinction of the disease followed by reintroduction from another region. However, as we do not possess more resolved data, it is impossible to reach a final conclusion with regards to this issue.

The patterns of rubella incidence in large populations are in contrast with those of measles. For the latter, the spectra are characterized by sharp and narrow peaks at 1 and 2 years (as opposed to the broad multi-annual peak and a narrow peak at 1 year observed for rubella) and thus correspond to much more regular dynamics. The transitions in measles behaviour owing to vaccination or change in birth rates are associated with transitions between the annual and biennial limit cycles of the deterministic model. In future work, it would be interesting to study stochastic measles dynamics for higher levels of seasonal forcing that correspond to chaos in the deterministic model [21].

Both measles and rubella are found to be close to the extinction boundary, and increasing the amplitude of seasonal forcing only extends the region of parameter space with frequent extinctions. Figure 9 illustrates typical stochastic trajectories from simulation in the susceptible–infected plane from which the spectra were computed. Interestingly, these patterns suggest that the mechanism accounting for high extinction probability is different for rubella and measles. For rubella, extinction occurs as a consequence of large stochastic fluctuations about a small (and globally less stable) annual limit cycle (figure 9*a*). For measles, extinctions are mainly due to the shape of a large (and globally more stable) biennial limit cycle, from which the system can be driven to extinction by even relatively small fluctuations (figure 9*b*).

We can use the framework developed to predict the effect on persistence of an effective reduction in *R*_{0} by vaccination or declining birth rates for rubella and measles. For measles in the biennial regime, either an increase or a decrease of *R*_{0} can lead to fluctuations around an annual cycle (rather than around a biennial cycle) that could result in lower extinction rates and thus higher persistence. For rubella, a reduction in *R*_{0} will lead to larger and more coherent oscillations that would unambiguously result in higher extinction probabilities, and thus lower persistence. Both these outcomes merit serious consideration in a public health context: vaccination against measles can make local elimination less likely [21], whereas vaccination against rubella is likely to increase local extinction, allowing the build-up of susceptible individuals in later age classes [5,7], potentially leading to an increase in the burden of CRS.

These conclusions point to a need for theoretical developments towards uncovering the mechanisms of stochastic extinctions in small populations based on the analysis of epidemic models (a thorough overview of studies in this area with the aim of understanding the persistence of measles is given in the recent work by Conlan *et al*. [57]). In the mathematical framework, we adopted in this paper, the approach to computation of the distribution of extinction times in an unforced stochastic epidemic model was proposed some time ago [58,59]. The disease persistence in stochastic epidemic models can also be studied using the so-called Wentzel–Kramers–Brillouin (WKB) approximation, but the method is only applicable to low-dimensional and unforced models [60,61]. Nevertheless, no analytical progress can be made along the same lines for the seasonally forced model we use here. We are aware of only one study [62] that addressed extinction probabilities in the periodic context using theoretical methods, but the method of Bacar & Ait Dads [62] has a limitation because it applies to the large population limit only. The development of the approaches to compute the time to extinction in seasonally forced models will be therefore a subject of further research.

Our focus has been on measles and rubella; however, the broad span of parameter space explored means that our results may shed light on the dynamics of other diseases whose dynamics can be described by a simple SIR formalism with seasonal forcing, for example, pertussis [31,41]. Although infection with pertussis does not confer permanent immunity, the SIR model has been shown to capture the qualitative patterns to some extent [41]. Taking the pertussis parameters before vaccination that are well established in independent data sources ([10,12]; 1/*ν* = 22 days and *R*_{0} = 17 for London and Ontario, Canada) from figure 5, we see that for pertussis the dominant period of stochastic fluctuations is 2–3 years. These periods are in agreement with the documented interepidemic periods [10,41]. We also find that the coherence of fluctuations is very low, which is compatible with the famously irregular dynamics of pertussis. The decrease of *R*_{0} due to vaccination would act to increase the dominant period and coherence which also agrees with the observation of the shift to more regular dynamics in the vaccine era [16].

To conclude, in this paper, we used a stochastic framework to explain the recurrent dynamics of rubella, particularly in comparison with measles. Our analysis revealed that while both rubella and measles are relatively close to their extinction boundary, the reasons for this are very different. Finally, our analysis showed that, for rubella, reducing *R*_{0} by vaccinating or a declining birth rate unambiguously result in higher extinction probabilities, whereas, for measles, outcomes can be more complicated; and both these facts have public health implications.

## Funding statement

G.R. gratefully acknowledges the financial support from the Portuguese Foundation for Science and Technology under contract no. SFRH/BPD/69137/2010. B.T.G. was supported by the RAPIDD programme of the Science and Technology Directorate, Department of Homeland Security and the Fogarty International Center, National Institutes of Health. J.E.M. and B.T.G. were supported by the Bill and Melinda Gates Foundation.

## Acknowledgements

We thank Kihei Terada for helpful discussion of rubella dynamics in Japan; Alan McKane for his careful reading and comments on the paper; and Ottar Bjørnstad for the help with the computation of the confidence intervals for data spectra.

- Received July 17, 2013.
- Accepted August 20, 2013.

© 2013 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.