Abstract

There has also been a shift over the same time period in the age group reporting the largest number of cases (aside from infants), from adolescents to 7—11 year olds. We use epidemiological modelling and a large case incidence dataset to explain the upsurge. We investigate several hypotheses for the upsurge in pertussis cases by fitting a suite of dynamic epidemiological models to incidence data from the National Notifiable Disease Surveillance System (NNDSS) between 1990—2009, as well as incidence data from a variety of sources from 1950—1989. We find that: 1) the best-fitting model is one in which vaccine efficacy and duration of protection of the acellular pertussis (aP) vaccine is lower than that of the whole-cell (wP) vaccine, (efficacy of the first three doses 80% [95% Cl: 78%, 82%] versus 90% [95% Cl: 87%, 94%]), 2) increasing the rate at which disease is reported to NNDSS is not sufficient to explain the upsurge and 3) 2010—2012 disease incidence is predicted well. In this study, we use all available U.S. surveillance data to: fit a set of mathematical models and determine which best explains these data and determine the epidemiological and vaccine-related parameter values of this model. We find evidence of a difference in efficacy and duration of protection between the two vaccine types, wP and aP (aP efficacy and duration lower than wP). Future refinement of the model presented here will allow for an exploration of alternative vaccination strategies such as different age-spacings, further booster doses, and cocooning.

Author Summary

There has also been a shift over the same time period in the age group reporting the largest number of cases (aside from infants), from adolescents to 7—11 year olds. We investigate several hypotheses for the upsurge in pertussis cases by fitting a suite of epidemiological models to incidence data from the National Notifiable Disease Surveillance System (NNDSS) between 1990—2009. We find that: 1) the best-fitting model is one in which the vaccine efficacy and duration of protection of the acellular pertussis vaccine is lower than that of the whole-cell vaccine, 2) increasing the rate at which disease is reported to NNDSS is not sufficient to explain the upsurge and 3) 2010—2012 disease incidence is predicted well. These results demonstrate that the resurgence in pertussis in the US. can be explained by past changes in vaccination policy. However, our findings suggest that the efficacy of the currently-used acellular vaccine is not much lower than that of the whole-cell vaccine, and booster doses may be sufficient to curtail epidemics while vaccine research continues.

Introduction

This upsurge occurred in the context of steadily rising reported disease from the early 1980s [4], despite the maintenance of high vaccination coverage levels (>90%) [5]. Pertussis remains a major cause of childhood mortality worldwide, responsible for 195,000 deaths in 2008 [6]. A vaccine for the disease was developed in 1942 in the U.S., and was included and administered as a killed whole-cell component of the diphtheria/tetanus/pertussis (DTP) combination vaccine.

Due to the possible reactogenicity of the whole-cell pertussis component of the DTP vaccine, acellular pertussis (aP) vaccine was developed as a replacement in 1991 (DTP became DTaP). This new vaccine was administered to children in 1992 [7,8] and then phased into the infant immunization schedule beginning in 1997 [9,10] , following approval by the U.S. Food and Drug Administration (FDA) and recommendations by the U.S. Advisory Committee on Immunization Practices (ACIP).

Even accounting for this steady upward trend, which some have attributed to improved surveillance and diagnostics [11,12], case numbers during 2004—2012 were notably high and many state public health departments reported outbreaks during these years. During the two most recent of these outbreaks, the worst-hit states were California (in 2010) [2,3] and Washington (in 2012) [13] with case counts unmatched since the mid-1940s. While the 2004—2005 outbreak was characterized by an increase in the number of cases among adolescents (10—20 year olds), as well as a less-pronounced increase across all ages, the most recent outbreaks have seen many more cases among 7—10 year olds and young adolescents (11—13 year olds) along with an overall increase in case numbers [2]. A number of explanations have been put forward for the increase in disease that has been observed, particularly for the outbreaks of the past 10 years [11,14—19]. These include: 1) The evolution of the circulating bacteria away from the targets of vaccine antigens; 2) A decline in vaccine coverage levels; 3) a change in vaccine efficacy (VE), and/or duration of protection due

Incidence of disease and vaccination coverage in the United States. The time varying incidence of disease in the US, summed over all ages, between 1950 and 2009, annotated for significant events relating to vaccination policy. The blue solid line illustrates the DTP3 (first three DTP doses) vaccine coverage level over time. to the shift to the acellular vaccine; or a specific lot of manufactured vaccine With lower efficacy; 4) The decline of natural boosting, through reduced exposure to naturally circulating pertussis bacteria, as a result of years of successful mass-vaccination; 5) An increase in disease reporting rates (e.g. due to improved diagnostics and surveillance) [20,21].

