Factors favouring the evolution of multidrug resistance in bacteria

The evolution of multidrug antibiotic resistance in commensal bacteria is an important public health concern. Commensal bacteria such as Escherichia coli, Streptococcus pneumoniae or Staphylococcus aureus, are also opportunistic pathogens causing a large fraction of the community-acquired and hospital-acquired bacterial infections. Multidrug resistance (MDR) makes these infections harder to treat with antibiotics and may thus cause substantial additional morbidity and mortality. Here, we develop an evolutionary epidemiology model to identify the factors favouring the evolution of MDR in commensal bacteria. The model describes the evolution of antibiotic resistance in a commensal bacterial species evolving in a host population subjected to multiple antibiotic treatments. We combine statistical analysis of a large number of simulations and mathematical analysis to understand the model behaviour. We find that MDR evolves more readily when it is less costly than expected from the combinations of single resistances (positive epistasis). MDR frequently evolves when bacteria are in contact with multiple drugs prescribed in the host population, even if individual hosts are only treated with a single drug at a time. MDR is favoured when the host population is structured in different classes that vary in their rates of antibiotic treatment. However, under most circumstances, recombination between loci involved in resistance does not meaningfully affect the equilibrium frequency of MDR. Together, these results suggest that MDR is a frequent evolutionary outcome in commensal bacteria that encounter the variety of antibiotics prescribed in the host population. A better characterization of the variability in antibiotic use across the host population (e.g. across age classes or geographical location) would help predict which MDR genotypes will most readily evolve.

EJ, 0000-0002-4568-283X; SL, 0000-0002-4236-828X; FD, 0000-0003-2497-833X; FB, 0000-0003-0591-2466 The evolution of multidrug antibiotic resistance in commensal bacteria is an important public health concern. Commensal bacteria such as Escherichia coli, Streptococcus pneumoniae or Staphylococcus aureus, are also opportunistic pathogens causing a large fraction of the community-acquired and hospitalacquired bacterial infections. Multidrug resistance (MDR) makes these infections harder to treat with antibiotics and may thus cause substantial additional morbidity and mortality. Here, we develop an evolutionary epidemiology model to identify the factors favouring the evolution of MDR in commensal bacteria. The model describes the evolution of antibiotic resistance in a commensal bacterial species evolving in a host population subjected to multiple antibiotic treatments. We combine statistical analysis of a large number of simulations and mathematical analysis to understand the model behaviour. We find that MDR evolves more readily when it is less costly than expected from the combinations of single resistances ( positive epistasis). MDR frequently evolves when bacteria are in contact with multiple drugs prescribed in the host population, even if individual hosts are only treated with a single drug at a time. MDR is favoured when the host population is structured in different classes that vary in their rates of antibiotic treatment. However, under most circumstances, recombination between loci involved in resistance does not meaningfully affect the equilibrium frequency of MDR. Together, these results suggest that MDR is a frequent evolutionary outcome in commensal bacteria that encounter the variety of antibiotics prescribed in the host population. A better characterization of the variability in antibiotic use across the host population (e.g. across age classes or geographical location) would help predict which MDR genotypes will most readily evolve.

Introduction
The emergence of antibiotic-resistant bacterial strains is a public health problem. Infections by antibiotic-resistant strains result in longer duration of hospital stay, poorer prognosis and additional costs compared with infections by sensitive strains [1][2][3][4][5][6][7][8][9]. Multidrug resistance (MDR) is even more critical, as it may be associated with an even higher risk of inappropriate therapy and a longer time to appropriate treatment than single resistances [10][11][12][13]. In spite of the clinical importance of MDR, there are much fewer models investigating the factors favouring the evolution of MDR than the evolution of single resistance.
It is useful to distinguish between the factors that explain the origin of MDR and those that explain their proliferation [14]. MDR strains can originate because a single biological mechanism confers resistance to multiple drugs, because several genes conferring resistances to several antibiotics are genetically linked together on the chromosome or a plasmid, or because multiple mutations conferring resistance to multiple antibiotics have evolved over the course of treatment in a host. This latter explanation may explain MDR in bacterial species causing chronic infections such as Mycobacterium tuberculosis [15]. Here, we focus on the evolutionary forces that promote the proliferation of MDR strains once they have originated. Three hypotheses have been proposed for MDR strain proliferation [14]: (i) When resistance is associated with a fitness cost, the overall cost of MDR may be lower than predicted by combining the costs of individual resistances, which would favour the evolution of MDR over single resistances. (ii) The use of multiple antibiotics in the host population may favour MDR. In particular, many bacterial species are mainly commensal, and thus primarily in contact with antibiotics not because they caused an infection themselves, but because their host is treated for infections caused by other species: this phenomenon is called 'bystander selection' [16]. These species would experience occasional antibiotic selection by the full diversity of antibiotics prescribed in the host population, which may select for MDR. (iii) If the host population is structured in different classes taking antibiotics at different rates, MDR may be overrepresented at the whole population level because MDR strains are favoured in the classes with high rates of antibiotic treatments.
Until recently, few theoretical studies had examined these three hypotheses for the proliferation of MDR. Lehtinen et al. [17] developed an epidemiological model to describe the evolution of MDR when the host population is structured (e.g. different age classes or geographical regions). This study found that MDR is over-represented when two criteria are met: firstly, variation in antibiotic consumption rates across host classes is correlated for different antibiotics ('aligned use'-classes with high consumption of one antibiotic also consume other antibiotics at a high rate); and secondly, the classes are strongly assortatively mixing: transmission within classes is much more frequent than transmission between classes. The structure of the host population may thus play a central role in the evolution of MDR but some questions remain unanswered. First of all, many previous models for the evolution of resistance do not explicitly consider hosts that are under treatment. These models assume that treatment either does not affect resistant bacteria at all, and that treatment instantaneously removes all sensitive bacteria from the host and is stopped. Yet, treated hosts form a protected niche in which resistant strains preferentially replicate, which can favour the coexistence of both resistant and sensitive strains [18]. Secondly, in a structured host population, which resistance genotypes evolve must depend on the interaction between the rate of transmission across different classes (mixing), and how the use of different antibiotics varies across classes. Thirdly, recombination between resistance loci should also affect the evolution of MDR. This is true both in the context of a single host population and when the host population is structured. In a structured host population, mixing and recombination of different genotypes from distinct classes may create maladapted genotypes that are distinct from those that would be favoured in the absence of recombination. It is not clear, in these more general scenarios of heterogeneity in antibiotic use, variable mixing and recombination, which genotypes will evolve. Lehtinen et al.'s study explored the role of non-aligned variation in antibiotic use, mixing and recombination, in separate extensions of their main model. But, to our knowledge, no study has examined the more general scenario with explicit treatment, arbitrary variations in antibiotic use and level of transmission between classes, and recombination.
Here, we develop a model examining the evolution of MDR in a structured host population that complements existing models in these aspects: we explicitly model the dynamics of untreated and treated individuals, in contrast with previous models often assuming that treatment instantaneously clears individuals colonized by a sensitive strain, and does not affect individuals colonized by a resistant strain. We model recombination between different strains, and finally, we extend this model to a host population structured in different classes and we examine how the variation in antibiotic use, the rate of transmission between classes and recombination affect MDR evolution.

