Geographic variation in pneumococcal vaccine efficacy estimated from dynamic modeling of epidemiological data post-PCV7

Although mean efficacy of multivalent pneumococcus vaccines has been intensively studied, variance in vaccine efficacy (VE) has been overlooked. Different net individual protection across settings can be driven by environmental conditions, local serotype and clonal composition, as well as by socio-demographic and genetic host factors. Understanding efficacy variation has implications for population-level effectiveness and other eco-evolutionary feedbacks. Here I show that realized VE can vary across epidemiological settings, by applying a multi-site-one-model approach to data post-vaccination. I analyse serotype prevalence dynamics following PCV7, in asymptomatic carriage in children attending day care in Portugal, Norway, France, Greece, Hungary and Hong-Kong. Model fitting to each dataset provides site-specific estimates for vaccine efficacy against acquisition, and pneumococcal transmission parameters. According to this model, variable serotype replacement across sites can be explained through variable PCV7 efficacy, ranging from 40% in Norway to 10% in Hong-Kong. While the details of how this effect is achieved remain to be determined, here I report three factors negatively associated with the VE readout, including initial prevalence of serotype 19F, daily mean temperature, and the Gini index. The study warrants more attention on local modulators of vaccine performance and calls for predictive frameworks within and across populations.


S1 Pneumococcus colonization model with vaccination
The proportions of children in different epidemiological classes change according to the equations: where the forces of infection are: λ V = β(I 1 V + I 1 V V + I 1 V N /2 + I 0 V + I 0 V V + I 0 V N /2); for vaccine serotypes and λ N = β(I 1 N + I 1 N N + I 1 V N /2 + I 0 N + I 0 N N + I 0 V N /2), for non-vaccine serotypes. The superscripts 0 and 1 denote non-vaccinated and vaccinated children respectively, assumed to homogeneously mix in a day-care setting. The serotype groups share similar transmission, clearance and competition rates.
This model can be reduced to a 10-dimensional system using Here ρ ∈ [0, 1] denotes the proportion of susceptibles vaccinated at birth, and w = 1 − V E ∈ (0, 1) denotes the reduced rate of acquisition of vaccine serotypes (VT) in vaccinated children.

S2.2 Age distribution of children in the different studies
The expected duration of DCC attendance in each setting, is approximately calculated from the relative frequencies of children in each age sub-group, multiplied by the remaining time until maximum age in that DCC setting.

Age-group
Mean % over all years of survey S3 Parameter estimates with recruitment rate µ based on mean age Table S2: Parameter estimates for pneumococcus transmission and PCV7 vaccine efficacy (V E) in day care settings, assuming a different birth/death rate µ=1/mean age in the SIS model. From Table 1 in the paper we can obtain µ as: 0.021 (Portugal), 0.025 (Norway), 0.043 (France), 0.021 (Greece), 0.018 (Hungary) and 0.021 (Hong Kong). Here, the means of the posterior distributions and the 95% credible intervals are given. The parameter k is set to a fixed value in 3 cases, equal to 0.05. The transmission rate per month is derived as: S4 Sensitivity of VE estimates with respect to µ and ρ S4.1 Keeping step-wise increasing vaccination coverage and varying µ I performed a sensitivity analysis for each country, varying susceptible recruitment rate µ in the range spanned by 1/mean age and 1/mean DCC attendance, using the same step-wise increasing vaccination coverages ( Table 1 in the paper) and visualized here in Figure S1.

Hong Kong
Year of survey ρ Figure S1: The increasing vaccine coverage rates in each setting. As extracted from the papers and reported in Table 1. These correspond to the updated ρ parameter values at the beginning of each simulation period corresponding to one year in the epidemiological dynamics with vaccination.  Figure S2: The relationship between assumed susceptible recruitment rate and estimated vaccine efficacies. For each country a range of five µ values was tested between 1/mean age and 1/mean DCC attendance. The estimated VE values from model fitting (same MCMC procedure explained in the paper) with step-wise increasing vaccination coverage are denoted with the colored circles, where the right-most one (filled) corresponds to the values reported in the paper.
The values of estimated vaccine efficacies are robust to changes in µ, as shown in Figure S2, and importantly, the ranking between countries is preserved. We observe that as the assumed µ increases, the estimate of VE decreases slightly for the same datasets under the original vaccination coverages. Yet, the change is not major, except for Hungary, where we cover the largest uncertainty in µ values. Furthermore the ranking is still preserved between Norway and Hungary displaying higher realized VE, than France and Greece, with Portugal and Hong-Kong displaying lowest realized protection, independently of µ.
S4.2 Assuming fixed constant coverage throughout the survey Next, I performed the same sensitivity analysis using a continuous vaccination assumption instead (as opposed to a step-wise increasing coverage), holding fixed the vaccination coverage per year, throughout the years of the survey. This is a drastic qualitative change in those model parameters which are assumed known. Even though it may be unrealistic, it is presented to complete the overall picture and clarify the importance of vaccination coverage. The range of variation for the constant ρ in this second sensitivity analysis was chosen relative to the median of the values in Table 1 of the paper, set typically from 0.5 × median to 1.5 × median, in each country. The medians are given in   Table 1 in the paper). These values were used to fit a different version of the model with constant coverage throughout the study period, starting from t = 0. This is qualitatively different, and somewhat an over-simplification, compared to the increasing step-wise coverage, adopted in the paper, which likely represents more accurately the actual pace of the intervention. However these values are used in the sensitivity analysis, presented for completion.