We fit the output of this model to 20 years of age-stratified incidence data from the National Notifiable Diseases Surveillance System (NNDSS). The model includes both vaccine and naturally-induced immunity (with both of these sources of protection waning over time), the current vaccination schedule administered to the population, an estimate of the underreporting of disease, and a decreased infectiousness of those Who have experienced greater than one infection.

Results

The most parsimonious (best fitting) model, Model 8, has a difference in efficacy and duration of protection between the whole-cell and acellular vaccines, and a difference between the duration of protection of whole-cell vaccine and natural infection (Table 1). Model 5, in which the duration of protection for whole-cell vaccine and natural infection is the same, has a slightly lower DIC value, but is very close to Model 8, and these two models are almost interchangeable in terms of their parameter values. Further details of the candidate models are provided in the Supporting Information. The model reproduces the trends in the NNDSS incidence data up to 2009; running the model forward over another three years shows a continued correspondence with observed incidence trends (Fig 2). A large proportion of forward simulations (>50%) show an upsurge in incidence around 1990, a phenomenon that is not recorded in the data. However, the mean of the model-predicted incidence remains much closer to the data, albeit with a large range of uncertainty. Table 1. Descriptions of the nested models that were fitted to the NNDSS incidence data. Model 1 Description Protection duration of whole cell vaccine same as natural infection; acellular vaccine same as whole-cell Protection duration of whole cell vaccine same as natural infection; different efficacy for acellular vaccine Protection duration of whole cell vaccine same as natural infection; different protection duration for acellular vaccine; Protection duration of whole cell vaccine different from natural infection; acellular vaccine same as whole-cell Protection duration of whole cell vaccine same as natural infection; protection duration and efficacy different for acellular vaccine Whole cell vaccine protection duration different from natural infection; different efficacy for acellular vaccine Whole cell vaccine protection duration different from natural infection; different protection duration for acellular vaccine Whole cell vaccine protection duration different from natural infection; protection duration and efficacy different for acellular vaccine The mean posterior values of the Deviance Information Criterion (DIC) of the models are given in the rightmost column. Fig 3 illustrates the age distribution of incident cases every two years from 1994—2012, for Model 8. Several features of the NNDSS data are reproduced by the model: 1) an overall rise in disease cases throughout the observed period; 2) an adolescent peak in the incidence curve for most of the years post-1995; 3) a gradual rise in cases among 5—10 year olds from 2006, forming a peak which is comparable in size with the adolescent peak by 2008.

This smaller peak moves forward one year at a time, growing in magnitude until 2006 when the two peaks are of similar size. By 2008, the cohort of individuals

Year

Incidence of disease in the United States, compared with modeled values. The time varying incidence of disease cases in the model and the US data (log scale) after 1990. The black dots are US disease incidence data, and the shaded regions represent the credible intervals (50% and 95%) obtained through model parameter estimation of model 8. The model has been run beyond the time over which it was trained to illustrate its continued correspondence with the 2010—2012 data (red crosses). Table 2. Parameter estimates for the best-fitting model, Model 8 (models outlined in Table 1). Parameter description Value Vaccine efficacies & waning

50 yrs protection Vaccine efficacy As acellular Epidemiological Parameters (with reference to naTve individuals) infection (with reference to primary-infected individuals) Year of reporting rate change None Mean reporting rate after change n/a trends in the data. The reporting rate is estimated at 6% [95% CI: 0.1%, 22%] with the wide credible intervals suggesting that a gradual change in reporting over time may not be precluded, even if a step change is not supported. A gradual change in reporting is consistent with the uptake of more sensitive diagnostic methods as well as greater awareness of pertussis in recent years.

We demonstrate that, in successive ages beyond when the final vaccine dose is administered, there is a decline in the 95% credibility envelope of the model-based VE estimates (gray-shaded region): the bottom of this envelope declines more rapidly than measured for the California study. The curve lying above the shaded region illustrates the results of a hypothetical VE study if it had been conducted in 1990, the pre-acellular era. It can be seen that VE declines more slowly in this case.

Discussion

