Within-host competition modulates pneumococcal antibiotic resistance in the pre-vaccination era

Background Antibiotic treatment is a key tool in the fight against pneumococcal infections. However, the ongoing emergence of antibiotic resistant strains and high frequencies of antibiotic resistance of pneumococci pose a major public health challenge. How and which ecological and evolutionary mechanisms help maintain the coexistence of strains susceptible and resistant to antibiotic treatment remains largely an open question. Methods/results Expanding on a Streptococcus pneumoniae modelling framework, we here explore how both between- and within-host mechanisms of transmission can sustain observed levels of pneumococcal resistance to antibiotics in the pre-vaccination era using a stochastic, individual-based model. Our framework considers that within-host competition for co-colonization between resistant and susceptible strains can arise via pre-existing immunity (immunological competition) or intrinsic fitness differences due to resistance costs (ecological competition). We find that beyond stochasticity, host-population structure or movement at the between-host level, competition at the within-host level can explain observed variation in resistance frequencies. Conclusion In a series of simulated scenarios informed by observed pneumococcal data in the European region, we demonstrate that ecological competition for co-colonization can explain much of the variation in co-existence observed at the country level in the pre-vaccination era. This work expands our understanding of how within-host pneumococcal competition facilitates the maintenance of antibiotic resistance in the pre-vaccination era. The demonstration of the effects of such underlying, often unmeasured competition-related components of pneumococcal dynamics improves our understanding of the mechanistic drivers for the emergence and maintenance of antibiotic resistance.


Introduction
The Gram-positive bacteria Streptococcus pneumoniae (or pneumococci) remains the most common cause of pneumonia in children and the elderly, and is likely responsible for over 500,000 yearly deaths in children aged 1 to 5 years [1] . Independently of age, most carriers are asymptomatic and only a small proportion may develop pneumonia, or other invasive disease such as meningitis or bacteraemia [2] . Vaccination remains the most successful public health strategy to date of reducing pneumococcal carriage and invasive disease [3,4] .
Pneumococci are characterized by a polysaccharide capsule comprising over 90 distinct antigenic types, termed serotypes [5,6] . Since the introduction of the pneumococcal conjugate vaccines, targeting seven (PCV7), ten (PCV10) or thirteen (PCV13) serotypes, there has been a rapid decrease in pneumococcal carriage and disease of the targeted serotypes across ages [7][8][9][10][11] . Concomitantly, increases in the frequencies of non-targeted serotypes have been recorded both in carriage and disease in a number of countries worldwide [12][13][14][15][16][17][18][19][20] . Several studies have explored the phenomenon of vaccine-induced serotype replacement (VISR), whereby the removal of competition from targeted serotypes allows non-targeted serotypes to fill the immunological and ecological niches now vacant due to the removal of cross-immunity or direct resource competition [21][22][23][24][25] . Particular changes in the genetic composition of pneumococci have also been observed in the post-vaccination era [26,27] , and have been proposed to reflect vaccine-induced metabolic shifts (VIMS) as a consequence of changes in competition on resource among pneumococci sharing metabolic niches [28][29][30] .
A second major route of combating pneumococcal carriage and disease is antibiotic treatment [31] . Similarly to vaccination, however, the effectiveness of antibiotic drugs is often impaired by the emergence and spread of escaping, antibiotic-resistant strains. Generally, vaccination has had a beneficial impact on the frequencies of antibiotic resistance among pneumococci, since vaccine types (serotypes targeted by vaccination, VT) are often with increased antibiotic resistance [32][33][34][35] . However, there are also reports of post-vaccination increases in resistance among non-VTs (NVT) [36][37][38][39][40][41][42][43][44] , likely due to changes in the herd-immunity and direct competition landscapes, although emergence dynamics and the serotypes involved have varied between regions.
Previously, we have put forward hypotheses regarding the possible drivers of post-vaccination increases in resistance among serotypes not targeted by vaccination. VIMS, for example, posits that the genetic profiles associated with resistance in vaccination-targeted serotypes can shift to non-targeted serotypes enjoying relaxed competition pressures in the post-vaccination era [45] . In a study by Obolski et al. [46] , we have shown that vaccination can reduce within-host competition between co-colonizing serotypes. This benefits antibiotic resistant strains disproportionately as vaccine-driven reductions in competition may relax the fitness cost of resistance. As reported in different countries [36][37][38][39][40][41][42][43][44] , vaccination can thus somewhat ameliorate the antibiotic resistance cost at the population-level and facilitate the increase in frequency of antibiotic-resistant non-vaccine serotypes.

