On the evolutionary ecology of multidrug resistance in bacteria

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.


Implicit modelling of balancing selection
In the main text, we show that for the purposes of understanding association between resistance determinants, models of competition between antibiotic resistance and sensitivity with host population structure and strain structure can both be simplified to a series of independent SIS models. In the case of strain structure, this requires the presence of balancing selection to maintain the strains differing in duration of carriage. We do not model this balancing selection explicitly: we assume coexistence of the strains and focus on modelling the competition between sensitivity and resistance within each strain.
To illustrate this simplification, we consider a model with n strains in which strain coexistence is maintained by scaling the transmission rate of strain i (β i ) by g(f i ), where f i is the frequency of strain i and g() is a decreasing function (this is the model in Lehtinen et al. [1]). For any strain i, with sensitive and resistant sub-strains is and ir, clearance rate µ i , cost of resistance c β (c β ≤ 1) and c µ (c µ ≤ 1) and population antibiotic consumption τ , the competition between sensitivity and resistance is captured by: This model predicts competitive exclusion within strain: when c β c µ (1 + τ µi ) > 1, strain i will be resistant (I is = 0, I ir > 0) and when c β c µ (1 + τ µi ) < 1, strain i will be sensitive (I is > 0, I ir = 0). Importantly, these expression depend only on strain i: the competition between sensitivity and resistance within strain i is independent of the other strains in the model. For our purposes, therefore, we can model each strain separately, replacing g(f i )U by a fixed number of hosts available to strain i, giving rise to the model presented in the main text.
2 Schematics of resistance profiles

Increasing threshold
Decreasing proneness

Increasing threshold
Decreasing proneness

Increasing threshold
Decreasing proneness

Model with strong assumptions
Competitive exclusion within strata Nested resistance profiles Intergroup transmission and/or recombination at duration locus Coexistence within strata is possible Nested resistance profiles

Recombination at resistance loci
Coexistence within strata is possible Non-nested resistance profiles observed

Increasing threshold
Decreasing proneness

Imperfectly correlated resistance proneness
Competitive exclusion within strata Non-nested resistance profiles observed A B C D Figure A: Schematic illustrating the effect of relaxing various model assumption on a set of resistance profiles from a system with five strata and four antibiotics. A: set of resistance profiles under the model with strong assumptions (main text Figure 2). B: Mixing between strata (integroup transmission or recomination at the duration locus) introduces some within-strata diversity in resistance profiles, but all observed resistance profiles have a nested structure. C: If the resistance proneness of a strata is not the same for all antibiotics, non-nested resistance profiles can arise. D: Recombination at the resistance loci breaks down associations between resistances and leads to the presence of non-nested resistance profiles and coexistence of multiple resistance profiles within strata.

