Interaction of Vaccination and Reduction of Antibiotic Use Drives Unexpected Increase of Pneumococcal Meningitis

Antibiotic-use policies may affect pneumococcal conjugate-vaccine effectiveness. The reported increase of pneumococcal meningitis from 2001 to 2009 in France, where a national campaign to reduce antibiotic use was implemented in parallel to the introduction of the 7-valent conjugate vaccine, provides unique data to assess these effects. We constructed a mechanistic pneumococcal transmission model and used likelihood to assess the ability of competing hypotheses to explain that increase. We find that a model integrating a fitness cost of penicillin resistance successfully explains the overall and age-stratified pattern of serotype replacement. By simulating counterfactual scenarios of public health interventions in France, we propose that this fitness cost caused a gradual and pernicious interaction between the two interventions by increasing the spread of nonvaccine, penicillin-susceptible strains. More generally, our results indicate that reductions of antibiotic use may counteract the benefits of conjugate vaccines introduced into countries with low vaccine-serotype coverages and high-resistance frequencies. Our findings highlight the key role of antibiotic use in vaccine-induced serotype replacement and suggest the need for more integrated approaches to control pneumococcal infections.

covariate in the model. Because carriage is typically asymptomatic, we supposed that all individuals were equally exposed to antibiotics, irrespective of their carriage status. We considered the following classes of antibiotics: beta-lactams (penicillins and cephalosporins) and macrolides, which represent 75-80% of antibiotic prescriptions in France 12 . Although we did not consider class-specific antibiotic exposure, we sought to reproduce those differences by computing a weighted average of each class's effect. Such a calculation might introduce bias if the relative frequency of each antibiotic class varies over time; as Fig. S2 shows, however, these relative frequencies varied little in France over 2001-2009, despite an overall reduction of antibiotic use 12 . The relative frequencies of each antibiotic class were 44% for penicillins, 28% for cephalosporins and 28% for macrolides. There are distinct mechanisms by which antibiotics can select for pneumococcal resistance in the population 13 . Here, we considered two mechanisms by which antibiotics favor pneumococcal resistance. First, antibiotics were supposed to enhance the clearance of carried pneumococci. This effect was modeled as an additional clearance rate under antibiotic treatment, (1 − ! ) for penicillin-susceptible strains and (1 − ! ) for penicillin-resistant strains, where ω represents the rate of antibiotic action, and σ S (respectively σ R ) the probability of non-decolonization of penicillin-susceptible (resp. penicillin-resistant) under antibiotic treatment. The two parameters σ S and σ R were first fixed for each antibiotic class, as follows. Penicillins are highly effective at clearing penicillin-susceptible pneumococci, and also partly effective at clearing penicillin-resistant pneumococci 14,15 . Therefore, we assumed that penicillin use cleared all susceptible pneumococci (that is, ! ! = 0) and 50% of resistant pneumococci ( ! ! = 0. 5). In contrast, cephalosporins and macrolides are less effective at clearing penicillin-susceptible pneumococci, and have little effect, if any, on penicillin-resistant pneumococci 14,15 ; therefore, we fixed ! ! = 0.5 and ! ! = 1 for cephalosporin use, and ! ! = 0.25 and ! ! = 1 for macrolides.
Second, we assumed that, under antibiotics, noncarriers had a lower risk of acquiring a susceptible strain (relative risk ϕ S <1) and a higher risk of acquiring a resistant strain (relative risk ϕ R >1). Again, these two parameters were fixed for each antibiotic class. Penicillins are narrowspectrum antibiotics that barely affect the nasopharyngeal commensal flora; it has been argued that penicillins do not increase the relative risk of acquisition of resistant strains 16 . Therefore, we fixed ! ! = 0.1 and ! ! = 1 for penicillin use. Cephalosporins and macrolides, however, impact the resistant ecology more markedly 17 ; therefore, for these two classes, we fixed 18 The final values (parameters without superscripts) were calculated as an average of the class-specific parameters, weighted by the relative frequency of the different antibiotic classes.
The fixed values for each antibiotic class and the final values used in the model are summarized in Model equations. We adopted the following notations. The population was divided into noncarriers S and carriers C, who were subdivided into four different categories, depending on vaccination status and antibiotic-exposure status: unexposed and unvaccinated (superscript UNV), exposed and unvaccinated (superscript ENV), unexposed and vaccinated (superscript UV), and exposed and vaccinated (superscript EV The total numbers of carriers for each type were: The strain-specific forces of colonization, λ X , were given by: where N is the population size, and β X the strain-specific transmission rates, assumed to be constant over time. The strain-specific infection rates were modeled as: where ρ X is the mean infection rate, the amplitude of seasonality, t the time (in years) and ϕ ρ the phase, which controls the time-location peak. This formulation led to the following deterministic differential equations; in the simulations, we used a stochastic analogue of these equations with Gillespie's τ-leap algorithm. Table S4 presents the complete list of model parameters.

Equations in the unexposed and unvaccinated population (superscript UNV)
Equations in the exposed and unvaccinated population (superscript ENV)

Equations for the number infected
Typical model outputs. In Fig. S3, we present typical model outputs for the population size, the number of vaccinees, and the number of individuals under antibiotic treatment. Because these outputs varied little across simulations, results of only one simulation are given.
Derivation of the VNV and SR models. As explained in the Main Text, we formulated two hypotheses regarding the cause of serotype replacement in PM, which were translated in two distinct models.
The Vaccine/NonVaccine (VNV) model integrated possible transmissibility and/or invasiveness differences between vaccine and nonvaccine serotypes. Therefore, for this model: In practice, we estimated (β V , β NV /β V ) and (ρ V , ρ NV /ρ V ) for this model.
The Susceptible/Resistant (SR) model integrated possible transmissibility and/or invasiveness differences between penicillin-susceptible and penicillin-resistant pneumococci.
Therefore, for this model: For the parameters β S and β V , i.e., absolute transmission rates, preliminary analyses showed a nonidentifiability, which reflected the lack of information in the data on absolute carriage levels in the population. To address this, we added a fifth component in the likelihood, 2 ), where c 0 is the initial carriage prevalence (Table S4). This distribution was chosen to add a baseline value for carriage prevalence in the population, but with a sufficiently large variance to allow variations over time, e.g., related to vaccination and reduced antibiotic use.
In practice, this solved the problem of identifiability for the two parameters by pushing their values to regions that resulted in carriage prevalence close to c 0 . For completeness, we also performed estimations without adding L C ; the other parameter estimates were identical (data not shown).
The likelihood was maximized with respect to the parameters using the maximum iterated filtering (MIF) algorithm 19 . This estimation procedure requires the specification of algorithmic parameters, which, once convergence has been reached, play no role in the results. For all estimations, we used the same initial variance multiplier, c 2 =9, and the same cooling factor,