1/17
In this study, we further explore the mathematical framework by Obolski et al., originally developed to study post-vaccination changes in the frequencies of antibiotic resistance. Here, it is expanded to study the co-existence of pneumococcal susceptible and resistant strains in the pre-vaccination era. We translate the original framework into an individual-based framework. This allows us to explore the effects of stochasticity, host population structure and mobility. We perform a series of sensitivity analyses based on the equilibrium states of stochastic simulations of the pre-vaccination era dynamics. We explore how the within-host competition parameters, defined by two different mechanisms, determine the levels of co-circulation of susceptible and resistant strains. We find that the model reproduces observed frequencies of resistance to penicillins and macrolides in Europe. In this region, direct (ecological) competition between susceptible and resistant strains explains more of the observed differences between countries than does immunological competition. Our results contribute to a better understanding of the processes driving the maintenance of antibiotic resistance in the absence of vaccination.

Epidemiological framework
We model the transmission dynamics of pneumococcal strains according to the epidemiological framework by Obolski et al. [46] . We define pneumococcal strains by the combination of their serotype and antibiotic resistance profile to two distinct antibiotics. Similarly to the nomenclature used in previous studies, each strain's genotype is conceptualised by a tuple , where determines the serotype and the antibiotic i, } { j i j resistance profile. We model two serotypes only, such that . For the entirety of our a, } i ∈ { b study, we define serotype as a vaccine type (VT) and serotype as a non-vaccine type a b (NVT). We also model two antibiotics and and consider only resistance to each x y independently, with , where refers to sensitive strains, to 00, 1, 0} j ∈ { 0 1 0 j = 0 1 j = 0 strains resistant only to antibiotic , and to strains resistant only to antibiotic . For y 0 j = 1 x instance, a strain defined by is of VT serotype and is sensitive to both , j 0 i = a = 0 a antibiotics; while a strain is of NVT serotype , sensitive to antibiotic and , j 1 x resistant to antibiotic . For simplicity, we do not consider double-resistance, i.e. in y 1 j = 1 our model. This is consistent with relatively low estimates of double resistance to penicillin and macrolides identified in Europe in pneumococci before the pneumococcal vaccine introduction [47] .
We keep the terminology of the original model formulation regarding the epidemiological parameters: is serotype-specific immunity; is the probability that a carried susceptible γ ψ strain ( ) will suppress host co-colonization by a resistant strain ( or ) due 0 j = 0 1 j = 0 0 j = 1 to the fitness cost of antibiotic resistance (details below); is the host life-span; is the /μ 1 /σ 1 host carriage duration; R0 is the basic reproduction number and is the transmission rate. β Individuals are born naive to all strains, with the size of the host population kept constant (deaths being replaced by newborns). Colonization results in carriage for an average 2/17 duration of days, from which recovery (clearance) may lead to complete serotype /σ 1 immunity if or partial life-long immunity if . Thus, if , hosts can be γ = 1 γ < 1 γ < 1 re-colonized by the same serotype throughout their lifetime, with probability . 1 − γ Co-colonization of up to two strains is allowed, unless and strains belong to the same γ = 1 serotype. Within-host strain competition interferes with co-colonization if , which ψ > 0 captures the degree to which intrinsic fitness differences (e.g. growth rates) between resistant and sensitive strains may allow a currently carried sensitive strain ( ) to 0 j = 0 suppress co-colonization by a resistant strain ( , ). We emphasize that models 1 j = 0 0 j = 1 ψ a form of competition between bacterial strains that is not mediated by immunity -from now onwards referred to ecological competition. The modelled processes related to an individual's colonization and epidemiological states are summarised in Figure 1A . The full model description is provided in the Supplementary Information . In this study we use a stochastic, individual-based model (IBM) implementation [48] of the original epidemiological framework [46] . We implement model events by sampling from a Poisson distribution with mean equal to the deterministic expectation of state transitions: e.g. at each time step, the number of deaths is obtained from a Poisson distribution with mean , where is the total population size, after which individuals are randomly selected for μ × N N a death event. The time unit is a day.

