Mechanisms that maintain coexistence of antibiotic sensitivity and resistance also promote high frequencies of multidrug resistance

Antibiotic resistance, and in particular, multidrug resistance are public health concerns. Yet, there has been little theoretical work on the evolutionary dynamics of multidrug resistance (MDR). Here, we present a generic model of MDR inspired by two pervasive trends in resistance dynamics. The first trend is the robust coexistence of antibiotic sensitivity and resistance in multiple bacterial species and for numerous antibiotics. The second is that resistance to different antibiotics tends to be concentrated on the same strains, giving rise to high MDR frequencies. We argue that these two observations are linked: mechanisms that maintain coexistence also promote high MDR frequencies. This argument is based on the recognition that, in many of the most plausible models of coexistence, the coexistence-maintaining mechanism is fundamentally similar: either strain or host population structure stratifies the pathogen population into sub-populations and introduces variation in the fitness effect of resistance between these sub-populations. We show that this model structure also gives rise to high MDR frequencies, because resistance against all antibiotics is concentrated in the sub-populations where the fitness advantage gained from resistance is high. We test predictions from this model on two pneumococcal datasets and find predicted trends are qualitatively consistent with those observed in data. This model provides a parsimonious explanation for the pervasiveness of high MDR frequencies and allows us to reconcile this trend with observed long-term stability in the prevalence of resistance.


Introduction
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 [18,5] and there is concern that the emergence of pan-resistant strains (pathogens resistant to all available antibiotics) will render some infections untreatable. For example, there are multiple reports of resistance in Neisseria gonorrhoeae leading to treatment failure with cefixime, the last remaining monotherapy option [32]. Emergence of pan-resistance has also been reported in other Gram-negative bacteria [15].
From the point of view of finding effective treatment options, multidrug resistance is particularly problematic because resistance to different antibiotics tends to be concentrated in the same strains: 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 [5]. The evolutionary processes that drive this trend of 'MDR over-representation' are not fully understood. Previously suggested explanations (discussed below) focus on specific mechanisms that might lead to different resistance determinants being found in the same strains, but do not consider the problem more broadly in terms of competition between strains with different resistance profiles.
Indeed, there is little theoretical work on the ecology of MDR strains in general. Generic models of competition between sensitivity and resistance generally focus on a single antibiotic [2,8]. Models of multidrug resistance, on the other hand, have either only included competition between a subset of the possible resistance profiles [7,26] or sought to answer specific, rather than general, questions about the behaviour of MDR strains (e.g. the impact of changes in antibiotic consumption rates) [34,17].
In this paper, we develop a generic model of multidrug resistance to explain the over-representation of MDR strains. Our approach builds on work exploring the mechanisms maintaining stable coexistence of antibiotic sensitivity and resistance in Streptococcus pneumoniae [8,29,6]. We identify structural similarities in many of the most plausible coexistence-maintaining mechanisms and show that these coexistence-maintaining mechanisms also give rise to MDR over-representation. More specifically, we argue that mechanisms that introduce variation in the fitness advantage gained from resistance within the pathogen population promote coexistence by creating niches for resistance and sensitivity (see below). Such mechanisms also promote MDR over-representation because all resistance determinants will tend to be found where the fitness advantage gained from resistance is the greatest. We therefore suggest that explaining stable coexistence and persistent MDR over-representation are in fact facets of the same problem. This model predicts patterns of resistance that are qualitatively consistent with those observed in S. pneumoniae datasets.
Before presenting the model, we will review i) evidence in support of, and ii) potential explanations for both stable coexistence and MDR over-representation.