Effect of resistance on clearance rate in multidrug model
In Sections 2.1.3 and 2.1.4 of the main text, when considering selection for resistance against antibiotic a, we formulate the model in a way that ignores the effect of treatment with, and resistance against, other antibiotics. Resistance against other antibiotics affects clearance rate: firstly, by determining whether the strain can be cleared through exposure to the other antibiotics and, secondly, through the effect of fitness cost of resistance (if applicable). Therefore, the absence/presence of resistance against other antibiotics will affect selection for resistance to antibiotic a.
To investigate the impact this effect has on our predictions about the association between resistance determinants, we consider a system with two antibiotics and competition between four resistance profiles: sensitive to both antibiotics (SS); sensitive to antibiotic 1, resistant to antibiotic 2 (SR); resistant to antibiotic 1, sensitive to antibiotic 2 (RS); resistant to both (RR). Within a particular stratum with total antibiotic consumption τ and clearance rate µ, the dynamics of these strains are captured by: where U + I SS + I SR + I RS + I RR = 1; β is the transmission rate; c β1 ≤ 1 and c β2 ≤ 1 are the cost of resistance to antibiotic 1 and 2 respectively on transmission and c µ1 ≤ 1 and c µ2 ≤ 1 are the costs on clearance; and γ 1 and γ 2 are the proportion of total antibiotic consumption corresponding to each antibiotic. Without loss of generality, we assume antibiotic 2 is prescribed at a higher rate than antibiotic 1 (0.5 > γ 1 = 1 − γ 2 ).
This system will give rise to incomplete linkage disequilibrium (D < 1) between the resistances (and non-nested resistance profiles) if there are strata (i.e. some values of resistance proneness τ µ ) where strain RS is selected for and other strata (i.e. other values of τ µ ) where strain SR is selected for. The basic reproductive +γ2τ and the basic reproductive number of SR is +γ1τ . SR out-competes RS when R0 SR > R0 RS . Re-arranging this gives: When the fitness cost of the more commonly prescribed antibiotic (antibiotic 2) is lower (c β1 c µ1 < c β2 c µ2 , c β1 ≤ c β2 and c µ1 ≤ c µ2 ), the inequality will hold independent of the value of µ and τ (as long as u > 0 and τ ≥ 0): the left-hand side is always positive, whereas the right-hand side is always negative. Thus, resistance to the more commonly prescribed antibiotic with lower fitness cost is always more beneficial than resistance to the other drug and SR is fitter than RS in all strata.
When the fitness cost of the more commonly prescribed antibiotics is greater (c β1 c µ1 > c β2 c µ2 , c β1 ≥ c β2 and c µ1 ≥ c µ2 ), the relative fitness of SR and RS depends on the antibiotic consumption rate and clearance rate and may therefore differ between strata. However, incomplete linkage disequilibrium will only arise if a change in τ µ leads from a stable equilibrium with SR at fixation to a stable equilibrium with RS at fixation (or, equivalently, vice versa).
We perform a stability analysis to determine the parameter space in which this transition occurs. In line with the result above, stability of the SR equilibrium requires: In addition, the following conditions must also be met: These two conditions are also required for the equilibrium with RS at fixation to be stable, in addition to: Therefore, the relative treatment rate (γ 1 ) and fitness cost parameters at which incomplete linkage disequilibrium may occur are constrained by: γ1cµ1cµ2c β1 . In addition, for incomplete linkage equilibrium to actually be observed, there must be at least one stratum with resistance proneness τ µ > c β1 cµ1−c β2 cµ2 cµ1cµ2(γ2c β2 −γ1c β1 ) and one stratum with resistance proneness τ µ < c β1 cµ1−c β2 cµ2 cµ1cµ2(γ2c β2 −γ1c β1 ) . This places considerable restrictions on the parameter range at which incomplete linkage equilibrium could be observed. For example, consider a system in which strata correspond to strains differing in clearance rate, with a uniform total antibiotic consumption rate across the strata. c µ2 = 0.9, γ 1 = 0.09, c β1 = c β2 = 1 and total antibiotic consumption rate τ of one prescription per year per person, incomplete linkage equilibrium will be observed only if there is at least one strain with clearance rate µ between 0.677 and 0.743 per month and another with clearance rate between 0.668 and 0.677 per month -or, equivalently, duration of carriage between 1.35 and 1.48 months and duration of carriage between 1.48 and 1.50 months (mean duration of carriage is 1 µ ). If these conditions are not met, incomplete linkage disequilibrium, though theoretically possible, will not be observed.
SI Figure B and SI Figure C further illustrate this sensitivity to parametrisation. In particular, the closer the rate at which the two drugs are consumed, the more restricted the parameter space in which incomplete linkage disequilibrium can arise.
Finally, it is worth noting that when the fitness cost of resistance affects transmission rate only, a similar analysis shows that incomplete linkage disequilibrium never arises. In summary, under specific conditions (the cost of resistance affecting clearance rate and more commonly used antibiotics having a greater fitness cost), incomplete linkage disequilibrium is theoretically possible. However, the narrowness of the parameter range in which this occurs suggests this mechanism leading to incomplete linkage disequilibrium and non-nested resistance profiles being observed is very rare.

Recombination rate
Allowing recombination at the resistance loci breaks down linkage disequilibrium and therefore decreases the association between resistance determinants: increasing the recombination rate parameter (r) in the model represented by Equation 10 in the main text decreases mean linkage disequilibrium (Figure 3 in main text). The r parameter is not directly interpretable as the frequency at which recombinant strains are transmitted: we assume recombination requires co-infection -how often a host infected by any particular strain transmits a recombinant strain therefore also depends on the frequency of other strains. Here, to make values of r more interpretable, we explore its relationship to the proportion of transmission events involving a recombinant strain: r x y I xy . (Note that this includes all recombination events, even those leading to transmission of an identical genotype). Mostowy et al. have estimated that the PMEN1 pneumococcal lineage undergoes recombination events (at any locus) at a rate of approximately 0.01 events per month [2]. Pneumococcal transmission rates, on the other hand, have been estimated to be around 1.3-4.3 per month [3], suggesting that transmission of recombinant strains, where recombination has occurred at any locus, accounts for 0.2 -0.8% of transmission. The upper bound of the r parameter range explored in the main text (r = 0.01), on the other hand, corresponds to transmission of strains where recombination has occurred specifically at one of the three resistance loci making up 1.7% of transmission events (SI Figure D). Some resistance loci in the pneumococcus are recombination hotspots, but even so, this recombination rate is unrealistically high. Values of linkage disequilibrium remain relatively high at this unrealistic rate of recombination (Figure 3 in main text).
It is worth noting, however, that recombination rates for mobile genetic element (plasmid or transposon) associated resistance determinants may be higher if they occur on different elements. Very high rates of recombination would also abolish coexistence -and this effect occurs more rapidly than the break down of linkage disequilibrium (SI Figure E). Therefore, when fitness variation in the effect of resistance acts as a mechanism of coexistence, it will also produce linkage disequilibrium, even in the presence of recombination.  . Very high rates of recombination eliminate coexistence before the association between resistance determinants. Note the difference in x-axis range: LD is only defined between loci with more than one allele so once coexistence is only observed for a single locus, mean LD can no longer be calculated.
In the multidrug model introduced in Section 2.1.1 (Equation 3), individuals are exposed to a single antibiotic at a time: the total antibiotic consumption rate is the sum of the consumption rates for individual antibiotics ( j τ j ). At low antibiotic consumption rates, this is approximately equivalent to the assumption that different antibiotics are prescribed entirely independently (if τ 1 and τ 2 are small, the probability of being prescribed both antibiotics simultaneously (τ 1 τ 2 ) is very small). However, as discussed in main text Table 1, the assumption of independent prescription is not realistic: aspects of clinical practice lead to correlated antibiotic exposure at the individual level. Conceptually, individual-level correlations can arise in two ways: either because sub-groups within the population have different antibiotic consumption profiles or because antibiotics are sometimes prescribed together as combination therapy. In this section, we explore the effect of modelling individual level correlations in antibiotic exposure within a stratum.