3/17
The original framework is expanded by including host-population structure and host movement within a meta-population formulation as in our previously published IBMs [49][50][51] . Although we model a closed population, we allow for random introductions of any strain under a fixed frequency ; i.e. with probability per time step a random host and strain /η 1 η are selected and colonization is attempted (which depends on the host's epidemiological state as in Figures 1A-B ). In summary, population structure is included by subdividing individual hosts into a spatially organised set of communities of population size , L × L c N forming a squared lattice of columns and rows. Individuals are assumed to mix L L homogeneously within each community. Local movement is considered, such that the force of infection of a strain in a particular community is dependent on its carriage levels within that community and in the neighbouring communities. Neighbouring communities are defined as the ones directly up, down, left, and right. Long distance host mobility is also considered and assumed to take place randomly across the meta-population with probability . That is, every time step a Poisson distributed number of individuals, with mean , ω ω × y i,j will attempt to transmit strain to a random individual in the meta-population. The , i j modelled processes related to population structure are summarised in Figure 1C . The model parametrization is given in the Model parameterization section of the Supplementary Information .

Antibiotic resistance and consumption data
Data source used was the European Centre for Disease Prevention and Control (ECDC) from the surveillance atlas (antimicrobial resistance) for the year 2005, vastly representative of a pre-PCV era across European countries (available at [52] ). The year 2005 is the earliest made available by the ECDC. The data variables used were (i) the percent of samples sensitive to penicillins and macrolides (independently), and (ii) consumption defined by daily doses per 1000 people per day. We considered solely countries which had a number of samples larger than 50, which included 19 countries. Resistance level was aggregated across countries (EU) using a logistic regression model with random-effect intercepts for each country. Raw data on resistance per country and aggregated EU is presented in Supplementary Figures S1-2 ; and for consumption in Supplementary Figure S3 .

Approximate Bayesian computation
To compare the model output to real-world estimates on the prevalence of antibiotic resistance, we employed approximate Bayesian computation (ABC). This method allowed us to approximate the posterior density of simulation parameters related to competition ( , ), ψ γ given priors on the parameters and data, when there is no explicit likelihood function [53] , as is the case for our IBM framework. The priors used for the parameters ( , ) were both ψ γ uniformly distributed in [0.1,0.9]. These priors were chosen for two reasons: First, we are not aware of any specific prior knowledge of these values -except that they are probably not close to zero or one. Second, the model run time is non-negligible, and is substantially increased when either of the competition parameters is close to zero or one. Hence, constrained by the relatively long computation time the model required, we sampled the prior a total of 6741 times (each with a unique combination of and . We used three data ψ ) γ points as 'targets' for the ABC: (i) the observed proportion of resistance to penicillins or macrolides aggregated at the EU level; (ii) the observed relative proportion of NVT in some EU-related epidemiological contexts ( Supplementary Figure S4 ); and (iii) the proportion of carriers with co-infection. The ABC distance function employed to our model was the Euclidean distance from the simulation results to targets (i-iii). The threshold for ε acceptance under the ABC framework was set to either 0.025 or 0.05 so that the closest 2.5 or 5% of simulations to the targets (i-iii) would not be rejected. Further details of the implementation of the ABC approach are given in the section Approximate Bayesian computation of the Supplementary Information .

Pneumococcal strain dynamics
We first assessed average model behaviour with regards to the levels of ecological and immunological competition, for which the parameters and were marginally varied When immunological competition was high and above , independently of the level γ .6 γ ≈ 0 of ecological competition , model dynamics were generally incompatible with observed ψ pneumococcal strain dynamics. In particular, carriage levels of VT and NVT serotypes tended to be similar ( Figure 2A ), there was no host co-infection ( Figure 2B ), and strains resistant to antibiotics dominated the population ( Figures 2C ). Such epidemiological dynamics are incompatible with what is observed for the pneumococcus, which ubiquitously include high co-infection [54,55] , reinfection [56][57][58] , and co-circulation of NVT and VT serotypes with dominance of VT [10,55,56,59,60] .
In contrast, when immunological competition was low to intermediate, the model better γ replicated known properties of pneumococcal strain dynamics, resulting in dominance of VT serotypes and strains sensitive to antibiotics, as well as allowing for host coinfection ( Figures 2A-C ). Generally, for (i.e. ), the carriage level of NVT serotypes was γ 0.35 .55] [ − 0 extremely low ( Figure 2A ); there was no co-infection ( Figure 2B ); and, while strains resistant to antibiotics did not dominate the population, their prevalence generally decreased for higher ecological competition ( Figure 2C ). In contrast, when was lower (i.e. ), ψ γ .35 < 0 carriage levels of NVT serotypes increased with ecological competition ( Figure 2A ); host coinfection was possible with varying levels of ecological competition ( Figure 2B ); and, both coinfection and sensitivity to antibiotics decreased with increasing ecological competition ( Figures 2B-C ). These general behaviours across the competition parameter space ( , ) ψ γ did not change substantially in a sensitivity analysis of host-population structure and mobility values ( Supplementary Figures S7-8 ). Hence, this, suggests that variation of demographic parameters, which often exists between different regions, does not play a significant role in changing the possible levels of co-existence between susceptible and resistant strains   Figures S1-2 ). While a wide range of and values appeared γ ψ compatible with these levels of resistance ( Figure 2C ), some of those combinations would nonetheless result in widely different levels of the NVT versus VT prevalence ( Figure 2A ) and levels of co-infection ( Figure 2B ). Using an ABC algorithm applied to a large number of numerical simulations exploring the model's parameter space, we next set out to identify which combinations of and would reproduce the observed levels of resistance while γ ψ respecting prior knowledge on the relative proportions of NVT versus VT and host co-infection (see Materials and Methods section for details).

