Resistance against different antibiotics appears on the same bacterial strains more often than expected by chance, leading to high frequencies of multidrug resistance. There are multiple explanations for this observation, but these tend to be specific to subsets of antibiotics and/or bacterial species, whereas the trend is pervasive. Here, we consider the question in terms of strain ecology: explaining why resistance to different antibiotics is often seen on the same strain requires an understanding of the competition between strains with different resistance profiles. This work builds on models originally proposed to explain another aspect of strain competition: the stable coexistence of antibiotic sensitivity and resistance observed in a number of bacterial species. We first identify a partial structural similarity in these models: either strain or host population structure stratifies the pathogen population into evolutionarily independent sub-populations and introduces variation in the fitness effect of resistance between these sub-populations, thus creating niches for sensitivity and resistance. We then generalise this unified underlying model to multidrug resistance and show that models with this structure predict high levels of association between resistance to different drugs and high multidrug resistance frequencies. We test predictions from this model in six bacterial datasets and find them to be qualitatively consistent with observed trends. The higher than expected frequencies of multidrug resistance are often interpreted as evidence that these strains are out-competing strains with lower resistance multiplicity. Our work provides an alternative explanation that is compatible with long-term stability in resistance frequencies.
Antibiotic resistance and, in particular, multidrug resistance (MDR) are public health threats. Multidrug resistant infections are associated with poorer clinical outcomes and higher cost of treatment than other infections [1,2] and there is concern that the emergence of pan-resistant strains (pathogens resistant to all available antibiotics) will render some infections untreatable [3]. From the point of view of finding effective treatment options, multidrug resistance is particularly problematic because resistance to different antibiotics tends to be concentrated on the same strains: positive correlations between resistance to different drugs have been found in multiple species (including Streptococcus pneumoniae, Neisseria gonorrhoeae, Staphylococcus aureus, Escherichia coli, Klebsiella pneumoniae, Pseudomonas aeruginosa and Mycobacterium tuberculosis) [2]. In other words, the frequency of MDR strains is higher than we would expect from the frequencies of individual resistance determinants if these were distributed randomly in the population (‘MDR over-representation’). Understanding the causes of this MDR over-representation is important for limiting the impact of resistance. A number of possible explanations have been suggested (Table 1) [2], but the extent to which these processes contribute to the trend remains uncertain. Many of the proposed mechanisms are specific to subsets of antibiotics and/or species. The pattern of MDR over-representation, on the other hand, is pervasive: correlations have been observed between resistance to antibiotics acting through different mechanisms, and between chromosomal and mobile genetic element (MGE) associated resistance determinants [2]. Explanations for MDR over-representation must therefore be either sufficiently general or sufficiently diverse to account for this pervasiveness. In this paper, we approach the problem of explaining MDR over-representation in terms of strain ecology: explaining why resistance to different antibiotics is often seen on the same strain requires an understanding of the competition between strains with different resistance profiles. For models of such competition to be credible, they must capture observed trends in resistance dynamics whilst being ecologically plausible. Developing models that fulfil these criteria has not been trivial: sensitive and resistant strains compete for the same hosts and simple models of competition therefore predict that the fitter strain will out-compete the other (‘competitive exclusion’) [19]. However, this is rarely observed: resistance frequencies have remained intermediate over long time periods in a number of species. For example, sustained intermediate resistance frequencies are observed in Europe for various antibiotics and numerous species, including E. coli, S. aureus and S. pneumoniae (European Centre for Disease Prevention and Control Surveillance Atlas, available at https: //atlas. ecdc. europa. eu). Stable coexistence is also observed in surveillance data from multiple other locations (Centre for Disease Dynamics, Economics and Policy, available at https: //resistancemap. cddep. org/AntibioticResistance. php). For further review of evidence for stable coexistence, see references [19,20]. Recent work has explored the role of i) host population structure [21–23], ii) pathogen strain structure [20,21] and iii) within-host dynamics [24] in maintaining the coexistence of antibiotic sensitivity and resistance. In this paper, we identify a structural similarity in the first two categories of model. In these models, coexistence arises through a combination of two factors. First, the presence of groups within the host or pathogen population in which the evolutionary dynamics of resistance are approximately independent from the other groups. Second, the presence of variation in the benefit gained from resistance between these groups, so that antibiotic resistance is selected for in some groups while sensitivity is selected for in others. We show that if this variation is correlated for different antibiotics, models with this structure also predict high levels of association between resistance to different antibiotics: all resistance determinants will tend to be found where the fitness benefit gained from resistance is the greatest. The observed high frequency of multi-drug resistance is therefore in line with ecologically plausible models of coexistence, making these models a parsimonious explanation for both trends.
In this section, we discuss competitive exclusion and previously proposed coexistence mechanisms in the context of multidrug resistance. We identify a structural similarity in plausible models of coexistence [20–22] and show that, in a multidrug context, models with this structure predict MDR over-representation. The model we present captures the dynamics of a bacterial species which is mostly carried asymptomatically (e. g. E. coli, S. aureus or S. pneumoniae), so the probability of a host being exposed to antibiotics does not depend on whether the host is infected with the pathogen [13]. Key results, however, are also applicable when this is not the case (see Discussion). In this section, we explore how introducing additional complexity to the simplified model affects our predictions about association between resistance determinants and nestedness.
In this paper, we approach the question of explaining observed patterns of association between resistance to different antibiotics (‘MDR over-representation’) in terms of understanding the competition between strains with different resistance profiles. We consider recent models of the coexistence of antibiotic sensitive and antibiotic resistant strains [20–23] in which coexistence is maintained by heterogeneity in the fitness effect of resistance, arising either from heterogeneity in the rate of antibiotic consumption and/or difference in duration of carriage. We present a generalised version of these types of models, in which competition between antibiotic sensitivity and resistance is simplified to a series of independent sub-models (strata). We show that this model structure also gives rise to MDR over-representation because resistance to all antibiotics will be selected for in the strata where the fitness benefit of resistance (‘resistance proneness’) is the highest. Therefore, our results suggest that two pervasive trends in resistance dynamics, the robust coexistence of antibiotic sensitive and resistant strains and the over-representation of multidrug resistance, can both be explained by heterogeneity in the fitness effect of resistance within the host or pathogen population. We first present a simplified model for conceptual insights and then explore how additional complexity affects predicted trends. Under the strong assumption of identical antibiotic prescription patterns in all strata and no recombination, this model predicts complete linkage disequilibrium (D′ = 1) between resistance to all antibiotics. Relaxing these assumption decreases the magnitude of linkage disequilibrium, giving rise to values of D′ similar to those observed in multiple bacterial datasets. High D′ is maintained even at unrealistically high recombination rates. A lower correlation in antibiotic consumption profiles between strata leads to lower values of D′. However, the effect is gradual and the magnitude of the decrease depends on whether the strata also differ in clearance rate. Thus, even in context where patterns of prescription differ considerably between host groups, we would still expect a degree of association between resistance determinants when variation in duration of carriage contributes to variation in the fitness effect of resistance. Although the model builds on work exploring the stable coexistence of antibiotic sensitivity and resistance and coexistence is robustly observed in multiple datasets, the prediction that variation in the fitness effect of resistance leads to MDR over-representation does not require coexistence to be stable. We would expect MDR over-representation in the presence of fitness variation, even when this variation is not enough to maintain stable coexistence: for all antibiotics, the increase of resistance frequencies towards fixation would occur most rapidly in the populations with the greatest selection pressure for resistance. Under these circumstances, fitness variation would give rise to transient MDR over-representation. Our results show that when variation in the fitness effect of resistance is present and when this variation is at least partially correlated for different antibiotics, it will give rise to MDR over-representation. The extent to which this mechanism accounts for observed patterns of MDR over-representation therefore depends on the extent to which this type of fitness variation is present in pathogen populations. It is not entirely straightforward to evaluate how common variation in the fitness effect of resistance is. Wide-spread coexistence of sensitivity and resistance is not direct evidence for the pervasiveness of fitness variation because coexistence may not always arise through this mechanism. Although the majority of mechanisms proposed to date [20–23] work through fitness variation, other mechanisms are also possible [19]. In particular, recent modelling suggests that co-infection with sensitive and resistant strains gives rise to frequency-dependent selection for resistance and thus promotes coexistence [24]. However, the magnitude of this effect depends on the nature of within-host competition [24], for which there is limited data. Thus while theoretically plausible, the extent to which this mechanism contributes in practice is still unclear. It is worth noting that different coexistence mechanisms are not mutually exclusive. If coexistence arises through a combination of fitness variation and other mechanisms, we would a priori still expect the fitness variation to give rise to MDR over-representation. In the work presented here, we consider fitness variation arising from heterogeneity in antibiotic consumption between host groups (hospitals vs communities, geographic regions, age classes) and from heterogeneity in duration of carriage between host groups (age classes) and between strains (pneumococcal serotypes). This is not an exhaustive list of possible sources of heterogeneity. For example, serotype does not fully account for heritable variation in pneumococcal duration of carriage [28], suggesting other genetic traits also play a role in determining carriage duration. In light of recent results suggesting wide-spread negative frequency-dependent selection in bacterial genomes [30,31], it is not implausible to suggest these duration of carriage loci may also be under frequency-dependent selection. If so, diversity at these loci would create another source of variation in the fitness effect of resistance and hence promote coexistence and MDR over-representation. More broadly, variation in the fitness effect of resistance may arise through different mechanisms for pathogens with a different ecology than modelled in this work. For example, we have modelled a pathogen that is mostly carried asymptomatically and therefore exposed primarily to antibiotics prescribed against other infections. For pathogens where antibiotics prescribed due to infection with the pathogen itself contribute to a significant proportion of antibiotic exposure, the presence of strains differing in invasiveness would give rise to between-strain variation in antibiotic exposure and heterogeneity in the fitness effect of resistance. For bacterial species able to multiply both in hosts and in the environment, the sort of structure and heterogeneity considered in this work may also arise from differences between environmental niches. This study does not fully address the role antibiotic prescription patterns in MDR over-representation: we highlight two important remaining questions. Firstly, in the modelling framework used in this study, the distribution of drug consumption within a stratum (i. e. a well-mixed population) does not have an impact on MDR over-representation (S1 Text Section 5). In other words, the presence of host groups consuming antibiotics at different rates only promotes MDR over-representation if there is very little transmission between these host groups: individual-level correlation in antibiotic exposure is not predicted to promote multi-drug resistance. We have not explored this result in detail—it may arise because the model predicts competitive exclusion within a stratum. Secondly, in contrast to the distribution of antibiotic consumption within a stratum, our results suggest that the distribution of antibiotic consumption between strata does matter: the prediction of MDR over-representation is sensitive to how correlated prescription profiles are and the extent of this sensitivity depends on whether variation in duration of carriage is also present. Relating these theoretical results to observed correlations in the antibiotic consumption between different host groups and to the extent of assortative mixing between these groups will provide additional insights into observed patterns of MDR (e. g. why the association between some drugs is higher than others). The fitness variation model playing a role in MDR over-representation does not preclude a potential role for other mechanisms in contributing to the trend (Table 1). This study does not address the relative extent to which the different possible mechanisms contribute to MDR over-representation. This is for two reasons. Firstly, it is unclear what the patterns of MDR predicted by alternative mechanisms of MDR over-representation are. Secondly, we do not have a full understanding of which host and pathogen characteristics are relevant in defining the strata so it is difficult to directly address whether these traits are predictors of MDR. One alternative strategy for establishing the extent to which the fitness variation model contributes to MDR over-representation would be to assess patterns of association between resistance determinants in a single strain circulating in a well-mixed host population (i. e. a single stratum). The fitness variation model predicts no MDR over-representation (as defined by D′ > 0) under these circumstances. Therefore, if linkage disequilibrium is observed under these conditions, this would indicate that fitness variation is not the only mechanism of MDR over-representation. Furthermore, the magnitude of linkage disequilibrium could inform the relative contribution of the fitness variation mechanism: observing similar levels of linkage disequilibrium within strata and within the whole population would suggest the fitness variation is not a necessary mechanism for generating MDR over-representation. From a public health perspective, the fitness variation model makes two concerning predictions. Firstly, we predict frequencies of pan-resistance will be high: in a perfectly nested set of resistance profiles, the frequency of pan-resistance is equal to the frequency of the rarest resistance. As a consequence, we would expect resistance arising in response to adoption of new antibiotics or increased usage of existing antibiotics to appear on already multidrug resistant lineages—an observation which has been made for the emergence of ciprofloxacin resistance in N. gonorrhoeae in the United States [32]. Secondly, our analysis has implications for the effectiveness of potential interventions against MDR. The variation in the fitness effect of resistance to different antibiotics need not be perfectly correlated for it to promote MDR over-representation. If the variation in fitness effect is maintained by multiple factors (e. g. differential antibiotic consumption between populations and variation in clearance rates), removing one of these factors (e. g. changing patterns of prescription so that consumption of different antibiotics is no longer correlated between host groups) may have limited impact on MDR over-representation. The fitness variation model provides an explanation for MDR over-representation that is consistent with long term stability in resistance frequencies. This is relevant when considering temporal trends in resistance frequencies and predicting the future burden of resistance: other explanations for MDR over-representation (e. g. cost epistasis, correlated antibiotic exposure at the individual level—see Table 1) often require MDR strains to have an overall fitness advantage over strains with lower resistance multiplicity. This would imply that the higher than expected frequency of MDR is evidence for MDR strains out-competing other strains and thus suggest that MDR strains will eventually take over. Conversely, in the model we present, MDR strains are not out-competing other strains: all resistance frequencies are at equilibrium and MDR over-representation arises from the distribution of resistance determinants. It is worth noting, however, that even in the context of the fitness variation model, on a very long time-scale, we might expect the frequency of resistance to rise if bacteria are able to evolve resistance mechanisms that carry a lower fitness cost. We show that previously proposed models in which coexistence of antibiotic sensitivity and resistance is maintained by heterogeneity in the fitness effect of resistance also predict high frequencies of multidrug resistance. The pervasive trends of coexistence and MDR over-representation can therefore be considered, at least partially, facets of the same phenomenon. We do not propose that the model we present fully explains observed patterns of association between resistance determinants. However, this effect should be considered when evaluating the role of antibiotic-specific MDR promoting mechanisms. From a public health point of view, the model we present is concerning because it predicts high frequencies of pan-resistance. On the other hand, heterogeneity in the fitness effect of resistance as an explanation for MDR over-representation allows reconciling this trend with long term stability in resistance frequencies.
The Maela pneumococcal dataset [33], collected from a refugee camp on the border of Thailand and Myanmar from 2007 to 2010, consisted of 2244 episodes of carriage, with associated antibiograms and carriage durations. Data were obtained from, and durations of carriage calculated by, Lees et al. [28] (S1 File). Data on antibiotic sensitivity was provided for ceftriaxone, chloramphenicol clindamycin, erythromycin, penicillin, co-trimoxazole (trimethoprim/sulfamethoxazole) and tetracycline. Ceftriaxone was excluded from the analysis because data was missing for a large proportion of isolates (44%). The Massachusetts pneumococcal dataset, collected as part of the SPARC (Streptococcus pneumoniae Antimicrobial Resistance in Children) project [34], was obtained from Croucher et al. (2013) [35] (data available from Croucher et al [35]). Croucher et al. reported minimum inhibitory concentrations (MICs) for penicillin, ceftriaxone, trimethprim, erithromycin, tetracycline and chloramphenicol. Tetracycline and chloramphenicol were excluded from the analysis because data was missing for a large proportion of isolates (47% and 67% respectively). Non-sensitivity was defined in accordance to pre-2008 Clinical and Laboratory Standards Institute breakpoints [36]. For both datasets, ‘resistance’ as used throughout the paper refers to non-sensitivity. The four hospital datasets were obtained from Chang et al. [2] (S2 File). All data were analysed anonymously. If the frequency of resistance to antibiotic a is pa and the frequency of resistance to antibiotic b is pb, the coefficient of linkage disequilibrium between resistance to antibiotics a and b is Dab = pab − papb, where pab is the frequency of resistance to both a and b. The normalised coefficient D a b ′ is given by: D a b ′ = D a b m i n (p a p b, (1 - p a) (1 - p b) ) if Dab < 0 and D a b ′ = D a b m i n (p a (1 - p b), (1 - p a) p b) if Dab > 0. In general the sign of D′ is arbitrary because it depends on which alleles are chosen for the calculation. We consistently calculate D′ using the frequency of resistance: positive D′ therefore means resistance to one antibiotic is associated with resistance to the other, while negative D′ means association between sensitivity and resistance. All described models were implemented in Wolfram Mathematica (version 11. 2. 0. 0). Modelling results are numerical solutions at t = 100000 (equilibrium is reached considerably earlier, see Fig H in S1 Text). For computing D′, numerical results for strain frequencies have been rounded to the nearest 10−10 to ensure strain frequencies for absent strains are zero (as opposed to zero within numerical error). The code is provided as a supporting file. To test the effect of relaxing the assumption that the pathogen dynamics can be divided into non-interacting sub-models, we include three additional models. First, we model the dynamics of resistance to three antibiotics (i. e. eight possible resistance profiles) spreading in a host population consisting of five host groups. The antibiotics make up different proportions of total antibiotic consumption (20,35 and 45% of total antibiotic consumption rate τ). The pathogen experiences a different clearance rate within each host class p (μp). In addition, sub-strain with resistance profile g experiences clearance from antibiotic exposure at rate τg which depends on its resistance status: τg = τ (ia0. 20 + ib0. 35 + ic0. 45), where ia = 1 if g is sensitive to antibiotic a and 0 otherwise. Resistance to each antibiotic decreases transmission rate by a factor of c. Uninfected hosts of class p (Up) are therefore infected at rate c n g β [ (1 - m) I g, p + m 4 ∑ x ∈ P ′ I g, x ], where ng is the number of antibiotics strain g is resistant to, m is a parameter that sets the extent of mixing between the classes and P′ is the set of population classes excluding p. The dynamics of strain g within population p are thus described by: d I g, p d t = c n g β [ (1 - m) I g, p + m 4 ∑ x ∈ P ′ I g, x ] U p - (τ g + μ p) I g, p (10) Second, we model the dynamics of resistance to three antibiotics in a single host population in pathogen with five strains differing in clearance rate (i. e. eight resistance profiles and five strains, giving a total of 40 possible sub-strains) with recombination at the duration of carriage locus. Strain i is cleared at rate μi and, as above, sub-strains with resistance profile g experience clearance from antibiotic exposure at rate τg which depends on its resistance status: τg = τ (ia0. 20 + ib0. 35 + ic0. 45). Resistance to each antibiotic decreases transmission rate by a factor of c. Balancing selection is modelled similarly to Lehtinen et al. [20], by scaling transmission rate of strain i by a factor ψi which depends on the strain’s prevalence: ψ i = (1 - [ ∑ x I x, i 1 - U - 1 5 ]) k, where k is a parameter setting the strength of balancing selection and U is the uninfected host class. Recombination at the duration of carriage locus is modelled by allowing hosts infected with strain i with resistance profile g to transmit strain j with resistance profile g at a rate r∑x Ix, j. Recombination therefore decreases the transmission of strain i with resistance profile g by ρg, i = rIg, i ∑x ∑y Ix, y and increases it by κg, i = r∑y ∑x Ig, y Ix, i. Note that the recombination rate parameter r captures the probability of co-infection, the probability of recombination occurring and the probability of transmitting the recombinant sub-strain. The dynamics of strain i with resistance profile g are described by: d I g, i d t = c n g ψ i β [ I g, i - ρ g, i + κ g, i ] U - (τ g + μ i) I g, i (11) The third model is the same as the one above, with the exception that recombination occurs at the resistance loci instead of the duration of carriage locus. It is therefore described by Eq (11), but the expressions for ρ and κ are different. We define resistance profile g a ′ as a resistance profile otherwise identical to g, but with the other allele at locus a (i. e. if g is sensitive to antibiotic a, g a ′ is resistant), Ng, a as the set of resistance profiles with the same allele at locus a as profile g and N g, a ′ as the set of resistance profiles with the different allele at locus a than profile g. Hosts infected with strain i with resistance profile g transmit a strain i with a resistance profile g a ′ at rate r ∑ j ∑ x ∈ N g, a ′ I x, j. Recombination can occur at any of the three resistance loci (we assume recombination rates are low enough to ignore the possibility of recombination occurring at multiple loci at the same time). Recombination therefore decreases the transmission of strain i with resistance profile g by ρg, i = 3rIg, i ∑x ∑y Ix, y and increases it by κ g, i = r (I g a ′, i ∑ y ∑ x ∈ N g, a I x, y + I g b ′, i ∑ y ∑ x ∈ N g, b I x, y + I g c ′, i ∑ y ∑ x ∈ N g, c I x, y). The parameter values for the results presented in Fig 3 are: c = 0. 95, β = 2, {μ1, ‥, μ5} = {1. 2,1. , 0. 8,0. 6,0. 4}, τ = 0. 12 and k = 5. To test the effect of relaxing the assumption that all host groups consume different types of antibiotics in identical proportions, we model the dynamics of resistance to two antibiotics in a population consisting of ten host groups. The dynamics within each host group are represented by Eq 3, with parameter values cβ = 0. 95 and cμ = 1 for both antibiotics, β = 2, and, unless otherwise stated μ = 1. There is no transmission between these host groups. For both drugs, five of these host groups consume the antibiotic at a rate which selects for resistance when μ = 1 (τhigh = 0. 075), and five at a rate which selects for sensitivity when μ = 1 (τlow = 0. 025). There are therefore six different ways in which the consumption rates of the two antibiotics can be combined: all populations consuming the first drug at a high rate also consume the second drug at high rate (Spearman’s rho: 1); four out of the five populations consuming the first drug at a high rate consume the second drug at a high rate (Spearman’s rho: 0. 6); etc. We run a simulation for each of these six possible configurations. To test the effect of additional variation in the resistance proneness of strata, we introduce variation in the clearance rate of these populations: the five host groups consuming the first drug at the high rate now have different clearance rates (evenly spaced between a maximum and minimum clearance rate), and similarly for the five host groups consuming the first drug at the low rate. We run a simulation for each of these possible ways the consumption of the second drug can be distributed among these host groups.
Antibiotic resistance is a serious public health concern, yet the ecology and evolution of drug resistance are not fully understood. This impacts our ability to design effective interventions to combat resistance. From a public health point of view, multidrug resistance is particularly problematic because resistance to different antibiotics is often seen on the same bacterial strains, which leads to high frequencies of multidrug resistance and limits treatment options. This work seeks to explain this trend in terms of strain ecology and the competition between strains with different resistance profiles. Building on recent work exploring why resistant bacteria are not out-competing sensitive bacteria, we show that models originally proposed to explain this observation also predict high multidrug resistance frequencies. These models are therefore a unifying explanation for two pervasive trends in resistance dynamics. In terms of public health, the implication of our results is that new resistances are likeliest to be found on already multidrug resistant strains and that changing patterns of prescription may not be enough to combat multidrug resistance.
antimicrobials
medicine and health sciences
pathology and laboratory medicine
pneumococcus
pathogens
drugs
microbiology
epistasis
antibiotic resistance
linkage disequilibrium
antibiotics
pharmacology
bacteria
bacterial pathogens
antimicrobial resistance
fitness epistasis
medical microbiology
streptococcus
microbial pathogens
natural selection
heredity
genetics
microbial control
biology and life sciences
evolutionary biology
evolutionary processes
organisms
On the evolutionary ecology of multidrug resistance in bacteria