These differences caused a shift in the age distribution of pertussis disease incidence: the adolescent peak in years prior to 2006 shifted to a younger age group (5—10 year olds) post-2006. This explanation for pertussis disease dynamics is in broad agreement with case-control and retrospective cohort studies showing that vaccine-in-duced immunity wanes faster than previously thought [22—24]. The average duration of protection we estimate for acellular vaccine is shorter than the essentially lifelong protection we estimate for both whole-cell vaccination and natural infection; however our estimate of the acellular protection duration is around 50 years, which is longer than might be expected from published VE studies [22,23]. Since our calculated waning rate is

Years since 5th vaccine dose

Case-control study results compared with modeled values. Vaccine Effectiveness (VE) as measured in the case-control study of Misegades et al [22] from 2010 in California (black crosses) compared with VE values generated by simulations of the case-control study using the model fitted to the incidence and these VE data (gray region). The gray shaded region represents the 95% credible interval of the model outputs. The dotted curve lying above the 2010 data and simulations was calculated by simulating a 1990 case-control study, and shows significantly slower waning of the VE value.

By simulating a prospective VE study, we find that even a small increased flow into the susceptible pool is sufficient to result in a strong signal of declining VE following the final vaccine dose (see Figs 4 and 5 and [22]). The efficacy values implied here for the full acellular vaccine regimen (i.e. 5 doses) are consistent with those estimated in the acellu-lar pertussis vaccine (APERT) trial [25] , lending further weight to the correspondence of our model with the broad range of data. Furthermore, VE values for individuals with mixed acellu-lar and whole-cell vaccine histories are consistent with a recent Australian report of disease rates and linked vaccination histories [26].

Our estimate of case reporting rates (approximately 6%) is in line with a recent multi-country analysis [20], but has wide 95% credible intervals. Hence, progressive changes in disease reporting, such as those expected by more sensitive diagnostic testing (e.g. PCR) or a greater awareness of pertussis among patients and physicians, is within the scope of our model. Further work with this model could investigate published reporting trends more Closely [l l,20,21,27,28] .

First, we estimate R0 to be in the range of 9—12 which is closer to the values found by previous mathematical models [fi,fl] than the often-quoted range of 12—17 [E]. Such values do not rule out elimination at high enough vaccination coverage levels. Second, previous recipients of vaccine whose protection has waned are estimated to be 32% as susceptible to infection as in-fection-na’ive individuals and, if infected, these individuals are estimated to be 17% as infectious as primary infections in immunologically naive individuals. These numbers suggest that adults and adolescents may be an important reservoir of infection. Third, the most parsimonious model with a good fit to the data suggests that the duration of immunity generated by the whole-cell vaccine is nearly equivalent to that induced by natural infection. Loss of immunity 6

Recovery from 'nfectio

Force of infection

Loss of Force of

y vaccine reinfectlon

Epidemiological model diagram. The compartmental model for pertussis infection and disease used, modified from the basic model of Aguas et a/ [44]. There are two infected compartments for primary (/1) and secondary (or higher) infection (/2). We assume that surveillance systems only capture those who are experiencing primary infection. Once an individual has recovered from primary or secondary infection, induced immunity may wane, whereupon they will reenter the susceptible state. Rates of flow between compartments are defined in Table 3.

We find a 'hypersensitive' boosting model [14] to be too unstable to mimic the steady post 1970s rise in disease; the incidence of disease for that model shows a large amplitude oscillatory behaviour. Population age-dependent contact patterns alone also appear to be insufficient to capture recent changes in U.S. incidence of disease [15]. Declines in vaccination coverage and the evolution of the circulating bacteria would lead to rising disease but not the changing age-distribution, although the recent discovery of pertactin-negative (or other genetic) variants may be playing a role in the vaccine efficacy decline we have found [32,33]. Changes in age-dependent disease reporting could mimic observed incidence patterns, of course, but such models would lack the simplicity of ours.