. General features of the model
We model the evolution of antibiotic resistance using a compartmental model describing the dynamics of different types of hosts. We assume that the overall host population size is stable (scaled to one without loss of generality) and we describe the dynamics of the frequency of hosts colonized by the bacterial strain i, under treatment m, denoted by X i,m . We define R Ã the set of all possible bacterial strains, and R the set of all possible colonization states of the hosts, R ¼ R Ã < ; where the empty set ; denotes the absence of colonization. We denote by T Ã the set of possible antibiotic treatments (m [ T Ã ) and by T ¼ T Ã < ; the set of all possible host treatment status, including no treatment. A treatment can be a combination of multiple drugs.
The following events change the frequencies of the different host types: colonization of uncolonized hosts by a bacterial strain (transmission), clearance of colonized hosts, antibiotic treatment and supercolonization. Supercolonization describes the colonization by a challenger strain of a host already colonized by a resident strain.
We assume that the bacterial species has a mainly commensal lifestyle. Therefore, we do not model the occasional infection by the focal species, but only its asymptomatic (commensal) carriage by hosts. This assumption corresponds to some common species (e.g. Escherichia coli, Streptococcus pneumoniae, Staphylococcus aureus). We also assume that de novo evolution of resistance is negligible: resistance is primary, i.e. caused by resistant strains already circulating in the population.
where β is the baseline transmission rate, χ i is the transmission cost of resistance of strain i (equal to 0 when the strain is fully sensitive) and ψ i,m is the effect of treatment m on the transmission rate (equal to 0 when the strain is resistant to treatment m). For a given genotype, the reduction in transmission due to the cost, 1 − χ i , is computed by multiplying the effects of individual resistances that the genotype carries and pairwise interactions between all pairs of resistances. This introduces cost epistasis, whereby the cost of carrying multiple resistance may not be equal to the multiplication of all individual resistances' costs. Negative cost epistasis in particular may favour the evolution of MDR. For example, if there are only two drugs, the reduction in transmission for the double resistant is 1 À where χ e is the transmission cost epistasis parameter. In a similar way, the effect of treatment in reducing transmission 1 − ψ i,m depends on the multiplicative effects of individual drugs present in the treatment and pairwise interactions for all pairs of drugs to which genotype i is sensitive [28,29].

Treatments
The rate of start of treatment with treatment m is denoted τ m , and the rate of treatment stop is denoted ν m . Thus, 1/ν m is the expected duration time of the treatment m. These rates do not depend on the colonization status of the host. This results from the assumption that the focal species is mainly a commensal: as a consequence, antibiotic treatment is most of the time linked to an infection unrelated to the focal species ('bystander' antibiotic selection, [16]).

Clearance
The strain may be cleared from a host, either naturally or by antibiotic treatment. We denote u i the natural rate of clearance of the strain i, and a i,m the rate of clearance of the strain i by antibiotic treatment m. We assume a i,m = 0 when the strain i is resistant to treatment m ( perfect resistance).

Supercolonization and within-host competition
Supercolonization corresponds to individuals already colonized by a strain i under treatment m that are colonized by a challenger strain j replacing the resident strain i. The outcome of supercolonization is decided by within-host competition between the resident and the challenger strain. Supercolonization occurs with a reduced probability compared to transmission to uncolonized hosts. The probability of successful supercolonization is denoted σ and is assumed to be the same for all pairs of (resident, challenger) strains. During supercolonization, recombination and within-host competition occur. Recombination happens with probability ρ and may produce new genotypes. Resistance to each drug is determined by a single diallelic locus (sensitive or resistant allele). These genotypes compete within the host and a single genotype eventually dominates. We assume that these processes of recombination and within-host competition occur on a fast timescale compared to epidemiological events, and are, therefore, instantaneously resolved upon supercolonization to give rise to a host colonized by a single genotype. We denote by η i,j,m→k the probability that a colonization by strains i and j within a host with treatment status m gives rise to a majority strain k. These probabilities verify P k[R h i,j,m!k ¼ 1. These probabilities are calculated by first considering whether recombination happens within the host (with probability ρ), and second deciding which of the genotypes emerges as the winner, which is affected by the strain's competitive ability ω k,m . We define the set ϱ i,j as the set of all genotypes that can be generated by recombination between genotypes i and j. Then, the probabilities defining supercolonization transitions are The first case corresponds to the situation where genotype k is different from both i and j, and cannot be generated by recombination. The second case corresponds to the situation where genotype k is different from both i and j and must, therefore, be generated by recombination. Recombination happens with probability ρ, and will give rise to genotype k only if k belongs to the set ϱ i,j of genotypes generated by recombination between i and j. The third case corresponds to the situation where k is one of i or j and can, therefore, be the winning strain in the event of recombination or not. The fractions describe the outcome of within-host competition, whereby a strain k present in the host will win the competition between all strains present with a probability equal to its within-host competitive ability ω k,m divided by the summed competitive abilities of all strains present. This model is a simple way to implement differential competitive abilities but does not explicitly describe withinhost evolutionary dynamics. In this model, recombination between i and j is assumed to always generate all intermediate Within-host competitive abilities ω k,m depend on the cost of resistances carried by the bacterial genotype (and the potential interactions between them), and on the effect of treatment on these genotypes (and the potential interaction between drugs). In mathematical terms, the within-host competitive ability . The term 1 − c i is calculated as the multiplication of all costs of individual resistances and pairwise epistasis between all pairs of resistances. The term 1 − ε i,m is calculated as the multiplication of all effects of individual drugs to which genotype i is sensitive and interactions for all pairs of drugs.
The instantaneous change in density of a strain i under treatment m resulting from all the possible events of supercolonization can be written: The first term composed of three sums, denoted by γ i,m (supercolonization gains), is the sum of all possible supercolonizations leading to strain i in a host of treatment m (all challenger strains k × treatments n × resident strains j ). Similarly, the second term, denoted by π i,m (supercolonization losses), is the sum of all possible supercolonizations removing the resident strain i under treatment m (all possible challenger strains k × treatments n).

Mathematical form of the model
We model the four following events: colonization of uncolonized (susceptible) hosts, treatment use, clearance of colonized hosts and supercolonization (during which recombination and within-host competition happen on a fast timescale). Hosts can be uncolonized and untreated (X ;,; ), uncolonized and treated (X ;,m ), colonized and untreated (X i,; ), and, finally, colonized and treated (X i,m ). The model can be described by the following set of ordinary differential equations: royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200105 Uncolonized, untreated individuals: (n n X ;,n À t n X ;,; ) Uncolonized, treated with drug m: Colonized with strain i, untreated: Colonized with strain i, treated with drug m:

Linear discriminant analysis of the simulation results
We investigated with simulations the model's behaviour as a function of its parameters. We ran 400 000 simulations of the model with randomly chosen sets of parameters and recorded the equilibrium values. Parameters were randomly chosen in plausible distributions (table 1). The rates of antibiotic clearance a i,m were set to 20 per month when clearance is successful (treatment m clears strain i) and 0 otherwise. To analyse the results of simulations, we used linear discriminant analysis (LDA; as encoded in the MASS package [30] in R [31]). LDA reduces the high-dimensional parameter space into a low dimensional space that best separates the outcome of the model in terms of a qualitative variable of interest (see [32] for an application of LDA to a similar topic). The qualitative variable that we used was the list of all strains present at a frequency greater than 1% in the system after 300 months, a time when the equilibrium of the system is reached (electronic supplementary material, figure S2). The basis of the reduced space (here, we chose to reduce space to two dimensions) is composed of linear combinations of these parameters. The impact of each parameter on the outcome of the model can thus be described by a vector of the contribution of this parameter on each dimension.
To check the robustness of the results, we bootstrapped the 400 000 simulations, and visually confirmed that the LDA gave similar results (electronic supplementary material, figure S3).

Basic reproductive number and invasion fitness
In addition to simulations, to gain analytical understanding of the model, we computed the basic reproductive number R 0,i of a strain i, as well as the invasion fitness r i,j of strain j invading a host population already colonized by the wild-type strain i equilibrium. The basic reproductive number R 0,i quantifies the number of newly colonized individuals per colonized individual in an otherwise fully susceptible host population. It is a relevant quantity to evaluate the spread of a transmissible disease in a fully susceptible (uncolonized) population. When R 0,i > 1, the strain i spreads in the population. However, R 0 might not be a very useful quantity to evaluate the spread of a new genotype (for example, a resistant strain) in a population where another genotype (for example, a sensitive strain) is already present [33]. We will see below that R 0,i is a relevant quantity to determine which strain becomes dominant only in an approximation of the above model. The invasion fitness r i,j is the growth rate of rare strain j invading a host population where another strain i is already present and at equilibrium [34]. When positive, the strain i will invade and reach a positive equilibrium. Furthermore, if r i,j > 0 and r j,i > 0, both strains i and j can invade the population of the other when rare, meaning that both strains can coexist [18,35].
All analyses and simulations, with the exception of the LDA, were done with Wolfram Research, Inc., Mathematica, v. 11.3, Champaign, IL (2018).

Results
Our aim is to describe the bacterial strain composition at the equilibrium as a function of the model parameters. To this end, we simulated the model under randomly chosen sets of parameter values, and complemented these simulations with an analysis of the strains' basic reproductive number and invasion fitness.
In the first part of the paper, we consider the dynamics of a bacterial species in a single, well-mixed population of hosts. Three antibiotic treatments are applied to the host population: a treatment by antibiotic 1 denoted by T 1 , a treatment by a distinct antibiotic 2 denoted by T 2 and combination therapy (denoted by T 1,2 ). Thus, the bacterial population is in contact with the two drugs. The bacterial population is composed of four strains: the fully sensitive strain S, the strain R 1 resistant to the drug 1 only, the resistant R 2 resistant to the drug 2 only and the double resistant strain R 1,2 resistant to both drugs. In the second part, we relax the assumption of a single wellmixed population and consider the evolution of resistance in two classes of hosts taking antibiotics at different rates.

Factors favouring the evolution of resistant and
multi-resistant strains 3.1.1. An advantage for multi-resistance over single resistances In the single well-mixed population model, the evolution of double resistance was the most common outcome (figure 1a). For randomly chosen parameter values (table 1), at equilibrium, 29% of the 400 000 simulations resulted in the double resistant strain dominating. Furthermore, coexistence between sensitive and multiple resistant bacterial strains was rarely maintained. A single strain was maintained in 80% of the simulations. In rarer cases, two strains coexisted at equilibrium (representing 18% of parameters), while coexistence between more than two strains represented 1.9% of parameters. This property of models of drug resistance evolution has been identified before in simpler models [36,37], and, for the most part, holds true in this more complicated model that includes an explicit description of treatment dynamics and supercolonization. Intuitively, several strains cannot easily coexist because sensitive and resistant strains compete to colonize the same hosts. These results can be understood by analysing a simplified model where treatment is implicit and there is no supercolonization. By implicit, we mean that treatment instantaneously clears sensitive strains (thus decolonizing the host) and has no effect at all on resistant strains. In such a simplified model, no coexistence is maintained at equilibrium and the winning strain is the one with the highest basic reproductive number or R 0,i (electronic supplementary material, sup. text).
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200105 In this simplified model, we find that the basic reproductive number of strain i is where δ i,m is an indicator variable defining the sensitivity status of strain i, equal to 1 when strain i is sensitive, and 0 when strain i is resistant to treatment m. The strain striking the best balance between the cost of resistance (that may both reduce the transmission rate b i,; and increase the clearance rate u i ) and the benefits of resisting treatments will be the eventual winner. The analysis of basic reproductive numbers captures reasonably well the domain of dominance of a strain. This was verified in additional simulations with chosen parameters (electronic supplementary material, figure S1). This was true even in the full simulated model, which differs by the consideration of explicit treatment, supercolonization and recombination (electronic supplementary material, figure S1B,C).
While the evolution of double resistance was the most frequent outcome, it was very rare that the two single resistant strains were present in the population at equilibrium. The examination of the basic reproductive numbers suggests that when both single resistances fare better than the sensitive strain, the double resistant strain should eventually dominate (electronic supplementary material, sup. text). Thus, interestingly, in a commensal bacterial population evolving in a host population subjected to multiple treatments, the evolution of double resistance may be a common outcome even in the absence of combination therapy. This confirms one of the hypotheses for the evolution of MDR in commensal bacteria (hypothesis (ii) in the introduction).

Linear discriminant analysis of the single population model
The LDA reveals which parameters are important to the evolution of resistant and multi-resistant strains ( figure 1b). The LDA identifies domains where a certain dynamical outcome predominates in the parameter space. In our analysis, the horizontal axis separated sensitivity from double resistance. The vertical axis separated resistance to drug 1 from resistance to drug 2. In addition, each parameter can be represented as a vector in the reduced space. The norm of the vectors quantifies the importance of the corresponding parameter for determining the final outcome, given the range of parameters sampled.   Figure 1. Results of the linear discriminant analysis. (a) The frequency (in %) of occurrence of each combination of strains among the 400 000 simulations. We considered that a strain was present when its prevalence at 300 months (a time when equilibrium was reached) was greater than 1%. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200105 In line with the results derived from the basic reproductive number, the best predictors of the success of a resistant strain were a high rate of the corresponding treatment (high τ m ), a long duration of this treatment (small ν m ) and a low cost of resistance in terms of transmission and within-host competition (figure 1b). Combination therapy (corresponding to the rate of treatment t T1,2 ), in particular, favoured the evolution of the double resistant strain as expected.
The costs of resistance had a strong influence on the outcome (i.e. had the largest norms) (figure 1b). Among these costs, for our choice of parameters, the cost on transmission χ i hindered resistant strains more than the cost on within-host competitive ability c i . This is because the cost on within-host competitive ability is only expressed upon supercolonization, which in our model occurred less frequently than transmission to an uncolonized host (this is modulated by the susceptibility to supercolonization σ, table 1). Furthermore, the double resistant strain was favoured when its cost was lower than expected from the two single resistances combined. This is evident on figure 1b from the vectors representing cost epistasis, be it on transmission (parameter χ e ) or on within-host competitive ability (parameter c e ). Negative cost epistasis is indeed one of the hypotheses for the evolution of MDR (hypothesis (i) in the introduction).
The LDA also shows that low transmission, low supercolonization and long carriage duration (low u i ) all favoured MDR (figure 1b,c). To understand why, we computed the invasion fitness of different strains. The invasion fitness of a mutant is its initial growth rate when invading a resident strain at equilibrium. When it is positive, the mutant will successfully establish in the population. The analytical expressions for invasion fitness of a strain given its resistance profile inform on the factors that favour sensitivity, single resistance of multiple resistance in the model (including explicit treatment dynamics and supercolonization, in contrast with the analysis based on basic reproductive numbers). Our main analytical result is the approximation of the invasion fitness under three assumptions: (i) the treatment rates are smaller than epidemiological rates (transmission, clearance), (ii) the probability of supercolonization σ is small and (iii) recombination is ignored.
Under these assumptions (details in electronic supplementary material, Mathematica notebook), the invasion fitness of a challenger strain j invading a resident strain i denoted r i,j is, to the first order: whereX i,m denotes the density of hosts colonized by the strain i under treatment m, when the resident strain i is at equilibrium. In this expression, the first two terms correspond to the growth rate of strain j in the population at equilibrium, in the absence of treatment and supercolonization. This growth rate results from a balance between colonization of untreated hosts and natural clearance. The third term represents colonization of treated uncolonized hosts. It is summed over all possible antibiotic treatments and modulated by D j,m ¼ ðb j,mX;,; þn m Þ=ðb j,;X;,; þ n m þ a j,m Þ.
The fourth term represents antibiotic clearance, again summed over all antibiotics, and again modulated by Δ j,m . The quantity Δ j,m is between zero and one in the biologically plausible case where the transmission rates are smaller in treated individuals (b j,m b j,; 8 m). For a perfectly sensitive strain a j,m → ∞, Δ j,m = 0: this strain cannot transmit to treated individuals and is completely cleared by antibiotics. For a perfectly resistant strain (a j,m = 0, b j,m ¼ b j,; 8 m [ T ), Δ j,m = 1: the strain can fully transmit to treated individuals and is not cleared by antibiotics. The last term in the equation represents the balance resulting from supercolonization events, where h i,j,;!j is the probability that supercolonization by a strain j of an untreated (the index ;) host colonized by a strain i, leads to a host colonized by strain j following within-host competition. Similarly, h j,i,;!i is the probability that supercolonization by strain i of an untreated host colonized by a strain j, leads to a host colonized by strain i following within-host competition. Under the assumptions of small treatment rates and small probability of supercolonization, the supercolonization of treated hosts that potentially favours resistant strains is a negligible phenomenon. We ignored recombination to obtain the equation. Recombination would potentially create other strains than the resident i and the challenger j, which would greatly complicate the analysis. Later on, we use simulations to investigate the impact of recombination of evolutionary dynamics of MDR.
To further facilitate the interpretation of the expression for the invasion fitness, we express it for the strain R 1,2 perfectly resistant to the two antibiotics (no effect of treatments on clearance or on transmission, a R1,2,m ¼ 0 and where r S,R1,2 is the invasion fitness of the double resistant strain R 1,2 invading the sensitive strain S at equilibrium. The double resistant strain invades the sensitive strain at equilibrium by multiplying in both untreated and treated individuals. The double resistant strain is not hindered by treatment. However, it may suffer a cost in the form of a reduced transmission rate compared to the sensitive strain (b R1,2,; , b S,; ). It also suffers a cost in supercolonization, as its within-host competitive ability is reduced in untreated hosts: thus, the supercolonization balance is negative for the double resistant strain. Analogously, the sensitive strain invading the doubly resistant strain has invasion fitness: where conversely, r R1,2,S is the invasion fitness of the sensitive strain S invading the double resistant strain R 1,2 at equilibrium. The sensitive strain grows only in untreated individuals and is cleared by treatment, but does not suffer a cost.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200105 These analytical expressions for invasion fitness delineated when a strain was able to invade another (electronic supplementary material, figure S1, for the example of coexistence of a sensitive and double resistant strain). It also explained the rare areas of coexistence of two strains observed in the simulations (figure 1). Indeed, when two strains can reciprocally invade the other, that is r i,j > 0 and r j,i > 0, coexistence between these two strains is possible. Sensitive and resistant strains use different strategies. The sensitive strain multiplies efficiently in the untreated hosts, while the resistant strains multiply less efficiently but both in untreated and treated individuals. In other words, the treatment structure means the environment of the pathogen is multidimensional and coexistence of sensitive and resistant strains can in theory occur [33]. In practice, coexistence was a rare outcome in our simulations.
The effects of transmission, supercolonization and carriage duration can be understood from the analysis of invasion fitness: Competition through supercolonization favours the sensitive strain because it transmits more and it is fitter than the resistant strains within untreated host. This benefit of supercolonization for the sensitive strain depends on the costs of resistance on transmission χ i and withinhost competitive ability c i : in additional simulations where we set these costs equal to zero, the benefit of supercolonization disappears (figure 2b dashed lines). -Longer carriage duration (smaller natural clearance rate) favours the resistant strain (figure 2c) [38].
The slope of the relationship between invasion fitness and a parameter reflects the strength of the effect of this parameter. These slopes (figure 2) are consistent with the LDA: supercolonization has a large effect compared to both transmission and carriage duration.  Figure 2. Invasion fitnesses of the sensitive strain S (blue lines) and the double resistant strain R 1,2 (red lines) invading the other strain at equilibrium. A positive value means that the invading strain can establish. In (a), the invasion fitness is shown with supercolonization ( plain lines, σ = 0.1) and without supercolonization (dashed lines, σ = 0). In (b), the invasion fitness is shown with within-host and transmission costs of resistance ( plain lines, , when varying the inverse treatment duration (ν m ), the expected fraction of treated individuals τ m /ν m is kept constant and equal to 0.05. Unless specified otherwise, parameter values in all panels are: royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200105 Lastly, and interestingly, a strategy consisting of more frequent treatments (high treatment rates τ m ) with shorter duration (high rate of treatment cessation ν m ) while keeping the fraction of treated individuals constant advantages the double resistant strain (figure 2d). This new result is an extension to multiple treatments of a previous result according to which more frequent, shorter-duration treatment favours resistant over sensitive strains. More frequent, shorter treatments disadvantages sensitive strains because more hosts are cleared of the sensitive strains [18].

Evolution of MDR in a structured host population
The structure of the host population may affect the evolution of MDR. Hosts of different ages, living in different locations, in the community or in the hospital, represent a heterogeneous environment in which multiple bacterial strains can multiply. Host structure favours coexistence of multiple strains when the different classes are treated at very different rates and when they are strongly isolated [18]. To study how structure in the host population influences the evolution of MDR, we developed a two-class extension of our model. We used 400 000 additional simulations of this extended model with randomly chosen parameter sets to examine how MDR evolves in these two host classes. The rate of treatment start, the rate of treatment stop and the duration of carriage vary across the two classes (table 2). The set of all classes is denoted by C. The pathogen moves across host classes by inter-class transmission, which we assume is typically smaller or equal to intra-class transmission. To control these rates of transmission, we denote by β i,m,y→x = α x,y β i,m the transmission rate of strain i colonizing a host of treatment m from class y to a host of class x. The parameter α x,y is the fraction of all transmission from hosts in class y that affect hosts in class x. Thus, 8 y [ C, P x[C a x,y ¼ 1. Two classes x and y are fully isolated if α x,y = α y,x = 0. The population is homogeneously mixed if 8 (x,y) [ C 2 , α x,y = 1/N C where N C is the number of classes. Here, we consider two classes A and B.
Our simulations show that aligned use of the two antibiotics across classes favours MDR. When the host population is structured, the rates of treatment by antibiotic 1 and 2 in the different classes are now vectors t T1 and t T2 whose elements are the values t T1;i and t T2;i in each class i of host. These vectors are perfectly aligned when t T1 is proportional to t T2 , in which case the classes differ only by the intensity of prescription of all antibiotics. This scenario promotes the evolution of MDR if the classes are sufficiently isolated, because MDR evolves in the classes with high use of all antibiotics [17]. By contrast, when the two vectors are orthogonal (t T1 Á t T2 ¼ 0), at most one antibiotic is used in each class. This scenario should instead promote the evolution of single resistances in each class, if the classes are sufficiently isolated. To verify this intuition, we computed the cosine of the angle between the antibiotic use vector as t T1 Á t T2 =(kt T1 kkt T2 k). This cosine varies from 0 (aligned uses) to 1 (orthogonal uses).
The effect of the alignment of antibiotic uses across classes on MDR was evident on the LDA of the 400 000 simulations of the two-classes model (electronic supplementary material, figure S4). Furthermore, the LDA also shows, just as in the single population model, that the evolution of double resistance in the two-classes model is favoured by high rates of treatments, combination therapy, low costs of resistance, negative cost epistasis, low supercolonization.
To further quantify MDR, we computed the linkage disequilibrium between the two resistance loci. Calling the frequencies of each bacterial genotype f S , f R1 , f R2 , f R1,2 , linkage disequilibrium is defined as LD ¼ f R1,2 f s À f R1 f R2 . LD varies between −0.25 (resistances are found in single resistant genotypes only) and +0.25 (resistances are found in double resistant genotypes only). LD declines with the level of polymorphism, and is 0 when there is no polymorphism at either locus. We also define the corrected linkage disequilibrium as The corrected LD varies from −1 (resistances are found in single resistant genotypes only) to +1 (resistances are found in double resistant genotypes only) and is unaffected by the level of polymorphism. We found that the mean linkage disequilibrium increased with the alignment of use of antibiotics 1 and 2 (figure 3a).
By focusing on the subset of simulations with very aligned uses (cosine of the angle greater than 0.95) and those with very orthogonal uses (cosine of the angle less than 0.05), we Table 2. Sampling distributions of the parameters in the two-classes model (other parameters as in the well-mixed single population model). We explore a wide range of proportion of inter-class transmission. For example, transmission between individuals from different regions or countries is likely to be small. The transmission rate between individuals of different age classes is around 10-50% that within age classes [39]. Sampling ranges for other parameters are justified as for the single population model. found that the emergence of positive or negative linkage disequilibrium depends on the classes being isolated one from another (inter-class transmission has to be much smaller than intra-class transmission) (figure 3b). However, the rate of recombination had little effect on linkage disequilibrium, both for aligned and orthogonal antibiotic uses (figure 3c). The build-up of positive or negative LD depends on different genotypes being locally favoured in one or the other class. For aligned uses, R 1,2 is favoured in one class and S is favoured in the other. For orthogonal use, R 1 is favoured in one class and R 2 is favoured in the other. In both cases, this results in LD emerging at the whole population level. Large inter-class transmission erodes polymorphism, which reduces LD. In the cases when polymorphism is maintained at high interclass transmission, it is maintained by mechanisms other than population heterogeneity among classes (i.e. treated individuals providing a niche for the resistant strain and/or supercolonization providing a frequency-dependent advantage to the sensitive strain), which do not generate LD. This results in the decline of corrected LD with inter-class transmission. Because LD is generated by heterogeneous selection across classes and limited mixing, there is little effect of recombination on LD: within each class the genetic diversity is small, and therefore recombination does not have an impact.

Discussion
We examined the factors favouring the evolution of MDR in commensal bacteria. We used an epidemiologicalevolutionary model describing transmission, natural clearance,  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200105 supercolonization, within-host competition and recombination. We investigated how bystander selection by antibiotic treatments selects for single or multiple resistance in commensal bacteria.
In the introduction, we recalled three non-mutually exclusive hypotheses for the proliferation of MDR strains: (i) negative cost epistasis, whereby the overall cost of MDR may be lower than predicted from the costs of individual resistances; (ii) the use of multiple antibiotics in the host population; and (iii) the structure in the host population. We confirmed the three hypotheses and showed that: -With the range of parameters chosen in our simulations, negative cost epistasis was the most potent factor selecting for MDR (figure 1b, parameters χ e and c e ); -Perhaps our most important result is that the application of several single drug treatments in the host population favoured MDR as much as combination therapy. Combination therapy of course favours MDR, but is comparatively rarer ( figure 1b, table 1). In other words, MDR may evolve even when drugs are never taken in combination in individual hosts (in agreement with previous results from a simpler model [17]). -In a structured host population, alignment in antibiotic use across classes favours MDR (figure 3a, electronic supplementary material, figure S3C).
Negative cost epistasis means that the overall cost of two resistances combined on a MDR genotype is lower than expected. This is most often called 'positive epistasis', as this corresponds to fitter than expected MDR genotypes. Positive epistasis between resistance mutations is common [40]. Several studies evaluated the costs of resistance and epistasis by measuring the fitness of sensitive, resistant and MDR strains in the absence of antibiotic in vitro [41,42]. For example, Trindade et al. [41] generated E. coli mutants resistant to the antibiotics nalidixic acid, rifampicin and streptomycin. Seventy-three percent of double mutants with significant epistasis had positive epistasis. Thus, MDR genotypes with these low-cost pairs of resistance mutations could readily evolve if the cost measured on the growth rate in vitro is a good proxy for the within-host and epidemiological costs of our model. Consistent with this idea, the pairs of resistance mutations that had the lowest in vitro costs were also those found most frequently in clinical isolates of M. tuberculosis [43]. We emphasize several assumptions of our model. First of all, we assumed that antibiotic treatment was given independently of the colonization status of the host (bystander exposure). We indeed consider commensal species that are most of the time carried asymptomatically and rarely cause infections. Treatment is given for another infection that is not modelled. However, for some drug-species combinations, the proportion of antibiotic exposure that is targeted to the species may be significant. For example, close to 20% of E. coli's exposure to quinolones is targeted rather than bystander exposure [16]. Indeed, urinary tract infections are commonly caused by E. coli and treated with quinolones. When antibiotic use is targeted, the resistance of the focal species to the first-line drug may result in treatment failure and use of a second drug. This sequential pattern of antibiotic use, not considered in the present model, could further promote MDR because single resistant strains would be more frequently exposed to a drug they are sensitive to. Second, we assumed for simplicity that only one strain could colonize a host at a time. Within-host competition between strains was assumed to be strong and to rapidly result in a single strain taking over the host. However, in some species with high prevalence like E. coli, multiple strains are frequently present, and can coexist over long time periods in a host [25,44]. In these species, a costly resistant strain in competition with a sensitive strain within a host could decline slowly or persist at low frequency. Less intense within-host competition could weaken the effect of within-host parameters on MDR evolution. It would also affect how transmission and supercolonization influence the evolution of resistance, as discussed below.

The role of transmission and supercolonization on MDR
Interestingly, higher transmission favoured the sensitive strain ( figure 1d). This opposes the frequent assumption that high transmission (e.g. poor hygiene) favours the proliferation of antibiotic resistance. A small number of empirical studies suggest that higher transmission favours resistance. Antibiotic resistance of several bacterial species was found to be associated with poorer infrastructure and poorer governance across 103 countries [45]. One interpretation of this correlation is that the dissemination of resistance is favoured by poor sanitation. More convincingly, the individual carriage of extended-spectrum beta-lactamase-producing E. coli (an enzyme conferring resistance to beta-lactams antibiotics) was strongly associated with an indicator of poor hygiene in a large cohort [46]. However, it is not at all obvious in theory why transmission would favour resistance. Larger transmission implies that both the resistant and the sensitive strain transmit more. Why would this preferentially favour the resistant strain? In a previous study, it was shown that transmission favours resistance because resistant strains, unlike sensitive strains, rely on transmitting to treated individuals for multiplication [18]. In the model analysed here, this first effect is present (figure 2a) but is overwhelmed by a second effect: the sensitive strain is favoured by transmission and supercolonization of already occupied hosts, because it enjoys a fitness advantage in competition with the resistant strain within the host (as in the model presented in [47]). In our simulations, that second effect is stronger than the first, and higher transmission, therefore, favours the sensitive strain. It is possible that in the contexts of the empirical studies mentioned above, supercolonization and within-host competition are less common than in the range of parameters explored in the simulations. The higher transmission would thus favour resistant strains as they are able to colonize empty treated hosts. Because of the two effects, larger transmission favours the invasion of a sensitive strain by a resistant, and reciprocally (figure 2a, both invasion fitnesses increase with transmission). Thus, large transmission should also promote the coexistence of multiple strains. This prediction was verified in the simulations: the probability that both resistances coexist with sensitivity (frequency of both resistances in [0.01, 0.99]) increases with transmission ( figure 4).
The effects of transmission and supercolonization on the dynamics of resistance and sensitivity are complex and would be interesting to elucidate further. For example, introducing the possibility of long-term coexistence between multiple strains within hosts could affect our results. In a theoretical study, it was shown that co-colonization with resistant and sensitive strains favours the resistant strain because the elimination of the sensitive strain upon antibiotic treatment releases the resistant strain [47]. The transmission rate, therefore, impacts the competition between sensitive and resistant strains in three ways: through the increased transmission of the resistant strain in treated hosts, through increased supercolonization (which favours the sensitive strain if it has an advantage in supercolonization) and through increased co-colonization (which may favour either sensitive or resistant strain). The overall effect of transmission on the frequency of resistance likely depends on the balance of these three effects.

The role of recombination
The rate of recombination between strains colonizing the same host had little effect on linkage disequilibrium, both in a single population (figure 1d), and in a structured population (figure 3c). This may seem surprising, as recombination is expected to break down LD. In a single population, the most common outcome is the loss of polymorphism and the dominance of a single strain: thus, recombination does not disrupt this equilibrium. In a structured population, polymorphism may be maintained at equilibrium. However, the maintenance of polymorphism implies that the different classes are strongly isolated (low inter-class transmission), and, effectively, distinct strains are rarely in contact in the same host. Thus, recombination also has little effect on the equilibrium.
However, recombination could affect the transient temporal dynamics of the model. In a scenario where the eventual equilibrium is coexistence between sensitive and double resistant strains, recombination accelerates the emergence of double resistance but results in a lower equilibrium frequency of double resistance (electronic supplementary material, figure S5). When the sensitive and two single resistant strains are initially present but the double resistant is initially absent, only recombination between the strains R 1 and R 2 can generate the double resistant strain R 12 . Thus, recombination will initially speed up the emergence of double resistance. However, at equilibrium, single resistant strains are constantly generated by recombination between the sensitive and double resistant, resulting in a lower equilibrium frequency of double resistance with increasing recombination. This results in the coexistence of all four possible strains, but only the sensitive and double resistant strains would coexist without recombination.
Previous work showed that the effects of recombination on the transient dynamics of MDR are potentially more complex. For example, when the MDR strain is less costly than expected (negative cost epistasis), the MDR strain is favoured but recombination slows down its spread by creating less fit single resistances [48]. Moreover, recombination is more frequent when the prevalence of the bacterial population in the host is higher. The larger transmission will result in higher prevalence and more recombination, which will affect the transient dynamics of resistance [48]. These processes should take place in our model as well.

Conclusion
The epidemiology and evolution of resistance to antibiotics depends on the complex ecology of bacteria [18,26,38,47]. Bacteria multiply in a diversity of hosts-already colonized by a strain or not, untreated or treated by a variety of drugs, belonging to different age classes, geographical locations, etc. We formulated an epidemiological-evolutionary simulation model that included these phenomena and elucidated its behaviour with statistical analysis of a large number of simulations and mathematical analysis. We generated a number of predictions, the most interesting of which are (i) multiple single treatments received by the host population as a whole can drive the evolution of MDR; (ii) in a structured host population, the alignment of treatments across classes favours MDR; (iii) transmission and supercolonization favour the evolution of sensitive strains; (iv) higher transmission results in a greater diversity of sensitive and resistant strains. How the alignment of uses of different antibiotics across classes drives MDR has not been empirically investigated so far. To do so, it would be necessary to identify the forms of population structure that are most relevant to bacterial antibioresistance. MDR should occur more frequently for the pairs of drugs with the more aligned uses across classes. More broadly, it would be interesting to link more closely large scale epidemiological datasets on resistance with these new theoretical results.