INTRODUCTION
Pertussis, also known as whooping cough, is a vaccine-preventable disease that has generated a great deal of concern recently due to its resurgence in many countries that maintain high vaccination coverage (Yih et al. Reference Yih, Lett, des Vignes, Garrison, Sipe and Marchant2000; Celentano et al. Reference Celentano, Massari, Paramatti, Salmaso and Tossi2005; Rohani & Drake, Reference Rohani and Drake2011; Jackson & Rohani, Reference Jackson and Rohani2013). Candidate explanations range from the vaccine-driven evolution of the aetiological agent Bordetella pertussis (Mooi et al. Reference Mooi, van der Maas and de Melker2013), to increased circulation of congeners (Pittet et al. Reference Pittet, Emonet, Schrenzel, Siegrist and Posfay-Barbe2014), changes in reporting and surveillance (Cherry, Reference Cherry2012) and loss of vaccine efficacy due to the switch from whole-cell pertussis (wP) vaccine to the less reactogenic acellular pertussis (aP) vaccines in the mid-1990s (Shapiro, Reference Shapiro2012). With so much uncertainty still surrounding basic issues in pertussis epidemiology, mathematical models are invaluable as tools for the synthesis of epidemiological data, the quantitative evaluation of different explanations and as aids in identifying alternative mechanisms that may be at play (Lavine & Rohani, Reference Lavine and Rohani2012; Blackwood et al. Reference Blackwood, Cummings, Broutin, Iamsirithaworn and Rohani2013). Recently, it has been suggested that resurgence may be the foreseeable consequence of decades of incomplete vaccination with an imperfect vaccine, an explanation that does not presuppose changes in the underlying transmission biology (Riolo et al. Reference Riolo, King and Rohani2013). Realistic age-structured transmission models robustly exhibit a so-called ‘honeymoon period’ – a prolonged interval of especially low disease incidence – following the inception of mass vaccination with an imperfect vaccine (McLean, Reference McLean1998). With parameters representative of pertussis transmission and immunity, such models indicate that this period can extend for decades, its conclusion being marked by a rebound in incidence perhaps comparable with the current resurgence (Mossong & Muller, Reference Mossong and Muller2003; Heffernan & Keeling, Reference Heffernan and Keeling2009).
It has been shown that the nature of vaccine-induced protection can leave distinct footprints in the transient disease incidence patterns following the roll out of a vaccination programme (Magpantay et al. Reference Magpantay, Riolo, Domenech de Cellès, King and Rohani2014). In particular, the mode of vaccine failure can determine the depth and duration of the honeymoon period, as well as characteristics of the resurgence. In this study, we exploited this insight to attempt to ascertain the nature and degree of protection provided by the aP vaccine by analysing the dynamics of pertussis during a period over which vaccine coverage was ramped up. To cover the range of plausible effects of the aP vaccine on immunity to pertussis, we considered the following aspects of vaccine protection:
Immunity against infection
The ideal vaccine protects against transmissible infection. Three, not mutually exclusive, modes by which vaccines might fail in this goal are: (a) Primary vaccine failure. A vaccine exhibits primary vaccine failure if it fails to provide any form of protection against infection to some fraction of vaccinated individuals. Primary vaccine failure is quantified by the fraction of vaccinated individuals the vaccine fails to protect. (b) Leakiness. A vaccine is said to be leaky when it reduces, but does not eliminate, the potential for infection (Halloran et al. Reference Halloran, Haber and Longini1992). The leakiness of a vaccine is measured by the probability of infection upon exposure for a vaccinated individual relative to the same probability for an unvaccinated individual. (c) Waning. Vaccine-induced protection is said to wane when it ceases after some time. Here, we quantify the speed of waning by the mean duration of protection or, equivalently, its reciprocal, the rate at which immunity is lost.
Immunity against transmission and disease
Even when a vaccine fails to protect against infection, it might still reduce the infection's transmissibility and/or the severity of disease symptoms. These effects are quantified by (a) relative infectiousness, which we model as the ratio of the transmission rate of a vaccinated person to that of an unvaccinated person, and (b) relative reporting probability, the ratio of the reporting probability of a vaccinated person relative to that of an unvaccinated person. The latter can reflect the degree to which infections in vaccinated individuals produce milder disease symptoms.
There are currently no known serologic correlates of protection to pertussis (Plotkin, Reference Plotkin2010) and therefore much uncertainty regarding the nature, degree, and duration of protection provided by pertussis vaccines. However, because immune memory shapes epidemiological dynamics, this uncertainty can be reduced by fitting mechanistic models to time series data (Lavine & Rohani, Reference Lavine and Rohani2012). We constructed stochastic transmission models that account for potential vaccine-induced protection against infection, transmission and disease. These models were confronted with time series incidence data from six regions of Italy during a period of dramatic change in the Italian national immunization coverage (1996–2009). Initially, national coverage with wP vaccine was approximately 30%. This was ramped up rapidly in the mid-1990s when the country switched to the aP vaccine. Coverage further increased in the first decade of the new century, when vaccines were made available free of charge, reaching an average of 95% by 2009 (Gonfiantini et al. Reference Gonfiantini, Carloni, Gesualdo, Agricola, Rizzuto, Iannazzo, Ciofi Degli Atti, Villani and Tozzi2014). Using recently developed statistical inference techniques (King et al. Reference King, Ionides, Bretó, Ellner, Ferrari, Kendall, Lavine, Nguyen, Reuman, Wearing and Wood2015b ; Ionides et al. Reference Ionides, Nguyen, Atchadé, Stoev and King2015), we estimated the values of model parameters needed to explain the dynamics of pertussis incidence over this period of abrupt change in vaccine coverage.
The evidence we describe below suggests that, in the absence of primary vaccine failure, vaccinated individuals whose protection against infection has failed are unlikely to be recorded as cases (possibly due to vaccine-induced protection against severe disease) but may be just as infectious as unvaccinated individuals. Under the assumption of zero primary vaccine failure, the best models point to substantial aP-induced protection against infection, concomitant reduction in pathogen circulation and considerable herd immunity. However, the data provided insufficient information to allow us to identify the mode of vaccine failure. Specifically, a range of models incorporating varying degrees of leakiness and rates of waning were roughly equally well-supported by the data, as measured by likelihood. We relaxed our assumption of zero primary failure, considering models with modest levels of aP primary vaccine failure. The data were incompatible with even 15% primary vaccine failure. Moreover, as the rate of primary failure varied, substantial differences in predicted age distribution of incidence appeared. This implies that age-specific incidence data of sufficiently high resolution – unavailable to us – contain the information needed to identify not only the rate of primary vaccine failure, but also where the aP vaccine lies along the leaky/waning spectrum. Finally, we examined model-predicted epidemiological futures under two extreme versions of the well-supported models, revealing that quite distinct future dynamical scenarios are compatible with the data in hand. This implies that similar studies directed to locations and periods with different dynamics may contain the information needed, again, to resolve the remaining issues. In sum, this work shows how, by confronting mechanistic models to disease dynamics data, one can reduce uncertainty and gain insight into key immunological determinants of epidemiology, and also locate the limits to sound inference and determine the nature of the data needed to expand them.
MATERIALS AND METHODS
Data
Pertussis notification data were obtained from the Italian Ministry of Health (Ministero della Salute, 2014a ). The data were available at the regional level, with monthly reports from the beginning of 1996 until the end of 2009. We used data from Lazio, Lombardia, Sardegna, Sicilia, Toscana and Umbria due to the separation of their major cities, and their geographic distribution from the North to the South of the country. We obtained 1990–2012 regional demographic data (population sizes, annual numbers of live births and deaths) from Eurostat (European Commission, 2014) and regional vaccine coverage data from the Ministry of Health (Ministero della Salute, 2014b ). The value of the coverage at each year was defined as the proportion of children born that year who received three or more doses of the combined diphtheria, tetanus and aP vaccine (DTP) by 24 months of age. Since the vaccine schedule in Italy prescribes that the three doses be taken by 11 months of age, the delay between birth and three doses of DTP should be less than 24 months, on average, and closer to 11 months. These data were only available from 2001 to 2012. Lacking data on coverage prior to 2001, we made the pragmatic assumption that coverage was at the national average of 30% in 1994 and linearly ramped up to the first recorded level in 2001. A plot of the notification data and vaccination coverage for the six regions is presented in Fig. 1. A summary of the features of each region is given in Table 1.
1996–2009 pertussis notification reports and vaccine coverage for the selected regions in Italy. Dashed lines indicate our assumption concerning the ramp-up in coverage prior to 2001.
Summary features of the six regions in Italy considered
All demographic values were calculated from 1996–2009 annual Eurostat data (European Commission, 2014).
Model
Building upon the standard susceptible-exposed-infected-recovered model (Keeling & Rohani, Reference Keeling and Rohani2008), we constructed a model of pertussis with eight compartments. It includes two compartments each of susceptible (S i ), exposed (E i ) and infected (I i ) compartments in order to distinguish individuals who were never vaccinated (i = 1) from those who were vaccinated (i = 2). All recoveries from infection go to the recovered (R) compartment. The vaccinated (V) compartment contained individuals who were vaccinated and still maintain some vaccine-derived protection against infection. The model schematic is shown in Fig. 2 and the system of equations are given in the supplementary material (Section S1).
Diagram of the Full Model of pertussis with eight compartments. The S 1, E 1 and I 1 compartments consist of the susceptible, exposed and infected individuals, respectively, who were never vaccinated. The S 2, E 2 and I 2 compartments are the corresponding compartments for those who were vaccinated. All individuals recovering from infection go to the R class. Vaccinated individuals who do not experience primary vaccine failure enter the V compartment and leave it via the leaky (blue) or waning route (red). For most of this paper we focus on the case when there is no primary vaccine failure (ε A = 0). We considered two restrictions of this model: in the Waning Model, ε L = 0, while in the Leaky Model, α = 0 (refer to Table 4). Both restricted models make the further assumption – suggested by fits of the Full Model to data – that η = 0 and θ = 1, i.e. that post-vaccine infections are perfectly in apparent and as transmissible as infections in naive individuals.
At every time step, a fraction p(t) (corresponding to the vaccine coverage) of newborns are vaccinated. The vaccine is assumed to ‘fail to take’ (McLean & Blower, Reference McLean and Blower1993) in a fraction, ε A , of these newborns: ε A is thus the amount of primary vaccine failure. The vaccinated individuals who do not experience primary vaccine failure are added into the V compartment while the remainder of the newborns go to the S 1 class. For most of this paper we assume that there is no primary vaccine failure (ε A = 0). We relax this assumption in the discussion.
Individuals can leave the V class either by being infected and going directly into the E 2 class (the leaky route) or by losing their vaccine protection to join the S 2 class (the waning route). The rate by which transitions through the leaky route occur is modulated by the leakiness parameter ε L , which is the relative probability of infection upon exposure of individuals in the V class relative to those in the susceptible classes (S 1 or S 2). The rate by which transitions occur along the waning route is controlled by the vaccine waning rate α. The other two vaccine parameters in the model measure the relative infectiousness θ and relative reporting probability η of the infected vaccinated individuals (I 2), relative to infected individuals who were either never vaccinated or experienced primary vaccine failure (I 1). These last two parameters are discussed further in the descriptions of the force of infection and reporting model.
The exposed (E 1 and E 2) and infected classes (I 1 and I 2) were each broken down into three subcompartments in order to yield Erlang-distributed latent and infectious periods (Lloyd, Reference Lloyd2001; Wearing et al. Reference Wearing, Rohani and Keeling2005; Keeling & Rohani, Reference Keeling and Rohani2008). The entire model was set up as a system of 16 stochastic difference equations and simulated using a time step of 0·01 year. Transitions between compartments were simulated using multinomial samples with rates as indicated in Fig. 2. Births were added to the S 1 or V compartments using the exact birth and assumed vaccination rates. Changes in the total population size were accounted for by rescaling all components proportionally at each step based on the known death rate and a calculated immigration rate. We refer to this model in the following as the Full Model.
The model was initialized in the beginning of 1994, 2 years before the monthly pertussis notification time series commences, coincident with increases in vaccine uptake according to national estimates (Gonfiantini et al. Reference Gonfiantini, Carloni, Gesualdo, Agricola, Rizzuto, Iannazzo, Ciofi Degli Atti, Villani and Tozzi2014), and the switch to aP vaccine (Rota et al. Reference Rota, D'Ancona, Massari, Mandolini, Giammancoc, Carbonari, Salmaso and Ciofi degli Atti2005; Gabutti & Rota, Reference Gabutti and Rota2012). Since it is widely considered that the wP vaccine induces immunity similar to that induced by infection, the model does not distinguish between vaccine-derived and infection-derived immunity prior to 1994 (Ryan et al. Reference Ryan, Murphy, Ryan, Nilsson, Shackley, Gothefors, Oymar, Miller, Storsaeter and Mills1998). For simplicity, infection- and wP-derived immunity was assumed to be perfect and lifelong, consistent with conclusions from other studies (Wearing & Rohani, Reference Wearing and Rohani2009; Blackwood et al. Reference Blackwood, Cummings, Broutin, Iamsirithaworn and Rohani2013). Though it is widely assumed that infection-derived immunity lasts longer than vaccine-derived immunity, over the 16-year period we examined, this assumption should have at most minor consequences. Individuals who have gained or lost infection-derived immunity prior to the start of the simulations were accounted for in the fitting of initial conditions for S 1, E 1, I 1 and R. The other components (V, S 2, E 2, I 2) were initialized at zero since these components become active only with the commencement of aP vaccination.
The transmission rate β(t) was assumed to be a periodic function of time, with a period of 1 year. The mean transmission rate is given by the parameter β 1 and other parameters control the amplitude of seasonality (β 2) and the peak timing in transmission (ϕ). The form of this function is given in Section S1 and alternative formulations are discussed in Section S8 in the supplementary material. Taking into account the relative infectiousness of the I 2 class, the force of infection experienced by susceptible individuals is,
$$\eqalign{ & \lambda (t) = \left( {\displaystyle{{\beta (t)({I_1} + \theta {I_2}) + {\beta _1}\iota} \over {N(t)}}} \right)\Delta Q, \cr & \quad \Delta Q \sim \Gamma \left( {{\rm shape =} \displaystyle{{\Delta t} \over {{\beta ^2}_{{\rm S}{\rm. D}{\rm.}}}}, \,{\rm scale =} \displaystyle{{{\beta ^2}_{{\rm S}{\rm. D}{\rm.}}} \over {\Delta t}}} \right).}$$
Here N(t) is the known population size and ι is the number of infected individuals who immigrate to the region per year (assumed to be at a constant rate). The force of infection allows for environmental stochasticity generated by drawing ∆Q from a gamma white noise process with intensity controlled by the parameter β s.d. and the stepsize ∆t. Here, ∆Q has unit mean and variance equal to β 2 s.d. /Δt.
Unvaccinated individuals were assumed to be reported as cases with probability ρ during their transition from the infected I 1 class to the recovered class R. Vaccinated individuals were reported as cases with probability ηρ (η represents relative reporting probability) in their transition from the I 2 class to R. To accommodate variability in case reporting and allow for overdispersion, we used a negative-binomial observation model. Specifically, if the number of transitions from I 1 to R in 1 month is denoted by C 1, and the number from I 2 to R by C 2, the number of cases reported is assumed to follow a negative binomial distribution with mean C = ρ (C 1 + η C 2) and variance C + τ 2 C 2.
The full set of model parameters and their allowed ranges are listed in Tables 2 and 3. To understand the role of stochasticity in this model, we considered both deterministic and stochastic versions of the model. The deterministic versions consisted of a set of difference equations with rescaling adjustments made to correct for population size at each step. In this case, the environmental stochasticity parameter β s.d. was set to zero and for any fixed set of parameters, the only source of variation in reports is the reporting model. The systems of equations involved in the deterministic model can be found in the supplementary material.
Description of demographic inputs and disease parameters that are fixed in the model
Description of fitted disease parameters, vaccine parameters and initial conditions
Special cases of the model
Because the parameters of the Full Model were not uniquely identifiable with the data at hand (see Results), we examined the extremes of the spectrum of well-supported models in order to see how model parameters varied across this spectrum. Specifically, the two Restricted Models we considered were: (1) the Leaky Model wherein the only vaccine parameter allowed to vary is the leakiness and vaccine protection is permanent, and (2) the Waning Model wherein the vaccine gives full protection which wanes with time. We further constrained these Restricted Models by fixing the two vaccine parameters that were well-identified. In particular, for infections among the vaccinated population, we fixed the relative transmissibility at 1 and the relative reporting probability at 0 (see Results). A summary of the models considered is given in Table 4. We considered both deterministic and stochastic forms of all of these models. The differences in the implementation of the two forms are covered in Section S1 in the supplementary material.
Summary of the differences between the Full Model and the two Restricted Models (Leaky- and Waning-Models)
In these models primary vaccine failure was set to zero. The values of the relative infectiousness and relative reporting probability for the Restricted Models were fixed at values determined from the results of fitting the Full Model (see Results).
Trajectory matching
Nine region-specific parameters (including three initial conditions) and four global vaccine parameters (three of these are fixed in the Restricted Models) were estimated by fitting the deterministic models using maximum likelihood (ML) estimation via trajectory matching (King et al. Reference King, Ionides, Bretó, Ellner, Ferrari, Kendall, Lavine, Nguyen, Reuman, Wearing and Wood2015b ). In this case, the negative binomial reporting model is the only source of variability in observations for any fixed set of parameters, thus likelihoods could be directly calculated and maximized using the Nelder–Mead algorithm. The search was initiated over 104 initial points generated using Latin hypercube sampling over the parameter ranges indicated in Table 3. The best fits from this initial run were used as starting points for subsequent iterations of trajectory matching to reveal the ML estimate.
Likelihood profiles over each of the vaccine parameters were computed for each region. A profile is derived by creating an array of fixed values for a parameter and maximizing the likelihood at each fixed value over the remaining model parameters. The array of values in our profiles consisted of 100 points spread out across the entire allowed parameter range (see Table 2). These points were derived by splitting the range into 100 equally-spaced sub-intervals on a logarithmic scale, collecting the best points from the initial search that fall within each interval of the array, then optimizing the likelihood over all other parameters using trajectory matching. These initial profiles were refined by employing simple continuation techniques to derive new starting points for further iterations (at least three) of the profiling procedure. The results across regions were aggregated by assuming independence of the monthly reports across the different regions and adding the log-likelihoods. More details on the profiling procedure is given in Section S2 in the supplementary material.
Maximization via iterated filtering
Ten region-specific parameters (including three initial conditions) and four global vaccine parameters (three of these are fixed in the Restricted Models) were estimated by fitting the stochastic model using the second-generation Iterated Filtering algorithm (IF2, Ionides et al. Reference Ionides, Nguyen, Atchadé, Stoev and King2015) implemented in the R package pomp (King et al. Reference King, Ionides, Bretó, Ellner, Ferrari, Kendall, Lavine, Nguyen, Reuman, Wearing and Wood2015b , in press). The only extra parameter estimated in this case is the environmental stochasticity β s.d. . For each region in Italy, an iterated filtering search was initialized from 2000 points generated using Latin hypercube sampling over the range of allowed parameter ranges given in Table 2. This initial run was conducted using 50 IF2 iterations and a hyperbolic cooling fraction of 0·25. The random walk standard deviation for estimated parameters was set to 0·05. Full details are available in the scripts that reproduce our results (http://dx.doi.org/10.5061/dryad.58q00).
Confidence intervals using the aggregated profiles over each vaccine parameter were also derived. These profiles were refined using multiple rounds of the profiling procedure (up to six) described in the previous section.
Model selection
We used the Akaike Information Criterion with small sample bias adjustment (AIC c ) to select the best parsimonious model of pertussis. This is discussed in Section S3. There were a total of 1008 data points (6 regions × 14 years × 12 months per year). A total of 58 free parameters for the deterministic form and 64 for the stochastic form of the Full Model were fitted to data from all six considered regions.
RESULTS
The results of estimating parameters for the Full Model are presented first, followed by the results of fitting simplified versions of this model (the Restricted Models). Each Restricted Model was derived from the Full Model by restricting attention to a single mode of vaccine failure and fixing two other model parameters (θ, η) at values determined by fits of the Full Model (see Table 4).
Full Model
The profiles over the relative infectiousness and relative reporting probability are shown in Fig. 3, for both the deterministic and stochastic versions of the Full Model. From this we can see immediately that the deterministic profile leads to much lower likelihoods and higher AIC c , indicating that deterministic models do not explain the data as well as the stochastic models (∆AIC c = 59). The aggregated profile over the relative infectiousness parameter (Fig. 3A) strongly suggests that the infectiousness of vaccinated infected individuals is equal to that of unvaccinated ones (θ = 1). The results (Fig. 3B) also favour a small relative reporting probability (η = 0). These results are also supported by the profiles over individual regions (refer to Section S4 for details), which yielded consistent estimates for both the deterministic and stochastic models. The maximum likelihood estimates (MLEs) for these parameters and the confidence intervals are given in Table 5.
Plots of the log-likelihood values (profiles) corresponding to the MLEs of the Full Model at each fixed value of (A) relative infectiousness and (B) relative reporting probability. The plots for the stochastic models are 40–60 units higher than those for the corresponding deterministic models. The 95% confidence intervals about the parameter values that yield the highest likelihood values are shaded and overlap for the different models. These intervals are also presented in Table 5. Abbreviation: MLEs, maximum likelihood estimates.
Maximum likelihood estimates of the vaccine parameters for the Full Model. Here the leakiness ε L and waning rate α were unidentifiable. The best AICc value is shown in boldface.
Abbreviation: AIC, Akaike Information Criterion.
Profiling over the two other vaccine parameters that were allowed to vary (leakiness and waning rate) yielded broad confidence intervals that spanned almost the entire range of allowed values. This was true for both the deterministic and stochastic versions of the model, and indicated that the Full Model has too many degrees of freedom for these parameters to be uniquely identified. In other words, there was a spectrum of models, differing in some details and not in others, that were all roughly equally well-supported by the data. To understand how these models differ, we examined the extreme cases of purely leaky and purely waning vaccine.
Restricted Models: leaky vs. waning immunity
The Restricted Models are constrained forms of the Full Model where the relative infectiousness parameter, θ, is set to one and the relative reporting probability, η, is set to zero. We derived a Leaky Model (ε L ∈ [0, 1] and α = 0) and Waning Model (ε L = 0 and α ⩾ 0) so that the two different routes of vaccine-derived immunity failure could be considered separately (see also Table 4). The profiles over the corresponding vaccine parameter for each model are shown in Fig. 4. As before, the stochastic profiles attain much higher likelihoods (and lower AIC c ) than do the deterministic profiles. Accordingly, we devote no further attention to the results of the deterministic models.
Plots of the log-likelihood values (profiles) corresponding to the MLEs of (A) the Leaky Model at each fixed value of the leakiness, and (B) the Waning Model at each fixed value of the waning rate. Both models have the relative infectiousness θ = 1 and relative reporting probability η = 0. The plots for the stochastic models are 50–80 units higher than those for the corresponding deterministic models. The 95% confidence intervals about the parameter values that yield the highest likelihood values are shaded and are also given in Table 5.
A summary of the results for the vaccine parameters are given in Tables 5 and 6. Further details on the profiles for each individual region are shown in Section S4. Additionally, the robustness of our findings against changes on the assumptions in vaccination coverage is demonstrated in Section S9.
Maximum likelihood estimates of the vaccine parameters for the Restricted Models. In these models the relative infectiousness is fixed at θ = 1 and the relative reporting probability is fixed at η = 0 following the results of the Full Model. The best AICc values are shown in boldface.
DISCUSSION
Failure of vaccine-induced immunity has been proposed as one potential explanation for pertussis’ resurgence as seen in several populations with high routine vaccine coverage. Existing assessments of the effectiveness of pertussis vaccines, however, vary across studies. Some have proposed that both the aP (Rohani et al. Reference Rohani, Zhong and King2010; Domenech de Cell ès et al. Reference Domenech de Cellès, Riolo, Magpantay, Rohani and King2014; Gambhir et al. Reference Gambhir, Clark, Cauchemez, Tartof, Swerdlow and Ferguson2015) and wP (Broutin et al. Reference Broutin, Guegan, Elguero, Simondon and Cazelles2005; Blackwood et al. Reference Blackwood, Cummings, Broutin, Iamsirithaworn and Rohani2013) vaccines confer substantial long-lasting protection from infection. Others, arguing from different lines of evidence, have maintained that the aP vaccine reduces disease severity without considerably reducing pathogen transmission (Warfel et al. Reference Warfel, Zimmerman and Merkel2013; Smallridge et al. Reference Smallridge, Rolin, Jacobs and Harvill2014). These different modes of vaccine failure have widely different implications for the control of the disease (Riolo & Rohani, Reference Riolo and Rohani2015; Althouse & Scarpino, Reference Althouse and Scarpino2015). In this paper, we attempted to better resolve the nature of the immunity induced by aP and how often and by what mode it fails. Specifically, we endeavoured to identify whether and to what extent aP protects against infection, transmission and disease as well as how durable and how leaky is its protection. Recognizing that these parameters have implications for transient disease dynamics (Riolo et al. Reference Riolo, King and Rohani2013; Magpantay et al. Reference Magpantay, Riolo, Domenech de Cellès, King and Rohani2014), we sought to estimate them on the basis of time series displaying the dynamics of pertussis incidence and vaccination coverage. Accordingly, we formulated the alternative hypotheses as process-based models of pertussis transmission, vaccination and immunity, and challenged them to explain a panel of regional incidence data from Italy reflecting the immediate aftermath of that country's switch to the aP vaccine. We then used likelihood-based methods to quantify their relative explanatory power and to infer the nature, magnitude and durability of vaccine-induced immunity.
In the models we entertained, vaccine effects are described by four vaccine parameters: (1) probability of protection (inversely related to leakiness), (2) rate of waning (inversely related to durability), (3) vaccine-induced reduction in transmissibility of infection and (4) reduced reporting probability (such as might be produced by reduction in disease severity). We fitted the model to 1996–2009 data from six regions in Italy distributed across the country. In all regions, pertussis reports were declining over most of this period, though the regions had different population sizes, birth rates, overall population growth (Table 1) and vaccine coverage (Fig. 1). The models were fit to each region separately, allowing us to seek a comprehensive estimate of the vaccine parameters under the assumption of independence of the dynamics between regions. Parameters were thus estimated using aggregated likelihood profiles.
Relative infectiousness and relative reporting probability
In all aggregated parameter profiles, the stochastic models consistently had higher likelihoods and lower AIC c than their deterministic counterparts, indicating that stochastic models provide better descriptions of the dynamics during the period of study. In the case of the Full Model, the majority of the individual region profiles agree in concluding that an infection in a vaccinated individual is as transmissible as an infection in a naive individual (θ = 1) while the relative reporting probability is very small (see Section S4). We argue that the very low relative reporting probability is best interpreted as indicating that the preponderance of infections in the vaccinated population are mild or asymptomatic. This suggests that even if the aP vaccine provides little protection against infection, it does benefit individuals directly by diminishing disease severity.
These findings conflict with the classical belief that asymptomatic pertussis infections contribute little to transmission (Schellekens et al. Reference Schellekens, von König and Gardner2005). However, they are compatible with the observation that in clinical cases, infections are most communicable during early stages of infection when symptoms tend to be mild and not clearly distinguishable from minor respiratory infections (Edwards & Decker, Reference Edwards and Decker2013). Furthermore, experimental infections of baboons with pertussis (Warfel et al. Reference Warfel, Zimmerman and Merkel2013) showed that aP-vaccinated animals could transmit the infection (suggesting relative infectiousness greater than zero) without displaying clinical disease symptoms (consistent with low reporting probability), a finding compatible with our results. Unfortunately, the small sample size in the Warfel et al. study precluded identification of differences in the relative transmissibility of the infections in vaccinated animals, and the experimental design afforded no measurement of the effects of aP vaccination on the durability or leakiness of vaccine protection.
Of the four vaccine parameters that we allowed to vary, three (relative infectiousness, leakiness and waning rate) affect the level of pathogen circulation in a vaccinated population. High levels of transmission are manifested in short inter-epidemic periods and a concentration of cases among the younger population (Anderson & May, Reference Anderson and May1991; Blackwood et al. Reference Blackwood, Cummings, Broutin, Iamsirithaworn and Rohani2012; Magpantay & Rohani, Reference Magpantay and Rohani2015). Recent studies have attributed observed epidemiological shifts – an increase in the inter-epidemic period, a pronounced drop in infant incidence and a rightward shift in the age distribution of cases – to a considerable decline in transmission associated with vaccination (Wearing & Rohani, Reference Wearing and Rohani2009; Rohani et al. Reference Rohani, Zhong and King2010; Domenech de Cellès et al. Reference Domenech de Cellès, Riolo, Magpantay, Rohani and King2014). Since the present study finds no evidence for vaccine-induced reduction of infection transmissibility, we propose that any decline in transmission must instead be due to decreased susceptibility of vaccinated individuals.
Vaccine impact
Many studies have considered different types of infection- and vaccine-derived immunity (Halloran et al. Reference Halloran, Haber and Longini1992; McLean & Blower, Reference McLean and Blower1993; Farrington, Reference Farrington2003; Gomes et al. Reference Gomes, White and Medley2004; Magpantay et al. Reference Magpantay, Riolo, Domenech de Cellès, King and Rohani2014). A quantity called the vaccine impact, denoted by ϕ, has been introduced as a measure of the effectiveness of imperfect vaccines (McLean & Blower, Reference McLean and Blower1993). ϕ is related to the basic reproduction number Rp in the presence of constant vaccination at coverage level of p:
$${R_p}\, = \,{R_0}({\rm 1} - \varphi p).$$
Here, R 0 is the basic reproduction number of the disease in the absence of vaccination. In Section S10, using the same reasoning in Magpantay et al. (Reference Magpantay, Riolo, Domenech de Cellès, King and Rohani2014), we derived the vaccine impact of our system to be
$$\varphi = \,({\rm 1} - {\varepsilon _A})[{\rm} ({\rm 1} - \theta ){\rm} + \theta ({\rm 1} - {\varepsilon _L})({\rm 1} - {\varepsilon _W})].$$
Here ε W is the probability that vaccine-derived immunity will wane within an individual's lifetime. It can be derived from the waning rate as follows. At the MLE of the Waning Model, the waning rate α = 0·10 year −1. Assuming an (exponentially distributed) lifespan of mean 75 years, we have
$$\eqalign{{\varepsilon _W} & = \displaystyle{{{\rm waning}\,{\rm rate}} \over {{\rm exit}\,{\rm rate}\,{\rm from}\,V\,{\rm compartment}}} \cr & \approx \displaystyle{{0.10} \over {0.10 + 1/75}} = 0.08.}$$
Hence we derive that the vaccine impact for the best Waning Model is 0·12; similarly, one computates that, for the best Leaky Model, ϕ = 0·82. The values of the vaccine impact for the Restricted Models at their MLEs are also presented in Table 7.
Vaccine parameters are calculated using (2)–(4). All values are based on the stochastic forms of the Restricted Models (where θ = 1 and η = 0 is fixed)
From the relation (2), we have ϕ = (1/p)(R 0−R p )/R 0. Thus vaccine impact measures the relative reduction in transmission within a population. This is a reflection of herd immunity, since we can derive an inverse relationship between vaccine impact and the vaccination coverage necessary for disease eradication. Furthermore, with θ = 1, as in the case of our Restricted Models, the vaccine impact measures the protection a vaccinated individual receives relative to an unvaccinated individual (Magpantay et al. Reference Magpantay, Riolo, Domenech de Cellès, King and Rohani2014).
Comparison of the best Leaky- and Waning-Models
The MLEs for the waning rate in the Waning Model and for leakiness in the Leaky Model both indicate a reduction in disease transmission. The two models were comparable in their capacity to explain the 1996–2009 Italian data. Nevertheless, they make very different predictions about the impact of the aP vaccine on both individual-level protection from infection and overall transmission. In particular, while the vaccine impact values for both the Leaky- and Waning-Models point to some vaccine protection, the Leaky Model implies a much higher level of protection for the vaccinated individual (ϕ = 0·82), relative to that estimated under the Waning Model (ϕ = 0·12). This contrast reveals that, though our analysis sheds some light on pertussis epidemiology, considerable uncertainty in regard to key parameters remains.
Simulations of epidemiological futures under the two contrasting models show that the alternative scenarios they represent are actually quite distinct at the population scale as well. For example, Fig. 5 shows simulations of the stochastic Leaky- and Waning-Models at the MLE estimated for Lazio. While no clear difference between the simulated pertussis incidence reports under the two models if vaccination is maintained at 98% (Fig. 5A and B), if vaccination coverage is allowed to gradually drop from 98% in 2010 to 80% in year 2070, the difference between the two models becomes quite apparent (Fig. 5C and D). In particular, under the Leaky Model, typical simulations yield future epidemics at intervals of 4–5 years, while under the Waning Model, epidemics occur about every 2 years. These simulations illustrate the different dynamics under these two, equally plausible models embodying distinct modes of vaccine protection. An implication is that, by examining data from countries that have experienced different kinds of dynamics, one might be able to better establish the values of ε L and ε W . That doing so should be a high priority has been recently demonstrated, for example, by work establishing that control efforts optimal for leaky vaccines can differ greatly from those tailored to waning vaccines (Riolo & Rohani, Reference Riolo and Rohani2015).
Sample simulations of the Leaky- and Waning-Models at the MLEs for Lazio. Abbreviation: MLEs, maximum likelihood estimates.
Sensitivity of results to assumptions on primary vaccine failure
To this point, we have been discussing models under the assumption of zero primary vaccine failure (ε A = 0). We now consider the effects higher values of primary vaccine failure have on our Restricted Model parameter estimates. We re-fitted the stochastic forms of the Restricted Models under the assumptions of 7.5% (ε A = 0·075) and 15% primary vaccine failure (ε A = 0·15). The new profiles over the leakiness for the Leaky Model and the waning rate for the Waning Model are shown in Fig. 6, along with the original profiles corresponding to ε A = 0.
Plots of the log-likelihood values (profiles) corresponding to the MLEs of the stochastic forms of (A) the Leaky Model at each fixed value of the leakiness, and (B) the Waning Model at each fixed value of the waning rate at different levels of primary vaccine failure. Since all the models here have the same numbers of parameters, the change in AIC c is the same as two times differences in likelihood. Abbreviations: AIC, Akaike Information Criterion; MLEs, maximum likelihood estimates.
One immediately observes that either assumption of non-zero primary vaccine failure results in lower likelihoods (and therefore worse AIC c ) relative to the original assumption. The very low likelihoods under the 15% assumption indicate that the data are unequivocal that this level of primary failure or higher is incompatible with the data. Increasing ε A also changes the estimates of the leakiness, ε L and waning rate, α. As shown in Table 7, the values of both ε L (under the Leaky Model) and α (under the Waning Model) shift to much lower values when primary vaccine failure is a possibility. These shifts inform intuition as to how the model is explaining the data. At high values of primary vaccine failure, we expect there to be more fully susceptible (S 1) individuals in the population. To maintain the observed level of incidence, therefore, transmission must be reduced. The Leaky Model accomplishes this by reducing the size of the S 2 class via decreased vaccine leakiness. Likewise, the Waning Model achieves this effect by decreasing the rate of waning.
Breaking down the results by region, we observed that varying the probability of primary vaccine failure affected the consistency of the patterns of parameter estimates. In particular, with ε A = 0, the MLEs for Sicilia indicated low transmission and high reporting probabilities, in contrast to results for the other regions and inconsistent with our understanding of pertussis as a highly infectious childhood disease (refer to Section S4). At ε A = 0·075, we found three regions (Lombardia, Sicilia and Toscana) for which MLEs indicated relatively low transmission values, and two regions (Sicilia and Toscana) with relatively high reporting probability. At ε A = 0·15, four regions (Lombardia, Sicilia, Toscana and Umbria) had MLEs indicating low transmission and two (Sicilia and Toscana) with high reporting probabilities. Thus, under the assumption of zero primary vaccine failure, the models not only achieve higher likelihoods, but tell a more consistent tale about model parameters than do those with appreciable primary vaccine failure.
Measures of transmission at the steady state levels of the different models are presented in Table 8, as calculated for Lazio region. We observe that the predicted mean age at first infection increases with increasing primary vaccine failure, due to changes in transmission rate in Lazio, which to explain the data must drop with increasing rates of primary vaccine failure. Though the assumption of zero primary failure leads to better AIC c values, the observed mean age of first infection for immunologically naive individuals agrees most closely with the estimates under the 7·5% primary vaccine failure assumption. In view of this discrepancy, we emphasize the need for caution in interpreting the results of the AIC c comparison: it would be a mistake to view such a comparison, based on these data alone, as decisive. In particular, it is likely that age-specific incidence data, affording a more highly resolved view of the case age distribution than is available to us at present, would lead to different conclusions regarding the rate of primary vaccine failure, the degree of leakiness and the rate of waning.
Measures of transmission for Lazio time-averaged over the years 2040–2070 if the vaccine coverage as well as demographic growth rates are assumed to be maintained at the final levels recorded in 2012
All values are based on the stochastic forms of the Restricted Models (where θ = 1 and η = 0 is fixed). Details on the calculation of the mean values and confidence intervals are outlined in Section S11.
The big picture
There have been few studies in which mechanistic stochastic models have been fitted to pertussis notification data (Blackwood et al. Reference Blackwood, Cummings, Broutin, Iamsirithaworn and Rohani2013; Lavine et al. Reference Lavine, King, Andreasen and Bjørnstad2013). These studies adduce evidence for considerable infection-derived and wP-induced herd immunity. Building upon this body of work, we constructed a general model of pertussis considering different aspects of vaccine efficacy and confronted it with incidence time series from the early phase of the aP vaccination effort in Italy. While we were able to estimate some parameters with reasonable precision, further work is needed to verify the generality of our results and to eliminate remaining uncertainties. These uncertainties persist in large part due to the unavailability of region-specific vaccine coverage data during the early phase of aP vaccination, because we lack information on the age distribution of cases, and because the period examined is brief relative to the duration of transients potentially excited by the ramp-up (Magpantay et al. Reference Magpantay, Riolo, Domenech de Cellès, King and Rohani2014).
We note that the MLE of the Waning Model at zero primary vaccine failure is inconsistent with some previous work supporting a much longer duration of aP-induced immunity (Rohani et al. Reference Rohani, Zhong and King2010; Gambhir et al. Reference Gambhir, Clark, Cauchemez, Tartof, Swerdlow and Ferguson2015). In contrast, the very different and equally well-supported Leaky Model at zero primary vaccine failure is consistent with the very considerable reduction in transmission observed in other studies of pertussis (Broutin et al. Reference Broutin, Viboud, Grenfell, Miller and Rohani2010; Domenech de Cellès et al. Reference Domenech de Cellès, Riolo, Magpantay, Rohani and King2014). The equivocation of the data on the question of whether the Leaky or Waning Model is better points to the difficulty of disentangling the different modes of vaccine failure and to the uncertainties that consequently remain with respect to aP vaccine impact on transmission. In particular, although vaccine impact was high (>70%) under all but one of the scenarios we examined, under the Waning Model at zero primary vaccine failure, it was only 12%.
Concluding remarks
Different modes of vaccine failure were built into mechanistic models of pertussis that also incorporated demographic and vaccine coverage data. To assess the contributions of these different modes, model parameters were estimated by fitting to regional pertussis notification data from 1996 to 2009, during the transition from low to high aP-vaccine coverage in Italy. Both deterministic and stochastic forms of the models were used. The stochastic forms always yielded better likelihood and AIC c values over the deterministic versions, indicating that stochastic models are better at explaining the data during this transition period and further highlighting the pitfalls of attempting to explain non-equilibrium stochastic dynamics with deterministic models (King et al. Reference King, Domenech de Cellès, Magpantay and Rohani2015a ).
The evidence described here favours models wherein infections in vaccinated individuals are unlikely to be reported – perhaps because of decreased disease severity – but are as transmissible as those in unvaccinated individuals. Estimates of the protection against infection provided by the vaccine yield a spectrum of equally viable alternative descriptions of the vaccine ranging from (1) a leaky vaccine with no primary vaccine failure and has the effect of reducing the probability of infection upon exposure to approximately 0·18 (0·14, 0·25) of what it would be for unvaccinated individuals, to (2) a waning vaccine with no primary vaccine failure and a mean protection of about 10 (5·9, 12·5) years. The range of models supported by the 1996–2009 data from Italy yield differ substantially in their predictions as to the impact of vaccination on reducing susceptibility of individuals to infection and on diminishing transmission at the population level. Moreover, they give very different predictions as to the nature of a pertussis resurgence should vaccine coverage fall below present high levels.
We varied our assumptions regarding the probability of primary vaccine failure, finding little evidence in these data for appreciable primary vaccine failure and noting that more highly resolved age-specific data might well improve estimates of this quantity. Our conclusions were limited by incomplete vaccine coverage data, lack of age stratification in the data and the brevity of the time series: future work, using higher-resolution data from other countries and time-periods, will be needed to test the generality of our conclusions and further reduce the uncertainty on key parameters.
In concluding, we note that though these data leave some questions unanswered, information from other sources contain some clues. In particular, we note that the estimated 10 years of immunity we obtained under the assumption of a perfectly protective but waning vaccine is in apparent contradiction with the pronounced aP-induced reductions in transmission seen in other studies (Carlsson & Trollfors, Reference Carlsson and Trollfors2009; Rohani et al. Reference Rohani, Zhong and King2010; Domenech de Cellès et al. Reference Domenech de Cellès, Riolo, Magpantay, Rohani and King2014). At the other extreme of the spectrum, our model with leaky but permanent vaccine-induced immunity predicts a considerable reduction in transmission (82%), in line with studies of pertussis in other regions and times. Clearly, further examination of the potential for aP to induce durable but leaky protection is warranted.
SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit http://dx.doi.org/S0031182015000979. Data and code are available at http://dx.doi.org/10.5061/dryad.58q00.
ACKNOWLEDGEMENTS
We thank the members of the Rohani and King Labs for their helpful comments on the implementation of this project and the interpretation of results and the two anonymous reviewers for their careful reading and useful comments.
FINANCIAL SUPPORT
This work is supported by the Research and Policy in Infectious Disease Dynamics programme (Science and Technology Directorate, Department of Homeland Security and the Fogarty International Center, National Institutes of Health) and by a research grant from the National Institutes of Health (grant number 1R01AI101155) and by MIDAS, National Institute of General Medical Sciences (grant number U54-GM111274).