Country
Median The values for estimated VE in these cases, assuming a fixed ρ throughout the survey period, are generally lower, in the range of 8-15% for Hong-Kong, Portugal, Greece, France and Hungary, and about 24% for Norway this time. These values are about half the values obtained with the step-wise increasing ρ(t). When ρ(t) changes over time, as assumed in the paper, low coverage values in the beginning of an intervention, increase VE parameter estimates for the same observations, an effect that is stronger than high vaccination coverage later on: early on processes feed back more strongly in the overall nonlinear dynamics of the system. That is why, for example, the VE estimate for Norway is about 40%, under the stepwise coverage assumption, and appears as 24% in the constant coverage model. From this analysis, we can conclude that the order-ranking of site-specific VE estimates is generally preserved, independently of slight changes in assumptions on µ and ρ. Clearly, a more realistic description of vaccination implementation in each country must take into account changes from year to year, thus the step-wise vaccination model, as adopted in the paper (and extended in S6.1), remains a more appropriate approximation, likely yielding better parameter accuracy.  Table 1. The range of µ varied between 1/mean age and 1/mean DCC attendance. Model-fitting using the same MCMC procedure was performed on the data for each combination of µ and ρ, and the resulting VE estimates are shown as a contour plot. The filled circle reflects the value of VE corresponding to µ = 1/meanDCC assumed in the paper, and ρ = median of the values in Table 1 Figure S4: Summary of the relationship between assumed vaccination coverage and estimated vaccine efficacies, under a constant coverage model. This figure summarizes figure S3. For each combination (ρ, µ), setting-specific estimates of VE were obtained, then averaged over µ. A) VE averages over µ for five values of vaccination coverage (colored circles). The higher the fixed coverage assumed, the lower the vaccine efficacy estimated, to explain the same dataset. B) The VE estimates assuming µ = µmean and ρ = median fixed for each country over the survey period. The ranking between settings is preserved, although the mean efficacy appears lower. The bias in estimates toward lower values is due to the reduction from a step-wise increasing coverage to a constant coverage assumption, which is an over-simplification.  Figure S5: The relationship between target serotype decline after vaccination and transmission intensity in a neutral model. a-b) Vaccine coverage ρ = 0.5. c-d) Vaccine coverage ρ = 0.25. Panels show prevalence ratio P V /P N of target vs. non-target serotype prevalence against time, and relative prevalence P V /(P V + P N ) of target serotypes following vaccination for the same initial conditions, for two different vaccine efficacies (green and blue lines), and R 0 . The line styles depict different values of R 0 : 5 (solid line), 2.5 (dashed), 1.5 (dotted). Other parameters: competition coefficient k = 0.05; monthly clearance rate γ = 0.7; monthly susceptible recruitment rate µ = 0.02.  Figure S6: 95% credible envelopes for model trajectories, accounting for sampling process over time. These envelopes capture simultaneously uncertainty around the posterior parameter values, and the variability introduced by the sampling process. Here the same sample size was assumed in each setting equal to the mean in Table 1, fixed throughout the years of the survey. 500 model simulations with random parameter combinations, and superimposed sampling of binomial proportions were considered to compute the 95% CI of model-predicted dynamics.     (Table S1) with regards to proportion of shared NVTs among the 4 top-most prevalent post-PCV7 (number/4). Assuming independence between pairwise comparisons, the association emerges to be significant (p-value for slope p < 0.05). B) Averaged NVT sharing across pairs for each setting (mean over 5 pair-wise comparisons) and deviation from mean VE (dashed line). A decreasing trend, although it must be taken with caution given the limitations of the present dataset, suggests that ecological niche overlap, even though too coarsely captured via NVT subset overlap, may be yet another factor influencing site-specific variation in realized vaccine protection against PCV7 pneumococcal serotypes.