Further data collection will help to narrow uncertainties in the model, for example, the range of R0 and, the infectiousness of those who have been previously infected, both relevant to control strategies. Alternative model structures capturing the development and waning of immunity and boosting of immunity may be contributing to discrepancies between the outputs of our best-fitting model and the available data, such as an overestimate of the incidence peak in adolescents after 2008 (Fig 3), and signs of a model underestimate of the overall incidence (Fig 2). Future work should attempt to incorporate the full range of possibilities for modeling boosting and immunity e. g. using the Table 3. State variables (population compartments) and model parameters for the model investigated. Symbol Description Value Population compartments (state variables) 8 Susceptible [1 Primary infection [2 Secondary (or greater) infection R Recovered from infection (whether primary or secondary) V Vaccinated D Reported cases of pertussis Model parameters A(t) Force of infection at time t (this subsumes the model Varies in value depending upon number parameter [3 (known as the transmission parameter) of individuals infected (Per capita per which is estimated by the model-fitting procedure) year) T Rate of recovery from primary or secondary infection 24 per year (i.e. duration of infection is A Rate of loss of immunity following natural infection Estimated here (see Table 2) F Rate of loss of immunity following vaccination Estimated here (see Table 2) 2 Relative rate of acquisition of secondary to primary Estimated here infection H Relative infectiousness of secondary infections Estimated here p1, p2 Rate at which infected individuals report disease Estimated here proposed model of Lavine et al. [14], as well as mixing matrix data that are as locally-specific as possible, following the work of Rohani et al. and Riolo et al. [15,34]. Clearly, the use of a mixing matrix from Great Britain is not ideal for the U.S. but was all that was available at the time of writing. In addition, fitting our model to data from the U.S. aggregates the incidence from a large number of more spatially localized epidemics, which may not be completely synchronized; this means that, while our model is convincing, it tends to describe phenomena such as the ‘cohort ef-fect’ and peaks in incidence more sharply than the data. Future work should focus on modeling local outbreaks and endemic pertussis with high spatial resolution, following Rohani et al. and Choisy et al. [35,36]. More intricate modeling frameworks, such as individual-based modeling, are also promising directions for future work, given their ability to represent human behavioural dynamics, which could be a major local contributor to vaccine uptake or refusal [37].

Jackson et al. recently reviewed such data and showed that global pertussis epidemiology is complicated and by no means follows the same patterns and explanations as the U.S. [38]. However, the authors of the review do note that the issue of the switch from the wP to the aP vaccine may be sufficiently specific to allow investigation of its impact upon global trends.

Pertussis epidemiology is complex and disease incidence is non-linearly related to vaccine coverage. With doubt around the efficacy and duration of protection of the acellular vaccine, modeling is an essential tool to help us better understand the changing epidemiology of pertussis.

Methods

Infected individuals can be either primary-infected or infected more than once (see Fig 5, Table 3 and ‘Model structure details’ below). The model also accounts for the annual aging of the population: the time-dependent birthrates and deathrates are used for the U.S. (obtained from [42]). The birth and deathrates prior to the analysis period of the study were fixed at the earliest values available, namely those of 1950. Mixing between age-groups is incorporated according to a European diary-based contact study (used in the absence of a suitable U.S. alternative [43]).

In both of these compartments, individuals have developed a protective immune response which decreases the probability of infection but eventually wanes. This vaccinated proportion may be thought of as coverage X vaccine e icacy and is equal to the protected proportion of the population with each dose. Annual vaccination coverage levels were obtained from the National Immunization Survey and its predecessor surveys [45].

We model the first three doses as a single removal of susceptible infants six months following birth (as an approximation of the building immunity due to the first three doses), the 18 month dose as two separate removals, each with one half of the overall coverage, at one and two years of age and, similarly, the final dose as two separate removals at four and five years of age. The 18 month and four to six year doses are represented as double doses to capture a fairly wide range of times at which vaccine is administered to children across the whole population. The whole-cell vaccine is first applied to the model population in the simulation year 1947 and a switch is made to the acellular vaccine for the final two doses in 1992 and then to the entire schedule in 1998. A Tdap booster vaccine is introduced to those aged 12 years in the model, in 2005, with coverage levels given by recent surveys [46—51].

This rate was permitted to change once in the time period modelled, to partially account for (i) a possible change in the sensitivity with which cases are found over time, and (ii) increased reporting associated with greater awareness of pertussis among patients and physicians. We also allow for some level of spatial and temporal heterogeneity in the reporting rate by modelling p as drawn from a Beta probability distribution (see Statistical details); this allows for a single average value of the reporting rate to be estimated for each of the two time periods, but for the uncertainty in this value to be quantified. Only those in the primary-infected compartment (I 1, Fig 5) were counted as disease cases; this is a simplification, but one that is justified by recent studies which show that second-ary-infected or vaccinated nonhuman primates had more limited disease when reinfected [52].

The fitting was performed using Markov Chain Monte Carlo (MCMC) methods, which find an ensemble of model parameters that fit the data and allow a credible interval for each of these parameters to be determined. The chains were run to 100,000 steps during equilibration and 50,000 to sample the posterior parameter distributions. Sensitivity of the parameter estimates was explored by performing MCMC separately using different subsets of the data.

These model components included possible changes in (i) duration of protection between whole-cell vaccination and natural infection; (ii) vaccine efficacy and (iii) duration of protection between the whole-cell and acellular vaccine. The candidate models were compared to determine the simplest and most biologically compelling explanation of the observed variations using the Deviance Information Criterion (DIC), a measure that accounts for both good-ness-of-f1t as well as model complexity [53] (the higher the value, the better the fit to the data).

In order to further constrain our parameter selections and values, we used our models to simulate this study and fitted the resulting VE values to those found by Misegades at each of the corresponding times. Misegades found that the VE decreased from 98% to 82% in 4 years following the final of five doses. As well as measuring our model-based VE over the same period, we used the same method to estimate VE during a period when the course of doses consisted entirely of the whole-cell vaccine; we selected 1990 for this purpose. Model structure details

The mixing matrix fi(i, j) represents the product of the contact rate (obtained from the ‘POLYMOD’ diary study of contact patterns in Great Britain (GB) [43]) and the transmission probability per contact between individuals in age groups 1' and j. Individuals in state 121- have an infectiousness value of 17 compared with those in state 11,-, and the relative susceptibility to infection of those in Rl- (the ‘recovered’ state) compared With those completely susceptible. Vaccination occurs with the DTP vaccine With Whole-cell pertussis component by shifting individuals into compartment Vll- and into V2,- for the DTaP and Tdap vaccines (the main text simplifies this flow into a single compartment V,- so as not to clutter the exposition). The equations corresponding to the model flow diagram (Fig 5, Table 3) are given below:

Statistical details

Pertussis case count data from NNDSS were aggregated into annual counts for each age group so that yl- (t) is the number of pertussis disease cases for age group i in year t. Our mathematical model outputs were also aggregated into annual counts for each age group i, so that x,- (t) was the model-derived case count for age group 1', during year t. These model-derived case counts are functions of the model structure and parameters, so that they might be better expressed as x,- (t|0, M), where 0 represents the parameter vector for model M. Since the pertussis cases counted by NNDSS may be a proportion of the true number (and the output of the natural history epidemiological model is designed to simulate the true number) we can say that the observed number of cases follows an observation model so that annual counts of pertussis cases (0195,- (t)) are governed by a binomial distribution:

The binomial log-likelihood of the data given these considerations is: Log Likelihood =

However, we find that it can be excessively difficult to explore the parameter space productively (i.e. moving from regions of poorer to relatively better fits to the data) since the binomial likelihood surface can take on high (i.e. better fitting) values in very small (and therefore very hard to locate) regions of parameter space. To attempt to remedy this problem (as well as to account for a more general degree of uncertainty in the disease reporting rate), we introduce into our observation model a further layer of uncertainty such that the proportion of individuals reporting disease is drawn from a Beta distribution, which has two parameters, Where the Beta distribution’s parameter values are time-dependent in the same way as we posited above i.e. : Using this structure for the reporting rate component of the observation model, we can recalculate the likelihood of a particular dataset given the model (removing the subscripts for simplicity of presentation): This integral allows us to compute an average likelihood allowing for the entire distribution of the reporting rate p, given its Beta distribution. Because the Beta and Binomial distributions are conjugate, the integral is straightforward. Where B(o<, fl) is the normalizing constant for the Beta probability distribution. The integrand can be seen to be the probability density function of the Beta distribution Beta(o</, fl, ) with the adjusted (or ‘renormalized’) parameter values, so that: So the integral above is the normalization constant of the Beta distribution (which is referred to as the Beta function, itself comprised of a set of gamma functions, see below) With the adjusted parameters. The likelihood therefore becomes: And for all of the data points, the log-likelihood becomes:

Supporting Information

81 Text includes a brief description of alternative model structures that have been investigated in the recent literature; a table outlining precisely which parameters are allowed to alter and which are fixed for each of the models we investigate; a series of plots outlining the best-fitting model’s fit to the age-dependent incidence of disease cases for each of the years between 1990 and 2009; and a model-based simulation of the case-control study outlined in the main article, showing the gradual effect of an aging cohort of individuals who have been entirely vaccinated with the DTaP rather than the DTP vaccine. (DOCX)

Acknowledgments

The authors wish to acknowledge Lara Misegades, Amanda Faulkner, Stacey Martin and Garrett Asay from CDC for invaluable and insightful discussions on the work and the manuscript.

Author Contributions

Performed the experiments: MG SC NMF. Analyzed the data: MG SC NMF. Contributed reagents/materials/analysis tools: MG SC NMF. Wrote the paper: MG TAC SC SYT DLS NMF.