Strain competition as determinants of antibiotic resistance
As illustrated in Figures 3A-B , reproduction of pre-vaccination resistance levels while respecting priors related to observed NVT versus VT prevalence and strain co-infection identified particular competition parameter regions. The latter are characterized by a combination of low immunological ( ) and low to intermediate ecological competition ( .25 γ < 0 ), which are compatible with observations of recurrent pneumococcal .2, ψ .7 ψ > 0 < 0 infection [56][57][58] , co-circulation of NVT and VT serotypes [10,55,56,59,60] as well as frequent co-infection [54,55] .
For penicillins, a restricted range for immunological competition (0.12 -0.18) and wide γ range of ecological competition (0.3 -0.7) were compatible with observed resistance ψ levels ( Figure 3C ). In contrast, resistance to macrolides was compatible with the model output under a wider range of (0.1 -0.25) and stricter and lower range of (0.22 -0.38). γ ψ

6/17
The intersection of the two compatible parameter regions for macroIides and penicillins was for and . Doubling the percentage of accepted simulations under the .175 γ ≈ 0 .375 ψ = 0 ABC framework ( ) or the target co-infection level (from 0.2 to 0.3) had no significant % ε = 5 impact on the identified parameter regions ( Supplementary Figures S9-10 ). In contrast, accepting higher margins of error (e.g. ) resulted in the ABC identifying substantially 0% ε = 1 wider regions of the parameter space as compatible, and thus making it difficult to infer relationships between the model parameterization and observed resistance levels. We further examined the associations between the competition parameters and resistance to penicillins and macrolides, performing an ABC which utilized each EU country's antibiotic resistance frequencies to each of the two mentioned antibiotic classes (see Supplementary Information for method details). For either of the antibiotic classes, ecological competition explained a higher percent of the variance of resistance than immunological competition ψ in a simple linear model (penicillins, of 0.9 versus 0.001; macrolides, of 0.71 versus γ R 2 R 2 0.52, respectively). Furthermore, for either of the antibiotic classes, a likelihood ratio test concluded that a bivariate model including both and did not have a significantly better γ ψ goodness of fit over a model containing alone. Hence, in our model, variations in ψ ecological competition could better explain differences of antibiotic resistance frequencies in the pre-vaccination era.