Coexistence of sensitivity and resistance
Stable coexistence of antibiotic sensitivity and resistance is common. Although there are specific examples of resistance reaching fixation (e.g. penicillin resistance in Staphylococcus aureus [30]), 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 Escherichia coli, S. aureus and S. pneumoniae (European Centre for Disease Prevention and Control Surveillance Atlas, available at https://atlas.ecdc.europa.eu). Thus, while coexistence of sensitive and resistant strains is not universal, it is very common across species and antibiotic classes.
In theory, if the dynamics of resistance were slow enough, apparently stable coexistence could arise from a slow increase of resistance towards fixation. This interpretation is implausible, however, because there is evidence for rapid changes in the prevalence of resistance in response to changes in antibiotic consumption. For example, the frequency of antibiotic resistance fluctuates with seasonal variation in antibiotic prescriptions in E. coli and S. aureus [37] and S. pneumoniae [13,3]. There is also evidence to suggest antibiotic resistance in N. gonorrhoeae is higher in the winter months, potentially relating to higher antibiotic consumption [23].
Stable coexistence of sensitivity and resistance is unexpected: simple models of competition between sensitive and resistant strains predict that the fitter strain would out-compete the other ('competitive exclusion') [8]. In species for which carriage is generally symptomatic, coexistence might be maintained by rates of drug prescription being dependent on the frequency of resistance (e.g. if prescriptions are discontinued once resistance reaches a certain threshold) [2]. However, for many species, drug consumption is unlikely to depend on the frequency of resistance: bacteria which are mostly carried asymptomatically (e.g. E. coli, S. aureus, S. pneumoniae) are primarily exposed to antibiotics prescribed against other infections.
There are a number of possible mechanisms which may maintain coexistence of sensitivity and resistance. We have previously suggested that bacterial duration of carriage affects the fitness advantage gained from resistance [29]. As a consequence, balancing selection maintaining variation in a pathogen population's clearance rate (e.g. serotype-specific acquired immunity in the pneumococcus) can maintain coexistence of sensitivity and resistance [29]: the strains differing in duration of carriage can be thought of as niches for sensitivity and resistance. Assortatively mixing sub-groups within the host population promote coexistence if the sub-groups differ in antibiotic consumption or clearance rates [6,26]. The structure of these mechanisms is very similar: heterogeneity in the host or pathogen population maintains variation in the fitness advantage of resistance, allowing sensitivity and resistance to coexist by occupying different niches. It is unclear whether coexistence can be fully explained by these mechanisms: a model incorporating all of the effects above did not fully replicate the extent of coexistence in S. pneumoniae [6]. The unexplained coexistence could arise from unmodelled sources of heterogeneity, or some as yet unidentified coexistence-promoting mechanisms. Nevertheless, heterogeneity in the fitness advantage gained from resistance is likely to play a major role in maintaining coexistence of sensitivity and resistance.

MDR over-representation
The frequency of multidrug resistance is higher than expected from the frequencies of resistance against individual antibiotics: isolates resistant to one antibiotic are more likely to also be resistant to other antibiotics. Positive correlations between resistance to different drugs have been found in multiple species, including S. pneumoniae, N. gonorrhoeae, S. aureus, Enterobacteriaceae and Mycobacterium tuberculosis [5]. This pattern is pervasive: these correlations have been observed between resistance to antibiotics acting through different mechanisms, and between chromosomal and mobile genetic element (MGE) associated resistance determinants. Explanations for MDR over-representation must therefore be either sufficiently general or sufficiently diverse to account for this pervasiveness.
The most straightforward explanation for MDR over-representation is that one resistance mechanism may confer resistance against multiple drugs. This is particularly relevant for antibiotics of the same class (e.g. β-lactamases and some penicillin-binding protein mutations conferring resistance against multiple β-lactams). However, shared mechanisms also exist for drugs of different classes: there are examples of efflux pumps acting on multiple drugs in numerous species [36] and evidence for clinical relevance of efflux pumps in multidrug resistance [35]. However, efflux pumps by themselves are not a sufficient explanation for MDR over-representation between drugs of different classes because excess MDR is also present where efflux pumps are not thought to be a major mechanism of resistance (e.g. between β-lactams and other classes of antibiotics in S. pneumoniae [5]).
If MDR over-representation is not simply due to multi-resistance mechanisms, it must also arise from multiple resistance genes accumulating in a single cell. Genetic linkage between resistance determinants has been suggested as a potential mechanism promoting this accumulation [5], particularly for resistance genes on mobile genetic elements (MGEs) [33]. Here, it is helpful to distinguish between resistance mechanisms where a particular allele of a gene confers resistance (e.g. target modification) and those where resistance is associated with the presence of a resistance gene (e.g. enzymes that break down the drug). For the former category, genetic linkage is not a plausible mechanism for generating MDR over-representation: a priori, linkage is no more likely to promote association between resistance determinants than association between resistance to one antibiotic and sensitivity to another.
Linkage seems like a more promising explanation for association between the second type of resistance determinants, particularly if these determinants are found on mobile elements. When resistance is associated with presence of a particular gene, absence of the gene from a particular MGE does not necessarily imply sensitivity to the antibiotic (the gene may be present on another element). As a consequence, spread of the MGE will spread resistance for resistance genes present on the element, but not spread sensitivity when resistance genes are absent. Linkage would thus asymmetrically favour association between resistance determinants. Even in this case, however, we would expect the effect of linkage to be transient: recombination and mutation eventually eliminate resistance determinants which do not confer a fitness advantage, even from mobile genetic elements. For example, in the PMEN1 pneumococcal lineage, there is evidence for loss of aminoglycoside resistance from the Tn916 transposon which encodes macrolide and sometimes tetracycline resistance [11].
In theory, accumulation of multiple resistance determinants in a cell could also arise from correlations in drug exposure -either because of use of combination therapy or through sequential drug exposure if the first choice of drug fails. In practice, however, it is unclear how significant a role these mechanisms play in MDR over-representation. Guidelines recommend single-agent antibiotic therapy over combination therapy under most circumstances [27], suggesting combination therapy is unlikely to represent a significant proportion of prescriptions for the majority of antibiotics. In the United Kingdom for example, monotherapy accounts for 98% of primary care prescriptions [12]. It is unclear whether treatment failure rates (below 20% for the most commonly prescribed antibiotics in UK primary care [12]) and patterns of subsequent drug prescription are enough to drive selection for multidrug resistance.
Finally, accumulation of multiple resistance determinants could also arise from 'cost epistasis' between resistance determinants (i.e. lower than expected fitness cost when both determinants are present). Indeed, there is evidence of cost epistasis between resistance determinants occurring in laboratory competition experiments in a number of species, including between streptomycin and rifampicin resistance in Pseudomonas aeruginosa [41] and in E. coli [1], and streptomycin and quinolone resistance in E. coli [38]. There is also in vitro evidence of cost epistasis between rifampicin and ofloxacin resistance in Mycobacterium smegmatis [4]. However, the extent to which epistasis plays a role in vivo remains unclear [43]. In particular, we would not, a priori, expect to observe cost epistasis between resistance determinants operating through entirely different mechanisms.
In summary, the mechanisms outlined above are likely to contribute to MDR over-representation for subsets of resistance determinants, but are not a satisfactory explanation for the trend overall: even in combination with each other, these mechanisms do not account for the pervasiveness of MDR over-representation. Figure 1: Competition between a sensitive and resistant sub-strain for a single strain circulating in a homogeneous host population (the 'triangle model'). U , I s and I r represent uninfected hosts, hosts infected with the sensitive sub-strain and hosts infected with the resistant sub-strain, respectively. β r and β s are transmission rates for resistant and sensitive sub-strains; µ r and µ s are clearance rates; τ is the population antibiotic consumption rate (we assume immediate clearance following antibiotic exposure). Resistance is associated with a fitness cost affecting the transmission (β r = β s c β , where c β < 1) and/or the clearance (µ r = µs cµ , where c µ < 1).

'Fitness variation' model of multidrug resistance
We present a generic model of competition between antibiotic sensitivity and resistance, which aims to account for both the stable coexistence of antibiotic sensitivity and resistance and the over-representation of MDR. The model captures the dynamics of a bacterial species which is mostly carried asymptomatically (e.g. E. coli, S. aureus or S. pneumoniae), where the probability of a host being exposed to antibiotics does not depend on whether the host is infected with the pathogen. Key results, however, are also applicable when this is not the case (see Discussion). The model, as presented here, makes very strong assumptions -we will discuss the effect of relaxing these assumptions in later sections of the paper. It is also worth noting that our aim is to explain the distribution of resistance determinants in the pathogen population rather than predicting their frequencies -the simplifications we make reflect this (see Supporting Information).
We start by considering a simple ('triangle') model of competition between resistant and sensitive variants of an otherwise genetically homogeneous pathogen (one strain) in a homogeneous host population ( Figure 1). To avoid ambiguity later in the paper, we will refer to the sensitive and resistant variants as 'sub-strains'. Uninfected hosts (U ) become infected with the resistant (I r ) or the sensitive (I s ) sub-strain at rate β r and β s ; infections are cleared at rate µ r and µ s ; the sensitive sub-strain experiences an additional clearance rate τ corresponding to the population antibiotic consumption rate (we assume immediate clearance following antibiotic exposure); and resistance is associated with a fitness cost affecting transmission (β r = β s c β , where c β < 1) and/or clearance (µ r = µs cµ , where c µ < 1). The dynamics of this model are described by: In this model, the fitness effect of resistance, defined as the ratio of the basic reproductive numbers (i.e. R 0 , the average number of new infections an infected host gives rise to in a fully susceptible population) of the resistant and sensitive sub-strains, depends only on population antibiotic consumption, the fitness cost of resistance, and the clearance rate of the strain: R0r and with c = c µ c β and µ = µ s . The model predicts competitive exclusion: when the fitness effect of resistance is above 1, the resistant sub-strain out-competes the sensitive sub-strain, and vice-versa when the fitness effect is below 1.
We generalise this model to heterogeneous host and strain populations by dividing the pathogen population into sub-populations (strata) within which competition between sensitivity and resistance is independent of other sub-populations. For competition within the strata to be considered independent, the frequency of resistance within one stratum must not affect the frequency of resistance within other strata. When this is the case, we can approximate the dynamics of resistance in the population as a whole by modelling each stratum as a separate triangle model.
In the case of a single strain circulating in a structured host population, the strata correspond to assortatively mixing sub-groups (e.g. age classes). We assume transmission between host sub-groups is either low enough to be ignored or high enough for the sub-groups to be modelled as a single population. In the case of multiple strains circulating in a homogeneous host population, the strata correspond to different strains (e.g. serotypes in the pneumococcus). We assume recombination between strains is low enough to be ignored. When both strain and population structure is present, each stratum corresponds to a particular strain circulating in a particular host sub-group.
While competition between sensitivity and resistance can be modelled the same way for host and strain strata, these two types of structure are not entirely identical. In the case of host population structure, the uninfected compartment is also structured: for example, if there are n age classes, there will be n uninfected compartments. In the case of strain structure, there is a single uninfected compartment. As a consequence, the strains are competing for the same hosts. If these strains differ in fitness, diversity must therefore be maintained through balancing selection, for example serotype-specific immunity in the case of pneumococcal serotypes. For the purposes of our MDR model, balancing selection maintaining strain structure effectively allows us to treat the uninfected compartment as if it were structured (see Supporting Information for further discussion of this point).
Because we aim to explain stable coexistence of sensitivity and resistance and the persistent overrepresentation of MDR, we are primarily interested in the model's steady state behaviour. We therefore only model genetic variation that is maintained by some form of balancing selection and do not consider transient genetic diversity during a selective sweep. Genetic factors like compensatory mutations may introduce between-strain variation in the fitness cost of resistance, but in the absence of evidence for balancing selection maintaining this diversity, we assume this variation will be transient and do not include it in the model.
We also assume that while host sub-groups can differ in the antibiotic consumption rate, the relative rate at which different antibiotics are prescribed remains the same in all sub-groups (antibiotic a accounts for proportion γ a of total antibiotic consumption, with a γ a = 1).
Under these assumptions, resistance to antibiotic a out-competes sensitivity in host sub-group p and strain s when: As before, there is no coexistence within the strata: below this threshold, sensitivity out-competes resistance.
We can separate Inequality 2 into stratum (i.e host sub-group and pathogen strain) and antibiotic related effects (left-hand and right-hand sides of Inequality 3, respectively): The ratio τp µps reflects how advantageous resistance is within a stratum ('resistance proneness' of a stratum), while the ratio 1 γa ( 1 ca − 1) reflects how advantageous resistance against a particular antibiotic needs to be for it to be selected for ('resistance threshold' of a drug). The strata can therefore be ranked in order of the fitness advantage gained from resistance, based on resistance proneness.
Finally, we extend this model to a multidrug context. Assuming no cost epistasis between resistance determinants (i.e. the fitness cost of resistance against antibiotics a and b is equal to c a c b ), the ranking of the strata is the same for all antibiotics. As a consequence, if resistance to an antibiotic with a high resistance threshold out-competes sensitivity in a stratum, resistance to antibiotics with lower thresholds will also be selected for (Inequality 3). Thus, rare resistance determinants will only be found in pathogens carrying more common resistance determinants, giving rise to MDR over-representation ( Figure 2). In other words, the resistance determinants will be nested within each other.

Duration of carriage is predictor of resistance multiplicity
The model presented above suggests that duration of carriage and antibiotic consumption rate within strata are predictive of resistance multiplicity. Fully testing this prediction is challenging, because we do not have a full understanding of which host and pathogen characteristics are relevant in defining the strata. To partially test the prediction of the model, we confirm that duration of carriage is indeed predictive of resistance multiplicity (Kendall rank correlation (τ ) 0.14 95% CI 0.11-0.16, see Figure 3), using a dataset of pneumococcal carriage episodes and associated durations of carriage [28] ('Maela dataset', n = 2244, see Methods). The direction of causality for this association is not entirely clear: the model we present suggests longer duration of carriage  Figure 2: Illustration of the fitness variation mode. Each triangle model represent competition between sensitivity and resistance within a strain circulating in a sub-population of a heterogeneous host and pathogen population. The height of the bars represent resistance proneness within each stratum ( τp µps ) and the arrows represent the resistance thresholds ( 1 γa ( 1 ca − 1)) for different drugs. If resistance to an antibiotic with a high threshold out-competes sensitivity given stratum, resistance to antibiotics with a lower threshold will also be selected for within this stratum. Resistance determinants will therefore be nested within each other.
leads to resistance, but resistance would also be expected to lead to longer duration of carriage through decreased clearance from antibiotic exposure. The association between duration of carriage and resistance multiplicity is also present at the serotype level (τ 0.27 95% CI 0.05-0.46, see Supporting Information), where differences in duration of carriage are thought to arise from the properties of serotype capsules [42] (rather than differences in antibiotic resistance). This suggests that resistance leading to longer duration of carriage is unlikely to fully account for the observed association.

Observed patterns of association between resistance determinants
The fitness advantage model also makes predictions about the patterns of association between resistance to different antibiotics. We test whether these predictions are consistent with pneumococcal data from two carriage studies: the Maela dataset discussed above and a dataset collected in Massachusetts (n = 603, with data on sensitivity to four antibiotics, see Methods for further details).

Pairwise association
The fitness advantage model predicts that rare resistance determinants will only be found in the presence of more common resistance determinants. Based on the model, we would therefore expect complete linkage disequilibrium (LD) between resistances: D = 1 for all pairs of resistances, where D is the normalised coefficient of linkage disequilibrium (see Methods). We observe high D between all pairs of antibiotics with the exception of chloramphenicol in the Maela dataset (Figure 4), for which estimates ranged from high positive LD to high negative LD. The low frequency of chloramphenicol resistance (6%) is a potential explanation for this variability: D is biased towards extreme values when allele frequencies are low. Furthermore, the fitness advantage variation model ignores the possibility that disadvantageous allele combinations may be reintroduced through recombination and mutation. If the rate at which alleles are re-introduced is independent of the allele frequency, we would expect the impact on LD to be greater when allele frequencies are low. The re-introduction of disadvantageous alleles may therefore contribute to the variability of observed LD for chloramphenicol resistance.
Complete LD (D = 1) was only observed between penicillin and ceftriaxone resistance (Massachusetts dataset). This likely reflects a shared mechanism, rather than the effect of variation in the fitness advantage of resistance as penicillin and ceftriaxone both target penicillin binding proteins.
In summary, we observe consistently high, but not complete, linkage disequilibrium between resistance determinants.

Nestedness
Because rarer resistance determinants are expected to appear only in the presence of more common ones, the fitness advantage model predicts that we will only observe a subset of the possible antibiogram profiles: all antibiograms of multiplicity m (i.e. with resistance against m antibiotics) will have the m most common resistance determinants. We call antibiograms which fulfil this criteria nested (see Figure 2). Nestedness is related, but not identical, to D : a set of resistance determinants in complete positive LD (D = 1 for all pairs) would also be perfectly nested, but this equivalence does not hold when LD is not complete. In other words, knowing D between all pairs of antibiotics does not fully inform the degree of nestedness of a set of antibiograms (see below). We quantify the degree of nestedness in the datasets as the proportion of observed antibiograms with a nested structure (while more complex metrics of nestedness have been proposed, these are less interpretablesee Methods). In the Maela dataset, 71% of isolates have a nested structure, 76% excluding chloramphenicol; in the Massachusetts dataset, 82% of isolates have a nested structure ( Figure 5). Similar results have also been reported in S. aureus: 78% of antibiograms in a Massachusetts collection spanning 14 years were nested (data for resistance to four antibiotics) [25].
We are interested in whether the fitness variation model is a better predictor of the overall structure of the antibiograms than the alternative MDR over-representation mechanisms. Ideally, therefore, we would compare the observed nestedness to that expected from these alternative mechanisms (shared resistance mechanisms, linkage, correlated drug exposure and epistasis). This is not possible, however, because it is unclear how much nestedness we would expect under these MDR mechanisms. We also considered comparing the observed nestedness against a null model. However, establishing an appropriate null hypothesis was not straightforward and although we considered a number of options, none of these were informative in our context (see Supporting Information).

Sensitivity of predictions to model assumptions
The fitness variation model makes a number of strong assumptions. We test the extent to which relaxing these assumptions affects the model's prediction about MDR over-representation and association between resistance determinants.
The model stratifies the pathogen and host populations into separate sub-populations within which competition between resistance and sensitivity is represented by the simple triangle model. This requires us to make two broad assumptions: first that it is possible to define non-mixing strata which can be considered To test the effect of relaxing the first assumption, we construct a model with three different antibiotics and five host sub-groups differing in clearance rate (see Methods). Allowing a small amount of transmission between these sub-groups gives rise to within-strata coexistence, but does not affect predictions about nestedness: mixing between model strata does not introduce a mechanism that would favour non-nested antibiograms ( Figure 6). Increasing the amount of transmission between host groups causes them to behave as a single population, thus abolishing the stratified population structure (see Supporting Information). Because of the analogy between strain and host population structure, the same results are likely to hold when mixing between strata arises from recombination at the duration of carriage locus.
These results have two implications. Firstly they suggest that modelling competition between sensitivity and resistance as non-mixing homogeneous strata is a robust approximation because the population appears to behave either as non-mixing sub-groups or fully mixing. Secondly, even when mixing introduces within-strata coexistence, predictions about nestedness remain valid. Therefore, the fitness variation model remains an explanation for MDR over-representation even when assumptions about mixing between, or structure within, strata are relaxed.
The effect of relaxing the second assumption (that the triangle model adequately captures competition within strata) depends on the way in which competition deviates from the triangle model. As discussed above, deviation from the triangle model arising from a small amount of mixing between strata do not affect predictions about nestedness. On the other hand, we would expect the re-introduction of disadvantageous resistance alleles through mutation or recombination to allow low frequencies of non-nested antibiograms to persist in the population. As discussed in the context of chloramphenicol resistance in the Maela dataset, the impact on nestedness may be particularly pronounced when allele frequencies are low. It is also possible that some other coexistence-maintaining mechanisms may introduce selection favouring non-nested antibiograms, but without plausible candidate mechanisms, we cannot test this.
In addition to the simplifications made in constructing the stratified model, we make assumptions that allow us to separate the variation in the fitness effect of resistance into strata-related (i.e. pathogen and host)  and antibiotic-related effects. This separability means that the ranking of the strata by fitness advantage of resistance will be the same for all antibiotics. The separation relies on two assumptions: i) the fitness cost of a particular resistance being the same in all strata (i.e. no variation in the fitness cost of resistance between strains) and ii) different types of antibiotics being consumed in the same proportions in all strata (i.e. variation in the rate at which host sub-groups consume antibiotics, but not in the types of antibiotics consumed by these groups). Violating these assumption may give rise to non-identical strata rankings and thus lead to deviations from nestedness. It is possible both of these assumptions represent oversimplifications of resistance dynamics. In the case of the first assumption, there is no direct evidence for stable variation in the fitness cost of resistance. However, the processes determining fitness cost are not fully understood and the fitness cost of resistance mutations is thought to depend on both genetic background and environment [31]. It is therefore difficult to rule out that stable between-strain variation might exist. Similarly, the assumption about antibiotic prescription rates is generally plausible, but might not hold in all contexts. For example, overall antibiotic prescription rates differ between regions of the United States, but relative prescription rates of different antibiotic classes are highly correlated (mean correlation in relative prescription rates 0.99, based on data reported by Drekonja et al. [14]). However, correlations between relative prescription rates may not be so high between all host sub-groups (e.g. children and adults).
To test the extent to which between-strata differences in the relative fitness advantage gained from resistance to different antibiotics affects patterns of nestedness, we simulated strata with increasingly uncorrelated resistance proneness for different antibiotics (see Methods). The proportion of nested antibiograms decreased with decreasing correlation (Figure 7). However, the effect was gradual: imperfectly correlated strata rankings gave rise to relatively high levels of nestedness (comparable to those seen in the data). Variation in the fitness advantage of resistance is therefore a mechanism promoting MDR over-representation even when this variation is not identical for all antibiotics.

