Abstract
Although pneumococcal conjugate vaccines (PCVs) have greatly reduced the incidence of invasive diseases caused by vaccine-targeted serotypes (VT) of Streptococcus pneumoniae, vaccine impact may be eroded by the increase in rates of disease caused by non-vaccine serotypes (NVT) – a phenomenon known as serotype replacement. Here, we investigated the effect of social contact patterns on the dynamics of vaccine impact and serotype replacement in carriage.
We developed a neutral, age-structured, susceptible–colonized (S–C) model incorporating VT-NVT co-colonization and childhood immunization with PCVs and verified it against real-world carriage data. Using empirically derived contact matrices from 34 countries, we assessed the impact of contact patterns of different age groups on the time-to-elimination, here defined as the time taken for the proportion of VT among circulating serotypes to drop to 5% of the pre-PCV level. Finally, we quantified the contribution of various parameters—such as vaccine efficacy, coverage, immunity waning, and population susceptibility—to the dynamics of VT elimination.
Our model recapitulated the observed prevalence of carriage of VTs observed in the real-world data and showed that varying the contact structure alone led to different time-to-elimination (range: 3.8 – 6 years). We found that higher total contact rate and assortativity in children under 5 were key factors in accelerating VT elimination. In addition, higher vaccine efficacy and coverage, and slower immunity waning led to shorter time-to-elimination.
These findings illuminate the mechanisms controlling the dynamics of vaccine impact and serotype replacement and may help predict the impact of the higher-valency PCVs in communities with different contact patterns.
Introduction
Streptococcus pneumoniae is one of the five leading pathogens for the estimated 7.7 million bacteria-associated deaths globally [1]. The first several generations of pneumococcal conjugate vaccines (PCVs) have reduced invasive pneumococcal disease (IPD) substantially in all age groups [2]. However, the reduction in rates of disease caused by vaccine-targeted serotypes of pneumococci (VT) was partially offset by an increase in rates of disease caused by non-vaccine-targeted serotypes (NVT) [3,4]. This phenomenon, known as “serotype replacement”, occurred because PCVs targeted a subset of the 100 identified serotypes [5], reducing the fitness of VT and changing the competitive balance between VT and NVT [6,7]. Nasopharyngeal carriage is a prerequisite for pneumococcal diseases, and the reduction in carriage in immunized children leads to indirect protection of unvaccinated children and adults [8]. Likewise, serotype replacement in carriage may erode the population-level impact of PCV and thus demands public health attention.
Observed serotype replacement in diseases was initially more pronounced in the UK than in the US, for which multiple possible explanations have been suggested: the distribution of risk factors, the vaccination schedule and coverage, and the pre-PCV composition of circulating serotypes [9]. While replacement in diseases is partial, replacement in carriage is almost complete, and it occurs faster in some populations than others [3,10,11]. However, the mechanisms driving such variation remain unclear. One potential determinant for the serotype replacement dynamics is the social contact structure in a population. Carriage studies have shown that social contact with preschool-age children is associated with higher prevalence of pneumococcal carriage [12,13]. While social contact structures are thought to be major drivers of infectious disease dynamics [14], there have not been studies investigating the effect of social contact structure on the dynamics of vaccine impact in pneumococcal carriage. Addressing this knowledge gap can elucidate the potential mechanisms controlling the dynamics of vaccine impact and serotype replacement and may help predict the impact of the higher-valency PCVs in communities with different contact patterns.
In this study, we developed a mathematical model parameterized with empirical data to simulate the dynamics of serotype changes after PCV introduction (Figure 1). Using inferred contact matrices from the literature, we interrogated the impact of social contact patterns on the trajectory of VT carriage decline and increases in NVT carriage, and quantified the effect of key parameters such as vaccine efficacy (VE) and population susceptibility (Figure 2). Our findings showed that variations in social contact structure alone could lead to different time-to-elimination. We found high association between the contact pattern features in children under 5 and time-to-elimination. More broadly, our findings highlight the need to consider social contact structure when assessing the impact of vaccines.
Results
Real-world parameter sets allow the model to reproduce observed VT-carrier prevalence in children
We formulated a deterministic, Susceptible–Colonized model that simulates the transmission of VT and NVT carriage before and after the introduction of PCVs. The model was an instance of neutral null models proposed by Lipsitch et al. for multistrain pathogens [15]. A key property of these models is the lack of a stable coexistence equilibrium, so that any initial level of coexistence will be maintained over time for identical strains. Hence, these models do not contain hidden mechanisms that artificially promote stable coexistence.
In general, the simulations using location-specific parameter sets (Table 2) captured the observed dynamics of VT-carrier prevalence in children in the post-PCV era in the UK, Alaska (US), and Massachusetts (US), and with some discrepancy, in France, where the observed VT-carrier prevalence declined more rapidly than in the simulation (Figure 3). The VT-carrier prevalence in the pre-vaccine era was higher in France (43.9%, 95% CI 38.4–49.4%) [16] and the UK (31.9%, 28.1–36.1%) [17] than in Alaska (20%, 15.7–24.7%) [18]. For Massachusetts, the VT-carrier prevalence was 9.7% half year after vaccine introduction [19]. In all locations, the rapid decline in VT-carriers immediately after vaccine introduction was followed by a slower decline as VT-carriers became less prevalent.
The time-to-elimination was predicted to be shortest in children aged 1-5
We defined time-to-elimination as the duration between vaccine introduction and the time point when replacement was considered complete (see Methods). Using contact matrices derived from census and survey data in 34 countries [20], our transmission model produced variable time-to-elimination, ranging from 3.8 to 6 years in newborns, which was a fully unvaccinated age population and thus reflected the indirect effect of PCV introduction. The time-to-elimination in adults was similar to that in age 0 (Figure 4). In contrast, the time-to-elimination was the shortest in children of age 1 and above until age 5 in most countries and until age 10– 11 in Ireland, the Netherlands, and the US. This finding corresponded well with the observation that PCV impact could be observed earlier in children than in adults [21], which is likely due to children of these ages having received the vaccine themselves and benefiting from both direct and indirect protections. They were also the age populations with the highest VT-carrier prevalence (Supplementary Figure 4) and contact rate (Figure 5A).
In the sensitivity analyses, we considered two scenarios: (1) a lower prevalence of carriers at age 0 due to the time lag from birth to first pneumococcal acquisition, and (2) a higher prevalence of carriers in all ages to simulate settings with higher pneumococcal burden (Supplementary Figure 4). The results remained similar (time-to-elimination range: 4.2–7.1 years, 4.4–6.9 years).
Time-to-elimination was highly dependent on contact patterns in children under 5
To delineate the effect of mixing patterns in different age groups, we looked at two age group-specific social contact features that may be important for respiratory infection transmission [22]: contact rate (total daily contacts) and assortativity (fraction of within-group contact).
Across countries, the contact rate increased from children under 5 (0–4y) to peak around school-age (5–9y) and teenage years (10–19y), and declined towards older age (65–84y) (Figure 5A). In general, assortativity tended to be the lowest in children under 5. The change of assortativity with age was more variable than that of contact rate across countries, and we noted four major patterns (Figure 5B). In some contact matrices (e.g., Australia, China, India, Japan, USA), the fraction of assortative contact increased with age, reaching the peak among teenagers, and either remained high (e.g., Australia, Japan, USA) or declined (e.g., China, India) in adulthood. In other contact matrices (e.g., Canada, South Africa), assortativity were similar from birth until teenage and then either remained (e.g., Canada) or declined towards older age (e.g., South Africa).
In our simulation, we assumed children under 5 had the highest prevalence of carriers based on a systematic review [23]. However, this age group bore lower contact rates than the other age groups (Figure 5A). In contrast, age groups with higher contact rates (5–39y) tended to have lower carriage prevalence (Figures 5A, Supplementary Figure 4).
A plot of simulated time-to-elimination against contact rate and assortativity revealed a strong negative correlation between the contact patterns and time-to-elimination in children under 5 but not in other age groups (Figures 5C, Supplementary Figure 8). Using a generalized linear model (GLM), we found that contact rate and assortativity in children under 5 explained most of the variability in the simulated time-to-elimination (R-squared of 0.95). Both features accelerated reduction of VT (Figure 5D): one standard deviation of increase in total contact rate and fraction of assortative contact shortened time-to-elimination by 5.2% (95%CI 3.7– 6.7%) and 7.7% (95%CI 6.3–9.2%) respectively. Hence, these two features of social contacts in children under 5 explained a large part of the variability in time-to-elimination, highlighting the key role of this age group in the transmission dynamics of pneumococcal carriage.
Higher vaccine efficacy and coverage and slower waning of vaccine immunity accelerate time-to-elimination
To investigate the effect of key parameters on time-to-elimination, we varied one parameter at a time and measured the simulated time-to-elimination using a subset of 7 contact matrices representative of all the contact matrices from [20]. The key parameters tested were VE against carriage acquisition, vaccine coverage in the target age group (1-year-old), waning rate of vaccine immunity, the initial proportion of VT- and NVT-carriers, and population susceptibility to carriage acquisition (with 3 levels considered: high, medium, or low). Table 1 shows the list of parameters in the model.
Among the key parameters studied, vaccine factors resulted in the most prominent changes in time-to-elimination. When vaccine coverage reached 90%, using a highly efficacious vaccine (VE=77%) led to a 1.9–2.6-year reduction in time-to-elimination compared with a less efficacious vaccine (VE=33%) (Figure 6A). At lower coverage (50%), VT elimination was slower (5.1–7.8 years vs. 3.8–6 years in 90% coverage), and the same increase in VE caused a greater reduction in time-to-elimination (Supplementary Figure 5). In addition to no waning, we tested various durations of vaccine-conferred immunity and found that rapid waning (immunity duration of 3 years) slowed elimination by 0.5–1.6 years compared with slow waning (immunity duration of 10 years) (Figure 5B).
Given a fixed pre-PCV total pneumococcal carriage, increasing the initial proportion of VT among colonizing serotypes (quantity F, see Methods) by changing initial VT:NVT:Co-carriers ratio slowed elimination slightly (Figures 6C, Supplementary Figure 6), while maintaining F led to constant time-to-elimination (Supplementary Figure 7).
To vary population susceptibility, we changed the age-specific susceptibility parameter β(i) by ±20%. As expected, the total pneumococcal carriage pre-PCV increased with population susceptibility. However, the predicted effect of this parameter was moderate: transitioning from low to high population susceptibility resulted in 0.6–1.8 years longer time-to-elimination. This result may be explained by the fact that, for higher total carriage, more circulating VT had to be replaced (Figure 6D).
In summary, of all the parameters tested, the vaccine parameters had the strongest impact on time-to-elimination, while the other parameters had a more moderate effect. This result highlights the need for accurate estimates of PCVs properties to predict the time scale of VT elimination in a target population.
Discussions
The main goal of this study was to assess the effect of social contact structure on the impact of PCVs. To do so, we designed a realistic model of pneumococcal carriage transmission, parameterized based on empirical data and the best available social contact matrices in different regions worldwide. Our results show that heterogeneity in contact structure alone may lead to a range of time-to-elimination and thus sensitively affect the impact of PCVs. In addition, they highlight the key role of contact features in children under 5 in VT elimination and provide new insights into the mechanisms of VT elimination. More broadly, these findings identify social contact structure as a new key variable affecting vaccine impact, with potential implications beyond PCVs.
Our model predicted a range of time-to-elimination (3.8–6 years) that is consistent with the literature [24–26], in support of the WHO’s recommendation that 5 years of post-PCV data are necessary to assess pneumococcal serotype replacement [27]. The modeled time-to-elimination in our study was the shortest among vaccinated age groups, who had the most social contacts, reflecting combined direct and indirect effectiveness, in contrast to age 0, who benefited from indirect effectiveness only (Figure 4). This finding aligns with the reported direct and indirect effects of PCV on carriage [10].
Different features of a contact matrix can have different effects on infectious disease dynamics. For example, assortativity, a well-studied feature measuring the extent of preferential mixing of individuals within the same demographic stratum, was shown to drive the spread of HIV infections differently in groups with different risks [28]. Another widely used feature is the number of social contacts, which was suggested as the main factor that explained the higher COVID rates among older adults in Italy [22]. We investigated these two features in our study and observed that the total contact rates in all countries followed a similar trend, with lower contact rates in extreme ages and a peak in the adult age groups (20–64y); contrastingly, there was a much higher variability in mean assortativity across countries. We found both features in children under 5 to be significant drivers for shorter time-to-elimination, indicating children under 5 as the key age group in driving the serotype replacement dynamics, despite having a lower contact rate than other age groups. This finding echoes the evidence that children under 5 are key for pneumococcal transmission [29].
Intuitively, higher total contact rates would speed up the transmission dynamics and shorten time-to-elimination. The effect of assortativity can be explained by the high carriage prevalence in this age group. Children under 5 had the highest carriage prevalence, so a more assortative contact in this age group would promote the within-group transmission. In general, infection spreads faster for a high-risk group with assortative mixing because contacts with low-risk groups slow down the transmission dynamic [28]. These findings point to the vital role of contact patterns in the high-prevalence groups in infection transmission and can be the basis of infection preventive strategies.
In addition to social contact structure, we found vaccine factors to be the most influential parameters in the serotype replacement dynamics. This finding is consistent with the epidemiological evidence that locations with high vaccine coverage saw rapid carriage replacement [24]. We also found that rapid waning led to longer time-to-elimination. Furthermore, initial VT:NVT:Co-carriers ratio and population susceptibility had slight to moderate effects. Given the same overall carriage, initial VT:NVT:Co-carriers ratio only affected the time-to-replacement if the proportion of VT among circulating serotypes, F, was changed: higher F led to a slightly longer time-to-elimination. The time-to-elimination was also longer in a more susceptible population whose overall carriage is higher. These results demonstrated that when the circulating VT burden is higher, it takes longer for replacement to be complete.
Our study has several limitations. When simulating the VT-carrier prevalence in children using location-specific parameter sets and considering the uncertainty in VE against colonization, our model captured the patterns in the observed data in the UK, Alaska (US), Massachusetts (US), but, with some discrepancy in the initial post-PCV era, in France. This discrepancy potentially stemmed from the partial uptake of PCV in the private market, reaching a vaccine coverage of over 20% in the year before vaccine introduction in our simulation [30]. In our model, we did not consider seasonal fluctuations in contact rates and assortativity, which could affect the transmission of infections [31,32]. How this biased the estimated time-to-elimination depends on whether holidays increase the total contacts considering the change in both inter-age and intra-age contacts. Most of the contact matrices used in this study came from high-income countries, limiting our findings’ generalizability. While there were published contact matrices for more countries, the ones used in our study offer the best age resolution to date. The variation of contact patterns across geographic locations and income settings is expected to be larger than observed in this study, and including them in future studies can help elucidate the phenomenon observed outside high-income settings. For instance, in high transmission settings, the pneumococcal reservoir may involve a wider age group, including school-age children [29]. Given the evidence of a higher extent of serotype replacement in indigenous children in Fiji [33] and in rural areas in Nigeria [34], future studies should explore further how contact patterns in sub-populations within the same country influence serotype replacement.
Despite these limitations, our study demonstrated how to combine contact matrices and mathematical modeling to unravel the dynamics between the host, the pathogen, and a public health intervention. The strengths of our study included using a neutral model without a hidden mechanism that promotes the co-existence of VT and NVT, which reduced bias, and using contact matrices of high age resolution, which allowed us to differentiate the transmission dynamics in every year of age. In addition, for parameters with variable estimates, we based our assumptions on non-linear models fitted to extracted data from observational studies. In conclusion, our findings demonstrate that, as for other vaccine-preventable diseases, social contact structure is a critical element for understanding the vaccine epidemiology of pneumococcus. Hence, we propose this element should be considered in future studies assessing the impact of PCVs and, more broadly, of other vaccines.
Methods
Data
a) Contact matrices
We used the inferred contact matrices Mij from [20]. The contacts, stratified by age yearly from 0 to 84, were derived from synthetic networks built using population census data and socio-demographic surveys in various settings, namely, household, school, workplace, and community. The overall contact matrix for a location is a weighted sum of the setting-specific contact matrices. Mij gives the total number of daily contacts between age groups i and j per person in age group i, we applied reciprocity correction on Mij and transformed it into , which gives the total annual contacts between age groups i and j per person in age group i and per person in age group j (density scale, as defined in [35]; see Supplementary Figure 2).
b) Carriage duration and prevalence
We extracted data about age-specific carriage duration and carriage prevalence from published studies identified through a scoping literature search (Supplementary Data 1, Supplementary Data 2).
Among the identified culture-based studies, we included the studies that reported median duration (n=8), because the duration of carriage has a left-skewed distribution, with few individuals showing lasting carriage. For the model of carriage duration with age, we used non-linear least square regression to estimate the parameters in the equation (Supplementary Figure 1): where a = 21 (standard error: 5.8), b = 62 (14.6), and c = 0.45 (0.4).
In the main analysis, we fixed the age-specific initial carriage prevalence based on [23] and the age-specific susceptibility parameter β(i) based on [47]. As sensitivity analyses, we used two other distributions of β(i) over age, considering a lower carriage prevalence in age 0 and a higher carriage prevalence in all ages, to reflect the observed data from the identified culture-based pre-PCV carriage prevalence studies (n=17) (Supplementary Figure 4). After calibrating β(i) for the assumed carriage prevalences, we re-simulated the time-to-elimination for all countries.
c) VT-carrier prevalence in children in the real world
We extracted the VT-carrier prevalence in children from the pre-to post-PCV era in 4 locations – France [16], UK [17], Alaska, US [18], and Massachusetts, US [19] (Data S3) – to verify our model’s ability to reproduce the decline in VT-carrier prevalence following PCV introduction. The observed data come from cross-sectional surveys among children attending daycare centers or primary care clinics. In all four included studies, the detection of S.pneumoniae was culture-based and the serotyping was either by traditional Quellung reaction or molecular methods. For point estimates of carriage reported without uncertainty, we calculated the standard error (SE) for proportion and indicated the uncertainty limits as 1.96×SE from the mean.
Model
a) Structure
We formulated a deterministic model that simulates the transmission of VT and NVT carriage based on the neutral null model proposed by [15]. Assuming a stable population (i.e., birth rate = death rate), susceptible individuals (S) become VT-carriers (CV) at the rate λV, or NVT-carriers (CN) at the rate λN. Mono-carriers CV (or CN) can be colonized by the other serotype at rate kN × λN (or kV × λV) and become co-carriers (CVN), and co-carriers return to mono-carriers CV (or CN) at rate c × kV × λV) (or c × kN × λN). The inter-serotype competition parameter, k, was assumed to be 0.5 based on published estimates [36]. When kV=0.5, VT is half as likely to colonize an individual already colonized by NVT. We further assumed this competition to be symmetrical (kV = kN) to ensure neutrality at initiation. The parameter c, representing the fraction of co-carriers returning to CV (or CN) upon re-infection with VT (or NVT), was fixed to 0.5 to ensure neutrality [37].
The vaccine was introduced at time tV and had a coverage of pV. Therefore, pV was zero before time tV and equal to pV starting from time tV.
In our age-structured model, individuals moved from one age to the next year of age at an aging rate δi=1 per year. The whole population of newborns was unvaccinated. As individuals moved from age 0 to age 1, a fraction (pV) of the population was vaccinated and partially protected from pneumococcal colonization (superscript “(V, 1)”). The rest (1– pV) of age 0 stayed unvaccinated as they reached age 1 (superscript “(N, 1)”).
For the dynamics of the vaccinated individuals, the rate of VT carriage acquisition λV was reduced by a factor ϵV, where ϵV represents the vaccine effectiveness against acquisition of VT carriage. Vaccine-conferred immunity was assumed to wane at a rate αV, so that 1/αV represents the average duration of vaccine protection.
The age-specific acquisition rate, λ(i), depends on β(i), the cumulative number of carriers in the contactee age groups, CC(j), and the per capita contact matrix, The carriage acquisition rates for VT and NVT were expressed as : where Here, q refers to the relative infectiousness with each serotype for co-carriers. Table 1 summarizes the parameters used in this study.
b) Outcome definition
In a neutral null model, one serotype is not assumed to have a fitness advantage over the other; therefore, co-carriers transmit either VT or NVT at equal probability. The relative infectiousness with each serotype for co-carriers, q, is set to 0.5, such that co-carriers are equally infectious as mono-carriers [15]. To ensure neutrality in the null model, we checked if F was stable over time in the model without an effective vaccine, as suggested by [15] (Supplementary Figure 3). F is given by: We defined time-to-elimination as the duration between vaccine introduction and the time when F dropped to 5% of its initial value in age 0, representing a fully unvaccinated population and reflecting the indirect effect of PCV introduction. As time-to-elimination in all ages were highly correlated in each country, the choice of age had a negligible effect on the analyses comparing countries.
Model assessment
To verify the model, we used location-specific parameter sets (Table 2) to simulate the VT-carrier prevalence in children from the pre-to post-PCV era in 4 locations (Table 2) and compared the simulated values to the observed ones qualitatively.
We calibrated the model for each location by estimating the parameters β(i). First, we performed a global search on 1000 values between –10 and 10, corresponding to β values between 0 and 1 on the logit-transformed scale. The values were sampled using Sobol’s sequence [38], a quasi-random sampling method, to ensure the global parameter space was searched thoroughly. The global search sought a set of β(i) that minimized the total squared difference in simulated versus observed pre-PCV VT-carrier prevalence on the logarithmic scale in all age groups. Here, the age groups were defined based on the observed prevalences as 0 y, 1–4 y, 5–17 y, 18–39 y, 40–59 y, and 60–84 y. The best five solutions from the global search were used as the starting value for a local search using the Subplex algorithm [39] until the total squared difference was minimized or could not be further reduced after a maximum of 1000 evaluations. Given the uncertainty around this parameter, we simulated the VT-carrier prevalence in children in each location using a range of VE against colonization [40].
Effect of contact features
To investigate the effect of social contact structure on the replacement dynamics, we first summarized the contact matrices using two age group-specific features – contact rate and assortativity – and then explored the relationship between these features and the time-to-elimination. Here, the age groups were defined as 0–4 y, 5–9 y, 10–19 y, 20–39 y, 40–59 y, and 60–84 y, to be consistent with the parameter value assignment (Table 1). We defined contact rate as the average total daily contacts in an age group and assortativity as the fraction of contacts from within the age group out of total contacts for each age in the age group.
We described the distribution of total contact and assortativity over ages in all 34 contact matrices (Figures 5A, B) and explored the association between time-to-elimination and these contact features in all age groups (Supplementary Figure 8).
Based on the strong negative correlation between the two contact features and time-to-elimination portrayed in children under 5 (Figure 5C), we performed a regression analysis on time-to-elimination with standardized contact and assortativity as covariates in this age group using a GLM with a log link (Figure 5D). We reported the effect estimates with 95%CI for both variables and assessed the goodness-of-fit with R-squared.
Effect of key parameters
To delineate the individual effect of the key parameters—VE, vaccine coverage, immunity waning, the initial proportions of VT and NVT carriers, and population susceptibility—on time to elimination, we varied them one at a time using a range of values and measured the time to elimination.
VE and coverage were considered key parameters because these contributed to the selective pressure that drives serotype replacement. We varied VE between 33–77% based on the observed efficacy with uncertainty in a community randomized trial [40], consistent with the findings of a systematic review [41]. Other than no waning, we tested a range of durations of vaccine-conferred immunity, ranging from 3 to 10 years [42,43]. Evidence suggests pre-PCV serotype distribution in carriage and diseases as important predictors of vaccine impact [44]; therefore, we tested a range of initial proportions of VT-carriers (fV(0)), NVT-carriers (fN(0)), and implicitly, co-carriers (1–fV(0)–fN(0)), either allowing the proportion of VT among colonizing serotypes (F) to fluctuate or be fixed at 0.65.
In these model experiments, the age-specific overall carriage remained constant. Lastly, to investigate the dynamics under different population susceptibilities, we changed the age-specific susceptibility parameter, β(i), by ±20% compared to the baseline value, which led to higher and lower overall carriage, respectively. In each simulation, we used 7 contact matrices from [20] to see if the effect of each key parameter differs by social contact structure.
Numerical implementation
All analyses were conducted in RStudio with R version 4.2.2 (R) and the non-linear model fitting was performed using the base package “stats” [45]. The transmission model was implemented using the package “pomp” version 4.6 [46]. All optimization procedures were implemented using the algorithms available in the package “nloptr” version 2.0.3 [47].
Data Availability
The code and data are available online at Edmond, the Open Data Repository from the Max Planck Society.
Code and data availability
The code and data are available from Edmond, the Open Data Repository from the Max Planck Society: https://doi.org/10.17617/3.RIGYAK.
Supplementary materials
To obtain the duration of carriage, we fitted a non-linear function of age to the extracted duration of carriage from published longitudinal carriage studies identified through a scoping literature search (Data 1). We first searched for observational studies on pneumococcal carriage duration on PubMed and further identified relevant studies from the references of the initially included studies. We selected culture-based studies to allow the inclusion of the maximum number of studies because the majority of the early studies relied on culture-based detection. Further, because the duration of carriage has a left-skewed distribution, with few individuals showing lasting carriage, we preferred median to mean reporting and included studies that reported median duration. A total of 8 studies were included [1–8]. For the model of carriage duration with age, we used the non-linear least square algorithm to estimate the parameters in the equation: where a = 21 (standard error: 5.8), - = 62 (14.6), and 6 = 0.45 (0.4).
Figure 1 shows the extracted data from observational studies as grey points and the modeled duration of carriage with a red curve.
Let M = (Mij) represent the original contact rate matrix calculated in [9]. By definition, where Eij represents the numer of daily contacts between age groups i and M j, and Ni the population size of age group i. Because of the necessary reciprocity of total contacts (i.e., Eij = Eji), the per capita contact matrix should be symmetric. To ensure this symmetry, one can calculate a contact matrix corrected for reciprocity based on the population structure in the study population [10]: As a result, the per capita contact matrix: is symmetric, as it should be.
Finally, we multiplied the per capita matrix by 365 to obtain a per capita annual contact matrix because all rates were per year in the simulations.
Figure 2 shows the pre- (A) and post-transformation (B) contact patterns of India, as an example, and the post-transformed matrices of the other six countries (C–H).
To check if our null model fulfils the neutrality criterion, we tracked proportion of VT among all colonizing serotypes (F) as proposed by [11]. The quantity F is given by: In a neutral null model, where one serotype is not assumed to have a fitness advantage over the other, any level of co-existence should be permitted and the proportions of the two serotypes should not converge to a global equilibrium (e.g., F = 50%) without explicit mechanism [11]. In other words, neutral models do not lead to a 50-50 coexistence by default. Specifically, ! should remain constant regardless of the initial proportions of VT and NVT under no intervention, where the intervention represents a mechanism that gives one serotype a fitness advantage over the other.
Figure 3 left panel shows that F remained stable when a null-impact vaccine was introduced for all contact matrices. The right panel shows that F declined after the introduction of a vaccine with VE=0.6 at 90% coverage, at different rates in different contact structures.
The transmission dynamic is described by the following system of ordinary differential equations:
Equations in (unvaccinated) newborns (i = 0)
Newborns are assumed to be non-carriers and, therefore, directly enter the S(N,0) compartment.
δ0 represents the aging rate, equal to , where a0 represents the age span of age group 0.
is the total population size (summed across all age groups).
μ is the per capita birth rate. To keep the population constant, this is calculated as μ−1 = a = ∑i ai, where ai is the age span of age groups i and a represents the assumed lifespan. For example, if one assumes a lifespan a = 80 years, then the birth rate equals per year.
In the model, newborns are born at the rate WN, then individuals age across the age groups, and all die exactly at the age of 84. This simplified demographic model is known as type-I mortality distribution [12].
Equations in infants receiving vaccination (i = 1)
Write v(t) = pV for (t ≥ tV), the time-varying vaccine coverage, where tV represents the time point of PCV introduction and pV the proportion of infants vaccinated after PCV introduction.
Vaccination is assumed to occur at age 0, that is, when the newborns age to the second age group. The other properties of the vaccine are its efficacy against acquisition of V-serotypes (denoted byϵV) and its rate of waning protection (denoted by αV, where 1/αV represents the average duration of vaccine protection).
Equations in older age groups (i = 2, …, A − 1)
In the older age groups, no more vaccination is assumed, and the dynamic is described by the following system of ordinary differential equations:
Initial conditions
The model was first run to simulate the pre-vaccine era, so all corresponding initial conditions were set to 0 = (X(V, i) = 0). For the initial conditions in the non-vaccinated groups, we defined the following parameters:
: initial prevalence of carriers in age group i
fV (0) : proportion of VT carriers / all carriers
fN (0): proportion of NVT carriers / all carriers
fV N (0)) = 1 − fV (0) − fN(0) : proportion of dual carriers/all carriers
With these parameters, the state variables were initialized as follows:
In Figure 4 left panels, the points show the observed carriage prevalence by age extracted from published observational studies (Data 2) conducted in high-income countries (HIC) (panel A,C) and low-income countries (LMIC) (panel E); the blue lines mark the assumed initial prevalence of carriers in the main analysis (A) and the sensitivity analyses (C, E). The right panels show the age-specific susceptibility parameter values used in the main analysis (B) and the sensitivity analyses (D, F).
In the main analysis, we assumed the initial overall carriage prevalence to be 50% in age 0-4, 20% in age 5-19, and 10% in age 20-84, based on a review [13], reflecting the carriage prevalence in the settings of high-income countries (panel A). We fixed the age-specific susceptibility parameter d(!) for low, medium, and high population susceptibility (panel B) when quantifying the influence of this parameter to time-to-elimination (results in main text Figure 6B).
As sensitivity analyses, we used different susceuptibility distributions over age considering two scenarios based on a scoping literature review.
We searched for observational studies on pneumococcal carriage prevalence on PubMed and further identified relevant studies from the references of the initially included studies. We included culture-based studies conducted in the pre-PCV era that reported either overall (including both VT and NVT) or VT carriage prevalence because VT-carriers accounted for most of the carriage in the pre-PCV era. After excluding publications that were based on the same sample, a total of 17 studies were included (HIC: n=11 [14–24], LMIC: n=6 [1,3,4,6,25,26]). We observed that the carriage prevalence at age 0 was lower that in children aged 1–4, and that carriage prevalence was higher at all ages in LMIC compared with HIC. The income group of countries were based on the World Development Indicators 2008 [27].
Based on these observations, we considered a lower prevalence of carriers at age 0 (10% instead of 50%) due to the time lag from birth to first pneumococcal acquisition (panel C) and optimized β(i) (panel D). The resulted time-to-elimination ranged from 4.2–7.1 years.
We then considered higher carriage prevalences for all ages (70% in age 0-4, 40% in age 5-19, and 20% in age 20-84) to mimic settings with higher pneumococcal burden (panel E), which is commonly observed in low-income countries (LIC) [13,28], and optimized β(i) (panel F). The resulted time-to-elimination ranged from 4.4–6.9 years.
The change in time-to-elimination due to vaccine efficacy was larger when vaccine coverage was lower. At 50% coverage, switching from a less efficacious vaccine (VE=33%) to a highly efficacious vaccine (77%) resulted in 2.5–3.5 years shorter time-to-elimination (left panel). When vaccine coverage reached 90%, using a highly efficacious vaccine led to a 1.9–2.6-year reduction in time-to-elimination compared with a less efficacious vaccine (right panel, same as main text Figure 6A).
We varied the initial proportion of VT-carriers (fV (0)) between 0.2 and 0.8, and the initial proportion of NVT-carriers (fN(0)) between 0.2 and 0.4, with the proportion of co-carriers (fCo(0)) implicitly determined (1 − fV(0) − fN(0)). The combinations of initial VT:NVT:Co-carriers ratio were restricted to two conditions: (1) fV (0) + fN(0) ≤ 1 because the sum of proportions cannot exceed one, and (2) fV(0) ≥ fN (0) because VT-carriers tended to be more prevalence than NVT carriers in the pre-PCV era (main text Table 2).
When fN(0) was fixed at 0.2, increasing fV (0) from 0.2 to 0.8 resulted in an increase in the initial proportion of VT among all colonizing serotypes (!) from 0.5 to 0.8, leading to a slightly longer time-to-elimination (left panel, same as main text Figure 6C).
This effect remained when fN (0) was fixed at 0.4 and fV (0) increased from 0.4 to 0.6, which resulted in an increase in initial ! from 0.5 to 0.6.
We varied the initial proportion of VT-carriers (fV (0)), of NVT-carriers (fN (0)), and implicitly, the initial proportion of co-carriers (1 − fV (0) – fN(0)), such that the initial proportion of VT among all colonizing serotypes (F) remained fixed at 0.65. The combinations of initial VT:NVT:Co-carriers ratio were restricted to two conditions: (1) fV (0) + fN (0) ≤ 1 because the sum of proportions cannot exceed one, and (2) fV (0) ≥ fN (0) because VT-carriers tended to be more prevalence than NVT carriers in the pre-PCV era (main text Table 2).
With initial F fixed at 0.65, the time-to-elimination remained constant regardless of the initial VT:NVT:Co-carriers ratio.
After simulating the transmission using the contact matrices from 34 countries in [1], we explored the relationship between contact features and time-to-elimination in each age group.
The time-to-elimination was measured for each age group. We then calculated the total number of daily contacts in an age group, divided by the age group width, such that this measure would not be inflated in the age groups of wider bands. For example, for age group 0–4y, the total number of daily contacts was divided by 5; for age group 20–39y, the total number of daily contacts was divided by 20. We defined assortativity as the average fraction of contacts from within the age group out of total contact for each age in the age group. We standardized both measures of contact features for easier comparison.
Figure 8 shows time-to-elimination’ correlation with standardized contact rate (x-axis) and standardized assortativity (colour scale) for all age groups, revealing a strong trend in children under 5. The R-squared value in each panel indicates the variability in time-to-elimination explained by standardized contact rate and assortativity as covariates in a generalized linear model (GLM) for each respective age group.