Discussion
How do ecological and evolutionary forces drive the coexistence of pathogen strains susceptible and resistant to treatments or prophylaxis? This question remains a conundrum of high public health concern. Examples of such coexistence include some of the most relevant human pathogens: antiviral resistance of the immunodeficiency virus (HIV) [61,62] , antimalarial resistance to malaria parasites [63] , drug and vaccine escape by Mycobacterium tuberculosis [64] , and antimicrobial resistance in Streptococcus pneumoniae [34,37,41,65,66] .
Specifically, the understanding of what allows resistant bacteria to persist in human populations with different epidemiological, demographic, antibiotic consumption and epidemiological control conditions remains elusive. Antibiotic consumption is considered a major driver of resistance at the population level for both commensal and pathogenic bacteria [65,67,68] . However, while antibiotic consumption undoubtedly selects for antibiotic resistance, a major fraction of the observed spatio-temporal variation in resistance prevalence remains to be explained. Indeed, assuming no substantial time lag between consumption and prevalence of resistance [69] , the ECDC data here explored shows that country-specific resistance is only mildly correlated with antibiotic consumption across countries, with consumption explaining 0.14 ( ) of the variation for macrolide resistance ( R 2 ) and 0.38 of penicillin resistance ( ) ( Supplementary Figure S3 ). .06 p = 0 .002 p = 0 Mathematical models have the potential to explore possible mechanistic drivers of resistance by explicitly considering factors modulating selective pressures. In this context, most pneumococcal modelling frameworks commonly overlook within-host interactions between phenotypically divergent strains (although for exceptions see [22,70] ). Neglecting these dynamics may result in missing or misinterpreting non-linear effects favouring the co-existence of antibiotic resistant and susceptible variants or the emergence of resistant strains. In this study, we expand a modelling framework previously used to explore within-host competition between pneumococcal strains as a driver of post-vaccination increases in antibiotic resistance [46] . Informed by available data for European countries, we expand the framework to a pre-vaccination era and explore how the combination of within-host ecological and immunological competition may drive observed population-level coexistence between resistant and susceptible pneumococcal strains. Other possible drivers of coexistence (e.g. age-dependent assortative mixing [68] ) were not explored in this study, although we have shown that high-level heterogeneities such as population structure and host mobility do not significantly affect model output ( Supplementary Figures S6-S7 ).
Our main results demonstrate the complex relationship between the frequency of NVT and VT pneumococci, co-infection, antibiotic resistance, and ecological and immunological competition factors. Based on data from European countries, our model indicates that pre-vaccination era antibiotic resistance levels are negatively associated with ecological competition. In our model, high ecological competition decreases antibiotic resistance levels, as this parameter governs the co-colonization competence of antibiotic resistant strains against susceptible strains, due to possible fitness costs incurred by acquisition of resistance [71] . The variation of ecological competition levels between susceptible and resistant strains across countries can be due to differences in both the bacterial and human populations. Different strains, with different genetic loads, and hence resistance fitness costs, can expand across countries due to chance and founder effects [42,72,73] . Additionally, differences in human population density and demography, environmental conditions, per capita hospital beds, and other factors can substantially affect ecological competition and its influence on co-existence of different strains [74][75][76] .

8/17
In contrast to the identified role of ecological competition in our model, we do not find that immunological competition significantly contributes to explaining the variation in observed resistance levels across European countries. In reality, immunological competition between strains is determined by a complex, largely unknown network of antigenic relationships between pneumococcal strains; with capsular polysaccharides believed to be the major determinant of strain-specific and cross-strain immunological responses [77] . As such, immunological competition levels would have to substantially differ between countries if they were a factor affecting rates of co-existence between resistant and susceptible strains. However, there is a lack of support for such substantial local differences in the antigenic determinants and overall prevalence levels of circulating antigenic groups [78,79] . Hence, the results of our model are compatible with the notion that population-level antigenic relationships and diversity are not dissimilar enough across countries to justify immunological competition as a major determinant of the co-existence of resistant and susceptible strains in the pre-vaccination era.
To conclude, we present a general yet explicit framework to model pneumococci, considering mechanisms of both between-host transmission and within-host competition as determinants of the pathogen's population dynamics. In a series of simulated scenarios, we demonstrate that ecological competition during co-colonization between resistant and susceptible pneumococcal strains can explain much of the variation of resistance frequencies observed at the country level in the pre-vaccination era. Our results demonstrate the framework's general potential to explore the drives of pneumococcal resistance dynamics. Elucidating the role of underlying competition-related components of pneumococcal dynamics improves the understanding of the mechanistic drivers for the emergence and maintenance of antibiotic resistance. Consequently, frameworks such as this can be used to project the influence of changes in vaccination or antibiotic treatment strategies on antibiotic resistance.

Conflict of interest
Authors declare no conflict of interest.

Disclaimer
The views and opinions of the authors expressed herein do not necessarily state or reflect those of ECDC. The accuracy of the authors' statistical analysis and the findings they report are not the responsibility of ECDC. ECDC is not responsible for conclusions or opinions drawn from the data provided. ECDC is not responsible for the correctness of the data and for data management, data merging and data collation after provision of the data. ECDC shall not be held liable for improper or incorrect use of the data.