Predicted and observed patterns of association
We present a model of multidrug resistance in which competition between antibiotic sensitivity and resistance is simplified to a series of independent sub-models. Mechanisms that introduce variation in the fitness effect of resistance between these sub-models promote both coexistence of antibiotic sensitivity and resistance, and higher than expected frequencies of multidrug resistance (MDR over-representation). Based on this model, Table 1: Table of discussed evolutionary processes that could explain MDR over-representation. The two last rows represent fitness variation processes that also predict coexistence.

Process
Could explain coexistence?
Compatible with unequal resistance frequencies? Plausibility Shared resistance mechanisms

No
Requires additional mechanisms Only applicable to a subset of antibiotics Linkage between resistance genes

No
Requires additional mechanisms In absence of selection maintaining LD, associations between alleles would eventually be lost. Time-scale at which this would occur is uncertain. Correlated drug exposure of individual host

No
Requires additional mechanisms Rates of combination therapy are generally low. Unclear how significant a role correlated exposure due to treatment failure might play and which antibiotics it affects.

Cost epistasis No Requires additional mechanisms
Little evidence of this outside laboratory context. Implausible when resistance mechanisms are unrelated and not plasmid-associated. Epistasis between resistance and carriage duration

Yes
Yes Correlation between duration of carriage and resistance.
Local adaptation to host sub-group Yes Yes Correlation between antibiotic consumption and resistance for various host sub-groups (e.g. countries, age classes).
we would expect all observed antibiograms to be 'nested', with rare resistance determinants appearing only in the presence of more common ones. Consistent with this fitness variation model, we find positive linkage disequilibrium between resistance determinants and a high degree of nestedness in two pneumococcal datasets. Similar levels of nestedness have also been reported in S. aureus antibiograms [25]. Contrary to the model's prediction, however, the observed linkage disequilibrium is not complete (i.e. D = 1) nor are all (i.e. 100%) antibiograms nested. The prediction of complete LD and perfect nestedness arises because we formulate the model so that the variation in the fitness effect of resistance is identical for all antibiotics. This requires us to make strong assumptions about relative prescription rates of different antibiotics (same in all host sub-groups) and the relative fitness cost of different resistances (same for all strains). Relaxing these assumptions gives rise to imperfectly correlated variation in the fitness effect of resistance for different antibiotics, which in turn generates imperfect nestedness. The patterns we observe are therefore consistent with variation in the fitness effect of resistance that is correlated, but not identical, for different antibiotics.