Sub-groups with different consumption profiles
The presence of sub-groups with difference antibiotic consumption profiles within a homogeneously mixing population does not affect resistance dynamics. To demonstrate this, we consider a population consisting of n sub-groups (note that these sub-groups do not constitute separate strata because these groups are not assortatively mixing). Overall in the population, antibiotics 1 and 2 are consumed at rates τ 1 and τ 2 , with sub-group n making up proportions p n,1 and p n,2 of the consumption of antibiotics 1 and 2 respectively (with n i p i,1 = 1 and n i p i,2 = 1). In a scenario with two sub-groups (a and b) for example, p a,1 = 1 and p b,2 = 1 would correspond to group a consuming antibiotic 1 only and group b consuming antibiotic 2 only; while p a,1 = 1 and p a,2 = 1 would correspond to group a consuming both antibiotics and group b consuming no antibiotics. With I SS = n i I i,SS , the dynamics of strain SS in subgroup a are described by: where β is the transmission rate and U a are uninfected hosts in sub-group a. As U = n i U i , the overall dynamics of strain SS are described by: thus recovering the equation from Section 2.1.1 of the main text (Equation 3). The same reasoning applies to all strains. Therefore, the model formulation presented in the main text captures any distribution of antibiotic consumption within a homogeneously mixing population, excluding combination therapy.

Combination therapy
The model formulation in the main text assumes antibiotics are always consumed separately. To investigate the effect of this assumption, we consider a two drug model which includes simultaneous exposure to both drugs (i.e. combination therapy). In this model, a proportion p (0 ≤ p ≤ 1) of the drug with the lower rate of consumption (labelled drug 1, without loss of generality) is given in combination with the drug with the higher prescription rate (labelled drug 2). The consumption rate of drug 1 by itself is therefore (1 − p)τ 1 , the consumption rate of drug 2 is τ 2 − pτ 1 and the consumption rate of combination therapy is pτ 1 . The SS strain is cleared by all three forms of treatment, the SR strain by drug 1 by itself and combination therapy, the RS strain by drug 2 and combination therapy and the RR strain is cleared by no form of treatment. The total clearance rate experienced by each strain (λ k in main text Equation 3) is therefore:  Figure F: The effect of simultaneous antibiotic prescription (combination therapy) on resistance dynamics among a homogeneously mixing population. Combination therapy increases the antibiotic consumption rate at which resistance to the antibiotic with the higher consumption rate becomes beneficial, but does not affect the total antibiotic consumption rate at which multidrug resistance becomes selected for. Total antibiotic consumption (x-axis) refers to the overall antibiotic consumption rate in the population (τ 1 + τ 2 ). τ 1 always makes up 30% of the total antibiotic consumption. Other parameters for these simulations are: β = 2, µ = 1, c µ1 = c µ2 = 1, c β1 = c β2 = 0.98.
This model can therefore be expressed in the same general form as the multidrug model in the main text (Equation 3) -the result about competitive exclusion therefore also applies when antibiotics are prescribed simultaneously. Furthermore, the proportion of antibiotics administered as combination therapy (p) only affects the clearance rate of the SS strain (λ k ). Inclusion of combination therapy therefore affects the competition between the SR and the SS strain, but not the SR and the RR strain. This is illustrated in Figure F.
Linkage disequilibrium is not defined when allele frequencies are either 0% or 100%, it is therefore not meaningful to discuss MDR over-representation in absence of coexistence. In this model, combination therapy therefore does not promote MDR over-representation within stratum. If within stratum coexistence arises through some other mechanism however, combination therapy might well contribute to promoting linkage disequilibrium within strata.