PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of interfacehomepageaboutsubmitalertseditorial board
 
J R Soc Interface. Nov 7, 2012; 9(76): 2959–2970.
Published online Jun 20, 2012. doi:  10.1098/rsif.2012.0432
PMCID: PMC3479927
Modelling the long-term dynamics of pre-vaccination pertussis
Ganna Rozhnova1,2* and Ana Nunes1
1Departamento de Física, Centro de Física da Matéria Condensada, Faculdade de Ciências da Universidade de Lisboa, 1649-003 Lisboa Codex, Portugal
2Kavli Institute for Theoretical Physics, University of California, Kohn Hall, Santa Barbara, CA 93106-4030, USA
*Author for correspondence (a_rozhnova/at/cii.fc.ul.pt).
Received May 27, 2012; Accepted May 28, 2012.
The dynamics of strongly immunizing childhood infections is still not well understood. Although reports of successful modelling of several data records can be found in the previous literature, the key determinants of the observed temporal patterns have not yet been clearly identified. In particular, different models of immunity waning and degree of protection applied to disease- and vaccine-induced immunity have been debated in the previous literature on pertussis. Here, we study the effect of disease-acquired immunity on the long-term patterns of pertussis prevalence. We compare five minimal models, all of which are stochastic, seasonally forced, well-mixed models of infection, based on susceptible–infective–recovered dynamics in a closed population. These models reflect different assumptions about the immune response of naive hosts, namely total permanent immunity, immunity waning, immunity waning together with immunity boosting, reinfection of recovered and repeat infection after partial immunity waning. The power spectra of the output prevalence time series characterize the long-term dynamics of the models. For epidemiological parameters consistent with published data, the power spectra show quantitative and even qualitative differences, which can be used to test their assumptions by comparison with ensembles of several-decades-long pre-vaccination data records. We illustrate this strategy on two publicly available historical datasets.
Keywords: pertussis, childhood infections, recurrent epidemics, stochastic fluctuations, power spectra
Childhood infections remain a public health concern as well as a modelling challenge, despite many decades of control measures and theoretical efforts. Large-scale vaccination programmes against some of these diseases started in the 1940s–1960s in developed countries and led to marked reductions of incidence levels, but in developing countries they are still a cause of significant levels of infant mortality [15]. The resurgence of pertussis reported in developed countries after decades of high vaccine coverage [610], as well as a recent upsurge of measles in eastern and southern Africa, prompted a renewed interest in assessing the efficacy of control policies for these childhood infections.
Such an assessment must rely on simulations based on mathematical models [1114], whose more recent versions highlight the importance of contact structure, stochasticity and seasonality as drivers of the observed temporal patterns of incidence [1522]. Concurrently, analytical tools to deal with these ingredients have been developed [2330].
Models of higher complexity involve many parameters, some of which are difficult to determine, leaving considerable room for data fitting [3134]. On the other hand, one of the features of childhood infections is the variability exhibited by different data records, even within those that correspond to a single disease in comparable social environments [4,17,35]. Therefore, successful modelling of particular datasets reported in the literature [17,36,37] has not closed the discussion about the key determinants of the observed dynamics. Measles is an exception: a parsimonious, stochastic, susceptible–infective–recovered (SIR)-based model has been shown to adequately describe disease incidence in large urban populations [15,38,39].
At the opposite extreme, pertussis keeps defying mathematical modelling, as illustrated by recent efforts in different directions [10,37,4043]. In particular, several proposals explore the hypothesis about disease-induced and vaccine-induced immunity [10,35,43], relying on assumptions about vaccine uptake and efficacy for the purpose of comparison with real data. However, we are still lacking a sound uncontroverted model for pre-vaccination pertussis.
Here, we have sought to contribute to the goal of using pre-vaccination data records to obtain information about the properties of disease-induced immunity. The paper focuses on exploring the influence of naive hosts' immune response in the long-term patterns of pertussis prevalence as given by the averaged power spectra of simulated time series corresponding to several decades. The power spectrum of the stochastic, seasonally forced, well-mixed SIR model is compared with those of four modifications of the model that reflect different assumptions about disease-induced immunity. These four different assumptions have been proposed in the literature in the framework of deterministic models that have all been shown to be compatible with available data. The five variants we consider deliberately avoid all the complications related to contact structure and spatial spread, as well as vaccine coverage and efficacy and waning of vaccine-induced immunity, in order to keep a relatively low number of free parameters in the model.
We show that the stochastic versions of the SIR model and its four variants have significantly different properties, which translate into quantitatively different prevalence and incidence power spectra. This opens the possibility of using the stochastic properties of long, well-resolved data records to constrain these and other variants of the model, with the power spectra of the pertussis incidence time series in large urban centres in the pre-vaccination era as the target long-term dynamics that the model should reproduce. We illustrate this strategy by applying it to two publicly available historical data records for pre-vaccination pertussis incidence.
2.1. Models
We consider five seasonally forced, stochastic compartmental models, summarized in figure 1, as candidates for the description of pre-vaccination pertussis. Their deterministic counterparts are described in the electronic supplementary material.
Figure 1.
Figure 1.
Five compartmental models for pertussis. Models (a) and (b) are the seasonally forced SIR and SIRS models, respectively, while models (ce) are extensions of the SIR(S) paradigm. Model (c) allows for boosting of immunity proportional to the force (more ...)
In all cases, the population includes three classes of individuals, the susceptible (S in (ad), S1 and S2 in (e)), the infectious (I in (ac), I1 and I2 in (d,e)) and the recovered (R in (a,b) and (d,e), R and W in (c)). Also, in all five cases, there is replenishment through births of naive susceptibles. The birth rate is kept constant and equals the death rate μ.
2.1.1. SIR model
Model (a) is the standard SIR model: recovery confers permanent immunity, and all infections are first infections. Infected individuals recover at a constant rate δ, and susceptible individuals are infected at a rate An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i1.jpg, where N is the population size, I is the number of infected individuals and β(t) is a periodic function of period 1 year that represents the variable contact rate associated with the school terms. To study the long-term temporal patterns of the fluctuations of this stochastic process, the particular form of β(t) is not crucial, and we shall take An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i2.jpg.
2.1.2. SIRS model
Model (b) is the standard SIRS model. It differs from the SIR model in that recovery does not confer permanent immunity. Instead, immunity wanes at a constant rate γ. The fraction of primary infections is given approximately by An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i3.jpg, where s* and i* are the equilibrium values for the susceptibles and infectives in the deterministic counterpart of the model (see the electronic supplementary material for details). Although there is no direct correspondence between primary infections and age groups in the model, it is reasonable to admit that a fraction of the infectives close to the fraction of non-primary infections, given by An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i4.jpg, are not school children, and that β1 should be reduced proportionally. This correction could be refined, but, as we shall see in §3, the properties of the stochastic fluctuations in this model are largely insensitive to the value of β1.
2.1.3. SIRWS model
Model (c) is a SIRS-type model allowing for immune boosting of recovered individuals upon re-exposure to infection. Infectives I recover at rate δ and susceptibles S get infected at rate λt, as in the SIR model. However, in addition to the S and I classes, there are two classes of recovered individuals denoted by R and W. The immunity of the former is not permanent and they move to the waning class W at a constant rate 2γ. The individuals in class W undergo two possible transitions: further immunity loss at rate 2γ to become susceptible S, or immunity boosting upon contact with infectious individuals to return to the recovered class R at a rate An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i5.jpg, where κ is the immunity boosting coefficient. We call this scheme the immune boosting model and denote it by SIRWS.
An age-structured version of this model with vaccination has been used to explain the recent re-emergence of pertussis cases, despite high vaccine coverage in Massachusetts, USA [10].
2.1.4. SIRIS model
Model (d), proposed by Aguas et al. [43], sets a scenario based on the SIRS model with a moderate rate γ of immunity loss in which recovered individuals are immune to severe disease but susceptible with reduced susceptibility to mild forms of the disease. The classes of individuals infected with severe and mild infections are denoted by I1 and I2, respectively. Recovery from both forms of the disease takes place at the same rate δ. Infectiousness of mild infections is reduced by a factor An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i6.jpg, and susceptibility of recovered individuals to mild infections is also reduced by a factor An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i7.jpg. Moreover, because mild infections typically occur in adults who are not affected by seasonal forcing, the force of infection is therefore taken to be An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i8.jpg for susceptible individuals S and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i9.jpg for recovered individuals R. We call this scheme the reinfection model and denote it by SIRIS.
The main feature of this model is that it exhibits a reinfection threshold [44], a value of infectiousness above which total disease incidence rises by one or more orders of magnitude, owing to high incidence levels of mild infections. These mild infections are sub-clinical and so the disease incidence/prevalence is associated with the class I1.
2.1.5. SIRSI model
Finally, model (e), proposed in the study of Wearing & Rohani [40], assumes an immune response that is a combination of (b) and (d), in the sense that recovered individuals are fully immune to the disease, but they lose this immunity at a certain rate γ to become susceptible to repeat infections, although with reduced susceptibility. The two classes of susceptible individuals are denoted by S1 for the naive susceptibles and S2 for the susceptibles generated by immunity waning, and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i10.jpg is the reduced susceptibility factor for repeat infections. The classes of individuals infected with first and repeat infections are denoted by I1 and I2, respectively. Recovery from both classes takes place at the same rate δ. The class of repeatedly infected individuals contributes with reduced infectiousness to the pool of infectives responsible for disease transmission, and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i11.jpg is the reduced infectiousness factor. Repeat infections play a similar role in this model to that of mild infections in the SIRIS model, and the force of infection is taken to be An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i12.jpg for S1 and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i13.jpg for S2, because repeat infections typically occur in adults whose contacts are not subject to school-term forcing. Similar to the SIRIS model, the disease incidence/prevalence is identified with the class I1. We call this scheme the repeat infection model and denote it by SIRSI.
We finish this description by pointing out that, for a fixed β0, the unforced (β1 = 0) deterministic counterparts of the five stochastic models introduced in this section have different densities of infectives in equilibrium. In figure 2, we plot these equilibrium densities as a function of β0. The range of An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i14.jpg in the plot corresponds to the basic reproductive ratio [45] An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i15.jpg. Throughout the main text, β0 is fixed so that An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i16.jpg. This and values of the remaining parameters will be justified later. For An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i17.jpg, the density of infectives in the deterministic versions of the five models oscillates with period 1 year around the equilibrium value of the unforced equations (see the electronic supplementary material). No other attractors show up for the parameter values that we will use in our analysis.
Figure 2.
Figure 2.
Equilibrium infective density of the unforced deterministic versions of the models as a function of β0 for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i18.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i19.jpg. The lines stand for the I/N in the SIR (solid black), SIRS (dashed black) and SIRWS (dot-dashed black) models, and for the I1/N in (more ...)
2.2. Data records
Historical data records for pre-vaccination pertussis incidence in Greater London, UK, in the 1946–1957 period (figure 3a) and in Ontario, Canada, in the 1914–1943 period (figure 4a) are available from the International Infectious Disease Data Archive (http://iidda.mcmaster.ca). The time series correspond to weekly data in the case of London, and monthly data in the case of Ontario. Both raw and detrended (but not rescaled) data are shown for Ontario, where significant changes in population size and reporting efficiency are apparent during most of that period. In order to analyse the populations involved in these two data records, we shall consider only the undetrended data for Ontario corresponding to the last 8 years of the period. Demographic dataset the population of Greater London close to 8 million and that of Ontario close to 3.5 million in that period. However, the average recorded rate of new cases in a month interval is similar for the two datasets, indicating a ratio of effective populations close to 1. Differences in surveillance coverage and/or in reporting efficiency could explain this. For London, estimates of pertussis reporting efficiency in the 5–25% range can be found in the literature [46]. The fact that the effective population for Ontario is approximately of the same size indicates that a higher reporting efficiency compensates for the smaller population.
Figure 3.
Figure 3.
(a) New weekly cases of pertussis in (Greater) London before vaccination. The dashed horizontal line is the average number of weekly cases. (b) Spectrum of the time series. (Online version in colour.)
Figure 4.
Figure 4.
(a,b) New monthly cases of pertussis in the Canadian province of Ontario before vaccination (a) and the detrended time series (b). The dashed horizontal line in (b) is the average number of monthly cases. (c) The solid black line is the spectrum of the (more ...)
It is known from several case notification datasets that the interepidemic periods for pertussis in the pre-vaccination era show significant multi-annual structure. The multi-annual periods range from 2 to 3 years, and annual outbreaks are also observed [17,40,47]. The power spectra computed from the data confirm these observations (figures 3b and and44c). In the spectra plots, the shaded region marks the range of frequencies, An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i20.jpg, corresponding to the interepidemic periods An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i21.jpg years. Because of limited length and resolution in time, the spectra have a short range and a low resolution in frequency that affects mostly the annual peak. The widths of the multi-annual peaks are comparable. For Ontario, the power spectra of the full time series and of the detrended and undetrended partial datasets show the same dominant frequency components. In what follows, we shall consider only the power spectrum of the undetrended data.
2.3. Parameter values
The values of the demographic and epidemiological parameters that are well established in independent data sources [17,45] are kept fixed. These are the birth/death rate μ (or the average lifespan 1/μ), the rate of recovery from infection δ (or average infectious period 1/δ) and the average contact rate β0. The value of the last corresponds to the basic reproductive ratio An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i22.jpg reported in the literature for the two datasets considered in this study [17]. If the reported value of R0 comes from average age at infection data, together with the assumption of SIR dynamics, then it has to be corrected for models with temporary immunity. In the electronic supplementary material, we explore the whole range of possible R0 values for the SIRS and SIRSI models and show that the analysis of §3 below is not affected by this correction. We also give for the SIRWS and SIRIS models the dependence of R0 computed from average age at infection data on the model parameters.
For the remaining parameters, we either use the estimates found in the literature for a particular model or explore a range of possible values if such estimates are absent. For the SIRS model, we take 1/γ to be equal to 20–40 years. This is above the accepted range of duration of naturally acquired immunity for pertussis [48], but the large upper limit allows for a comparison with the SIR model. On the other hand, as we shall see, taking 1/γ smaller than 20 years would worsen the model's predictions. For the SIRIS model, we take 1/γ = 20 years [43] and the relative infectiousness and relative susceptibility of mild infections An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i23.jpg. This analysis can be extrapolated to other values of 1/γ, using the fact that this model has the SIRS model as one of its limiting cases. For the SIRWS model, we take the value of the boosting coefficient An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i24.jpg and 1/γ = 10 years considered in Lavine et al. [10] and, in addition, we study the behaviour of the model for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i25.jpg years and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i26.jpg. For the SIRSI model, we set 1/γ = 20 years and the relative infectiousness of and relative susceptibility to repeat infections An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i27.jpg [40]. We discuss in the electronic supplementary material the behaviour of this model for other values of the duration of immunity.
For illustration purposes, we use only three values of the amplitude of seasonal forcing β1 = 0, 0.05, 0.1 and discuss the predictions of the models for intermediate values of β1 whenever necessary. A summary of the parameter values is given in table 1.
Table 1.
Table 1.
The epidemiological description of the parameters and the range of their values for pertussis.
2.4. Stochastic simulations
To study the behaviour of the five candidates for the description of pre-vaccination pertussis, we use stochastic simulations of the processes described in figure 1. The simulations are based on a modification of Gillespie's algorithm [49], which accounts for the explicit time dependence in the contact rate [50,51]. In §§3.1–3.5, each simulation run starts from a random initial condition and the prevalence of the disease is recorded with a time step of 0.05 years for 450 years after 50 years of transient period. From each simulation run, the power spectrum of the prevalence time series is computed with the use of the discrete Fourier transform. The final spectra are averages of 103 simulations. The population size used in the simulations in §§3.1–3.5 is 5 × 105. Changing population size and/or computing the incidence in a given interval instead of the prevalence of the disease may change somewhat the amplitude of the multi-annual and annual peaks but does not influence their position.
We investigate the dynamics of the stochastic, seasonally forced, SIR, SIRS, SIRWS, SIRIS and SIRSI models by comparing the power spectra of long time series of the five models in the relevant parameter ranges. We first describe the performance of the SIR model and then assess whether each of its four extensions improves or worsens the results. In all cases, the spectra correspond to the fluctuations around the equilibrium of the unforced deterministic versions of each model (β1 = 0) or around the associated period 1 year stable attractor (An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i31.jpg). Numerical and simulational results for the parameter ranges we have explored never showed evidence of other attractors. A direct quantitative comparison of these power spectra with the spectra shown in figures 3 and and44 based on the amplitude of the multi-annual peaks and their precise location is undermined by the low resolution in frequency and poor statistics of the two datasets, which correspond to much shorter periods than those that can be used in the numerical simulations. That is why our first criterion in making the comparison between the models' predictions and the spectra computed from the data will be the position of the multi-annual peak. In order to be compatible with data records for pre-vaccination pertussis, a model's spectrum should have the multi-annual peak located in the range of frequencies An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i32.jpg (figures 3 and and4).4). This region is the shaded region shown in all plots. All simulation spectrum plots are shown in lin-log scale, and a dashed line indicates the dominant frequency of the stochastic multi-annual peak.
3.1. SIR model
The model's power spectrum computed from stochastic simulations for different values of the forcing amplitude β1 is shown in figure 5. The forced spectra exhibit two peaks: a dominant annual peak due to the deterministic limit cycle and a broad subdominant multi-annual stochastic peak similar to that of the unforced case. The contribution of annual epidemics in the time series increases as β1 increases, resulting in an enhanced annual peak. The multi-annual stochastic peak reflects the presence of noisy oscillations in the incidence time series with that dominant frequency. The mechanisms by which internal noise excites these resonant fluctuations has been treated in several papers [26,27,30,41]. They can be described analytically using van Kampen's expansion of the master equation of the underlying stochastic process around the attractor of the deterministic system [52].
Figure 5.
Figure 5.
The average power spectra from simulations of the SIR model for different amplitudes of seasonal forcing, β1. Parameters: green line, β1 = 0; orange line, 0.05; black line, 0.1. (Online version in colour.)
For the SIR model, the dominant period of stochastic fluctuations is 2.7 years. The stochastic peak lies in the shaded region, and its frequency and shape are largely independent of β1. Therefore, in this study, the dominant frequency of the main stochastic peak in the seasonally forced models can be computed from the unforced model whose analysis is easier (see the electronic supplementary material for details on the analytical computation of the power spectra). In §§3.2–3.5, we will discuss how the SIRS, SIRWS, SIRIS and SIRSI extensions of the SIR model are reflected in the position and amplitude of the stochastic peak.
3.2. SIRS model
In the SIRS model, the rate of immunity waning γ is a free parameter. In the case of lifelong immunity, γ = 0, the SIRS model reduces to the SIR model considered in the previous section. Figure 6 shows the model's spectrum for different values of the forcing amplitude β1 and two typical values of γ. In the unforced SIRS model, β1 = 0, the spectrum has a pronounced stochastic peak. However, it is situated outside the region of interest both for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i33.jpg years and for 1/γ = 20 years. For fixed γ and increasing β1, the frequency of the stochastic peak is slightly shifted towards the target region, but it remains outside this region even for large values of γ and, at the same time, the power associated with the annual peak becomes very large. Note that the dominant frequency of the stochastic peak would be even further away from the shaded region if β1 were decreased to take into account that non-primary infections do not undergo seasonal forcing. We conclude that the SIRS model's spectrum is compatible with the London and Ontario data only in the limit of An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i34.jpg when it approaches the SIR model. The performance of the SIRS model for other values of R0 is discussed in the electronic supplementary material.
Figure 6.
Figure 6.
The average power spectra from simulations of the SIRS model for different amplitudes of seasonal forcing, β1, and two values of the immunity waning rate, γ. Parameters: black line, 1/γ = 20 years; orange line, 40 years. (Online (more ...)
3.3. SIRIS model
Next, we analyse the power spectra from simulations of the SIRIS model. This model has three free parameters apart from β1: the rate of immunity loss γ, and the relative infectiousness η and relative susceptibility σ. For the limiting case of η = 0 and σ = 0, the model reduces to the SIRS model whose multi-annual peaks are located outside the shaded region. The results for fixed An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i35.jpg, 1/γ = 20 years and different values of η and σ are shown in figure 7. The spectrum for small values of η and σ is reminiscent of the SIRS model. A comparison of the spectra in figure 7ac shows that the stochastic and deterministic peaks are the highest for the smallest values of both η and σ. This is because it is easier to generate incidence oscillations in this parameter range. Increasing η and σ results in a flattening of the spectrum occurring through a gradual disappearance of, first, stochastic and, then, deterministic peaks. Moreover, varying the amplitude of seasonal forcing does not change this picture. For example, for the smallest values of η and σ used in figure 7, we observe that the position and the shape of the stochastic peaks agree perfectly for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i36.jpg and that they are always outside the region of interest (this result is given in the electronic supplementary material). The parameter region with higher values of η and σ is even less interesting because of the flattening of the spectrum. Thus for all relevant ranges of η, σ and β1, the stochastic multi-annual peak in this model is situated outside the region of interest except when it approaches the limit of the SIR model.
Figure 7.
Figure 7.
The average power spectra from simulations of the SIRIS model for β1 = 0.05, 1/γ = 20 years and different values of the relative infectiousness, η, and relative susceptibility, σ, of mild to severe infections. The spectra (more ...)
3.4. SIRWS model
Here, we show the power spectra computed from simulations of the SIRWS model for boosting coefficient κ = 20 and different values of γ and β1 (figure 8). In the absence of seasonality, β1 = 0, the spectra have a dominant multi-annual peak situated in the target region for the whole range of An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i37.jpg years. When β1 increases, a deterministic annual peak appears in the spectrum, while the dominant frequency of the multi-annual peak stays unchanged. Similar to the SIR model, the higher the amplitude of seasonal forcing β1, the higher the annual peak. However, unlike in that model, this phenomenon is accompanied by an overall complication of the structure of the spectrum for higher values of β1. For example, the power spectrum for 1/γ = 10 years and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i38.jpg (black line in figure 8c) has several secondary multi-annual peaks, one of which is situated near but outside the shaded region. Thus, the SIRWS model predicts a rather broad range of frequencies for the stochastic fluctuations.
Figure 8.
Figure 8.
The average power spectra from simulations of the SIRWS model for different amplitudes of seasonal forcing, β1, and two values of the immunity waning rate, γ. Parameters: black line, 1/γ = 10 years; orange line, 40 years. The immune (more ...)
Complementary results for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i39.jpg are given in the electronic supplementary material. For An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i40.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i41.jpg years, the stochastic multi-annual peaks of all power spectra lie in the shaded region. The same happens for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i42.jpg and 1/γ = 10 years but not for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i43.jpg. For An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i44.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i45.jpg, the SIRWS model exhibits a SIRS-like behaviour, and thus the multi-annual peaks of the SIRWS model are outside the region of interest.
3.5. SIRSI model
Finally, we present the results for the SIRSI model. Figure 9 shows the spectra for different values of An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i46.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i47.jpg and 1/γ = 20 years. As for SIR, the simulation spectra of the SIRSI model have a high stochastic multi-annual peak in addition to the annual peak. The amplitude of the deterministic peaks increases as β1 increases while the stochastic peaks remain almost unchanged. The increase, however, is smaller in the SIRSI model than in the SIR model. This means that, notwithstanding the similarity of their simulation spectra, the contribution of annual epidemics is smaller in the former than in the latter, corresponding to a qualitative difference in the time series of the two models.
Figure 9.
Figure 9.
The power spectrum from simulations of the SIRSI model. Parameters: green line, β1 = 0; orange line, 0.05; black line, 0.1; 1/γ = 20 years. The values of η and σ were chosen from a large set of possible values explored (more ...)
For the three sets of parameter values of figure 9, the dominant stochastic multi-annual peaks lie in the shaded region, and their frequencies and shapes are largely independent of β1. To explore the consequences of varying η and σ for fixed γ, we considered in the interval (0,1) nine equally spaced values and constructed a grid of 81 points in the (η, σ) space. As the position and the shape of the stochastic peak in this model is independent of β1, we can calculate the stochastic peak's frequency from the unforced model (see the electronic supplementary material for more details). For example, for 1/γ = 20 years used in figure 9, the number of spectra for the SIRSI model with the stochastic peak within the shaded region is 62 (for An external file that holds a picture, illustration, etc.
Object name is rsif20120432-i48.jpg years, this number varies from 43 to 78). We thus conclude that the SIRSI model's spectrum illustrated in figure 9 is robust with respect to variation of all free parameters, exhibiting in the accepted parameter ranges some variability in the amplitude of the stochastic peak. The robustness of these results with respect to variations of R0 too is illustrated in the electronic supplementary material.
3.6. Comparison with real data
The results for the power spectra of stochastic simulations show that, out of the four extensions of the SIR model under consideration, only SIRWS and SIRSI have frequency values for the stochastic multi-annual peak that preserve or improve the compatibility with data for pre-vaccination pertussis of the SIR model. Here we compare the other features of spectra of the SIR model and of the SIRWS and SIRSI extensions with the spectra for the London and Ontario datasets shown in figures 3 and and44.
We determine the population size N for each model, as a first step towards achieving this. As discussed before, both datasets have a similar average rate of new cases per month and thus correspond to similar effective populations. However, the deterministic counterparts of the models, and consequently the stochastic models too, predict different densities of infectives in the steady state. To resolve this, for each of the three models, we make simulations with the same length and sampling time as those of the London or Ontario datasets and record the incidence of the disease in the sampling interval. The transient period in all simulations is taken to be 50 years, and each simulation starts from a random initial condition. We calibrate the effective population size N for a given model by imposing that the incidence averaged over 103 different simulation runs is the same as the time averaged incidence in the corresponding dataset, given by the dashed line in figures 3 and and4.4. In figure 10, we show typical incidence time series for β1 = 0.05 and the effective population size N for each of the three models from which individual power spectra are computed.
Figure 10.
Figure 10.
Incidence of the disease in a simulation of the (a,b) SIRSI, (c,d) SIR and (e,f) SIRWS models. The horizontal dashed line is the average number of weekly (London) or monthly (Ontario) cases (compare these with the dashed lines in figures figures (more ...)
The second step is to determine for each model the level of seasonality that corresponds to the dataset. We estimate the value of the amplitude of the seasonal forcing β1 as the value for which the amplitude of the annual deterministic peak in the power spectrum of incidence time series averaged over 103 different simulation runs is equal to the amplitude of the annual peak in the spectrum computed from the data.
Using the values of N and of β1 determined in this way, we produce from the stochastic models time series mimicking the London or Ontario datasets, and average their power spectra over 103 simulation runs. The power spectra from simulations (dashed red (black) lines) and the spectra computed from the datasets (black solid lines) are shown in figure 11. The power spectra obtained with the SIR model have stochastic multi-annual peaks with amplitudes that are similar to, though on average larger than, those of the datasets' power spectra. For the SIRWS extension, the average amplitude of the stochastic peak is even larger, while the SIRSI extension gives a better agreement with this feature of the data with respect to the SIR model.
Figure 11.
Figure 11.
The solid black line is the power spectrum of the London or Ontario datasets as indicated in each panel. The dashed red (black) line is the average power spectrum of the incidence time series computed from 103 simulation runs with the same length and (more ...)
The results of figure 11 thus suggest that the SIRSI model is the one that best reproduces the stochastic properties of both datasets. However, as shown by the four examples plotted in the figure for each case (dashed grey lines), the poor statistics of the short-term samplings translate into a considerable variability in the power spectra computed from each time series. The two historical data records considered here illustrate the principle that the incidence power spectra can be used to discriminate the three variants of the model, but a final conclusion on their performance would require the analysis of longer time series.
Different assumptions about disease- and vaccine-induced immunity have been proposed in the recent literature to model the dynamics of pertussis in the presence of mass vaccination. Our motivation has been to study separately the influence of disease-acquired immunity only, in order to reduce the number of free parameters, and to consider explicitly the stochastic properties of the incidence time series as part of the model's output, in order to look for further constraints.
We have used the basic SIR model and four variants that reflect different immune responses. We characterized the dynamic patterns of each of these five models through the average power spectra of an ensemble of long-term prevalence time series obtained from stochastic simulations. We have found that the power spectra of the five alternative models show quantitative and even qualitative differences. One of them (SIRS) fails to reproduce a stochastic peak in the frequency range found in real data for pertussis, and another one (SIRIS) is too stiff to produce stochastic fluctuations with the observed amplitudes.
Assessing the performance of the other two variants, SIRWS and SIRSI, with respect to the SIR model requires a quantitative comparison of the simulated power spectra with real data; for that purpose, the effective population size and the forcing amplitude of each data record must be determined. We illustrated this procedure by considering two publicly available historical data records for pre-vaccination pertussis incidence. On the basis of these relatively short time series, we found that the output of the SIRSI model agrees with the main aspects of the phenomenology of the disease in a broad parameter range. It reproduces the multi-annual interepidemic periods and the amplitude of fluctuations found in the data with a better agreement than SIR for the latter, while the SIRWS variant tends to produce stochastic fluctuations even larger than SIR.
Considerations about the recovery profile and its influence on the model's behaviour concur in favouring models such as SIRSI, with a much flatter spectrum than SIRWS. For unforced models, it is known [29] that the effect of changing from an exponentially distributed recovery profile to a unimodal gamma distributed recovery profile with the same average is an enhancement of the fluctuations, together with a small displacement of the main stochastic peak towards higher frequencies. We have checked (results not shown) that this holds too for the SIRSI and SIRWS models in the parameter region explored in the previous section. Therefore, under more realistic recovery profiles, the output of the SIRWS model would be further away from the target power spectra, as would be that of the SIR model.
However, the short length of the time series used in our illustrative study prevents any definite conclusion about the performance of the SIR model alternatives. Indeed, the ensemble of the power spectra obtained from simulated time series of this length exhibits enough variability to accommodate the two samples of real data. An extension of this analysis to a larger set of data could help to further elucidate the dynamics of disease-acquired immunity, setting a solid ground for more complex models to deal with vaccine-acquired immunity.
Acknowledgments
The authors thank Gabriela Gomes for helpful discussions. Financial support from the Portuguese Foundation for Science and Technology (FCT) under contracts no. POCTI/ISFL/2/261, SFRH/BD/32164/2006 and SFRH/BPD/69137/2010 is gratefully acknowledged. G.R. was also supported by the Calouste Gulbenkian Foundation Program ‘Stimulus for Research’ and by the National Science Foundation under grant no. NSF PHY05-51164.
1. Trottier H., Philippe P., Roy R. 2006. Stochastic modeling of empirical time series of childhood infectious diseases data before and after mass vaccination. Emerg. Themes Epidemiol. 3, 9 (doi:10.1186/1742-7622-3-9) doi: 10.1186/1742-7622-3-9. [PMC free article] [PubMed] [Cross Ref]
2. Yorke J. A., London W. P. 1973. Recurrent outbreaks of measles, chickenpox and mumps. II. Systematic differences in contact rates and stochastic effects. Am. J. Epidemiol. 98, 469–482. [PubMed]
3. Wallinga J., Lévy-Bruhl D., Gay N. J., Wachmann C. H. 2001. Estimation of measles reproduction ratios and prospects for elimination of measles by vaccination in some Western European countries. Epidemiol. Infect. 127, 281–295 (doi:10.1017/S095026880100601X) doi: 10.1017/S095026880100601X. [PubMed] [Cross Ref]
4. Broutin H., Guégan J. F., Elguero E., Simondon F., Cazelles B. 2005. Large-scale comparative analysis of pertussis population dynamics: periodicity, synchrony, and impact of vaccination. Am. J. Epidemiol. 161, 1159–1167 (doi:10.1093/aje/kwi141) doi: 10.1093/aje/kwi141. [PubMed] [Cross Ref]
5. Broutin H., Mantilla-Beniers N. B., Simondon F., Aaby P., Grenfell B. T., Guégan J. F., Rohani P. 2005. Epidemiological impact of vaccination on the dynamics of two childhood diseases in rural Senegal. Microbes Infect. 7, 593–599 (doi:10.1016/j.micinf.2004.12.018) doi: 10.1016/j.micinf.2004.12.018. [PubMed] [Cross Ref]
6. Centers for Disease Control and Prevention 2002. Pertussis: United States, 1997–2000 Morb. Mortal Wkly Rep. 51, 73.–76. [PubMed]
7. Skowronski D. M., De Serres G., MacDonald D., Wu W., Shaw C., Macnabb J., Champagne S., Patrick D. M., Halperin S. A. 2002. The changing age and seasonal profile of pertussis in Canada. J. Infect. Dis. 185, 1448–1453 (doi:10.1086/340280) doi: 10.1086/340280. [PubMed] [Cross Ref]
8. de Melker H. E., Schellekens J. F., Neppelenbroek S. E., Mooi F. R., Rümke H. C., Conyn-van Spaendonck M. A. 2000. Reemergence of pertussis in the highly vaccinated population of the Netherlands: observations on surveillance data. Emerg. Infect. Dis. 6, 348–357 (doi:10.3201/eid0604.000404) doi: 10.3201/eid0604.000404. [PMC free article] [PubMed] [Cross Ref]
9. Crowcroft N. S., Pebody R. G. 2006. Recent developments in pertussis. Lancet 367, 1926–1936 (doi:10.1016/S0140-6736(06)68848-X) doi: 10.1016/S0140-6736(06)68848-X. [PubMed] [Cross Ref]
10. Lavine J. S., King A. A., Bjørnstad O. N. 2011. Natural immune boosting in pertussis dynamics and the potential for long-term vaccine failure. Proc. Natl Acad. Sci. USA 108, 7259–7264 (doi:10.1073/pnas.1014394108) doi: 10.1073/pnas.1014394108. [PubMed] [Cross Ref]
11. Soper H. 1929. The interpretation of periodicity in disease prevalence. J. R. Stat. Soc. 92, 34–73 (doi:10.2307/2341437) doi: 10.2307/2341437. [Cross Ref]
12. London W. P., Yorke J. A. 1973. Recurrent outbreaks of measles, chickenpox and mumps. I. Seasonal variation in contact rates. Am. J. Epidemiol. 98, 453–468. [PubMed]
13. Fine P. E., Clarkson J. A. 1982. Measles in England and Wales—I: an analysis of factors underlying seasonal patterns. Int. J. Epidemiol. 11, 5–14 (doi:10.1093/ije/11.1.5) doi: 10.1093/ije/11.1.5. [PubMed] [Cross Ref]
14. Fine P. E., Clarkson J. A. 1986. Seasonal influences on pertussis. Int. J. Epidemiol. 15, 237–247 (doi:10.1093/ije/15.2.237) doi: 10.1093/ije/15.2.237. [PubMed] [Cross Ref]
15. Grenfell B. T., Bjørnstad O. N., Kappey J. 2001. Travelling waves and spatial hierarchies in measles epidemics. Nature 414, 716–723 (doi:10.1038/414716a) doi: 10.1038/414716a. [PubMed] [Cross Ref]
16. Earn D. J. D., Rohani P., Bolker B. M., Grenfell B. T. 2000. A simple model for complex dynamical transitions in epidemics. Science 287, 667–670 (doi:10.1126/science.287.5453.667) doi: 10.1126/science.287.5453.667. [PubMed] [Cross Ref]
17. Bauch C. T., Earn D. J. D. 2003. Transients and attractors in epidemics. Proc. R. Soc. Lond. B 270, 1573–1578 (doi:10.1098/rspb.2003.2410) doi: 10.1098/rspb.2003.2410. [PMC free article] [PubMed] [Cross Ref]
18. Keeling M. J., Rohani P., Grenfell B. T. 2001. Seasonally forced disease dynamics explored as switching between attractors. Physica D 148, 317–335 (doi:10.1016/S0167-2789(00)00187-1) doi: 10.1016/S0167-2789(00)00187-1. [Cross Ref]
19. Wearing H. J., Rohani P., Keeling M. J. 2005. Appropriate modeling of infectious diseases. PLoS Med. 2, e239 (doi:10.1371/journal.pmed.0020239) doi: 10.1371/journal.pmed.0020239. [PMC free article] [PubMed] [Cross Ref]
20. Simões M., Telo da Gama M. M., Nunes A. 2008. Stochastic fluctuations in epidemics on networks. J. R. Soc. Interface 5, 555–566 (doi:10.1098/rsif.2007.1206) doi: 10.1098/rsif.2007.1206. [PMC free article] [PubMed] [Cross Ref]
21. Edmunds W. J., Kafatos G., Wallinga J., Mossong J. R. 2006. Mixing patterns and the spread of close-contact infectious diseases. Emerg. Themes Epidemiol. 3, 10 (doi:10.1186/1742-7622-3-10) doi: 10.1186/1742-7622-3-10. [PMC free article] [PubMed] [Cross Ref]
22. Read J. M., Eames K. T., Edmunds W. J. 2008. Dynamic social networks and the implications for the spread of infectious disease. J. R. Soc. Interface 5, 1001–1007 (doi:10.1098/rsif.2008.0013) doi: 10.1098/rsif.2008.0013. [PMC free article] [PubMed] [Cross Ref]
23. Hethcote H. W. 1997. An age-structured model for pertussis transmission. Math. Biosci. 145, 89–136 (doi:10.1016/S0025-5564(97)00014-X) doi: 10.1016/S0025-5564(97)00014-X. [PubMed] [Cross Ref]
24. Allen L. J., Driessche P. 2006. Stochastic epidemic models with a backward bifurcation. Math. Biosci. Eng. 3, 445–458 (doi:10.3934/mbe.2006.3.445) doi: 10.3934/mbe.2006.3.445. [PubMed] [Cross Ref]
25. Keeling M. J., Eames K. T. 2005. Networks and epidemic models. J. R. Soc. Interface 2, 295–307 (doi:10.1098/rsif.2005.0051) doi: 10.1098/rsif.2005.0051. [PMC free article] [PubMed] [Cross Ref]
26. Alonso D., McKane A. J., Pascual M. 2007. Stochastic amplification in epidemics. J. R. Soc. Interface 4, 575–582 (doi:10.1098/rsif.2006.0192) doi: 10.1098/rsif.2006.0192. [PMC free article] [PubMed] [Cross Ref]
27. Rozhnova G., Nunes A. 2009. Fluctuations and oscillations in a simple epidemic model. Phys. Rev. E 79, 041922 (doi:10.1103/PhysRevE.79.041922) doi: 10.1103/PhysRevE.79.041922. [PubMed] [Cross Ref]
28. Volz E., Meyers L. A. 2009. Epidemic thresholds in dynamic contact networks. J. R. Soc. Interface 6, 233–241 (doi:10.1098/rsif.2008.0218) doi: 10.1098/rsif.2008.0218. [PMC free article] [PubMed] [Cross Ref]
29. Black A. J., McKane A. J., Nunes A., Parisi A. 2009. Stochastic fluctuations in the susceptible–infective–recovered model with distributed infectious periods. Phys. Rev. E 80, 021922 (doi:10.1103/PhysRevE.80.021922) doi: 10.1103/PhysRevE.80.021922. [PubMed] [Cross Ref]
30. Rozhnova G., Nunes A. 2010. Stochastic effects in a seasonally forced epidemic model. Phys. Rev. E 82, 041906 (doi:10.1103/PhysRevE.82.041906) doi: 10.1103/PhysRevE.82.041906. [PubMed] [Cross Ref]
31. Metcalf C. J., Bjørnstad O. N., Grenfell B. T., Andreasen V. 2009. Seasonality and comparative dynamics of six childhood infections in pre-vaccination Copenhagen. Proc. R. Soc. B 276, 4111–4118 (doi:10.1098/rspb.2009.1058) doi: 10.1098/rspb.2009.1058. [PMC free article] [PubMed] [Cross Ref]
32. Grassly N. C., Fraser C. 2006. Seasonal infectious disease epidemiology. Proc. R. Soc. B 273, 2541–2550 (doi:10.1098/rspb.2006.3604) doi: 10.1098/rspb.2006.3604. [PMC free article] [PubMed] [Cross Ref]
33. Meyers L. A., Pourbohloul B., Newman M. E., Skowronski D. M., Brunham R. C. 2005. Network theory and SARS: predicting outbreak diversity. J. Theor. Biol. 232, 71–81 (doi:10.1016/j.jtbi.2004.07.026) doi: 10.1016/j.jtbi.2004.07.026. [PubMed] [Cross Ref]
34. Eubank S., Guclu H., Kumar V. S., Marathe M. V., Srinivasan A., Toroczkai Z., Wang N. 2004. Modelling disease outbreaks in realistic urban social networks. Nature 429, 180–184 (doi:10.1038/nature02541) doi: 10.1038/nature02541. [PubMed] [Cross Ref]
35. Broutin H., Viboud C., Grenfell B. T., Miller M. A., Rohani P. 2010. Impact of vaccination and birth rate on the epidemiology of pertussis: a comparative study in 64 countries. Proc. R. Soc. B 277, 3239–3245 (doi:10.1098/rspb.2010.0994) doi: 10.1098/rspb.2010.0994. [PMC free article] [PubMed] [Cross Ref]
36. Rohani P., Keeling M. J., Grenfell B. T. 2002. The interplay between determinism and stochasticity in childhood diseases. Am. Nat. 159, 469–481 (doi:10.1086/339467) doi: 10.1086/339467. [PubMed] [Cross Ref]
37. Nguyen H. T., Rohani P. 2008. Noise, nonlinearity and seasonality: the epidemics of whooping cough revisited. J. R. Soc. Interface 5, 403–413 (doi:10.1098/rsif.2007.1168) doi: 10.1098/rsif.2007.1168. [PMC free article] [PubMed] [Cross Ref]
38. Ferrari M. J., Grais R. F., Bharti N., Conlan A. J., Bjørnstad O. N., Wolfson L. J., Guerin P. J., Djibo A., Grenfell B. T. 2008. The dynamics of measles in sub-Saharan Africa. Nature 451, 679–684 (doi:10.1038/nature06509) doi: 10.1038/nature06509. [PubMed] [Cross Ref]
39. Mantilla-Beniers N. B., Bjørnstad O. N., Grenfell B. T., Rohani P. 2010. Decreasing stochasticity through enhanced seasonality in measles epidemics. J. R. Soc. Interface 7, 727–739 (doi:10.1098/rsif.2009.0317) doi: 10.1098/rsif.2009.0317. [PMC free article] [PubMed] [Cross Ref]
40. Wearing H. J., Rohani P. 2009. Estimating the duration of pertussis immunity using epidemiological signatures. PLoS Pathog. 5, e1000647 (doi:10.1371/journal.ppat.1000647) doi: 10.1371/journal.ppat.1000647. [PMC free article] [PubMed] [Cross Ref]
41. Black A. J., McKane A. J. 2010. Stochasticity in staged models of epidemics: quantifying the dynamics of whooping cough. J. R. Soc. Interface 7, 1219–1227 (doi:10.1098/rsif.2009.0514) doi: 10.1098/rsif.2009.0514. [PMC free article] [PubMed] [Cross Ref]
42. Rohani P., Zhong X., King A. A. 2010. Contact network structure explains the changing epidemiology of pertussis. Science 330, 982–985 (doi:10.1126/science.1194134) doi: 10.1126/science.1194134. [PubMed] [Cross Ref]
43. Aguas R., Gonalves G., Gomes M. G. 2006. Pertussis: increasing disease as a consequence of reducing transmission. Lancet Infect. Dis. 6, 112–117 (doi:10.1016/S1473-3099(06)70384-X) doi: 10.1016/S1473-3099(06)70384-X. [PubMed] [Cross Ref]
44. Gomes M. G., White L. J., Medley G. F. 2004. Infection, reinfection, and vaccination under suboptimal immune protection: epidemiological perspectives. J. Theor. Biol. 228, 539–549 (doi:10.1016/j.jtbi.2004.02.015) doi: 10.1016/j.jtbi.2004.02.015. [PubMed] [Cross Ref]
45. Anderson R. A., May R. M. 1991. Infectious diseases of humans. Oxford, UK: Oxford University Press.
46. Clarkson J. A., Fine P. E. M. 1985. The efficiency of measles and pertussis notification in England and Wales. Int. J. Epidemiol. 14, 153–168 (doi:10.1093/ije/14.1.153) doi: 10.1093/ije/14.1.153. [PubMed] [Cross Ref]
47. Rohani P., Earn D. J. D., Grenfell B. T. 1999. Opposite patterns of synchrony in sympatric disease metapopulations. Science 286, 968–971 (doi:10.1126/science.286.5441.968) doi: 10.1126/science.286.5441.968. [PubMed] [Cross Ref]
48. Wendelboe A. M., Van Rie A., Salmaso S., Englund J. A. 2005. Duration of immunity against pertussis after natural infection or vaccination. Pediatr. Infect. Dis. J. 24, S58–S61 (doi:10.1097/01.inf.0000160914.59160.41) doi: 10.1097/01.inf.0000160914.59160.41. [PubMed] [Cross Ref]
49. Gillespie D. T. 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22, 403–434 (doi:10.1016/0021-9991(76)90041-3) doi: 10.1016/0021-9991(76)90041-3. [Cross Ref]
50. Anderson D. F. 2007. A modified next reaction method for simulating chemical systems with time dependent propensities and delays. J. Chem. Phys. 127, 214107 (doi:10.1063/1.2799998) doi: 10.1063/1.2799998. [PubMed] [Cross Ref]
51. Lu T., Volfson D., Tsimring L., Hasty J. 2004. Cellular growth and division in the Gillespie algorithm. Syst. Biol. 1, 121–128 (doi:10.1049/sb:20045016) doi: 10.1049/sb:20045016. [PubMed] [Cross Ref]
52. van Kampen N. G. 1992. Stochastic processes in physics and chemistry. Amsterdam, The Netherlands: Elsevier.
Articles from Journal of the Royal Society Interface are provided here courtesy of
The Royal Society