Alternative explanations for MDR over-representation
We present the fitness variation model as an alternative to a model in which MDR over-representation arises from a number of different mechanisms, each affecting a subset of antibiotics (summarised in Table 1). We do not mean to suggest these mechanisms do not contribute to MDR over-representation. Indeed, mechanisms alternative to fitness variation are a plausible explanation for some of the associations we observe in our data. For example, it is likely a shared mechanism of resistance accounts for the strong association between penicillin and ceftriaxone resistance in the Massachusetts dataset. Linkage may contribute to the association between MGE-associated resistance determinants (chloramphenicol, tetracycline, erythromycin and clindamycin): as discussed in the Introduction, we expect the effect of linkage to be transient, but the time-scale at which these associations would be broken down in absence of selection is unclear.
Nevertheless, we also observe associations which do not seem attributable to any of the alternative explanations (e.g. the association between penicillin and the other resistance determinants). Furthermore, the fitness variation model is an attractive explanation for the MDR over-representation trend as a whole.
Firstly, the fitness variation model is compatible with unequal resistance frequencies. Other possible MDR over-representation mechanisms promote MDR strains, either because these have an advantage over singly resistant strains (epistasis, correlated drug exposure) or because resistant determinants are genetically linked (shared mechanisms, linkage). As a consequence, these mechanisms also predict similar resistance frequencies against different antibiotics. Unlike the fitness variation model, these explanations therefore need additional mechanisms to account for observed differences in resistance frequencies.
Secondly, the fitness variation model is a parsimonious explanation for the pervasiveness of MDR overrepresentation. Each of the alternative explanations predict MDR over-representation only for subsets of antibiotics, yet the trend is ubiquitous. This pervasiveness is compatible with the fitness variation modelwhereas it is unclear if the other possible mechanisms are enough to explain MDR over-representation for all resistance determinants.