Comparison of National Reference Center for Pneumococci (NRCP) versus Medicalization
Program of Information System (MPIS) data. To assess the exhaustiveness of the NRCP data, we compared annual PM numbers notified to NRCP with those reported to the MPIS (Table S1).
This comparison did not reveal any trend in the notification to NRCP; however, notification differed systematically between odd and even years. Some Regional Pneumococcus Observatories do not send their pneumococcal isolates to NRCP on even years, thereby explaining this pattern. Therefore, we fixed two notification probabilities from Vaccine coverage by birth cohort (Table S2).

Text S3. Detailed results and supplementary analyses
Parameter estimates. For the parameter estimates for the VNV and SR models, see Table S5 and To what extent those results, obtained in children, can be compared to our estimate in a general population is not straightforward. Indeed, carriage duration decreases with age, thereby progressively limiting competition among pneumococci. Therefore, because of the lack of age-structure, the numerical value of our competition parameter may be difficult to interpret. This interpretation is further complicated by the existence of immune mechanisms responding to carriage, co-carriage in the nasopharynx, etc., which were beyond the scope of this study.
In the VNV and SR models, vaccine effectiveness against carriage was estimated to be high (>75%), a value higher than previous estimates in 50%-60% range 11 , but should be interpreted with caution. Indeed, in any vaccination model, the proportion of children protected by the vaccine is given by e*p, where p is the proportion of vaccinated children and e the vaccine effectiveness 22 . Therefore, any misspecification of the vaccine coverage p is compensated by a change in the estimated vaccine effectiveness e, which might explain the high value we obtained.

Sensitivity analyses.
To test the robustness of our results, we conducted a series of sensitivity analyses. In particular, we sought to determine whether the misspecification of antibiotic-related parameters (see Text S1) might have impacted our main results on fitness costs; i.e., parameters β R /β S and ρ R /ρ S . To do so, we varied the values of several pairs of parameters and performed the estimations as before.
Fig . S5 presents the sensitivity-analysis results for β R /β S , when varying ϕ S and ϕ R .
Increasing ϕ R values shifted the competitive balance in favor of resistant strains, an effect counterbalanced by higher estimated fitness costs; the opposite was true for ϕ S . More generally, the estimated fitness costs increased for higher values of parameters favoring resistant strains (e.g., σ R and ϕ R ) and decreased for higher values of parameters favoring susceptible strains (e.g., σ S and ϕ S ). The fitness costs also increased with longer durations of carriage 1/γ. These findings are consistent with theory, and were reported previously 23 .  Taking antibiotics

Fig. S9. Qualitative comparison of yearly age-stratified incidences between the data (top graph) and simulations from an age-structured model (bottom graph)
. We extended our model to include the following age structure: [0, 3) y, [3,16) y, and ≥16 y. Durations of carriage were assumed to decrease with age: 2 months in [0, 3) y, 1 month in [3,16) y, and 2 weeks in ≥16 y 27 . The contact matrix, C ij , was calculated using reported daily contacts in Great Britain 28 , converted to units of contacts per year. We then defined the transmission matrix, β ij =q i ×C ij , where q i =(1.7×10 -2 , 4.3×10 -3 , 4.3×10 -3 ) represents age-specific probabilities of pneumococcal acquisition given exposure, calibrated to reach 50% carriage prevalence in [0, 3) y, 30% in [3,16) y, and 10% in ≥16 y 27 . Although we did not investigate further, decreasing susceptibility with age might be explained by acquired immunity or differences in the nature of contacts (e.g., riskier contacts for younger children). Consistent with our goal to allow a qualitative comparison with the French data, we did not attempt to estimate the other age-specific parameters. In particular, age-specific mean infection rates were fixed at arbitrary values, and, for each age group, we represent the yearly incidence (per 100, 000) divided by the total incidence (per 100, 000) in 2001.
Therefore, we insist that this comparison is qualitative, not quantitative. For the deterministic simulation, β R /β S (fitness cost on transmission) was fixed at 0.945, all other parameters were fixed to values estimated (Table S5) or fixed (Table S4)     Parameter unit, if any, is indicated between parentheses. *We were unable to fit the VNV model with 7 estimated parameters. Therefore, we estimated 4 more parameters (antibiotic-related parameters) for this model. Notably, the estimates of σ S and σ R were unrealistic in this case (σ S >