Coexistence and MDR over-representation are linked by theory
The model we present is generic. We identify diversity in host antibiotic consumption rates and pathogen clearance rates as mechanisms which introduce variation in the fitness advantage of resistance, but do not elaborate on how this diversity arises. Antibiotic consumption rate is age-dependent [21] which, in combination with age assortative mixing, would give rise to the type of fitness variation described in the model. For the pneumococcus, clearance rate is known to depend on both serotype and host age [22], which would introduce further variation in the fitness effect of resistance. Furthermore, serotype does not fully account for heritable variation in pneumococcal duration of carriage [28], suggesting there may be also be other genetic factors giving rise to stable variation in duration of carriage. This is particularly plausible in light of recent results suggesting wide-spread negative frequency-dependent selection in bacterial genomes [9,20].
Potential sources of variation in the fitness advantage of resistance have not been extensively studied in pathogens other than S. pneumoniae. However, the wide-spread coexistence of sensitivity and resistance suggests that the presence of such variation in other bacterial species is plausible. This variation may arise from sources other than those included in the model we present. For example, if a pathogen has strains differing in invasiveness, this variation could lead to between-strain variation in antibiotic exposure and thus in the fitness effect of resistance.
The observation that variation in the fitness advantage of resistance promotes both coexistence and MDR over-representation highlights the importance of identifying coexistence-maintaining mechanisms. In S. pneumoniae, a model incorporating various sources of fitness variation does not fully explain the extent of observed coexistence [6]. Establishing whether the remaining coexistence arises from unidentified sources of fitness variation (rather than entirely different coexistence-promoting mechanisms) would clarify how applicable the fitness variation model is. Furthermore, identifying these sources of fitness variation would not only inform our understanding of individual resistance dynamics, but also provide further insights into the mechanisms of MDR over-representation.
Although the model builds on work exploring the stable coexistence of antibiotic sensitivity and resistance, the prediction that variation in the fitness advantage 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.

Public health implications
From a public health perspective, the fitness variation model makes two concerning predictions. Firstly, the nestedness prediction implies frequencies of pan-resistance will be high: in a perfectly nested set of antibiograms, the frequency of pan-resistance is equal to the frequency of the rarest resistance. As a consequence, we would expect resistance to last-line antibiotics to be found on isolates which are already resistant to all other available antibiotics -an observation which has been made for the emergence of ciprofloxacin resistance in N. gonorrhoeae in the United States [19]. This effect may be mitigated by the impact of mutation when allele frequencies are low. Indeed, we observe no pan-resistant isolates in the Maela dataset, where the frequency of the rarest resistance (chloramphenicol) is low (6%).
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) may have limited impact on MDR.
On the other hand, the fitness variation model provides an explanation for MDR over-representation that is consistent with long term stability in resistance frequencies. As discussed above, alternative explanations for MDR over-representation 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. The equilibrium frequencies are determined by the number of strata within which resistance is advantageous and thus depend on the antibiotic consumption rate, duration of carriage and fitness cost of resistance. It is worth noting, however, that even in the context of the fitness varition 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.

Conclusion
We present a conceptual model of multidrug resistance where variation in the fitness effect of resistance in a structured host and pathogen population i) maintains coexistence of resistance and sensitivity and ii) gives rise to high frequencies of multidrug resistance. Coexistence and multidrug resistance can therefore be considered different facets of the same phenomenon, highlighting the importance of identifying the mechanisms which maintain stable coexistence. While the model we present predicts high frequencies of pan-resistance, it allows us to reconciles MDR over-representation with long term stability in resistance frequencies.

Datasets
The Maela dataset [39], 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]. 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 dataset, collected as part of the SPARC (Streptococcus pneumoniae Antimicrobial Resistance in Children) project [16], were obtained from Croucher et al. (2013) [10]. 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 Clinical and Laboratory Standards Institute breakpoints [24].
For both datasets, 'resistance' as used throughout the paper refers to non-sensitivity.

Linkage disequilibrium
If the frequency of resistance to antibiotic a is p a and the frequency of resistance to antibiotic b is p b , the coefficient of linkage disequilibrium between resistance to antibiotics a and b is D ab = p ab − p a p b , where p ab is the frequency of resistance to both a and b. The normalised coefficient D ab is given by: if D ab > 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.
Confidence intervals for D were generated by bootstrapping: the set of m antibiograms was randomly sampled m times with replacement to generate a random set of the same size as the original set and D was computed on this set. We repeated this 10000 times to generate a distribution and computed the 95% confidence interval from this distribution.
We use D as a measure of linkage disequilibrium (as opposed to the other commonly used metric r 2 ) because the fitness effect model makes prediction specifically about D . The coefficient of determination, r 2 , measures the extent to which isolates that are resistant to one antibiotic are also resistant to another antibiotic while D captures the extent to which resistance to two antibiotics will be found in the same isolates, given the observed resistance frequencies. Therefore, r 2 is affected by how similar resistance frequencies are and by the distribution of the resistant determinants, while D is only affected by the later.

Measures of nestedness
We quantify nestedness in terms of the proportion of nested antibiograms. We chose this measure for two reasons. Firstly, it is invariant to inverting absences and presences (the proportion of nested antibiograms is the same if we quantify nestedness with respect to sensitivity, instead of resistance). Secondly, it is easily interpretable. There are a number of other possible metrics (reviewed in reference [40]), some of which also fulfil the first criteria. These metrics seek not only to quantify whether antibiograms are nested, but also how nested they are. However, there is no intuitive scale for degree of nestedness, which makes these metrics more difficult to interpret than simply reporting the proportion of nested antibiograms.

Effect of mixing between strata
To test the effect of relaxing the assumption that the pathogen dynamics can be divided into non-interacting sub-models, we model the dynamics of resistance to three antibiotics (i.e. eight strains in total) spreading in a host population consisting of five host classes. The antibiotics are all prescribed at the same rate τ , but incur different fitness costs. Strain g experiences a different clearance rate within each host class p (µ p ) and an additional clearance rate τ g for each antibiotic it is not resistant to (i.e. τ g = (3 − n g )τ , where n g is the number of antibiotics strain g is resistant to). Resistance to each antibiotic decreases transmission rate by a factor of c. Uninfected hosts of class p (U p ) are therefore infected at rate c ng β (1 − 4m)I p g + m x∈P I x g , where 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: For the results presented in the main text, m = 0.001.

Effect of imperfectly correlated strata rankings
To test the effect of relaxing the assumptions that lead to the ranking of the strata by fitness advantage of resistance being the same for all antibiotics (i.e. antibiotics consumed in the same proportions within all strata and the cost of resistance being the same within all strata), we created simulated resistance data similar to the Maela dataset (2244 isolates, 6 antibiotics). We randomly assigned a resistance proneness (i.e. τp µps ) to each isolate. We then set resistance thresholds for each antibiotic. Isolates with a resistance proneness above this threshold were resistant. (The resistance thresholds were chosen to give the same resistance frequencies as observed in the Maela data). We then calculated the proportion of resistant antibiograms for this simulated dataset.
To generated imperfectly correlated rankings, we replicated the resistance proneness vector for each antibiotic to generate a resistance proneness matrix and then randomly reassigned the resistance proneness of a proportion of the entries in this matrix (proportion labelled as 'noise in rankings' in Figure 7). Resistance was then assigned and nestedness calculated as above. By increasing the proportion of reassigned entries, we generate increasingly uncorrelated fitness advantage rankings.