Understanding the evolution of multiple drug resistance in structured populations

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: June 16, 2021 (This version)
Accepted Manuscript published: June 1, 2021 (Go to version)
Accepted: May 28, 2021
Received: December 10, 2020

1. Part of Collection
Evolutionary Medicine: A Special Issue

Edited by George H Perry et al.
Further reading

Article
Figures and data
Abstract
Introduction
Results
Discussion
Materials and methods
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

The evolution of multidrug resistance (MDR) is a pressing public health concern. Yet many aspects, such as the role played by population structure, remain poorly understood. Here, we argue that studying MDR evolution by focusing upon the dynamical equations for linkage disequilibrium (LD) can greatly simplify the calculations, generate more insight, and provide a unified framework for understanding the role of population structure. We demonstrate how a general epidemiological model of MDR evolution can be recast in terms of the LD equations. These equations reveal how the different forces generating and propagating LD operate in a dynamical setting at both the population and metapopulation levels. We then apply these insights to show how the LD perspective: (i) explains equilibrium patterns of MDR, (ii) provides a simple interpretative framework for transient evolutionary dynamics, and (iii) can be used to assess the consequences of different drug prescription strategies for MDR evolution.

Introduction

Antibiotic resistance is one of the biggest current public health problems, with antibiotic resistant infections responsible for tens of thousands of deaths annually (O’Neill, 2015). Of particular concern is the evolution of multidrug resistant (MDR) pathogens, that is, pathogens resistant to multiple classes of antibiotics. Despite its importance, understanding the evolution of MDR remains an ongoing challenge, as it is typically not captured by our understanding of the evolution of single drug resistance (for which there is a large body of theory; e.g., Blanquart, 2019; Bonhoeffer et al., 1997; Lipsitch et al., 2000; Bergstrom et al., 2004; Austin and Anderson, 1999). For instance, suppose we have two drugs, $A$ and $B$ , and that a fraction $f_{A B}$ of infections caused by the pathogen of interest are resistant to both drugs. To understand MDR evolution, we need to understand what determines the frequency $f_{A B}$ . If $f_{A}$ and $f_{B}$ are the frequency of infections resistant to drugs $A$ and $B$ , and $D$ denotes any non-random association between resistance to drugs $A$ and $B$ , then

f_{A B} = f_{A} f_{B} + D .

If $D = 0$ , then the evolution of resistance to each drug is independent, and so multiple drugs do not qualitatively alter the evolutionary dynamics of single drug resistance. However, whenever $D \neq 0$ , understanding the fitness costs and benefits of resistance to each drug in isolation is insufficient to understand the evolution of MDR, because doing so will not tell us what factors govern the propagation of $D$ , which in turn will affect $f_{A}$ and $f_{B}$ . Thus the challenge of understanding MDR evolution can be recast as understanding the dynamics of $D$ . The quantity $D$ is referred to as linkage disequilibrium (LD), and it has been extensively studied in population genetics (e.g. Lewontin, 1964; Felsenstein, 1965; Ohta, 1982a; Barton, 1995; Rice, 2004; Slatkin, 2008), particularly as it relates to population structure (Ohta, 1982b; Slatkin, 1975; Li and Nei, 1974; Nei and Li, 1973; Lenormand and Otto, 2000; Martin et al., 2006). However, there has been little attempt to apply these insights to MDR evolution; often the dynamics of doubly resistant infections are neglected to simplify the analysis of single drug resistance (e.g. Bergstrom et al., 2004; Bonhoeffer et al., 1997; Beardmore et al., 2017).

Here, we consider a simple epidemiological model of a primarily asymptomatically carried pathogen (e.g. Staphylococcus spp. or Enterococcus spp.) in a structured host population. We show how this model relates to general dynamical equations for LD (Day and Gandon, 2012), in turn revealing the role of population structure in MDR evolution. We then use these equations to show how analyzing problems from the LD perspective: (i) reveals the evolutionary logic underlying patterns of MDR at equilibrium, which we use to build on a recent paper on MDR evolution (Lehtinen et al., 2019); (ii) provides a framework for understanding transient evolutionary dynamics; and (iii) provides insight on the consequences different drug prescription strategies have on MDR, which we apply to a hospital-community setting.

Results

In what follows we will introduce and analyze a model of MDR evolution. We will highlight the most important aspects here while providing more extensive details in the Materials and methods 'Model derivation'. All notation used is summarized in Table 1.

Table 1

Notation used in main text.

In all cases, a quantity indexed with a superscript $x$ is the population $x$ quantity, whereas the absence of a superscript $x$ implies the quantity is for the metapopulation.

Symbol	Description
$I_{i j}^{x}$	Density of $i j$ -infections in population $x$ , where $i = A$ (resp. $i = a$ ) if infection is resistant (resp. sensitive) to drug $A$ and $j = B$ (resp. $j = b$ ) if infection is resistant (resp. sensitive) to drug $B$ .
$I^{x}$	Density of total infections in population $x$ .
$f_{d}^{x}$ , ${\bar{f}}_{d}$	Frequency of infections resistant to drug $d$ in population $x$ and the metapopulation, respectively.
$D^{x}$ , $\bar{D}$ , $D_{M}$	Linkage disequilibrium (LD) in population $x$ , average LD across populations and metapopulation LD, respectively.
$m^{x \to y}$	Per-capita rate at which hosts migrate from population $x$ to $y$ .
$r^{x}$ , $\bar{r}$	Per-capita growth rate of sensitive infections in population $x$ (or ‘baseline’ per-capita growth rate) and average across populations, respectively.
$s_{d}^{x}$ , ${\bar{s}}_{d}$	Additive selection coefficient for resistance to drug $d$ in population $x$ and average selection across populations, respectively.
$s_{E}^{x}$ , ${\bar{s}}_{E}$	Epistasis in fitness across drug resistance loci in population $x$ and average across populations, respectively.
$ϕ μ_{i j}^{x}$ , $ϕ ρ_{i j}^{x}$	Net change in $i j$ -infections in population $x$ due to mutation or recombination, respectively.
$μ_{i}^{x}$ , ${\bar{μ}}_{i}$	Per-capita rate at which mutations generate allele $i$ in population $x$ and average across populations, respectively.
$ρ_{i}^{x}$ , ${\bar{ρ}}_{i}$	Per-capita rate at which recombination leads to gain of allele $i$ in population $x$ and average across populations, respectively.
$s^{x}$ , $\bar{s}$	Average selection for drug resistance in population $x$ and average across populations, respectively.
$cov (X, Y)$	Covariance between the variables $X$ and $Y$ , that is, $cov (X, Y) = 𝔼 [X Y] - 𝔼 [X] 𝔼 [Y]$ , where $𝔼 [X]$ denotes the expectation of quantity $X$ .
$coskew (X, Y, Z)$	Coskewness between the quantities $X$ , $Y$ , $Z$ , that is, $coskew (X, Y, Z) = 𝔼 [(X - 𝔼 [X]) (Y - 𝔼 [Y]) (Z - 𝔼 [Z])]$ .

Consider an asymptomatically carried pathogen in a metapopulation consisting of $N$ host populations in which two drugs $d$ are prescribed, specifically, drug $d = A$ and drug $d = B$ . Focus upon population $x$ . Let $S^{x}$ and $I_{i j}^{x}$ denote the density of susceptible hosts and $i j$ -infections, respectively, at time $t$ , where $i$ indicates if the infection is resistant ( $i = A$ ) or not ( $i = a$ ) to drug $A$ and $j$ indicates if the infection is resistant ( $j = B$ ) or not ( $j = b$ ) to drug $B$ . Susceptible hosts contract $i j$ -infections at a per-capita rate $β_{i j}^{x} I_{i j}^{x}$ , where $β_{i j}^{x}$ is a rate constant, while $i j$ -infections are naturally cleared at a per-capita rate $α_{i j}^{x}$ . Hosts are treated with drugs $A$ , $B$ , or both in combination at per-capita rates $τ_{A}^{x}$ , $τ_{B}^{x}$ , and $τ_{A B}^{x}$ , respectively. Treatment is instantaneous and resistance is complete, that is, if the host that receives treatment is infected by a strain sensitive to the drug, the infection is cleared instantaneously, whereas if the host that receives treatment is infected by a strain resistant to the drug, treatment has no effect. Hosts move from population $x$ to $y$ at a per-capita rate $m^{x \to y}$ . Transmission between infected hosts leads to superinfection with probability $σ$ in which either strain is equally likely to instantaneously outcompete the other (Nowak and May, 1994; Alizon, 2013). We therefore do not allow for prolonged coinfection (Materials and methods 'Model derivation'). Finally, individual infections acquire allele $i$ through either mutation or recombination (during superinfection) at per-capita rates $μ_{i}^{x}$ and $ρ_{i}^{x}$ , respectively (note that $ρ_{i}^{x}$ depends upon infection densities, see Materials and methods Equation (13)).

From these epidemiological assumptions, the change in $i j$ -infections in population $x$ can be written as the sum of four processes

\frac{d I_{i j}^{x}}{dt} = (\overset{p e r - c a p i t a g r o w t h}{\overset{⏞}{r^{x} + 1_{A} s_{A}^{x} + 1_{B} s_{B}^{x} + 1_{A} 1_{B} s_{E}^{x}}}) I_{i j}^{x} + \overset{mutation}{\overset{⏞}{ϕ μ_{i j}^{x}}} + \overset{recombination}{\overset{⏞}{ϕ ρ_{i j}^{x}}} + \overset{migration}{\overset{⏞}{\sum_{y = 1}^{N} (m^{y \to x} I_{i j}^{y} - m^{x \to y} I_{i j}^{x})}},

where $𝟏_{d}$ is equal to 1 if the $i j$ -infection is resistant to drug $d$ and 0 otherwise (e.g., if $i j = A B$ , then the per-capita growth is $r^{x} + s_{A}^{x} + s_{B}^{x} + s_{E}^{x}$ ) and $ϕ μ_{i j}^{x}$ and $ϕ ρ_{i j}^{x}$ denote the net change in $i j$ -infections due to mutation and recombination (Figure 1; Materials and methods Equations (11) and (13)). To faciliate comparison with previous results, we have broken the per-capita growth term into four components: the ‘baseline’ per-capita growth rate, $r^{x}$ , the (additive) selection coefficients for resistance to drugs $A$ and $B$ , $s_{A}^{x}$ and $s_{B}^{x}$ , and any epistatic interactions, $s_{E}^{x}$ . These latter terms have the standard interpretation. If $s_{A}^{x} > 0$ (resp. $s_{B}^{x} > 0$ ), then resistance to drug $A$ (resp. $B$ ) is selected for. If $s_{E}^{x} > 0$ , there is positive epistasis, and the per-capita growth rate of doubly-resistant infections is greater than would be expected by consideration of the per-capita growth rate of singly-resistant infections. Thus although Equation (2) is derived from a specific model, the partitioning is very general and applies to many epidemiological scenarios. We stress that any of the terms $s_{d}^{x}$ , $s_{E}^{x}$ , $ϕ μ_{i j}^{x}$ , and $ϕ ρ_{i j}^{x}$ may themselves depend upon population densities (see Figure 1 for a concrete example). Note that this partitioning is not arbitrary, particularly as it applies to the selection coefficients and epistasis. The additive selection coefficients and epistasis are defined in terms of their effect upon fitness. In continuous time models, fitness is per-capita growth. Thus, the selection coefficient for allele $k$ measures the additive contribution of allele $k$ to fitness, while epistasis measures the excess of the fitness of strain $A B$ over its value if fitness were additive across the two loci (e.g. Felsenstein, 1965; Karlin, 1975; Rice, 2004; Kouyos et al., 2009) (see also Box 1). Because epistasis is defined in terms of fitness, how costs of resistance are modeled will typically have implications for whether epistasis occurs or not; for example, multiplicative costs will generate epistasis (Box 1; Materials and methods 'Equilibrium analysis of metapopulation consisting of independent populations'). We will return to this point in the examples.

Figure 1

Download asset Open asset

Schematic of the dynamics of system (2).

The metapopulation consists of $N$ connected populations. Each population has four possible types of infections, linked by one-step mutation or recombination (blue and red arrows), whose per-capita rates are independent of genetic background. The ‘baseline’ per-capita growth rate of sensitive infections is $r^{x}$ , the additive selection coefficients for drug $A$ and $B$ resistance are $s_{A}^{x}$ and $s_{B}^{x}$ , respectively, while $s_{E}^{x}$ denotes any epistatic interactions. In the inset, we compute these quantities for the specific model introduced in the main text, using the notation that $Δ z_{d}^{x}$ and $Δ z_{E}^{x}$ are the contribution of trait $z$ to the additive selection coefficient (for resistance to drug $d$ ) and to epistasis, respectively, in population $x$ (e.g., $Δ β_{A}^{x} = β_{A b}^{x} - β_{a b}^{x}$ and $Δ β_{E}^{x} = β_{A B}^{x} - β_{a b}^{x} - Δ β_{A}^{x} - Δ β_{B}^{x}$ ).

Box 1.

Costs of resistance, epistasis, and multidrug resistance.

The spread of multidrug resistance (MDR) is driven by selection acting on each drug resistance locus, but also on the linkage disequilibrium (LD), which can be produced by epistasis in fitness. Epistasis measures the interaction between resistance alleles (mutations) at different loci and is defined in terms of the per-capita growth rates of different genotypes as: $s_{E}^{x} \equiv r_{A B}^{x} + r_{a b}^{x} - r_{A b}^{x} - r_{a B}^{x}$ .

Selection at each locus, e.g. $s_{A}^{x} = r_{A b}^{x} - r_{a b}^{x}$ , depends on the effects of the mutations on the phenotypic traits of the pathogen. However, non-additive interactions among these mutations can create epistasis (see inset in Figure 1). To better see how these non-additive effects can emerge, consider the costs of drug resistance on pathogen transmission. Let $c_{β_{d}}^{x}$ denote the parameter controlling the cost of resistance to drug $d$ in population $x$ . Then using the notation of Figure 1.

	Transmission rates			Epistasis
	$β_{A b}^{x}$	$β_{a B}^{x}$	$β_{A B}^{x}$	$Δ β_{E}^{x}$
Additive	$β_{a b}^{x} - c_{β_{A}}^{x}$	$β_{a b}^{x} - c_{β_{B}}^{x}$	$β_{a b}^{x} - c_{β_{A}}^{x} - c_{β_{B}}^{x}$	0
Multiplicative	$β_{a b}^{x} (1 - c_{β_{A}}^{x})$	$β_{a b}^{x} (1 - c_{β_{B}}^{x})$	$β_{a b}^{x} (1 - c_{β_{A}}^{x}) (1 - c_{β_{B}}^{x})$	$β_{a b}^{x} c_{β_{A}}^{x} c_{β_{B}}^{x}$

Hence, only multiplicative costs generate non-additive interactions between loci on transmission, $Δ β_{E}^{x}$ , which leads to epistasis (inset of Figure 1); and, in turn epistasis produces LD which affects MDR evolution.

Of course, the magnitude of drug resistance costs and the interaction between these costs at multiple loci need not be additive nor multiplicative. Subject to appropriate constraints on the choice of costs (e.g. $0 \leq β_{i j}^{x} \leq β_{a b}^{x}$ ), our general framework can account for any pattern of epistasis (see Figure 1). The important point is that since epistasis is a key contributor to multilocus evolution, understanding when it occurs and what is producing it (in this case, assumptions about the cost of resistance) can provide valuable insight into the evolutionary dynamics of MDR.

While system (2) contains all the information necessary to analyze MDR evolution, as currently written it is particularly opaque for providing insight. Therefore, we would like to transform it to a form which brings to the forefront the different factors that promote or impede MDR evolution; the way to do this is by focusing upon the dynamical equations for linkage disequilibrium (LD) (Day and Gandon, 2012; Slatkin, 2008). However, the inclusion of multiple populations means that doing so is not as simple as Equation (1) would suggest since there are different scales at which LD and MDR can be measured. As the scale which is of most interest will depend upon the specifics of the problem, in what follows we will consider MDR evolution at both the population- and metapopulation-level.

Population-level multidrug resistance

To understand MDR evolution in a given population, say $x$ , we need to understand the dynamics of the frequency of infections resistant to drug $A$ and $B$ , $f_{A}^{x}$ and $f_{B}^{x}$ , and the dynamics of population LD, $D^{x}$ . First, consider the dynamics of $f_{A}^{x}$ (mutatis mutandis $f_{B}^{x}$ ). Using Equation (2), it is straightforward to compute

\begin{array}{ll} \frac{d f_{A}^{x}}{d t} & = \underset{d i r e c t s e l e c t i o n}{\underset{⏟}{s_{A}^{x} f_{A}^{x} (1 - f_{A}^{x})}} + \underset{i n d i r e c t s e l e c t i o n}{\underset{⏟}{s_{B}^{x} D^{x}}} + \underset{epistasis}{\underset{⏟}{s_{E}^{x} f_{A}^{x} (1 - f_{A}^{x}) \frac{f_{A B}^{x}}{f_{A}^{x}}}} + \\ \underset{m u t a t i o n a n d r e c o m b i n a t i o n}{\underset{⏟}{(μ_{A}^{x} + ρ_{A}^{x}) (1 - f_{A}^{x}) - (μ_{a}^{x} + ρ_{a}^{x}) f_{A}^{x}}} - \underset{migration}{\underset{⏟}{\sum_{y = 1}^{N} m^{y \to x} \frac{I^{y}}{I^{x}} (f_{A}^{x} - f_{A}^{y})}} . \end{array}

where $I^{x}$ is the total density of infections in population $x$ and $f_{A B}^{x} = D^{x} + f_{A}^{x} f_{B}^{x}$ is the frequency of doubly-resistant infections. A related formulation to Equation (3) can be found in Day and Gandon, 2012 (see also Rice, 2004).

Equation (3) is partitioned into recognizable quantities. First, if resistance to drug $A$ is selectively advantageous, $s_{A}^{x} > 0$ , then drug $A$ resistance will increase due to direct selection whose strength is dictated by the genetic variance at the locus, $f_{A}^{x} (1 - f_{A}^{x})$ (Fisher, 1930). Second, if doubly-resistant infections are over-represented in the population, $D^{x} > 0$ , and resistance to drug $B$ is selected for, $s_{B}^{x} > 0$ , then drug $A$ resistance will increase due to indirect selection upon resistance to drug $B$ . Third, if epistasis is positive, $s_{E}^{x} > 0$ , and there is genetic variance at the locus, drug $A$ resistance will increase due to the disproportionate growth of doubly-resistant infections. Fourth, mutation and recombination will increase drug $A$ resistance when there is a mutation or recombination bias toward gain of drug $A$ resistance, $μ_{A}^{x} > μ_{a}^{x}$ or $ρ_{A}^{x} > ρ_{a}^{x}$ , and the frequency of infections sensitive to drug $A$ exceeds the frequency of infections resistant to drug $A$ , $1 - f_{A}^{x} > f_{A}^{x}$ . Finally, migration acts to reduce differences between populations.

It follows that drug $B$ treatment alters the predicted dynamics of resistance to drug $A$ via two main effects: (i) the influence of epistasis and (ii) indirect selection on resistance to drug $B$ mediated through the presence of LD ( $D^{x} \neq 0$ ). Thus, consider the dynamics of $D^{x}$ ,

\begin{array}{ll} \frac{d D^{x}}{dt} = \underset{selection}{\underset{⏟}{(s_{A}^{x} - s^{x} + s_{B}^{x} - s^{x}) D^{x}}} - \underset{m u t a t i o n a n d r e c o m b i n a t i o n}{\underset{⏟}{(μ^{x} + ρ^{x}) D^{x}}} + \underset{epistasis}{\underset{⏟}{s_{E}^{x} f_{A B}^{x} f_{a b}^{x}}} - \underset{migration}{\underset{⏟}{\sum_{y = 1}^{N} m^{y \to x} \frac{I^{y}}{I^{x}} (D^{x} - D^{y} - (f_{A}^{x} - f_{A}^{y}) (f_{B}^{x} - f_{B}^{y}))}}, \end{array}

where $s^{x} = f_{A}^{x} s_{A}^{x} + f_{B}^{x} s_{B}^{x} + f_{A B}^{x} s_{E}^{x}$ is the average selection for resistance, $f_{a b}^{x} = 1 - f_{A}^{x} - f_{B}^{x} + f_{A B}^{x}$ is the frequency of doubly-sensitive infections, and $μ^{x}$ and $ρ^{x}$ are the total per-capita rates of mutation and recombination, respectively (e.g. $μ^{x} = μ_{a}^{x} + μ_{A}^{x} + μ_{b}^{x} + μ_{B}^{x}$ ; Materials and methods 'Model derivation').

Equation (4) is partitioned into four key processes. First, excess selection for resistance to drug $A$ (resp. $B$ ), $s_{A}^{x} - s^{x}$ , can cause pre-existing LD ( $D^{x} \neq 0$ ) to increase or decrease. For example, if $s_{A}^{x} > s^{x}$ and $D^{x} > 0$ then LD will increase. This is because drug $A$ resistant infections are fitter than the average resistant infection and so will increase in frequency. If $D^{x} > 0$ , it is more likely this increase will occur in doubly-resistant infections, thereby increasing $D^{x}$ . Second, mutation and recombination removes any LD present at a rate proportional to the LD (Rice, 2004; Slatkin, 2008). Third, epistasis generates same-sign LD, that is, positive epistasis, $s_{E}^{x} > 0$ , leads to MDR over-representation, $D^{x} > 0$ (Felsenstein, 1965; Lewontin and Kojima, 1960; Lewontin, 1964). Positive epistasis could occur if double-resistance costs are less than expected (Trindade et al., 2009; MacLean et al., 2010; Hall and MacLean, 2011) or drugs are prescribed in combination (Bretscher et al., 2004; Day and Gandon, 2012).

Migration is the final term of Equation (4) and reveals how the metapopulation structure affects population LD. Like epistasis, migration does not require preexisting LD to operate on LD (Li and Nei, 1974; Slatkin, 1975; Feldman and Christiansen, 1974; Ohta, 1982a; Ohta, 1982b). In particular, LD in population $x$ will be generated whenever the frequencies of resistance to drugs $A$ and $B$ differ between population $x$ and any other connected population, say $y$ . If both types of resistance are more common in one population than the other $(f_{A}^{x} - f_{A}^{y}) (f_{B}^{x} - f_{B}^{y}) > 0$ , then migration will generate positive LD in both populations, $D^{x} > 0$ and $D^{y} > 0$ . If instead drug $A$ resistance is more prevalent in one population, while drug $B$ resistance is more prevalent in the other, migration will generate negative LD in both populations.

Notice the presence of the multiplier $I^{y} / I^{x}$ in the final term of Equation (4). If the populations have roughly the same density of infections, then this term is unimportant. However, when one population, say $y$ , has much fewer total infections than population $x$ , $I^{y} ≪ I^{x}$ , the term $I^{x} / I^{y}$ will be very large, whereas $I^{y} / I^{x}$ will be very small. Consequently, the ability of migration to propagate LD will be greater in population $y$ than $x$ , and so all else being equal we would predict the population with a lower density of infections will have a greater magnitude of LD than the population with a higher density of infections.

The next insight shows the importance of also taking into account Equation (3). In particular, if we only inspected the migration term of Equation (4) we might conclude that as the per-capita migration rate, $m^{y \to x}$ , increases, so too will the ability of migration to propagate LD. However, the magnitude of population LD is actually maximized at intermediate migration rates (Figure 2). The reason is because the quantity $m^{y \to x}$ has two effects. On the one hand, it directly multiplies the migration term in Equation (4) thereby magnifying migration’s potential role in LD build-up, while on the other hand, it also balances infection frequencies between populations (Equation (3)), which in turn will reduce the magnitude of $(f_{A}^{x} - f_{A}^{y}) (f_{B}^{x} - f_{B}^{y})$ in Equation (4). These conflicting forces mean the magnitude of population LD tends to be maximized when migration is neither too infrequent nor too frequent (Figure 2).

Figure 2

Download asset Open asset

The effect of migration upon LD at equilibrium depends upon the scale at which LD is measured.

Here, we show equilibrium LD in a metapopulation consisting of four populations. Two scenarios are shown. In the first scenario (solid lines), drug $A$ and drug $B$ are both prescribed in the same two populations while the other two populations receive no drugs, thus $cov (τ_{A}, τ_{B}) > 0$ ; this yields $cov (f_{A}, f_{B}) > 0$ and so positive population, average, and metapopulation LD, that is, $D^{x}, \bar{D}, D_{M} > 0$ . In the second scenario (dashed lines), drug $A$ is prescribed in two populations and drug $B$ is prescribed in the other two populations, thus $cov (τ_{A}, τ_{B}) < 0$ ; this yields $cov (f_{A}, f_{B}) < 0$ and negative population, average, and metapopulation LD, i.e., $D^{x}, \bar{D}, D_{M} < 0$ . Because we assume identical treatment rates and costs of resistance for either drug, in the second scenario all the populations have the same LD, whereas in the first scenario, since the drugs are prescribed unequally across populations, the LD observed in each of the two pairs of populations diverge. Specifically, populations experiencing greater selection due to increased drug prescription also have greater LD; this follows from the first term in Equation (4).

Metapopulation-level multidrug resistance

Now what happens to LD and MDR evolution at the metapopulation-level? Here we will use $\bar{X}$ to denote the metapopulation average of quantity $X^{x}$ , e.g., ${\bar{f}}_{A}$ is the average drug $A$ resistance in the metapopulation (see Materials and methods 'Metapopulation LD and MDR' for further details). Using this notation, then analogously to the population case, metapopulation LD is defined as $D_{M} \equiv {\bar{f}}_{A B} - {\bar{f}}_{A} {\bar{f}}_{B}$ . A more informative, but mathematically equivalent, description of metapopulation LD, however, is to define it in terms of the population variables as

D_{M} \equiv \bar{D} + cov (f_{A}, f_{B}),

that is, $D_{M}$ is the sum of the average population LD, $\bar{D}$ , and the spatial covariance between the frequencies of resistance to drugs $A$ and $B$ . Equation (5) shows that even if there is no population LD, that is, $D^{x} = 0$ and so $\bar{D} = 0$ , there may still be metapopulation LD; likewise, there may be population LD, $D^{x} \neq 0$ , but no metapopulation LD, $D_{M} = 0$ (Nei and Li, 1973; Ohta, 1982a; Feldman and Christiansen, 1974).

With this in mind, the change in frequency of infections resistant to drug $A$ (mutatis mutandis drug $B$ ) can be written

\begin{array}{ll} \frac{d {\bar{f}}_{A}}{d t} = \overset{d i r e c t s e l e c t i o n}{\overset{⏞}{{\bar{s}}_{A} {\bar{f}}_{A} (1 - {\bar{f}}_{A})}} + \overset{i n d i r e c t s e l e c t i o n}{\overset{⏞}{{\bar{s}}_{B} D_{M}}} + \overset{e p i s t a s i s}{\overset{⏞}{{\bar{s}}_{E} {\bar{f}}_{A} (1 - {\bar{f}}_{A}) \frac{{\bar{f}}_{A B}}{{\bar{f}}_{A}}}} \\ + \underset{m u t a t i o n a n d r e c o m b i n a t i o n}{\underset{⏟}{({\bar{μ}}_{A} + {\bar{ρ}}_{A}) (1 - {\bar{f}}_{A}) - ({\bar{μ}}_{a} + {\bar{ρ}}_{a}) {\bar{f}}_{A}}} + \underset{h e t e r o g e n e i t y i n^{'} b a s e l i n e^{'} g r o w t h}{\underset{⏟}{cov (r, f_{A})}} + \underset{h e t e r o g e n e i t y i n i n d i r e c t s e l e c t i o n}{\underset{⏟}{{\bar{f}}_{B} cov (s_{B}, \frac{f_{A B}}{f_{B}})}} . \end{array}

The first four terms in Equation (6) are the metapopulation-level analogues of the first four terms in Equation (3) and so share the same interpretation. The last two terms, however, arise due to spatial heterogeneity in ‘baseline’ growth and selection and so are the consequence of population structure. As these terms are zero in the absence of spatial heterogeneities, they will be our focus here.

First, spatial heterogeneity arises through differences in the ‘baseline’ per-capita growth (i.e. $r^{x} \neq r^{y}$ ) coupled with differences in the frequencies of drug $A$ resistant infections (i.e. $f_{A}^{x} \neq f_{A}^{y}$ ). This is the spatial covariance between ‘baseline’ per-capita growth and the frequency of drug $A$ resistant infections, $cov (r, f_{A})$ . In particular, more productive populations (larger $r^{x}$ ) will have a disproportionate effect on the change in drug $A$ resistance. For example, if more productive populations also have a greater frequency of drug $A$ resistance, then heterogeneity increases the population frequency of drug $A$ resistance. Heterogeneity in baseline growth could arise through a variety of mechanisms, such as availability of susceptible hosts, treatment rates differences, or pathogen traits (e.g. transmissibility and duration of carriage).

Second, spatial heterogeneity arises through differences in indirect selection for resistance to drug $B$ (i.e. $s_{B}^{x} \neq s_{B}^{y}$ ) coupled with differences in the probability that drug $B$ resistant infections are also doubly-resistant (i.e. $f_{A B}^{x} / f_{B}^{x} \neq f_{A B}^{y} / f_{B}^{y}$ ). This is the spatial covariance between selection on resistance to drug $B$ and the conditional probability that a drug $B$ resistant infection is doubly-resistant, $cov (s_{B}, f_{A B} / f_{B})$ . In particular, populations experiencing greater selection for resistance to one drug will have a disproportionate effect on the change in frequency of infections resistant to the other drug, whenever populations differ in frequency of doubly-resistant infections. As an example, if populations experiencing stronger selection for drug $B$ resistance also have a greater probability of drug $B$ -resistant infections being doubly-resistant, heterogeneity in indirect selection increases the frequency of drug $A$ resistance in the metapopulation.

Next, the dynamics of metapopulation LD can be written as

\begin{array}{ll} \frac{d D_{M}}{d t} = \overset{selection}{\overset{⏞}{({\bar{s}}_{A} - \bar{s} + {\bar{s}}_{B} - \bar{s}) D_{M}}} - \overset{m u t a t i o n a n d r e c o m b i n a t i o n}{\overset{⏞}{(\bar{μ} + \bar{ρ}) D_{M}}} + \overset{epistasis}{\overset{⏞}{{\bar{s}}_{E} {\bar{f}}_{a b} {\bar{f}}_{A B}}} \\ + \underset{h e t e r o g e n e i t y i n^{'} b a s e l i n e^{'} g r o w t h}{\underset{⏟}{cov (r, D) + coskew (r, f_{A}, f_{B})}} + \sum_{d \in {A, B}} \underset{h e t e r o g e n e i t y i n r e s i s t a n c e s e l e c t i o n}{\underset{⏟}{(1 - {\bar{f}}_{d}) {\bar{f}}_{d} cov (s_{d}, \frac{f_{A B}}{f_{d}})}}, \end{array}

where $coskew (r, f_{A}, f_{B})$ is the spatial coskewness between $r$ , $f_{A}$ , and $f_{B}$ and we have assumed population differences in mutation and recombination are negligible (see Materials and methods 'Metapopulation LD and MDR'). The first three terms in Equation (7) are the metapopulation level analogues of the first three terms of Equation (4) and so share the same interpretation. The last two terms, however, arise due to spatial heterogeneity in ‘baseline’ growth and selection and so will be our focus here.

First, spatial heterogeneity arises through spatial differences in the ‘baseline’ per-capita growth (i.e. $r^{x} \neq r^{y}$ ) coupled with spatial heterogeneities in LD (i.e. $D^{x} \neq D^{y}$ ) or resistance frequencies (the coskewness term). The logic of the first term is clear: when population LD differs, more productive populations will disproportionately contribute to metapopulation LD. For the second term, when populations covary in frequency of resistance to drug $A$ and $B$ , more productive populations will disproportionately contribute to the covariance, $cov (f_{A}, f_{B})$ and so disproportionately contribute to metapopulation LD (through the second term in Equation (5)).

Second, spatial heterogeneity arises through differences in selection for resistance ( $s_{d}^{x} \neq s_{d}^{y}$ ) coupled with differences in the proportion of drug $d$ resistant infections that are doubly-resistant ( $f_{A B}^{x} / f_{d}^{x} \neq f_{A B}^{y} / f_{d}^{y}$ ). The logic here is that populations experiencing stronger selection for resistance are more likely to see an increase in resistant infections. If this increase occurs disproportionately in doubly-resistant infections, then from Equation (1) metapopulation LD will increase, whereas if this increase occurs disproportionately in singly-resistant infections, metapopulation LD will decrease. The magnitude of this effect is scaled by ${\bar{f}}_{d} (1 - {\bar{f}}_{d})$ since selection cannot operate without genetic variation. In the absence of population LD, then $f_{A B}^{x} / f_{A}^{x} = f_{B}^{x}$ and $f_{A B}^{x} / f_{B}^{x} = f_{A}^{x}$ , and so if populations experiencing stronger selection for resistance to one drug also have a greater frequency of infections resistant to the other drug, metapopulation LD will increase. This could occur if, for example, some populations experience greater treatment rates.

As a final note, observe that in contrast to Equation (4), in Equation (7) the per-capita migration rates $m^{y \to x}$ are nowhere to be found. The reason for this is intuitive: as migration does not affect the total density of infecteds, nor the resistance status of an infection, it will not change the quantities ${\bar{f}}_{A B}$ , ${\bar{f}}_{A}$ , or ${\bar{f}}_{B}$ , and so cannot change metapopulation LD. As a consequence, migration only affects metapopulation LD indirectly by reducing differences in infection frequency between populations, thereby dampening the magnitude (and hence the effect) of $cov (r, D)$ , $cov (r, f_{d})$ , and $cov (s_{ℓ}, f_{A B} / f_{d})$ in Equation (7). It follows that, all else being equal, the magnitude of $D_{M}$ is a decreasing function of the per-capita migration rate, and so is maximized when migration is infrequent (Figure 2).

Modeling the dynamics of LD: why bother?

To this point, we have focused upon developing the LD perspective to provide a conceptual understanding of MDR evolution in structured populations. However, framing the LD perspective in terms of general quantities has meant this conceptual understanding is somewhat abstract. What we now wish to demonstrate, through the consideration of three scenarios, is how the LD perspective can be used to tackle practical problems. In the first scenario, we show how the LD perspective provides additional insight into a recent paper on the effect of spatial structure on equilibrium patterns of MDR. In the second scenario, we show how the LD perspective allows for an understanding of transient dynamics, and we apply this understanding to patterns of MDR observed in Streptococcus pneumoniae. In the third scenario, we show how the LD perspective generates practical insight into designing drug prescription strategies across populations, with a focus upon a hospital-community setting.

LD perspective explains equilibrium patterns of MDR

Understanding the patterns of MDR in structured populations was first tackled in an important paper by Lehtinen et al., 2019. The paper by Lehtinen et al., 2019 (see also Jacopin et al., 2020) focused upon MDR evolution in a metapopulation consisting of independent host populations (so migration is restricted, $m^{x \to y} \approx 0$ ). For example, each population could represent a different Streptococcus pneumoniae serotype maintained by serotype-specific host immunity (Henriques-Normark and Tuomanen, 2013; Cobey and Lipsitch, 2012; Lehtinen et al., 2017). Lehtinen et al., 2019 found that at equilibrium, population differences could lead to MDR over-representation ( $D_{M} > 0$ ), and that populations with a longer duration of pathogen carriage were more likely to exhibit MDR, a result they attributed to an increased likelihood of antibiotic exposure per carriage episode. Here, we show how employing the LD perspective: (i) reveals the evolutionary logic behind what populations differences can maintain metapopulation LD at equilibrium and (ii) using these insights allows us to build upon the results of Lehtinen et al., 2019 to understand how epidemiological factors other than duration of carriage can play an important role. For simplicity, we will assume that costs are additive (see Box 1), and so there is no epistasis (i.e. $s_{E}^{x} = 0$ ), but as this differs from Lehtinen et al., 2019 who use multiplicative costs, we discuss this assumption in more depth in Materials and methods 'Equilibrium analysis of metapopulation consisting of independent populations'.

To maintain metapopulation LD at equilibrium, there needs to be at minimum some mechanism maintaining metapopulation resistance diversity, otherwise $D_{M} = 0$ . There are variety of ways in which this could occur (Lipsitch et al., 2009; Colijn et al., 2010; Davies et al., 2019; Lehtinen et al., 2017; Jacopin et al., 2020; Krieger et al., 2020), but Lehtinen et al., 2017, Lehtinen et al., 2019 assume it is due to some variation among populations in the conditions favoring resistance evolution. This mechanism maintains diversity at the scale of the metapopulation but leads to the fixation or the extinction of drug resistance locally. Thus $D^{x} = 0$ , and it follows from Equation (5) that $D_{M} = cov (f_{A}, f_{B})$ . Therefore, in order for metapopulation LD to exist, $f_{A}$ and $f_{B}$ must covary across populations. Specifically, whenever $f_{A}^{x}$ and $f_{B}^{x}$ (or their dynamical equations, Equation 3), are uncorrelated, the metapopulation will be in linkage equilibrium. From Equation (3) we see that if the additive selection coefficients, $s_{A}^{x}$ and $s_{B}^{x}$ , are uncorrelated, then so too are the dynamics of $f_{A}^{x}$ and $f_{B}^{x}$ , and so $cov (f_{A}, f_{B}) = 0$ . Hence only when population differences generate correlations between the selection coefficients will they generate LD.

Using this insight, why are populations with a longer duration of carriage associated with MDR (Lehtinen et al., 2019)? And should we expect associations between MDR and any other population attributes? Our primary focus is whether (and how) the selection coefficients are correlated. Letting $Δ z_{k}^{x}$ be the contribution of trait $z$ to the additive selection coefficient for resistance to drug $d$ in population $x$ (e.g. $Δ β_{A}^{x} = β_{A b}^{x} - β_{a b}^{x}$ ), then it is straightforward to compute (see Materials and methods 'Equilibrium analysis of metapopulation consisting of independent populations'),

\begin{array}{ll} s_{A}^{x} = Δ β_{A}^{x} S^{x} - Δ α_{A}^{x} + τ_{A}^{x}, \\ s_{B}^{x} = Δ β_{B}^{x} S^{x} - Δ α_{B}^{x} + τ_{B}^{x} . \end{array}

where we have used slightly different notation from Lehtinen et al., 2019. Now, consider a scenario in which both the treatment rates and the parameters controlling the (additive) costs of resistance are uncorrelated (i.e. $Δ β_{d}^{x} = Δ β_{d}$ , $Δ α_{d}^{x} = Δ α_{d}$ and $τ_{d}^{x} = τ_{d}$ ); this is one of the scenarios presented in Figure 4 of Lehtinen et al., 2019, with the key difference that they considered ‘multiplicative’ rather than ‘additive’ costs. From Equation (8), the only remaining source of correlation is susceptible density, $S^{x}$ , which plays a role whenever there are explicit transmission costs, $Δ β_{d} < 0$ . Although Equation (8) always holds, in keeping with Lehtinen et al., 2019 if we focus upon the equilibrium case, $S^{x}$ will be determined by pathogen traits such as transmission and duration of carriage, such that ‘fitter’ populations (i.e. those in which pathogens are more transmissible or have longer duration of carriage) will more substantially deplete susceptibles. By reducing $S^{x}$ , ‘fitter’ populations lower the transmission costs for resistance to either drug, and so double-resistance is more likely to be selectively advantageous, even when treatment rates are uncorrelated. In turn, this over-representation of doubly-resistant infections will generate metapopulation LD.

Thus, when costs are ‘additive’ (Box 1), although variation in duration of carriage can lead to MDR evolution and LD through its effect upon susceptible density (Figure 3a), it is neither necessary (the same pattern can be produced by variation in transmissibility; Figure 3b) nor sufficient (variation in duration of carriage has no effect without explicit transmission costs, Figure 3c). More broadly, if there are more than two drugs, then provided that there are explicit transmission costs for resistance to each drug, susceptible density will generate a correlation between all the selection coefficients, which in turn will yield the pattern of ‘nestedness’ observed by Lehtinen et al., 2019. What is critical for this effect to be prominent, however, is (i) the existence of population differences in susceptible density, and (ii) the costs of resistance (i.e. $Δ β_{d}$ ), are large enough so as to ensure a strong correlation amongst selection coefficients.

Figure 3

Download asset Open asset

Duration of carriage is one of many potential explanations for MDR over-representation at equilibrium.

When costs are additive and there is no epistasis, variation in duration of carriage across independent populations can lead to linkage disequilibrium (subplot a), but it is neither necessary (b), nor sufficient (c). We simulate 1000 populations (blue bars), each consisting of 20 independent populations in which treatment rates for each population are randomly chosen to be either $τ_{\max} = 0.075$ or $τ_{\min} = 0.025$ with equal probability while simultaneously satisfying $cov (τ_{A}, τ_{B}) = 0$ . The solid blue line is the mean LD across the simulations for each scenario. In subplot a, duration of carriage varies across populations and there are transmission resistance costs; in subplot b, transmission varies and there are transmission resistance costs; while in subplot c, duration of carriage varies and there are no transmission costs. These simulations diverge slightly from those of Lehtinen et al., 2019 in that their model always includes epistasis (see Materials and methods 'Equilibrium analysis of metapopulation consisting of independent populations'), whereas here we only consider non-epistatic scenarios.

LD perspective explains transient patterns of MDR

The predictions of Lehtinen et al., 2019 were used to explain the patterns of MDR observed in surveillance data. One of these data sets was a surveillance study that documented both the serotype as well as antibiotic resistance to a number of different drugs in S. pneumoniae infections sampled in Maela, Northern Thailand (Turner et al., 2012; Lehtinen et al., 2019). Although the prediction of positive (metapopulation) LD was met for most drug combinations (Lehtinen et al., 2019), inspection of the data set reveals significant serotype LD (Figure 4). This is notable because, as we have detailed above, at equilibrium the simplest version of the model used in the previous section will result in each serotype being in linkage equilibrium, $D^{x} = 0$ . How can we reconcile these conflicting observations? Although there are various possible explanations (e.g. an additional mechanism capable of maintaining diversity within-serotype), here we focus upon relaxing the assumption that the metapopulation is at equilibrium. That is, we are interested in whether long-term transient dynamics unfolding over months and years could plausibly suggest an alternative explanation for the observed serotype LD.

Figure 4

Download asset Open asset

Linkage disequilibrium for different drug pairs in *Streptococcus pneumoniae*.

Data is from the Maela surveillance data set of Lehtinen et al., 2019; Turner et al., 2012. The light red circles are the observed serotype LD, $D^{x}$ , the dark red circles are the average LD across serotypes, $\bar{D}$ , while the blue circles are the metapopulation LD, $D_{M}$ . We have restricted the data to serotypes involving 100 or more samples (serotypes 14, 6A/C, 6B, 15B/C, 19F, 23F). The drugs considered are: A = chloramphenicol, B = clindamycin, C = erythromycin, D = penicillin, E = sulphatrimethoprim, and F = tetracycline.

To do so, consider a metapopulation consisting of independent serotypes, differing in their transmissibility and duration of carriage (as in the model of Lehtinen et al., 2017; Lehtinen et al., 2019). Assume that there is no epistasis and that the additive selection coefficients take the form of Equation (8), where the parameters $Δ β_{d}^{x}$ , $Δ α_{d}^{x}$ and $τ_{d}^{x}$ do not depend upon serotype $x$ (Materials and methods 'Transient dynamics and MDR in streptococcus pneumoniae'). Suppose that initially the metapopulation is treated exclusively with drug $A$ at sufficiently high rates such that resistance to drug $A$ goes to fixation in each serotype, that is, ${\bar{f}}_{A} \to 1$ and ${\bar{f}}_{B} \to 0$ , and so $D_{M} = 0$ . Now suppose at time $t = 1000$ months that drug $B$ is ‘discovered’ and subsequently prescribed at a high rate in the metapopulation, while owing to its reduced efficacy, prescription of drug $A$ is reduced. Although the treatment rates do not vary by serotype, serotype differences in transmissibility and duration of carriage mean that the changes to treatment rates will differentially affect serotype density, which in turn will differentially affect the serotype-specific availability of susceptible hosts, $S^{x}$ . Since the serotype-specific selection coefficients, $s_{A}^{x}$ and $s_{B}^{x}$ , and baseline per-capita growth, $r^{x}$ , directly depend upon $S^{x}$ , the variation in $S^{x}$ introduces heterogeneity in Equation (7), which in turn generates metapopulation LD. Because the selection coefficients are positively correlated (due to the shared dependence upon $S^{x}$ ), the metapopulation LD generated will be positive, that is, $D_{M} > 0$ (Figure 5a,d). From Equation (7), once metapopulation LD is generated, it will be amplified by directional selection (first term of Equation 7) which is initially positive since resistance to drug $B$ is favored; this leads to a rapid build up of $D_{M}$ (Figure 5a,d). However, this initial increase in $D_{M}$ is transient; for this particular choice of parameter values, at equilibrium $D_{M} \to 0$ . Crucially, however, the changes to $D_{M}$ can unfold over a very long time (here the time units are months), such that surveillance data would detect little change in the metapopulation dynamics and so suggest a population roughly in equilibrium.

Figure 5

Download asset Open asset

Transient dynamics coupled with epistasis can explain patterns of serotype LD in *Streptococcus pneumoniae*.

In all simulations, serotypes differ based upon duration of carriage and transmissibility. At $t = 0$ , the pathogen is sensitive to both drugs; however, as hosts are initially treated with drug $A$ at a rate of $τ_{A} = 0.12$ per month, resistance to drug $A$ emerges and fixes in all serotypes. At $t = 1000$ (months), drug $B$ is introduced, and drug $A$ prescription reduced, $(τ_{A}, τ_{B}) = (0.07, 0.1)$ (note that the drugs are never prescribed in combination, $τ_{A B} = 0$ ). In the first column, there is no epistasis, thus although metapopulation LD builds up, serotype LD does not. In the second column, there is negative epistasis, which generates negative serotype LD. In the third column, there is positive epistasis which produces positive serotype LD. The thin lines denote the within-serotype dynamics, while the thick lines denote the metapopulation dynamics. In all cases, at equilibrium both the serotypes and the metapopulation will be in linkage equilibrium, however, transient LD can occur on sufficiently long timescales so as to appear permanent (see Materials and methods 'Transient dynamics and MDR in streptococcus pneumoniae' for more details).

Although this scenario can lead to considerable (transient) metapopulation LD, there is still nothing generating serotype LD. Our analysis of Equation (4) revealed two possible (deterministic) mechanisms capable of generating population (serotype) LD. First, migration between populations can lead to metapopulation LD spilling over into population LD. In this example, ‘migration’ between serotypes would correspond to serotype ‘switching’ (Croucher et al., 2015), whereby infections exchange serotypes through recombination. However, this is an unlikely explanation as the rate of serotype switching would have to be unrealistically large for serotype LD to be substantially altered. The second term in Equation (4) capable of generating serotype LD is epistasis, which will generate same sign serotype LD (Felsenstein, 1965; Lewontin, 1964; Lewontin and Kojima, 1960). Indeed, in the model considered, negative epistasis, $s_{E}^{x} < 0$ , generates transient negative serotype LD (Figure 5b,e), while positive epistasis generates transient positive serotype LD (Figure 5c,f; Materials and methods 'Transient dynamics and MDR in streptococcus pneumoniae'). Notably, although negative epistasis produces negative serotype LD (and so $\bar{D} < 0$ ), at the scale of the metapopulation this effect is swamped by the positive covariance in frequency of resistance and so metapopulation LD is positive, $D_{M} > 0$ . (Figure 5e).

Thus, transient dynamics coupled with epistasis could provide a potential explanation for the significant within-serotype LD observed in S. pneumoniae (Figure 4). More generally, the potential complexity of competing selective pressures associated with multilocus dynamics can lead to prolonged, but transient, polymorphisms and LD, and so surveillance data showing limited temporal change in resistance frequency should be treated cautiously and not assumed to be due to a stable equilibrium.

LD perspective helps identify drug prescription strategies limiting the evolution of MDR

Understanding the evolutionary consequences of different antibiotic prescription strategies across populations can have practical relevance for public health. The populations of interest could correspond to physically distinct groups such as a hospital and its broader community, or different geographical regions (e.g. countries). From a public health perspective, when considering different prescription strategies, a variety of factors must be considered, but in general the goal is to successfully treat as many people as possible, thereby reducing the total burden (Bonhoeffer et al., 1997; Abel zur Wiesch et al., 2014). In this circumstance, the LD in the metapopulation and/or populations can provide important information about the likelihood of treatment success. In particular, for a given population frequency of drug $A$ and drug $B$ resistance, negative LD (MDR under-representation) increases the likelihood that if treatment with one drug fails (due to resistance), treatment with the other drug will succeed. On the other hand, positive LD (MDR over-representation) increases the likelihood of treatment failure, since a greater proportion of resistant infections are doubly-resistant and so cannot be successfully treated with either drug. Equations (4) and (7) show that to generate negative LD, drugs should be deployed in a population specific fashion, that is, drug $A$ should be restricted to some populations and drug $B$ restricted to the remaining populations (see also Lehtinen et al., 2019; Day and Gandon, 2012; Jacopin et al., 2020). Doing so will create a negative covariance in selection, such that resistance to drug $A$ (resp. drug $B$ ) will be favored in some populations and disfavored in the others. This negative covariance in selection will give rise to negative LD and MDR under-representation (Figure 2).

As an application of this principle, consider two populations connected by migration, corresponding to a ‘community’ and a much smaller ‘hospital’. Drug prescription occurs at a fixed (total) rate in each population, while the prescription rate is much higher in the hospital (see Materials and methods 'Contrasting drug prescription strategies in a hospital-community setting'). Consider three antibiotic prescription strategies: (i) drugs can be randomly prescribed to individuals (mixing); (ii) drugs can be prescribed exclusively in combination; or (iii) prescription of drug $A$ and $B$ can be asynchronously rotated between the hospital and community, that is, if the hospital uses drug $A$ then the community uses drug $B$ , and vice versa (cycling). As both drugs are prescribed at higher rates in the hospital than the community, both mixing and combination generate a positive covariance in selection across populations, producing positive LD and MDR over-representation (see Equation (4); Figure 6). Thus, over the short- and long-term, mixing and combination produce similar results: doubly-resistant infections are favored, while singly-resistant infections are disfavored (Figure 6). Now consider cycling. When drugs are rotated rapidly between populations, infections in either population are likely to be exposed to both drugs. Because prescription rates are higher in the hospital, this effectively creates a positive covariance in selection (i.e. cycling behaves like mixing) and so when resistance emerges, infections tend to be doubly-resistant (MDR over-representation). When drugs are rotated less frequently, infections are more likely to be exposed to a single drug, creating a negative covariance in selection across populations. In this circumstance, although single resistance can emerge at lower treatment rates then when rotations are more frequent, the negative LD produced by the negative covariance in selection inhibits the emergence of double-resistance (MDR under-representation; Figure 6).

Figure 6

Download asset Open asset

Different antibiotic prescription strategies generate different patterns of LD at equilibrium.

Here, we focus upon a population divided into a community and a hospital. Individuals enter the hospital at a fixed rate and spend a fifth of the time in the hospital that it takes to naturally clear a sensitive infection. The hospital/community size split corresponds to 20 beds per 1000 people, while individuals in the hospital receive antibiotics at 15x the rate they do in the community. We integrate system (2) until equilibrium is reached; the final state of the system is what is shown. For cycling, we compute the average state over the last two rotations (i.e., over the last period, $T$ ; in this case $T = 100$ ). In panel a, we show the metapopulation LD, $D_{M}$ for the three treatment scenarios (combination, mixing, cycling). Combination and mixing generate identical LD in this example. In panel b we show the frequency of infections in the metapopulation resistant to drug $d$ , ${\bar{f}}_{d}$ (for our choice of parameters, ${\bar{f}}_{A} = {\bar{f}}_{B}$ ; Materials and methods 'Contrasting drug prescription strategies in a hospital-community setting'), and doubly-resistant, ${\bar{f}}_{A B}$ , for each scenario. Note that for mixing and combination treatments (solid curves), ${\bar{f}}_{A} = {\bar{f}}_{B} = {\bar{f}}_{A B}$ , whereas cycling (dashed curves) leads to singly-resistant infections at low treatment rates (see Materials and methods 'Contrasting drug prescription strategies in a hospital-community setting').

These results emphasize an important trade-off: delaying the evolution of MDR (e.g. by decreasing time between rotations) promotes the evolution of single drug resistance, whereas delaying the evolution of single drug resistance promotes the evolution of MDR. This is logical: when we maintain a constant treatment rate per individual, decreasing selection for the ‘generalist’ strategy (MDR) necessarily increases selection for the ‘specialist’ strategy (single drug resistance) (Wilson and Yoshimura, 1994). Thus, cycling can either be the best, or worst, option for single drug resistance (Beardmore et al., 2017), but critically, this has concomitant effects for MDR (see also Figure 6). Indeed, mixing, combination and cycling have been exhaustively compared in the context of single drug resistance (e.g. Bonhoeffer et al., 1997; Lipsitch et al., 2000; Bergstrom et al., 2004; Beardmore and Pena-Miller, 2010; Beardmore et al., 2017; Abel zur Wiesch et al., 2014; Tepekule et al., 2017); yet these studies largely ignored the consequences for MDR evolution. Our analysis suggests controlling for single drug resistance will have important consequences for MDR, and so it should not be considered in isolation. More generally, whether it is optimal to either delay single drug resistance or prevent MDR will depend upon what metric is used to evaluate what constitutes a ‘success’ or ‘failure’.

Discussion

The evolution of multidrug-resistant pathogens is a pressing health concern and is a topic which is increasingly gaining attention from evolutionary biologists and mathematical modellers alike. However, the typical process in studying the problem of MDR is to introduce a model of the form of (2), and then either proceed to a numerical analysis of these equations or simplify the model further by neglecting the dynamics of double resistant infections (Bergstrom et al., 2004; Bonhoeffer et al., 1997; Beardmore et al., 2017). This is because models of MDR evolution rapidly become intractable, a problem which is particularly acute when incorporating aspects of population structure. Here, we have argued that a more insightful and simplifying approach is the ‘linkage disequilibrium perspective’: after specifying the model of interest, as in (2), it is desirable to transform the model into the form of Equations (3), (4), (6), and (7), which brings to the forefront the role played by linkage disequilibrium for MDR evolution in structured populations. The LD perspective is particularly useful for analyzing and understanding transient evolutionary dynamics (Figure 5), which cannot be understood by, for example, invasion analysis.

Our analysis emphasizes that metapopulation structure alone can generate and maintain LD (and so MDR), even in the absence of epistasis (Nei and Li, 1973; Ohta, 1982a; Li and Nei, 1974; Slatkin, 1975). Since in natural populations metapopulation structure is often hidden (e.g. Rosen et al., 2015), patterns of MDR should not be assumed to be due to epistasis, even if no structure is readily apparent. Moreover, caution must be taken when measuring LD (and MDR) at a particular scale, as doing so can lead to erroneous conclusions: even if the metapopulation is in linkage equilibrium, $D_{M} = 0$ , the populations need not be, $D^{x} \neq 0$ , and vice versa (Equation (5)), while in more extreme cases, population and metapopulation LD can be of opposite sign (Figure 5b,e). These are not merely esoteric points; the presence or absence of LD (and MDR), and its source (epistasis or metapopulation structure) is critically important. For example, when MDR is due to metapopulation structure rather than epistasis, prescribing drugs across different populations so as to create a negative covariance in selection can reduce the prevalence of MDR (Figure 6; Day and Gandon, 2012; Jacopin et al., 2020; Lehtinen et al., 2019), while distinguishing between population and metapopulation LD can provide additional insight toward evaluating hypotheses (Figures 4 and 5).

Our analysis assumed that the evolutionary dynamics were deterministic, thus neglecting the influence of stochasticity. However, it is widely appreciated in population genetics that stochasticity can play an important role in multilocus dynamics. For example, LD can be generated through genetic drift (Hill and Robertson, 1966; Barton, 1995; Lenormand and Otto, 2000; Otto and Barton, 2001; Keightley and Otto, 2006; Martin et al., 2006), which in turn can interfere with the strength of selection (Hill and Robertson, 1966; Neher and Shraiman, 2011; Slatkin, 2008). Similarly, the (random) genetic background a rare mutation finds itself upon is critically important for its success (Kouyos et al., 2006; Gillespie, 2000; Neher, 2013), and in finite populations this alone can generate LD. However, little has been done to relate these results to evolutionary epidemiology, or to understand how epidemiological feedbacks can influence their predictions. The little work to date has relied upon complex simulations (e.g. Althaus and Bonhoeffer, 2005; Kouyos et al., 2009), which necessarily sacrifice general insight for specificity. Thus, the role of stochasticity in the evolution of MDR remains an area in which further investigation is warranted.

Understanding the evolution of MDR is a research topic of pressing concern. Here, we have argued that using the linkage disequilibrium perspective leaves us better equipped to determine what factors are responsible for generating MDR, and their generality. Moreover, taking such an approach leads to a more straightforward comparison with existing models and results.

Materials and methods

Here, we provide more comprehensive details on the analysis presented in the main text. We start by deriving the general epidemiological model for the dynamics of the different strains which are characterised by their multilocus genotype (Materials and methods 'Model derivation'). We then convert this model into an equivalent system which tracks the dynamics of allele frequencies at each locus, and the LD at the population level (Materials and methods 'Population LD and MDR'), before considering the set of equations for the dynamics of allele frequencies at each locus and the LD at the metapopulation level (Materials and methods 'Metapopulation LD and MDR').

We conclude by providing a detailed mathematical analysis of the three examples presented in the main text: (1) using the LD perspective to explain equilibrium patterns of MDR (Materials and methods 'Equilibrium analysis of metapopulation consisting of independent populations'); (2) using the LD perspective to explain transient patterns of MDR (Materials and methods 'Transient dynamics and MDR in streptococcus pneumoniae'); and (3) applying the LD perspective to identify drug prescription strategies limiting MDR evolution (Materials and methods 'Contrasting drug prescription strategies in a hospital-community setting').

Model derivation

Request a detailed protocol

Our focus is on an asymptomatically carried bacteria species in a metapopulation consisting of $N$ populations. Focus upon an arbitrarily chosen population $x$ . Let $S^{x}$ and $I_{i j}^{x}$ denote the density of susceptible hosts and $i j$ -infections, respectively, at time $t$ , where $i$ indicates if the infection is resistant ( $i = A$ ) or not ( $i = a$ ) to drug $A$ and $j$ indicates if the infection is resistant ( $j = B$ ) or not ( $j = b$ ) to drug $B$ . Susceptible hosts contract $i j$ -infections at a per-capita rate $β_{i j}^{x} I_{i j}^{x}$ , where $β_{i j}^{x}$ is a rate constant, while $i j$ -infections are naturally cleared at a per-capita rate $α_{i j}^{x}$ . Hosts in population $x$ are treated with antibiotics $A$ , $B$ , or both in combination, at per-capita rates $τ_{A}^{x}$ , $τ_{B}^{x}$ , and $τ_{A B}^{x}$ , respectively. Hosts move from population $x$ to population $y$ at a per-capita rate $m^{x \to y}$ .

The resistance profile of an infection changes through two processes. First, there may be de novo mutation, and so let $μ_{i}^{x}$ be the per-capita rate at which an infection in population $x$ acquires allele $i$ through mutation. Second, a $i j$ -infection may be super-infected by a $k ℓ$ -strain (Day and Gandon, 2012); in this circumstance recombination may occur. Specifically, $k ℓ$ -strains are transmitted to $i j$ -infections at rate $β_{k ℓ}^{x} I_{k ℓ}^{x} I_{i j}^{x}$ , whereupon with probability $σ$ super-infection occurs. In the event of super-infection, with probability $1 - ρ$ , recombination does not occur, in which case with equal probability the $i j$ -infection either remains unchanged or becomes a $k ℓ$ -infection. With probability $ρ$ , recombination does occur, in which case with equal probability the $i j$ -infection becomes either an $i ℓ$ - or $k j$ -infection. Because our focus is upon the role of population structure, we do not allow for co-infection or within-host competitive differences based upon resistance profiles (e.g. Davies et al., 2019) but these are straightforward extensions. Moreover, at this stage, we do not make any further specification of the dynamics of uninfected hosts, be they susceptible or recovered, as doing so is not essential for a qualitative understanding of MDR evolution.

Rather than immediately writing down the set of differential equations corresponding to these epidemiological assumptions, we instead group the terms based upon the four biological processes that are occurring. In particular, the change in $I_{i j}^{X}$ can be written as the sum of:

The net change due to mutation, denoted $ϕ μ_{i j}^{x}$ . As an example, focus upon the change in $A b$ -infections in population $x$ due to mutation, $ϕ μ_{A b}^{x}$ . These infections can increase through mutation in one of two ways: (i) $a b$ -infections acquiring allele $A$ at rate $μ_{A}^{x} I_{a b}^{x}$ or (ii) $A B$ -infections acquiring allele $b$ at rate $μ_{b}^{x} I_{A B}^{x}$ . On the other hand, $I_{A b}^{x}$ infections are lost due to mutation whenever they (i) acquire allele $a$ at a per-capita rate $μ_{a}^{x}$ , or (ii) acquire allele $B$ at a per-capita rate $μ_{B}^{x}$ . Combining this information gives the change in $A b$ -infections in population $x$ as
$ϕ μ_{A b}^{x} = μ_{A}^{x} I_{a b}^{x} + μ_{b}^{x} I_{A B}^{x} - (μ_{a}^{x} + μ_{B}^{x}) I_{A b}^{x},$
which is mathematically equivalent to
$ϕ μ_{A b}^{x} = μ_{A}^{x} (I_{a b}^{x} + I_{A b}^{x}) + μ_{b}^{x} (I_{A b}^{x} + I_{A B}^{x}) - μ^{x} I_{A b}^{x},$
where $μ^{x} \equiv μ_{a}^{x} + μ_{A}^{x} + μ_{b}^{x} + μ_{B}^{x}$ is the per-capita mutation rate in population $x$ . The only difference between the two formulations is interpretation: Equation (9) shows only mutations which lead to a change in state, whereas Equation (10) shows all possible mutations, even those which do not. This is why the per-capita loss term, $μ^{x}$ , in (10) can be considered the total per-capita mutation rate in population $x$ . More generally, we can write $ϕ μ_{i j}^{x}$ as
$ϕ μ_{i j}^{x} \equiv μ_{i}^{x} (I_{a j}^{x} + I_{A j}^{x}) + μ_{j}^{x} (I_{i b}^{x} + I_{i B}^{x}) - μ^{x} I_{i j}^{x} .$
The net change due to recombination, denoted $ϕ ρ_{i j}^{x}$ . Specifically, let $ρ_{i}^{x}$ be the per-capita rate at which infections gain allele $i$ through recombination. For example, consider $ρ_{A}^{x}$ . In particular, $i j$ -infections are challenged by strains carrying allele $A$ at rate $(β_{A b}^{x} I_{A b}^{x} + β_{A B}^{x} I_{A B}^{x}) I_{i j}^{x}$ . With probability $σ$ , a superinfection event occurs. Given an superinfection event, with probability $ρ$ recombination happens, in which case with probability $1 / 2$ the recombinant strain $A j$ will replace the $i j$ -infection. Thus
$ρ_{A}^{x} = ρ \frac{σ}{2} (β_{A b}^{x} I_{A b}^{x} + β_{A B}^{x} I_{A B}^{x}),$
and $i j$ -infections acquire allele $A$ in population $x$ at rate $ρ_{A}^{x} I_{i j}^{x}$ . Therefore, the change in $i j$ -infections in population $x$ due to recombination is
$ϕ ρ_{i j}^{x} \equiv ρ_{i}^{x} (I_{a j}^{x} + I_{A j}^{x}) + ρ_{j}^{x} (I_{i b}^{x} + I_{i B}^{x}) - ρ^{x} I_{i j}^{x}$
where $ρ^{x}$ is the per-capita rate of recombination in population $x$ , that is,
$ρ^{x} \equiv ρ σ \sum_{k ℓ} β_{k ℓ}^{x} I_{k ℓ}^{x} = ρ_{a}^{x} + ρ_{A}^{x} + ρ_{b}^{x} + ρ_{B}^{x} .$
The net change due to host migration between populations,
$- \sum_{y = 1}^{N} m^{x \to y} I_{i j}^{x} + \sum_{y = 1}^{N} m^{y \to x} I_{i j}^{y} .$
The net change due to per-capita growth,
$r_{i j}^{x} \equiv β_{i j}^{x} S^{x} - α_{i j}^{x} - 𝟏_{a} (i) τ_{A}^{x} - 𝟏_{b} (j) τ_{B}^{x} - (1 - 𝟏_{A} (i) 𝟏_{B} (j)) τ_{A B}^{x} - (1 - ρ) \frac{σ}{2} \sum_{k ℓ} (β_{k ℓ}^{x} - β_{i j}^{x}) I_{k ℓ}^{x},$
where $𝟏_{i} (j)$ is an indicator variable and is equal to 1 if $i = j$ and 0 otherwise.

With these four processes in hand, the dynamics of infection densities are given by the system of $4 N$ differential equations.

\frac{d I_{i j}^{x}}{dt} = ϕ μ_{i j}^{x} + ϕ ρ_{i j}^{x} - \sum_{y = 1}^{N} (m^{x \to y} I_{i j}^{x} - m^{y \to x} I_{i j}^{y}) + r_{i j}^{x} I_{i j}^{x}, x = 1, 2, \dots, N, i \in {a, A}, j \in {b, B} .

Population LD and MDR

Request a detailed protocol

In what follows, we provide more details for the calculations of population LD and MDR. First, we define the following frequencies of infections in population $x$ as

f_{A}^{x} = \frac{\sum_{j} I_{A j}^{x}}{I^{x}}, f_{B}^{x} = \frac{\sum_{i} I_{i B}^{x}}{I^{x}}, and f_{i j}^{x} = \frac{I_{i j}^{x}}{I^{x}},

where $I^{x} = \sum_{i j} I_{i j}^{x}$ is the total density of infections in population $x$ . Using these definitions, the standard measure of linkage disequilibrium in population $x$ is

D^{x} = f_{A B}^{x} - f_{A}^{x} f_{B}^{x},

which is mathematically equivalent to

D^{x} = f_{A B}^{x} f_{a b}^{x} - f_{A b}^{x} f_{a B}^{x} .

The three dynamical equations of interest for studying MDR in population $x$ are

\begin{array}{ll} \frac{d f_{A}^{x}}{d t} = s_{A}^{x} f_{A}^{x} (1 - f_{A}^{x}) + s_{B}^{x} D^{x} + s_{E}^{x} f_{A}^{x} (1 - f_{A}^{x}) \frac{f_{A B}^{x}}{f_{A}^{x}} + (μ_{A}^{x} + ρ_{A}^{x}) (1 - f_{A}^{x}) - (μ_{a}^{x} + ρ_{a}^{x}) f_{A}^{x} - \sum_{y = 1}^{N} m^{y \to x} \frac{I^{y}}{I^{x}} (f_{A}^{x} - f_{A}^{y}), \\ \frac{d f_{B}^{x}}{d t} = s_{B}^{x} f_{B}^{x} (1 - f_{B}^{x}) + s_{A}^{x} D^{x} + s_{E}^{x} f_{B}^{x} (1 - f_{B}^{x}) \frac{f_{A B}^{x}}{f_{B}^{x}} + (μ_{B}^{x} + ρ_{B}^{x}) (1 - f_{B}^{x}) - (μ_{b}^{x} + ρ_{b}^{x}) f_{B}^{x} - \sum_{y = 1}^{N} m^{y \to x} \frac{I^{y}}{I^{x}} (f_{B}^{x} - f_{B}^{y}), \\ \frac{d D^{x}}{dt} = (s_{A}^{x} - s^{x} + s_{B}^{x} - s^{x}) D^{x} - (μ^{x} + ρ^{x}) D^{x} + s_{E}^{x} f_{A B}^{x} f_{a b}^{x} - \sum_{y = 1}^{N} m^{y \to x} \frac{I^{y}}{I^{x}} (D^{x} - D^{y} - (f_{A}^{x} - f_{A}^{y}) (f_{B}^{x} - f_{B}^{y})) . \end{array}

System (19) contains a number of quantities that we now define in more detail. First, the (additive) selection coefficient for resistance to drugs $A$ and $B$ in population $x$ are defined as

s_{A}^{x} = r_{A b}^{x} - r_{a b}^{x} and s_{B}^{x} = r_{a B}^{x} - r_{a b}^{x},

respectively, while epistasis in population $x$ is $s_{E}^{x} = r_{A B}^{x} + r_{a b}^{x} - r_{A b}^{x} - r_{a B}^{x}$ . It follows that we can write each of the per-capita growth rates, $r_{i j}^{x}$ , as

r_{i j}^{x} = r^{x} + 𝟏_{A} (i) s_{A}^{x} + 𝟏_{B} (j) s_{B}^{x} + 𝟏_{A} (i) 𝟏_{B} (j) s_{E}^{x} .

This is why $r_{a b}^{x} = r^{x}$ can be thought of as ‘baseline’ per-capita growth. We define the average selection for resistance in population $x$ as

s^{x} = s_{A}^{x} f_{A}^{x} + s_{B}^{x} f_{B}^{x} + s_{E}^{x} f_{A B}^{x} .

Note that the average per-capita growth rate in population $x$ is therefore $r^{x} + s^{x}$ , that is, average per-capita growth rate is the sum of the ‘baseline’ per-capita growth rate and the average selection for resistance.

Metapopulation LD and MDR

Next, consider metapopulation (or total) LD and MDR. First, let $p^{x} = I^{x} / \sum_{j = 1}^{N} I^{j}$ be the fraction of total infections in population $x$ . Then the metapopulation quantities equivalent to Equations (16) are

{\bar{f}}_{i} = \sum_{x = 1}^{N} p^{x} f_{i}^{x} and {\bar{f}}_{i j} = \sum_{x = 1}^{N} p^{x} f_{i j}^{x} .

The standard measure of linkage disequilibrium at the level of the metapopulation is

D_{M} = {\bar{f}}_{A B} - {\bar{f}}_{A} {\bar{f}}_{B} .

which in terms of the population level variables is

D_{M} \equiv \sum_{x = 1}^{N} p^{x} D^{x} + \sum_{x = 1}^{N} p^{x} f_{A}^{x} f_{B}^{x} - (\sum_{x = 1}^{N} p^{x} f_{A}^{x}) (\sum_{x = 1}^{N} p^{x} f_{B}^{x}) = \bar{D} + cov (f_{A}, f_{B})

where $\bar{D}$ is the average population LD and $cov (f_{A}, f_{B})$ is the spatial covariance between frequency of resistance to drug $A$ and frequency of resistance to drug $B$ .

Using these variables, the three dynamical equations for studying metapopulation MDR are

\begin{array}{ll} \frac{d {\bar{f}}_{A}}{d t} = {\bar{s}}_{A} {\bar{f}}_{A} (1 - {\bar{f}}_{A}) + {\bar{s}}_{B} D_{M} + {\bar{s}}_{E} {\bar{f}}_{A} (1 - {\bar{f}}_{A}) \frac{{\bar{f}}_{A B}}{{\bar{f}}_{A}} \\ + ({\bar{μ}}_{A} + {\bar{ρ}}_{A}) (1 - {\bar{f}}_{A}) - ({\bar{μ}}_{a} + {\bar{ρ}}_{a}) {\bar{f}}_{A} + cov (r, f_{A}) + {\bar{f}}_{B} cov (s_{B}, \frac{f_{A B}}{f_{B}}), \\ \frac{d {\bar{f}}_{B}}{d t} = {\bar{s}}_{B} {\bar{f}}_{B} (1 - {\bar{f}}_{B}) + {\bar{s}}_{A} D_{M} + {\bar{s}}_{E} {\bar{f}}_{B} (1 - {\bar{f}}_{B}) \frac{{\bar{f}}_{A B}}{{\bar{f}}_{B}} \\ + ({\bar{μ}}_{B} + ρ_{B}) (1 - {\bar{f}}_{B}) - ({\bar{μ}}_{b} + {\bar{ρ}}_{b}) {\bar{f}}_{B} + cov (r, f_{B}) + {\bar{f}}_{A} cov (s_{A}, \frac{f_{A B}}{f_{A}}), \\ \frac{d D_{M}}{d t} = ({\bar{s}}_{A} - \bar{s} + {\bar{s}}_{B} - \bar{s}) D_{M} - (\bar{μ} + \bar{ρ}) D_{M} + {\bar{s}}_{E} {\bar{f}}_{a b} {\bar{f}}_{A B} + cov (r, D) + coskew (r, f_{A}, f_{B}) \\ + \sum_{d \in {A, B}} (1 - {\bar{f}}_{d}) {\bar{f}}_{d} cov (s_{d}, \frac{f_{A B}}{f_{d}}) + (1 - {\bar{f}}_{A}) Λ_{A a} - {\bar{f}}_{A} Λ_{a A} + (1 - {\bar{f}}_{B}) Λ_{B b} - {\bar{f}}_{B} Λ_{b B} . \end{array}

Note that in the equation $d D_{M} / d t$ , there are terms involving $Λ_{i j}$ which we chose to neglect in Equation (7) given in the main text. These terms are

Λ_{A a} = cov (μ_{A} + ρ_{A}, \frac{f_{a B}}{1 - f_{A}}) and Λ_{a A} = cov (μ_{a} + ρ_{a}, \frac{f_{A B}}{f_{A}}),

while

Λ_{B b} = cov (μ_{B} + ρ_{B}, \frac{f_{a B}}{1 - f_{B}}) and Λ_{b B} = cov (μ_{b} + ρ_{b}, \frac{f_{A B}}{f_{B}}) .

Thus the expression

(1 - {\bar{f}}_{A}) Λ_{A a} - {\bar{f}}_{A} Λ_{a A} + (1 - {\bar{f}}_{B}) Λ_{B b} - {\bar{f}}_{B} Λ_{b B}

in the equation $d D_{M} / d t$ is the effect upon $D_{M}$ of spatial heterogeneity in mutation and recombination rates ( $μ_{i}^{x} \neq μ_{i}^{y}$ and/or $ρ_{i}^{x} \neq ρ_{i}^{y}$ ) coupled with differences in the proportion of infections with allele $i$ (e.g. $i = A$ or $i = a$ ) that are resistant to the other drug ( $j = B$ ). In particular, populations in which infections are more likely to acquire resistance through mutation/recombination disproportionately affect metapopulation LD through an increase in doubly-resistant infections. However, these terms are likely to be quite small because they require that substantial differences in mutation/recombination rates exist between populations. Since these terms are unlikely to be a significant contributor to the dynamics of $D_{M}$ , we ignore them in the main text.

There remains a number of other quantities in system (Equation 26) that we now define in more detail. First, the probability that an infection resistant to drug $d$ is found in population $x$ is

p^{x} \frac{f_{d}^{x}}{{\bar{f}}_{d}} .

For example, if we apply our variable definitions, it is straightforward to show that

p^{x} \frac{f_{A}^{x}}{{\bar{f}}_{A}} = \frac{I_{A b}^{x} + I_{A B}^{x}}{\sum_{y = 1}^{N} (I_{A b}^{y} + I_{A B}^{y})} .

Next, to compute the metapopulation-level selection coefficients, and mutation/recombination rates, we need to compute the weighted average of the population quantities, where the weights are the probability that an infection of a particular type is in population $x$ (calculated above). Applying this logic, the metapopulation-level selection coefficients and epistasis are

{\bar{s}}_{i} = \sum_{x = 1}^{N} p^{x} \frac{f_{i}^{x}}{{\bar{f}}_{i}} s_{i}^{x} and {\bar{s}}_{E} = \sum_{x = 1}^{N} p^{x} \frac{f_{A B}^{x}}{{\bar{f}}_{A B}} s_{E}^{x} .

The average selection for resistance in the metapopulation is

\bar{s} = {\bar{s}}_{A} {\bar{f}}_{A} + {\bar{s}}_{B} {\bar{f}}_{B} + {\bar{s}}_{E} {\bar{f}}_{A B} .

The per-capita mutation and recombination rates follow similarly. Recall that $μ_{ℓ}$ and $ρ_{ℓ}$ are the per-capita rates at which infections gain allele $ℓ$ . Thus, for example,

{\bar{μ}}_{A} = \sum_{x = 1}^{N} p^{x} \frac{1 - f_{A}^{x}}{1 - {\bar{f}}_{A}} μ_{A}^{x} and {\bar{μ}}_{a} = \sum_{x = 1}^{N} p^{x} \frac{f_{A}^{x}}{{\bar{f}}_{A}} μ_{a}^{x} .

Similar calculations can be made to arrive at ${\bar{μ}}_{B}$ , ${\bar{μ}}_{b}$ , and the various ${\bar{ρ}}_{ℓ}$ . The total per-capita mutation and recombination rates are

\bar{μ} = {\bar{μ}}_{a} + {\bar{μ}}_{A} + {\bar{μ}}_{b} + {\bar{μ}}_{B} and \bar{ρ} = {\bar{ρ}}_{a} + {\bar{ρ}}_{A} + {\bar{ρ}}_{b} + {\bar{ρ}}_{B} .

Covariance and coskewness

Request a detailed protocol

Finally, we also use a number of covariance terms and a coskewness terms. Let $𝔼 [c]$ denote the expectation of the quantity $c$ . Then applying the definition of covariance, we have

\begin{array}{ll} cov (f_{A}, f_{B}) & = E [f_{A} f_{B}] - E [f_{A}] E [f_{B}] \\ = \sum_{x = 1}^{N} p^{x} f_{A}^{x} f_{B}^{x} - (\sum_{x = 1}^{N} p^{x} f_{A}^{x}) (\sum_{x = 1}^{N} p^{x} f_{B}^{x}) \end{array}

Following the same procedure, we can calculate $cov (r, f_{A})$ and $cov (r, D_{M})$ . When the covariance involves quantities that also specifically depend upon particular allele(s), the only difference is that when computing the expectation the probability used is the probability that an allele $ℓ$ is in population $x$ . For example,

\begin{array}{ll} cov (s_{A}, \frac{f_{A B}}{f_{A}}) & = E [s_{A} \frac{f_{A B}}{f_{A}}] - E [s_{A}] E [\frac{f_{A B}}{f_{A}}] \\ = \sum_{x = 1}^{N} p^{x} \frac{f_{A}^{x}}{{\bar{f}}_{A}} s_{A}^{x} \frac{f_{A B}^{x}}{f_{A}^{x}} - (\sum_{x = 1}^{N} p^{x} \frac{f_{A}^{x}}{{\bar{f}}_{A}} s_{A}^{x}) (\sum_{x = 1}^{N} p^{x} \frac{f_{A}^{x}}{{\bar{f}}_{A}} \frac{f_{A B}^{x}}{f_{A}^{x}}) \\ = \sum_{x = 1}^{N} p^{x} \frac{s_{A}^{x} f_{A B}^{x}}{{\bar{f}}_{A}} - (\sum_{x = 1}^{N} p^{x} \frac{f_{A}^{x}}{{\bar{f}}_{A}} s_{A}^{x}) (\sum_{x = 1}^{N} p^{x} \frac{f_{A B}^{x}}{{\bar{f}}_{A}}) \\ = \sum_{x = 1}^{N} p^{x} \frac{f_{A B}^{x}}{{\bar{f}}_{A}} (s_{A}^{x} - {\bar{s}}_{A}) . \end{array}

The covariance terms involving the recombination and mutation rates follow similarly, with the appropriate exchanges of variables. Finally, we have the coskewness term, which can be calculated as

\begin{array}{ll} coskew (r, f_{A}, f_{B}) & = E [(r - E [r]) (p_{A} - E [f_{A}]) (f_{B} - E [f_{B}])] \\ = c o v (r, f_{A} f_{B}) - {\bar{f}}_{B} cov (r, f_{A}) - {\bar{f}}_{A} cov (r, f_{B}) . \end{array}

Specific examples

Equilibrium analysis of metapopulation consisting of independent populations

Request a detailed protocol

This is a version of one of the models presented in Lehtinen et al., 2019. The metapopulation consists of $N$ populations. The populations are independent (i.e, there is no migration between populations), and each population is assumed to be of a fixed size of unity, so $S^{x} = 1 - \sum_{i j} I_{i j}^{x}$ . Resistance is gained and lost through unbiased mutation occurring at rate μ and there is no recombination. Therefore

\frac{d I_{i j}^{x}}{d t} = (β_{i j}^{x} S^{x} - α_{i j}^{x} - 𝟏_{a} (i) τ_{A}^{x} - 𝟏_{b} (j) τ_{B}^{x} - (1 - 𝟏_{A} (i) 𝟏_{B} (j)) τ_{A B}^{x}) I_{i j}^{x} + μ (\sum_{ℓ} (I_{ℓ j}^{x} + I_{i ℓ}^{x}) - 4 I_{i j}^{x}) .

Let $Δ z_{d}^{x}$ and $Δ z_{E}^{x}$ denote the contribution of parameter $z$ to the additive selection coefficient (for drug $d$ -resistance) and epistasis, respectively, in population $x$ . Specifically,

\begin{array}{ll} Δ β_{A}^{x} = β_{A b}^{x} - β_{a b}^{x}, Δ β_{B}^{x} = β_{a B}^{x} - β_{a b}^{x}, Δ β_{E}^{x} = β_{A B}^{x} + β_{a b}^{x} - β_{A b}^{x} - β_{a B}^{x} \\ Δ α_{A}^{x} = α_{A b}^{x} - α_{a b}^{x}, Δ α_{B}^{x} = α_{a B}^{x} - α_{a b}^{x}, Δ α_{E}^{x} = α_{A B}^{x} + α_{a b}^{x} - α_{A b}^{x} - α_{a B}^{x} . \end{array}

Then if we let $r_{i j}^{x}$ denote the per-capita growth term of an $i j$ -infection in subpopulation $x$ (the first term in brackets in Equation (36)), we can partition this as

r_{i j}^{x} = r^{x} + 𝟏_{A} (i) s_{A}^{x} + 𝟏_{B} (j) s_{B}^{x} + 𝟏_{A} (i) 𝟏_{B} (j) s_{E}^{x}

where

\begin{aligned} r^{x} & = β_{a b}^{x} S^{x} - α_{a b}^{x} - τ_{A}^{x} - τ_{B}^{x} - τ_{A B}^{x} \\ s_{A}^{x} & = Δ β_{A}^{x} S^{x} - Δ α_{A}^{x} + τ_{A}^{x} \\ s_{B}^{x} & = Δ β_{B}^{x} S^{x} - Δ α_{B}^{x} + τ_{B}^{x} \\ s_{E}^{x} & = Δ β_{E}^{x} S^{x} - Δ α_{E}^{x} + τ_{A B}^{x} \end{aligned}

This notation and formulation differs from that of Lehtinen et al., 2017; Lehtinen et al., 2019 in that they assumed costs were multiplicative, that is,

β_{a b}^{x} = β^{x}, β_{A b}^{x} = β^{x} c_{β_{A}}^{x}, β_{a B}^{x} = β^{x} c_{β_{B}}^{x}, β_{A B}^{x} = β^{x} c_{β_{A}}^{x} c_{β_{B}}^{x}

and

α_{a b}^{x} = α^{x}, α_{A b}^{x} = \frac{α^{x}}{c_{α_{A}}^{x}}, α_{a B}^{x} = \frac{α^{x}}{c_{α_{B}}^{x}}, α_{A B}^{x} = \frac{α^{x}}{c_{α_{A}}^{x} c_{α_{B}}^{x}}

where $0 \leq c_{β_{ℓ}}^{x} \leq 1$ and $0 \leq c_{α_{ℓ}}^{x} \leq 1$ (note the slightly different notation used for multiplicative costs in 1). There are two consequences of multiplicative costs . First, multiplicative costs produce epistasis. For the model of Lehtinen et al., 2017; Lehtinen et al., 2019:

s_{E}^{x} = β^{x} (1 - c_{β_{A}}^{x}) (1 - c_{β_{B}}^{x}) S^{x} - α^{x} \frac{(1 - c_{α_{A}}^{x}) (1 - c_{α_{B}}^{x})}{c_{α_{A}}^{x} c_{α_{B}}^{x}} + τ_{A B}^{x} .

Thus, in this model, there exists epistasis whenever there is a cost of resistance or drugs are prescribed in combination, $τ_{A B}$ . More specifically, transmission costs and combination treatment will produce positive epistasis, while duration of carriage costs will produce negative epistasis. Inclusion of epistasis (through multiplicative costs) is not necessarily a problem, and for epidemiological reasons multiplicative costs may be preferable. Indeed, because epistasis plays a central role in multilocus dynamics, it is valuable to recognize if/when epistasis is occurring. However, our analysis in the main text focused upon how population variation in susceptible densities can create correlations in the selection coefficients, favoring MDR, and so we excluded the possibility of epistasis.

The second consequence of multiplicative costs is that they have implications for the cost of resistance within a population. For the above model, and assuming $c_{α_{A}}^{x} = c_{α_{A}}$ , $c_{β_{A}}^{x} = c_{β_{A}}$ , the selection coefficient for resistance to drug $A$ can be written

s_{A}^{x} = τ_{A}^{x} - α^{x} \frac{1 - c_{α_{A}}}{c_{α_{A}}} - β^{x} \frac{1 - c_{β_{A}}}{c_{β_{A}}} S^{x} .

From Equation (43), we see that the costs of resistance depend upon the population’s epidemiological parameters (by assumption). Specifically, ignoring concomitant effects upon $S^{x}$ , populations that are more transmissible (larger $β^{x}$ ) with a shorter duration of carriage (larger $α^{x}$ ) pay higher costs of resistance due to how the cost parameters interact with the epidemiological parameters. For example, if there were no costs to transmission ( $c_{β_{A}} = 1$ ), then Equation (43) predicts that populations with longer duration of carriage (smaller $α^{x}$ ) are more likely to become resistant, because they pay disproportionately lower costs. Indeed, if the populations represent serotype, than this is an example of epistasis between serotype and resistance, which favors LD between duration of carriage and resistance.

To put this in a biological context, if (for example) the populations correspond to the different capsular serotypes of S. pneumoniae, it is possible that the differences between capsules interact with the mechanism of resistance so as to make resistance more costly for more transmissible capsular serotypes or those capsular serotypes associated with longer duration of carriage (multiplicative costs), but it is also possible that no interaction occurs between the capsule differences and the mechanism of resistance (additive costs) or that resistance is less costly for more transmissible serotypes or for capsular serotypes associated with a shorter duration of carriage.

Irrespective of whether the costs are multiplicative or additive, we would attach the constraints that $0 \leq β_{i j}^{x} \leq β_{a b}^{x}$ and $0 \leq α_{a b}^{x} \leq α_{i j}^{x}$ , that is, carriage of one or more resistance alleles will never increase transmissibility or decrease duration of carriage, respectively; otherwise there are no costs. As an example, if we were interested in attaching additive resistance costs to transmission, if we let $c_{β_{d}}^{x}$ and $c_{β_{A B}}^{x}$ be the (additive) transmission cost of resistance to drug $d$ and epistatic transmission cost, respectively, in population $x$ , so that

β_{i j}^{x} = β_{a b}^{x} - 𝟏_{A} (i) c_{β_{A}}^{x} - 𝟏_{B} (i) c_{β_{B}}^{x} - 𝟏_{A B} (i j) c_{β_{A B}}^{x},

then we could attach the constraints

β_{i j}^{x} = \min (0, β_{a b}^{x} - 𝟏_{A} (i) c_{β_{A}}^{x} - 𝟏_{B} (i) c_{β_{B}}^{x} - 𝟏_{A B} (i j) c_{β_{A B}}^{x}),

or simply assume that

c_{β_{A}}^{x} + c_{β_{B}}^{x} + | c_{β_{A B}}^{x} | \leq β_{a b}^{x} .

There are many possible ways of implementing the constraints. One point to note is that even with additive costs, the choice of constraints can also potentially create epistasis; this could occur in the case of Equation (45).

In Figure 3, we consider three scenarios; whenever possible we choose parameter values to agree with those of Figure 4 in Lehtinen et al., 2019. In each scenario, we assume there are 20 independent populations, that the per-capita mutation rate is $μ = 10^{- 10}$ , and there is no epistasis, $s_{E}^{x} = 0$ . In subplot 3a, we set $β_{a b}^{x} = 2$ , while duration of carriage varies among populations from $α_{a b}^{x} = 0.25$ to $α_{a b}^{x} = 1.75$ . In subplot 3b we set $α_{a b}^{x} = 0.5$ , while transmission varies among populations from $β_{a b}^{x} = 1$ to $β_{a b}^{x} = 3$ . In both subplots 3a and 3b, $Δ α_{A}^{x} = Δ α_{B}^{x} = 0$ , while $Δ β_{A}^{x} = Δ β_{B}^{x} = - 0.1$ . Finally in subplot 3 c, $Δ β_{A}^{x} = Δ β_{B}^{x} = 0$ , while duration of carriage varies among populations from $α_{a b}^{x} = 0.25$ to $α_{a b}^{x} = 1.75$ , with $Δ α_{A}^{x} = Δ α_{B}^{x} = 0.05$ .

Transient dynamics and MDR in Streptococcus pneumoniae

Request a detailed protocol

Here we use a variant of the model originally proposed by Lehtinen et al., 2017, Lehtinen et al., 2019 in which the populations represent different serotypes. Resistance is gained and lost through unbiased mutation at a per-capita rate µ and there is no recombination of resistance loci.

Applying these assumptions and using the notation presented with our model from the main text, this yields

\begin{array}{ll} \frac{d I_{i j}^{x}}{d t} = (β_{i j}^{x} ν (I, x) S - α_{i j}^{x} - 1_{a} (i) τ_{A} - 1_{b} (j) τ_{B} - (1 - 1_{A} (i) 1_{B} (j)) τ_{A B}) I_{i j}^{x} + μ (\sum_{ℓ} (I_{ℓ j}^{x} + I_{i ℓ}^{x}) - 4 I_{i j}^{x}) \end{array}

where

ν (I, x) = {(1 - [\frac{\sum_{i j} I_{i j}^{x}}{\sum_{k = 1}^{N} \sum_{i j} I_{i j}^{k}} - \frac{1}{N}])}^{ω}

is a balancing function intended to mimic the stabilizing effect adaptive host immunity has upon serotype diversity ( $ω$ controls the strength of this effect; see Lehtinen et al., 2017; Lehtinen et al., 2019). The treatment rates in Equation (47) are assumed to be independent of serotype. Note that although we could mechanistically model the susceptible hosts available to each subpopulation, $S^{x}$ , it is more straightforward and computationally simpler to use the phenomenological model given above in which $S^{x} = ν (I, x) S$ . The function $ν (I, x)$ ensures that $S^{x}$ has the two properties we are interested in: (i) there can be variation across populations of available susceptibles, and (ii) this variation is linked to population attributes (e.g. transmissibility and duration of carriage). The primary conclusions of our analysis would hold if a mechanistic model for $S^{x}$ were used instead.

If we let $r_{i j}^{x}$ denote the per-capita growth term of an $i j$ -infection belonging to serotype $x$ (the first term in brackets in Equation (47)), we can partition this as

r_{i j}^{x} = r^{x} + 𝟏_{A} (i) s_{A}^{x} + 𝟏_{B} (j) s_{B}^{x} + 𝟏_{A} (i) 𝟏_{B} (j) s_{E}^{x}

where if we use the notation introduced in Equation (37),

\begin{aligned} r^{x} & = β_{a b}^{x} ν (I, x) S - α_{a b}^{x} - τ_{A} - τ_{B} - τ_{A B} \\ s_{A}^{x} & = Δ β_{A}^{x} ν (I, x) S - Δ α_{A}^{x} + τ_{A} \\ s_{B}^{x} & = Δ β_{B}^{x} ν (I, x) S - Δ α_{B}^{x} + τ_{B} \\ s_{E}^{x} & = Δ β_{E}^{x} ν (I, x) S - Δ α_{E}^{x} + τ_{A B} \end{aligned}

For simplicity, we keep total population size constant, and so set $S = 1 - \sum_{x = 1}^{N} \sum_{i j} I_{i j}^{x}$ .

The simulations in Figure 5 assume the metapopulation is initially treated at per-capita rates $(τ_{A}, τ_{B}, τ_{A B}) = (0.12, 0, 0)$ , until $t = 1000$ when these rates switch to $(τ_{A}, τ_{B}, τ_{A B}) = (0.07, 0.1, 0)$ . Other parameters values used are $N = 12$ , $ω = 3$ , $Δ β_{A}^{x} = Δ β_{B}^{x} = - 0.2$ , $Δ α_{A}^{x} = Δ α_{B}^{x} = 0.05$ and $μ = 10^{- 8}$ . Finally, because Streptococcus serotypes differ based upon duration of carriage and transmissibility, and there is evidence of a positive correlation between the two (Weinberger et al., 2009; Zafar et al., 2017), $α_{a b}^{x}$ was chosen to assume evenly spaced parameter values from $α_{a b}^{x} = 0.2$ to $α_{a b}^{x} = 0.7$ , while $β_{a b}^{x}$ was chosen to assume evenly spaced parameter values from $β_{a b}^{x} = 3.25$ to $β_{a b}^{x} = 3$ . Cost epistasis, when it is present, is assumed to solely effect transmissibility (i.e. $Δ α_{E}^{x} = 0$ ). When there is positive epistasis, $Δ β_{E}^{x} = 0.065$ , whereas for negative epistasis, $Δ β_{E}^{x} = - 0.015$ .

Contrasting drug prescription strategies in a hospital-community setting

Request a detailed protocol

When we model the hospital and community, we use Equation (2) and assume the susceptible host density is controlled by

\begin{array}{ll} \frac{d S^{x}}{d t} = θ^{x} - d S^{x} - m^{x \to y} S^{x} + m^{y \to x} S^{y} - \sum_{i j} β_{i j}^{x} I_{i j}^{x} S^{x} \\ + \sum_{i j} (α_{i j}^{x} - d) I_{i j}^{x} + \sum_{i j} (1_{a} (i) τ_{A}^{x} + 1_{b} (j) τ_{B}^{x} + (1 - 1_{A} (i) 1_{B} (j)) τ_{A B}^{x}) I_{i j}^{x} \end{array}

where $θ^{x}$ is the influx of new hosts and $d$ is the background mortality rate.

In the hospital/community model, we assume population $C$ is the ‘community’ and population $H$ is the ‘hospital’. Therefore, we let $θ^{H} = 0$ , and $m^{C \to H} = m / \sum_{i j} I_{i j}^{C}$ be the rate at which individuals are admitted to the hospital, which is independent of population size. Individuals exit the hospital at a constant rate $m^{H \to C}$ , so they spend on average $1 / m^{H \to C}$ time units in hospital (assuming background mortality is low). The specification of the migration rates in this way allows us to ensure the ‘community’ is always much larger than the ‘hospital’. In Figure 6, we assumed that the total prescription rate per population, $τ^{C}$ or $τ^{H}$ , was the same for each strategy, and that drugs are prescribed at 15× the rate in the hospital versus the community, that is, $τ^{H} = 15 τ^{C}$ . For ‘mixing’, this means $(τ_{A}^{x}, τ_{B}^{x}, τ_{A B}^{x}) = (τ^{x} / 2, τ^{x} / 2, 0)$ , whereas for ‘combination’, this means $(τ_{A}^{x}, τ_{B}^{x}, τ_{A B}^{x}) = (0, 0, τ^{x})$ . Finally, for ‘cycling’ drug $A$ and $B$ were rotated from hospital to community every 50 time units so that either $(τ_{A}^{x}, τ_{B}^{x}, τ_{A B}^{x}) = (τ^{x}, 0, 0)$ or $(τ_{A}^{x}, τ_{B}^{x}, τ_{A B}^{x}) = (0, τ^{x}, 0)$ depending on the rotation. Therefore in Figure 6 the period is of length $T = 100$ . In Figure 6, we numerically integrate the system until $t = 10^{4}$ ; the final state at $t = 10^{4}$ is then what is plotted for the combination and mixing scenarios, whereas for cycling we plot the average state across the last two rotations (i.e. the average state over the final period, from $t = 9, 900$ to $t = 10^{4}$ ). Figure 7 shows how changing the length of time between drug rotations affects the evolution of single- and multi-drug resistance. Specifically, when drug rotations are frequent (with period of $T = 1$ ), cycling behaves like mixing and so positive LD is produced (Figure 7a). As drug rotations become less frequent (period of $T = 24$ and $T = 160$ ), cycling generates a negative covariance in selection, which in turn produces negative LD (Figure 7a). Thus, when drug rotations are more frequent, single-drug resistance is delayed and emerges at higher treatment rates, but the evolution of MDR occurs at lower treatment rates (Figure 7b). When drug rotations are infrequent, single-drug resistance emerges at lower treatment rates, but MDR evolution is delayed, emerging at higher treatment rates (Figure 7b). Parameters used in Figure 6 and Figure 7 were $β_{a b}^{x} = 2$ , $Δ β_{A}^{x} = Δ β_{B}^{x} = - 0.4$ , $α_{a b}^{x} = 0.1$ , $Δ α_{A}^{x} = Δ α_{B}^{x} = 0.02$ , $d = 0.01$ , $θ^{C} = 0.2$ , $θ^{H} = 0$ , $m^{H \to C} = 0.5$ , $m = 0.2$ , $μ = 10^{- 7}$ , $σ = 0$ .

Figure 7

Download asset Open asset

Time between drug rotations affects the evolution of both single- and multi-drug resistance.

When drugs are rotated every 0.5 time units (period of $T = 1$ ; magenta curves), cycling behaves like mixing and positive LD is generated. As we increase the time between rotations (period of $T = 24$ and $T = 160$ ), a negative covariance in selection is generated, producing negative LD (dashed and solid blue curves). In panel a we show the metapopulation LD, $D_{M}$ , while in panel b, we show the frequency of resistance in the metapopulation. When drugs are rotated frequently, single drug resistance emerges at higher treatment rates but MDR emerges at lower treatment rates, as compared to when drugs are rotated infrequently. Thus, there is a trade-off (indicated by the arrows in panel b) associated with time between drug rotations: we can delay single drug resistance but promote MDR (frequent drug rotations), or delay MDR but promote single drug resistance (infrequent drug rotations). In all cases, we integrated the system until $t = 10^{4}$ , then averaged the system state over the final two rotations (i.e. over a single period). The remaining parameter values are provided in Materials and methods 'Contrasting drug prescription strategies in a hospital-community setting'.

Data availability

All data used was from a previously published study (Lehtinen et al. 2019 PLoS Pathogens and Turner et al. 2012 PLoS ONE); this data has been uploaded by those authors to a public repository, we downloaded it from there and have provided the details.

The following previously published data sets were used

(2019) PLoS Pathogens
S1 File. Resistance profiles and duration of carriage estimates for the Maela dataset.
https://doi.org/10.1371/journal.ppat.1007763.s002

References

(2014) Cycling empirical antibiotic therapy in hospitals: meta-analysis and models
PLOS Pathogens 10:e1004225.
https://doi.org/10.1371/journal.ppat.1004225
- PubMed
- Google Scholar
1. Alizon S
(2013) Co-infection and super-infection models in evolutionary epidemiology
Interface Focus 3:20130031.
https://doi.org/10.1098/rsfs.2013.0031
- PubMed
- Google Scholar
1. Althaus CL
2. Bonhoeffer S
(2005) Stochastic interplay between mutation and recombination during the acquisition of drug resistance mutations in human immunodeficiency virus type 1
Journal of Virology 79:13572–13578.
https://doi.org/10.1128/JVI.79.21.13572-13578.2005
- PubMed
- Google Scholar
1. Austin DJ
2. Anderson RM
(1999) Studies of antibiotic resistance within the patient, hospitals and the community using simple mathematical models
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 354:721–738.
https://doi.org/10.1098/rstb.1999.0425
- PubMed
- Google Scholar
1. Barton NH
(1995)
Linkage and the limits to natural selection

Genetics 140:821–841.
- PubMed
- Google Scholar
(2017) Antibiotic Cycling and Antibiotic Mixing: Which One Best Mitigates Antibiotic Resistance?
Molecular Biology and Evolution 34:802–817.
https://doi.org/10.1093/molbev/msw292
- PubMed
- Google Scholar
1. Beardmore RE
2. Pena-Miller R
(2010) Antibiotic cycling versus mixing: the difficulty of using mathematical models to definitively quantify their relative merits
Mathematical Biosciences and Engineering : MBE 7:923–933.
https://doi.org/10.3934/mbe.2010.7.923
- PubMed
- Google Scholar
(2004) Ecological theory suggests that antimicrobial cycling will not reduce antimicrobial resistance in hospitals
PNAS 101:13285–13290.
https://doi.org/10.1073/pnas.0402298101
- PubMed
- Google Scholar
1. Blanquart F
(2019) Evolutionary epidemiology models to predict the dynamics of antibiotic resistance
Evolutionary Applications 12:365–383.
https://doi.org/10.1111/eva.12753
- PubMed
- Google Scholar
(1997) Evaluating treatment protocols to prevent antibiotic resistance
PNAS 94:12106–12111.
https://doi.org/10.1073/pnas.94.22.12106
- PubMed
- Google Scholar
(2004) Recombination in HIV and the evolution of drug resistance: for better or for worse?
BioEssays : News and Reviews in Molecular, Cellular and Developmental Biology 26:180–188.
https://doi.org/10.1002/bies.10386
- PubMed
- Google Scholar
1. Cobey S
2. Lipsitch M
(2012) Niche and neutral effects of acquired immunity permit coexistence of pneumococcal serotypes
Science 335:1376–1380.
https://doi.org/10.1126/science.1215947
- PubMed
- Google Scholar
1. Colijn C
2. Cohen T
3. Fraser C
4. Hanage W
5. Goldstein E
6. Givon-Lavi N
7. Dagan R
8. Lipsitch M
(2010) What is the mechanism for persistent coexistence of drug-susceptible and drug-resistant strains of Streptococcus pneumoniae?
Journal of the Royal Society, Interface 7:905–919.
https://doi.org/10.1098/rsif.2009.0400
- PubMed
- Google Scholar
(2015) Selective and genetic constraints on pneumococcal serotype switching
PLOS Genetics 11:e1005095.
https://doi.org/10.1371/journal.pgen.1005095
- PubMed
- Google Scholar
1. Davies NG
2. Flasche S
3. Jit M
4. Atkins KE
(2019) Within-host dynamics shape antibiotic resistance in commensal Bacteria
Nature Ecology & Evolution 3:440–449.
https://doi.org/10.1038/s41559-018-0786-x
- PubMed
- Google Scholar
1. Day T
2. Gandon S
(2012) The evolutionary epidemiology of multilocus drug resistance
Evolution; International Journal of Organic Evolution 66:1582–1597.
https://doi.org/10.1111/j.1558-5646.2011.01533.x
- PubMed
- Google Scholar
1. Feldman MW
2. Christiansen FB
(1974) The effect of population subdivision on two loci without selection
Genetical Research 24:151–162.
https://doi.org/10.1017/S0016672300015184
- PubMed
- Google Scholar
1. Felsenstein J
(1965)
The effect of linkage on directional selection

Genetics 52:349–363.
- PubMed
- Google Scholar
Book
1. Fisher RA
(1930) The Genetical Theory of Natural Selection
Clarendon Press.
https://doi.org/10.5962/bhl.title.27468
- Google Scholar
1. Gillespie JH
(2000)
Genetic drift in an infinite population. The pseudohitchhiking model

Genetics 155:909–913.
- PubMed
- Google Scholar
1. Hall AR
2. MacLean RC
(2011) Epistasis buffers the fitness effects of rifampicin- resistance mutations in Pseudomonas aeruginosa
Evolution; International Journal of Organic Evolution 65:2370–2379.
https://doi.org/10.1111/j.1558-5646.2011.01302.x
- PubMed
- Google Scholar
1. Henriques-Normark B
2. Tuomanen EI
(2013) The pneumococcus: epidemiology, microbiology, and pathogenesis
Cold Spring Harbor Perspectives in Medicine 3:a010215.
https://doi.org/10.1101/cshperspect.a010215
- PubMed
- Google Scholar
1. Hill WG
2. Robertson A
(1966) The effect of linkage on limits to artificial selection
Genetical Research 8:269–294.
https://doi.org/10.1017/S0016672300010156
- PubMed
- Google Scholar
(2020) Factors favouring the evolution of multidrug resistance in Bacteria
Journal of the Royal Society Interface 17:20200105–20200114.
https://doi.org/10.1098/rsif.2020.0105
- Google Scholar
1. Karlin S
(1975) General two-locus selection models: some objectives, results and interpretations
Theoretical Population Biology 7:364–398.
https://doi.org/10.1016/0040-5809(75)90025-8
- PubMed
- Google Scholar
1. Keightley PD
2. Otto SP
(2006) Interference among deleterious mutations favours sex and recombination in finite populations
Nature 443:89–92.
https://doi.org/10.1038/nature05049
- PubMed
- Google Scholar
(2006) Effect of varying epistasis on the evolution of recombination
Genetics 173:589–597.
https://doi.org/10.1534/genetics.105.053108
- PubMed
- Google Scholar
(2009) Recombination and drug resistance in HIV: population dynamics and stochasticity
Epidemics 1:58–69.
https://doi.org/10.1016/j.epidem.2008.11.001
- PubMed
- Google Scholar
(2020) Population structure across scales facilitates coexistence and spatial heterogeneity of antibiotic-resistant infections
PLOS Computational Biology 16:e1008010.
https://doi.org/10.1371/journal.pcbi.1008010
- PubMed
- Google Scholar
1. Lehtinen S
2. Blanquart F
3. Croucher NJ
4. Turner P
5. Lipsitch M
6. Fraser C
7. C. Fc
(2017) Evolution of antibiotic resistance is linked to any genetic mechanism affecting bacterial duration of carriage
PNAS 114:1075–1080.
https://doi.org/10.1073/pnas.1617849114
- PubMed
- Google Scholar
(2019) On the evolutionary ecology of multidrug resistance in Bacteria
PLOS Pathogens 15:e1007763.
https://doi.org/10.1371/journal.ppat.1007763
- PubMed
- Google Scholar
1. Lenormand T
2. Otto SP
(2000)
The evolution of recombination in a heterogeneous environment

Genetics 156:423–438.
- PubMed
- Google Scholar
1. Lewontin RC
(1964) The interaction of selection and linkage. I. general considerations; Heterotic models
Genetics 49:49–67.
https://doi.org/10.1093/genetics/49.1.49
- PubMed
- Google Scholar
1. Lewontin RC
2. Kojima K
(1960) The evolutionary dynamics of complex polymorphisms
Evolution 14:458–472.
https://doi.org/10.1111/j.1558-5646.1960.tb03113.x
- Google Scholar
1. Li WH
2. Nei M
(1974) Stable linkage disequilibrium without epistasis in subdivided populations
Theoretical Population Biology 6:173–183.
https://doi.org/10.1016/0040-5809(74)90022-7
- PubMed
- Google Scholar
(2000) The epidemiology of antibiotic resistance in hospitals: paradoxes and prescriptions
PNAS 97:1938–1943.
https://doi.org/10.1073/pnas.97.4.1938
- PubMed
- Google Scholar
1. Lipsitch M
2. Colijn C
3. Cohen T
4. Hanage WP
5. Fraser C
(2009) No coexistence for free: neutral null models for multistrain pathogens
Epidemics 1:2–13.
https://doi.org/10.1016/j.epidem.2008.07.001
- PubMed
- Google Scholar
(2010) The population genetics of antibiotic resistance: integrating molecular mechanisms and treatment contexts
Nature Reviews. Genetics 11:405–414.
https://doi.org/10.1038/nrg2778
- PubMed
- Google Scholar
(2006) Selection for recombination in structured populations
Genetics 172:593–609.
https://doi.org/10.1534/genetics.104.039982
- PubMed
- Google Scholar
1. Neher RA
(2013) Genetic draft, selective interference, and population genetics of rapid adaptation
Annual Review of Ecology, Evolution, and Systematics 44:195–215.
https://doi.org/10.1146/annurev-ecolsys-110512-135920
- Google Scholar
1. Neher RA
2. Shraiman BI
(2011) Statistical genetics and evolution of quantitative traits
Reviews of Modern Physics 83:1283–1300.
https://doi.org/10.1103/RevModPhys.83.1283
- Google Scholar
1. Nei M
2. Li WH
(1973) Linkage disequilibrium in subdivided populations
Genetics 75:213–219.
https://doi.org/10.1093/genetics/75.1.213
- PubMed
- Google Scholar
1. Nowak MA
2. May RM
(1994) Superinfection and the evolution of parasite virulence
Proceedings. Biological Sciences 255:81–89.
https://doi.org/10.1098/rspb.1994.0012
- PubMed
- Google Scholar
1. Ohta T
(1982a) Linkage disequilibrium due to random genetic drift in finite subdivided populations
PNAS 79:1940–1944.
https://doi.org/10.1073/pnas.79.6.1940
- PubMed
- Google Scholar
1. Ohta T
(1982b) Linkage disequilibrium with the island model
Genetics 101:139–155.
https://doi.org/10.1093/genetics/101.1.139
- PubMed
- Google Scholar
1. Otto SP
2. Barton NH
(2001) Selection for recombination in small populations
Evolution; International Journal of Organic Evolution 55:1921–1931.
https://doi.org/10.1111/j.0014-3820.2001.tb01310.x
- PubMed
- Google Scholar
Website
1. O’Neill J
(2015) Tackling a global health crisis: initial steps
The Reviewon Antimicrobial Resistance. London: Wellcome Trust and Government Ofthe United Kingdom. Accessed September 6, 2018.

https://amr-review.org/sites/default/files/RARJ3003_Global_health_crisis_report_20.03.15_OUTLINED.pdf
Book
1. Rice SH
(2004)
Evolutionary Theory: Mathematical and Conceptual Foundations

Sunderland: Sinauer Associates.
- Google Scholar
1. Rosen MJ
2. Davison M
3. Bhaya D
4. Fisher DS
(2015) Microbial diversity. Fine-scale diversity and extensive recombination in a quasisexual bacterial population occupying a broad niche
Science 348:1019–1023.
https://doi.org/10.1126/science.aaa4456
- PubMed
- Google Scholar
1. Slatkin M
(1975) Gene flow and selection in a two-locus system
Genetics 81:787–802.
https://doi.org/10.1093/genetics/81.4.787
- PubMed
- Google Scholar
1. Slatkin M
(2008) Linkage disequilibrium--understanding the evolutionary past and mapping the medical future
Nature reviews. Genetics 9:477–485.
https://doi.org/10.1038/nrg2361
- PubMed
- Google Scholar
(2017) Modeling antibiotic treatment in hospitals: a systematic approach shows benefits of combination therapy over cycling, mixing, and mono-drug therapies
PLOS Computational Biology 13:e1005745.
https://doi.org/10.1371/journal.pcbi.1005745
- PubMed
- Google Scholar
1. Trindade S
2. Sousa A
3. Xavier KB
4. Dionisio F
5. Ferreira MG
6. Gordo I
(2009) Positive epistasis drives the acquisition of multidrug resistance
PLOS Genetics 5:e1000578.
https://doi.org/10.1371/journal.pgen.1000578
- PubMed
- Google Scholar
1. Turner P
2. Turner C
3. Jankhot A
4. Helen N
5. Lee SJ
6. Day NP
7. White NJ
8. Nosten F
9. Goldblatt D
(2012) A longitudinal study of Streptococcus pneumoniae carriage in a cohort of infants and their mothers on the Thailand-Myanmar border
PLOS ONE 7:e38271.
https://doi.org/10.1371/journal.pone.0038271
- PubMed
- Google Scholar
1. Weinberger DM
2. Trzciński K
3. Lu YJ
4. Bogaert D
5. Brandes A
6. Galagan J
7. Anderson PW
8. Malley R
9. Lipsitch M
(2009) Pneumococcal capsular polysaccharide structure predicts serotype prevalence
PLOS Pathogens 5:e1000476.
https://doi.org/10.1371/journal.ppat.1000476
- PubMed
- Google Scholar
1. Wilson DS
2. Yoshimura J
(1994) On the coexistence of specialists and generalists
The American Naturalist 144:692–707.
https://doi.org/10.1086/285702
- Google Scholar
1. Zafar MA
2. Hamaguchi S
3. Zangari T
4. Cammer M
5. Weiser JN
(2017) Capsule type and amount affect shedding and transmission of Streptococcus pneumoniae
mBio 8:e00989-17.
https://doi.org/10.1128/mBio.00989-17
- PubMed
- Google Scholar

Decision letter

George H Perry

Senior and Reviewing Editor; Pennsylvania State University, United States
Sonja Lehtinen

Reviewer; ETH Zurich, Switzerland
Stephen M Kissler

Reviewer; Harvard TH Chan School of Public Health, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Acceptance summary:

This paper addresses the important question of multidrug resistance evolution, which is of both theoretical and applied interest. The authors' efforts to carefully distinguish population and metapopulation linkage disequilibrium and to develop a framework to rigorously analyze the relationship between the two represent an advance in our understanding of microbial population dynamics.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting your work entitled "Understanding the evolution of multiple drug resistance in structured populations" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by a Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Brian J Arnold (Reviewer #2).

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

We have appended to this message a summary of the two reviewers' comments, including both major and minor concerns. Briefly, there were concerns about the extent of the claims made about the model explaining features of data, the clarity of the model development and notation, and links between the results and previous literature. Please see below for details.

Reviewer 1:

MDR is more common than expected by random chance, especially in particular bacterial species. Understanding how this linkage between resistance elements arises is important for public health but also for generally making sense of genomic data, e.g. classifying genomic patterns as outcomes of specific biological processes. Here, McLeod and Gandon show that spatial heterogeneity across subdivided populations can create LD between resistance elements, even in the absence of epistasis.

The novelty of this manuscript appears to be an extension of Day and Gandon, (2012) to a metapopulation model as well as comparing output of the model to a dataset on pneumococcus, which has sort of served as a model for the ecology and evolution of MDR. One particular result of interest was the observation of simulated transient LD even without epistasis. Though this particular observation is not necessarily new (Martin et al., 2006; Day and Gandon, 2012), casting this metapopulation model in an epidemiological framework appears to be new to my knowledge.

However, the intuition this model provides for D_tot > 0 at equilibrium is intriguing: variability in susceptibility density across populations can create covariance in selection coefficients for the two alleles, in turn creating correlations between subpopulation allele frequencies, driving D_tot > 0. This insight is extremely general and has nothing to do with the biology of a particular species.

While this isn't a major concern, I wanted to point out what I thought was confusing language that made it take longer for me to understand the results.

1) I found the number of ways to quantify LD a little confusing and done with ambiguous language. D is first used in it's standard form D = P_ab – P_a*P_b (Equation 1), which made it confusing when I then encountered D within the definition of D_tot (Equation 5), where D is the weighted average of D = P_ab – P_a*P_b for each of the subpopulations. I think it would be clearer if D perhaps consistently referred to a subpopulation, the average D across subpopulations were instead indicated with a summation in Equation 5.

Also with Equation 5, it took me a while to figure out what cov(p_a, p_b) actually referred to, as the definition "spatial covariance between resistance to drugs A and B", along with being primed by historical literature, made me think of this covariance as LD (i.e. how frequently alleles are found within the same genome). However, this is referring to covariance of *allele frequencies across subpopulations* and has nothing to do with whether or not they're found in the same genome. More specific/clear language will help readers distinguish between the different covariances studied here: between alleles across genomes (LD) and between allele frequencies across subpopulations (cov(p_a, p_b)).

2) Also, I'm left not knowing the ultimate conclusion of applying the model to the pneumococcal data. On line 226 you claim that transient dynamics and epistasis can explain patterns in pneumococcus, whereas on line 301 you claim that variation in S^X across subpopulations can generate LD. Some short conclusion regarding these two claims would be useful.

Citations:

Day & Gandon, (2012). The evolutionary epidemiology of multilocus drug resistance. Evolution. https://doi.org/10.1111/j.1558-5646.2011.01533.x

Martin, Otto & Lenormand, (2006). Selection for recombination in structured populations. Genetics, 172(1), 593-609. https://doi.org/10.1534/genetics.104.039982

1) For the case of additive selection, you provide one particular example of parameters that give rise to transient positive LD (Figure 4), but it would seem that many situations would also give rise to negative LD. Is there a symmetry to the model such that if you were to "integrate" across all possible starting conditions, the expectation would be zero? Or is there an asymmetry that makes transient LD more likely to be positive? I assume the latter if you add epistasis, but it's unclear with additive selection. This would help map intuition from your model to results such as Figure 3, in which it seems there's a skew towards positive LD. Is epistasis absolutely needed to explain this skew?

2) One interesting implication of this work is what to expect when bacterial "populations" are either not defined at all or even incorrectly defined. Since Dtot is very positive under many scenarios (even with negative epistasis! Figure 4), you will observe an excess of positive LD between resistance elements (given these starting conditions in Figure 4, that is).

Reviewer 2:

This is an interesting and important topic and the paper promised to be a new take on it, from the point of view of dynamics of linkage disequilibrium.

My major concerns are (1) clarity of the mathematics, notation and model development and (2) with the claims that the model explains data, without a wide range of simulations or a direct comparison of model to data (for example see below re line 239 onwards).

The model in Equation (2) looks linear and the reader has to see figure 1 to note that there is actually a nonlinear term here (the susceptible * infectious term, and a product of two infectious terms) because the s_A etc depend on the prevalence.

The model derivation is cumbersome and far from straightforward; I was unable to directly derive (3) from (2). Instead, for the first three terms of (3) I obtained

dpA/dt = (1-pA)(pA sA + pAB (sAB – sA – sB) ) + sB D

instead of the expression in the paper which is

(1-pA) (pA sA + sAB pAB) + sB D

This would seem to carry through to the other equations. However my derivation attempts were simply using sA, sAB, sB as given in the text, and putting in the rate of growth of the AB variant (r + sAB) as in (2) for example, without substituting in the expressions from Figure 1 for sA, sB, and sAB. This may work it out (?) However this left me a bit sceptical.

Bottom of page 4: if doubly-resistant infections are overrepresented in the population, D^X > 0. < ----> middle page 5 "since D^X > 0" Which is it?

Page 10: "…, drug B resistance is often more likely to occur in an infection with a genetic background resistant to drug A." -- more likely than what?

Line 239++ : This section makes the claim:

"Thus transient dynamics coupled with epistasis can explain the significant within-serotype LD observed in Streptococcus pneumoniae".

This is a big claim, based on a qualitative paragraph with no direct comparison to data, in which there is an underlying assumption of fixation of resistance to drug A in many serotypes. To my knowledge we often do not see fixation of resistance in Streptococcus pneumoniae or in the Maela data in particular, but rather long-term coexistence. I am not convinced that the dynamics of this model explain the significant LD in the data as claimed.

Equation (2) – there is a note above that rho_l^X depends on infection densities; the nature of this dependence should be specified.

Page 12: why, if there is no mechanism maintaining within-population diversity, does this necessarily mean D^X = 0 (line 282) ?

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Understanding the evolution of multiple drug resistance in structured populations" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and George Perry as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Sonja Lehtinen (Reviewer #1); Stephen M Kissler (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

While the reviewers of the previous version of your paper were not available to reconsider the revised manuscript, the reviewers of this version of the paper did take the previous comments and your responses into consideration when composing their own reviews.

Collectively we see much potential in your manuscript, but we request the following points to be addressed at minimum (see below reviews for further details and more suggestions):

1. Expanded explanation of the assumptions behind the LD framework. In particular:

a) Is epistasis necessarily defined in terms of an additive expectation on growth rate? (If this is a standard result in population genetics, a reference to an accessible explanation would be enough).

b) The s coefficients are dependent on variable densities – does this matter for the partitioning and interpretation of equations 3 and 4?

2. If the authors disagree with Reviewer 1's questions about the interpretation that variation in susceptible density explains the effect in example 1 and in Lehtinen et al., 2019, more explanation is needed for the following points:

a) Why variation in clearance rate affects selection for resistance when there is a multiplicative cost on clearance rate, but no cost on transmission rate (for the re-interpretation of Lehtinen specifically)?

b) What is happening in example 2, where the effect is currently explained in terms of serotype-specific susceptible density, but the model does not include serotype-specific susceptibles?

3. For examples 2 and 3, please either explain why additive transmission costs are reasonable (given the concerns highlighted below) or change the modeling of cost.

Reviewer #1 (Recommendations for the authors):

My major scientific recommendations are covered by my review. There are some additional less fundamental points that I am happy to share with the authors if they are interested but won't include here as this review is already quite long.

In terms of presentation, I thought everything was mostly very clear, with the exception of the section on serotype dynamics. Here it would be more helpful to make clear the shift from a structured susceptible population to a shared susceptible population.

Reviewer #2 (Recommendations for the authors):

I was asked to review this manuscript after it had already gone through a round of revisions. My comments reflect both my own initial reading of the manuscript and my assessment of the authors' responses to the previous reviews.

While I did not have access the authors' initial submission, it appears that the authors have sufficiently addressed the previous reviewers' comments.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Understanding the evolution of multiple drug resistance in structured populations" for further consideration by eLife. Your revised article has been evaluated by George Perry (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed as outlined below in edited comments from Reviewer Sonja Lehtinen, as outlined in points 2-5 below. In addition, I have reviewed the discussion with the eLife editorial office about the article structure and have considered this question myself, ultimately coming to point 1.

1. Please do create a methods section for this paper. There is no need to remove anything from your current main text. Rather, in the new methods section please first overview the methods even if they are already described in more detail elsewhere in the paper (some repetition is ok, and I agree with the level of detail you have provided in the Results section, given the type of paper) in one or more subheadings, and then import all of the sub-headings and the full content of the appendix. Since we don't have page limitations and since this is important to the paper, I think it should be included in the methods section of the main text of the paper, rather than as an appendix.

2. The point the authors are making about density of susceptibles and that this point is not limited to additive costs is understood. However, the section on equilibrium patterns of MDR still should be updated because the difference between additive and multiplicative costs is not explicitly flagged in the main text. (I think the discussion in the Sup Mat is very helpful, but most readers will not read the Supplement). Specifically, when considering equilibrium patterns of MDR, the paragraph starting line 283 is only true for additive costs. Without reading the supplement, this paragraph comes across as a general result. If costs are multiplicative, variation in duration of carriage is both necessary (i.e. variation in transmission rate does not produce LD at equilibrium) and sufficient (explicit transmission costs are not needed). Thus, overall, it is true to say that variation in duration of carriage is not necessary, but not that it is not sufficient. It would be helpful to explicitly make this difference clear somewhere in this section. Please also rephrase the sentence lines 270-273 to avoid implying that the costs in Lehtinen et al., are additive.

3. The main text would benefit from a discussion of additive costs. Specifically:

a. The authors make a good point that cost parameters are subject to constraints however they are modelled. The reviewer's point was that these constraints are qualitatively different for additive transmission costs: reasonable additive costs for single resistance can lead to negative transmission rate for dual resistance (as in the example in my review) in the absence of cost epistasis. This suggests that for cost parameters where this is the case, assuming no epistasis cannot be appropriate. Although the specific parameter values the authors use don't give rise to this problem, the authors also present more general results for a model with additive costs and no epistasis, so it would be helpful to highlight this constraint somewhere.

b. The authors' careful consideration of how the specification of costs affects interpretation is indeed very useful and interesting. In addition to this conceptual point, the authors say that they are aiming to highlight mechanisms – e.g. variation in transmission rate – that could plausibly give rise to patterns of MDR. In the case of equilibrium patterns, variation in transmission rate only gives rise to LD when transmission costs are additive. The plausibility of variation in transmission rate as an explanation for equilibrium patterns of MDR therefore depends on the plausibility of additive transmission costs. Therefore, if the authors want to suggest that variation in transmission rate is a plausible explanation for patterns of MDR, it is necessary to include a discussion of whether/how additive transmission costs might arise.

4. Add "per carriage episode" to the summary of the biological interpretation of the duration of carriage effect in Lehtinen et al., 2017/2019 (lines 244-245), to make it clear that they were not suggesting that (host) populations with longer durations of pathogen carriage have greater antibiotic exposure. (This point was missed in the original review).

5. One final observation, which is not an essential revision and the authors should only address it if they think it will improve the paper. Points 2 and 3 are both related to the result that a model with a multiplicative transmission cost predicts, at equilibrium, LD between duration of carriage and resistance, but not transmission rate and resistance. Would the authors explain this in terms of variation in the density of susceptibles giving rise to the LD with duration of carriage, but this effect being offset by the epistatic interaction between transmission rate and cost of resistance leading to no LD between transmission rate and resistance? This might be interesting to explicitly discuss in the manuscript.

https://doi.org/10.7554/eLife.65645.sa1

Author response

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

Essential revisions:

Reviewer 1:

MDR is more common than expected by random chance, especially in particular bacterial species. Understanding how this linkage between resistance elements arises is important for public health but also for generally making sense of genomic data, e.g. classifying genomic patterns as outcomes of specific biological processes. Here, McLeod and Gandon show that spatial heterogeneity across subdivided populations can create LD between resistance elements, even in the absence of epistasis.

The novelty of this manuscript appears to be an extension of Day and Gandon, (2012) to a metapopulation model as well as comparing output of the model to a dataset on pneumococcus, which has sort of served as a model for the ecology and evolution of MDR. One particular result of interest was the observation of simulated transient LD even without epistasis. Though this particular observation is not necessarily new (Martin et al., 2006; Day and Gandon, 2012), casting this metapopulation model in an epidemiological framework appears to be new to my knowledge.

However, the intuition this model provides for D^tot > 0 at equilibrium is intriguing: variability in susceptibility density across populations can create covariance in selection coefficients for the two alleles, in turn creating correlations between subpopulation allele frequencies, driving D_tot > 0. This insight is extremely general and has nothing to do with the biology of a particular species.

Yes, one of the main insights we gain from our analysis is a clarification and a generalisation of the results discussed recently by Lehtinen et al., (2019) about the conditions favouring MDR at equilibrium. We believe our Equation (5) provides a new and useful way to understand the evolution of MDR. We would also like to highlight that the dynamical equations we derive are very general and can be used to understand, and perhaps limit, the spread of MDR. For instance, the final example regarding the impact of cycling antibiotics between a hospital and community is new (previous studies did not study the influence of cycling drugs in a metapopulation on MDR evolution) and opens up new perspectives on the management of MDR resistance.

In our revisions, we extended and rewrote our discussions of the examples and we show more explicitly how equations 3 to 7 can help understand the complex evolutionary dynamics of drug resistance.

While this isn't a major concern, I wanted to point out what I thought was confusing language that made it take longer for me to understand the results.

1) I found the number of ways to quantify LD a little confusing and done with ambiguous language. D is first used in its standard form D = P_ab – P_a*P_b (Equation 1), which made it confusing when I then encountered D within the definition of D_tot (Equation 5), where D is the weighted average of D = P_ab – P_a*P_b for each of the subpopulations. I think it would be clearer if D perhaps consistently referred to a subpopulation, the average D across subpopulations were instead indicated with a summation in Equation 5.

Also with Equation 5, it took me a while to figure out what cov(p_a, p_b) actually referred to, as the definition "spatial covariance between resistance to drugs A and B", along with being primed by historical literature, made me think of this covariance as LD (i.e. how frequently alleles are found within the same genome). However, this is referring to covariance of *allele frequencies across subpopulations* and has nothing to do with whether or not they're found in the same genome. More specific/clear language will help readers distinguish between the different covariances studied here: between alleles across genomes (LD) and between allele frequencies across subpopulations (cov(p_a, p_b)).

We agree and we have rewritten many sections of the manuscript with this comment in mind. Additionally, we changed the notation to try improve the readability. We believe this rewording has significantly increased the clarity of our arguments.

2) Also, I'm left not knowing the ultimate conclusion of applying the model to the pneumococcal data. On line 226 you claim that transient dynamics and epistasis can explain patterns in pneumococcus, whereas on line 301 you claim that variation in S^X across subpopulations can generate LD. Some short conclusion regarding these two claims would be useful.

We have rewritten both of these examples in light of the reviewers’ comments, and hope the conclusions now are clearer. Specifically, variation in S^x (density of susceptibles) across populations can generate metapopulation LD, and in fact can lead to MDR overrepresentation in the metapopulation. However, this example is focused upon an equilibrium analysis. Consequently, although metapopulation LD exists at equilibrium, because nothing maintains diversity within-populations (i.e., within-serotypes), there will be no LD within serotypes at equilibrium. This would be true with or without epistasis; i.e., with epistasis, we would expect at equilibrium that $D^x = 0$ even if $D_M \not = 0$.

In the next example (the transient + epistasis example), when we consider transient dynamics we are specifically focused upon mechanisms that can generate (transient) LD within-serotype, i.e., $D^x \not=0$. As we detail in this example, there are two possibilities: epistasis or migration (serotype recombination or serotype switching), however migration (serotype recombination) is an unlikely explanation because the rate at which this occurs would have to be unrealistically high. Therefore, we focus upon exploring the potential role played by epistasis, and show that it can (transiently) generate LD within-serotype. Of course, since there is nothing maintaining within-host diversity within-serotype, at equilibrium $D^x = 0$ as in the previous example. We believe that the way these examples have been rewritten makes it clearer what we are trying to show.

1) For the case of additive selection, you provide one particular example of parameters that give rise to transient positive LD (Figure 4), but it would seem that many situations would also give rise to negative LD. Is there a symmetry to the model such that if you were to "integrate" across all possible starting conditions, the expectation would be zero?

In this example, the shared dependence of the additive selection coefficients upon S^x means that there is a tendency to produce positive (metapopulation) LD, both transiently and at equilibrium. In order for this model to produce negative LD, we would need there to be some (strong) factor generating negative LD to counteract this effect. One possibility would be negative epistasis, but this needs to be very strong, as it needs to overcome the positive covariance between additive selection coefficients induced by the shared dependence on S^x.

In terms of the initial conditions, we are starting from a state in which there is no LD, and then changing the selective conditions faced by the metapopulation. This then leads to alleles segregating within the population in response to the new selective pressures, which in turn generates LD.

If we were to instead manipulate the initial conditions (i.e., set up negative LD or positive LD) without changing the selective pressures, this is a slightly different problem and is related to the Hill Robertson effect in which LD can interfere with directional selection. In this case, HR theory predicts negative LD would persist longer in the population than positive LD by weakening directional selection. We now discuss HR effects in more detail in the conclusion.

Or is there an asymmetry that makes transient LD more likely to be positive? I assume the latter if you add epistasis, but it's unclear with additive selection. This would help map intuition from your model to results such as Figure 3, in which it seems there's a skew towards positive LD. Is epistasis absolutely needed to explain this skew?

We detailed the emergence of LD at the metapopulation and at the population scale in the revised version. At the metapopulation scale, the positive LD is due to the heterogeneity in baseline growth rate among serotypes (our equation 7 shows how this leads to D_M>0). At the scale of serotypes, we find that, indeed, epistasis is required to generate LD (see figure 5e and 5f).

2) One interesting implication of this work is what to expect when bacterial "populations" are either not defined at all or even incorrectly defined. Since Dtot is very positive under many scenarios (even with negative epistasis! Figure 4), you will observe an excess of positive LD between resistance elements (given these starting conditions in Figure 4, that is).

Yes. We like this comment regarding the fact that we often lack a proper description of the metapopulation structure of microbial communities. We now develop this idea further in the conclusions section. Thanks.

Reviewer 2:

This is an interesting and important topic and the paper promised to be a new take on it, from the point of view of dynamics of linkage disequilibrium.

My concerns are (1) clarity of the mathematics, notation and model development

We reworded several sections and changed the notation to improve the readability of our work. As pointed out above, we extended our discussion of the results in the examples section. We believe this analysis helps to show how our analytical work (equations 3 to 7) can help understand the complex dynamics.

and (2) with the claims that the model explains data, without a wide range of simulations or a direct comparison of model to data (for example see below re line 239 onwards).

We claim that our analysis can help to understand observed patterns. For instance, we provide a general understanding on the conditions promoting MDR. This analysis clarifies the link between heterogeneities among mechanisms of clearance and the evolution of drug resistance.

But we do not try to fit quantitatively observed patterns. This would be feasible in principle (parameter estimation) but, in the present work, our main objective was to present a qualitative understanding of the interplay between selection, migration and recombination.

The model in Equation (2) looks linear and the reader has to see figure 1 to note that there is actually a nonlinear term here (the susceptible * infectious term, and a product of two infectious terms) because the s_A etc depend on the prevalence.

We clarified that (line 66-67).

The model derivation is cumbersome and far from straightforward; I was unable to directly derive (3) from (2). Instead, for the first three terms of (3) I obtained

dpA/dt = (1-pA)(pA sA + pAB (sAB – sA – sB) ) + sB D

instead of the expression in the paper which is

(1-pA) (pA sA + sAB pAB) + sB D

This would seem to carry through to the other equations. However my derivation attempts were simply using sA, sAB, sB as given in the text, and putting in the rate of growth of the AB variant (r + sAB) as in (2) for example, without substituting in the expressions from Figure 1 for sA, sB, and sAB. This may work it out (?) However this left me a bit sceptical.

We have changed the notation to avoid this type of confusion. s_AB was meant to describe epistasis in fitness (not selection on AB genotypes; which is s_A + s_B + s_AB). We now use the notation s_E for epistasis in fitness, and we are very explicit in the main text and in the appendix how we derive the different equations.

Bottom of page 4: if doubly-resistant infections are overrepresented in the population, D^X > 0. < ----> middle page 5 "since D^X > 0" Which is it?

We have clarified the wording in this paragraph.

Page 10: "…, drug B resistance is often more likely to occur in an infection with a genetic background resistant to drug A." -- more likely than what?

We have rewritten this entire example and so this sentence no longer appears.

Line 239++ : This section makes the claim:

"Thus transient dynamics coupled with epistasis can explain the significant within-serotype LD observed in Streptococcus pneumoniae".

This is a big claim, based on a qualitative paragraph with no direct comparison to data, in which there is an underlying assumption of fixation of resistance to drug A in many serotypes. To my knowledge we often do not see fixation of resistance in Streptococcus pneumoniae or in the Maela data in particular, but rather long-term coexistence. I am not convinced that the dynamics of this model explain the significant LD in the data as claimed.

It was not our intent to claim that this simple example can fully explain the data. Our intention was to demonstrate that our dynamical equations help understand what we observe in the simulations. This is key because, at first sight, multilocus models may seem overwhelmingly complex. We show with this example that our analysis helps to understand both the short-term and long-term dynamics of drug resistance.

Note as well, that for the parameter values we have chosen (taken when possible from the Lehtinen et al., 2019 paper for Streptococcus Pneumonaie), it is straightforward to generate scenarios (as in Figure 5) in which there is a considerable amount of transient polymorphism with or without epistasis in figures 5a, b and c. The point here is that things which look like they are in equilibrium based upon surveillance data need not be.

Moreover, the existing explanation for the Maela data presented in Lehtinen et al. 2019 makes a very clear and testable prediction (as we discuss in the text): there should be no LD within-serotype. Yet this prediction is not met in the empirical data (see Figure 4). This would suggest a number of possibilities, including that either: (1) the populations are in equilibrium, but the coexistence mechanism proposed by Lehtinen et al., is not operating, or (2) the coexistence mechanism proposed by Lehtinen et al. is in effect, but the population is not in equilibrium.

In our paper we focus upon examining whether relaxing the conditions in (2), i.e., that the population is in equilibrium, could plausibly explain the patterns of LD in combination with the coexistence mechanism proposed by Lehtinen et al. We show that it could. We stress that this need not be the only explanation (e.g., condition (1) may be at play), and more work would be needed to definitively answer this question.

However, we recognize that the language we previously used in the manuscript may have been misleading about what we were and were not trying to do, and what our analysis does and does not say. As such, we have carefully gone through the manuscript to ensure that the language is not misleading about any conclusions of our analysis, and we have been more specific about what our objectives were for this example.

Equation (2) – there is a note above that rho_l^X depends on infection densities; the nature of this dependence should be specified.

We have added a reference on line 52-53 indicating the equation in the Sup. Mat. which specifies the full dependencies.

Page 12: why, if there is no mechanism maintaining within-population diversity, does this necessarily mean D^X = 0 (line 282) ?

This sentence specifically refers to equilibrium in a deterministic model consisting of independent populations. If nothing maintains within-population diversity, then some combination of alleles will be fixed in population $x$, e.g., $f^x_Ab = 1$ (using the notation where “f_ij” refers to the frequency of genotype ij). Since $D^x = f_AB^x f_ab^x – f_Ab^x f_aB^x$, when $f^x_Ab = 1 $and so $f_AB^x = f_ab^x = f_aB^x = 0$, then it follows that $D^x = 0$.

[Editors’ note: what follows is the authors’ response to the second round of review.]

Essential Revisions:

While the reviewers of the previous version of your paper were not available to reconsider the revised manuscript, the reviewers of this version of the paper did take the previous comments and your responses into consideration when composing their own reviews.

Collectively we see much potential in your manuscript, but we request the following points to be addressed at minimum (see below reviews for further details and more suggestions):

1. Expanded explanation of the assumptions behind the LD framework. In particular:

a) Is epistasis necessarily defined in terms of an additive expectation on growth rate? (If this is a standard result in population genetics, a reference to an accessible explanation would be enough).

Yes, for continuous time models epistasis is defined in terms of additive expectations on growth rate. This is because in continuous time models, fitness is per-capita growth rate, r. For discrete time models, epistasis is multiplicative (see e.g., Sup. Info. of Kouyos et al., 2009 Epidemics). We have added a paragraph stating this on lines 70-78.

But to see why it is a necessity, denote per-capita growth of genotype (strain) jk ∈ {ab,Ab,aB,AB} in population x as $r_{jk}^{x}$ . Consider evolution restricted to locus A when allele b is fixed, and let $p_{A}^{x}$ be the frequency of allele A in population x. Then the evolution of allele A is governed by

$\frac{{d p}_{A}^{x}}{d t} = (r_{Ab}^{x} - r_{ab}^{x}) p_{A}^{x} (1 - p_{A}^{x})$ Therefore allele A will increase in frequency whenever $r_{Ab}^{x} > r_{ab}^{x}$ ; thus $s_{A}^{x} = r_{Ab}^{x} - r_{ab}^{x}$ is the selection coefficient for allele A. Now epistasis is defined as any interaction between loci affecting fitness. Since $s_{A}^{x}$ and $s_{B}^{x}$ are the additive contribution to baseline fitness, $r_{ab}^{x}$ , by strains carrying alleles A or B, the fitness of strain AB is the ‘baseline’ fitness, plus the change due to carrying allele A and B ( $s_{A}^{x}$ and $s_{B}^{x}$ , respectively), plus any remaining fitness unaccounted for due to interactions between loci, $s_{E}^{x}$ (epistasis): $r_{AB}^{x} = r_{ab}^{x} + s_{A}^{x} = s_{B}^{x} + s_{E}^{x}$ Rearranging gives $s_{E}^{x} = r_{AB}^{x} - r_{ab}^{x} - s_{A}^{x} - s_{B}^{x} = r_{AB}^{x} + r_{ab}^{x} - r_{Ab}^{x} - r_{aB}^{x}$ Thus epistasis in fitness is defined in terms of the per-capita growth rates.

There are many different ways we can model the costs of carrying resistance. All are a priori valid but what we want to stress here is the consequences of the way we model these costs on epistasis and consequently on LD.

For a generic epidemiological model, let β^xS^x and µ^x be the per-capita transmission and mortality rate in population x, where S^x is the available density of susceptibles. Let $c_{β_{k}}^{x} S^{x}$ and $c_{μ_{k}}^{x}$ be the decrease in transmission and increase in clearance due to carriage of the resistance allele k. Finally, let $s_{E}^{x}$ denote any epistatic interactions.

Then the per-capita growth of strain jk is

$r_{jk}^{x} = (β^{x} - 1_{A} (j) c_{β_{A}}^{x} - 1_{B} (k) c_{β_{B}}^{x}) S^{x} - 1_{a} (j) τ_{B}^{x} - 1_{b} (k) τ_{B}^{x} - μ^{x} - 1_{A} (j) c_{μ_{A}}^{x} - 1_{B} (k) c_{μ_{B}}^{x} + 1_{A} (j) 1_{B} (k) s_{E}^{x},$ where $1_{𝓁} (𝓏)$ is 1 if $1 = 𝓏$ and zero otherwise. Therefore $s_{k}^{x} = τ_{k}^{x} - c_{β_{k}}^{x} S^{x} - c_{μ_{k}}^{x}, k \in {A, B}$ For example, in Lehtinen et al., they de_ne per-capita growth of strain k in population

i, resistant to the set of antibiotics R_k and sensitive to the set S_k as

$r_{k}^{x} = β^{x} \prod_{{jϵR}_{k}}^{} c_{β_{j}} S^{x} - μ^{x} \prod_{{jϵR}_{k}}^{} c_{μ_{j}} - \sum_{{jϵS}_{k}}^{} τ_{j}^{x}$ We have slightly modified their notation here, but the gist is the same; we have also ignored the possibility of combination treatment (but see Sup. Mat. 4.1 for an example in which it is included). Therefore in the case of two antibiotics, A and B,

we have $c_{β_{k}}^{x} = (1 - c_{β_{k}}) S^{x}$ and $c_{μ_{k}}^{x} = (c_{μ_{k}} - 1)$ and so

$s_{k}^{x} = τ_{k}^{x} - (1 - c_{β_{k}}) β^{x} S^{x} - (c_{μ_{k}} - 1) μ^{x}, k \in {A, B}$ and $s_{E}^{x} = (1 - c_{β_{A}}) (1 - c_{β_{B}}) β^{x} S^{x} - (c_{μ_{A}} - 1) (c_{μ_{B}} - 1) μ^{x}$ Thus there exists epistasis in fitness whenever there is a cost of resistance. More specifically, transmission costs will produce positive epistasis, while duration of carriage costs will produce negative epistasis. We have added further clarification on lines 716-738.

b) The s coefficients are dependent on variable densities – does this matter for the partitioning and interpretation of equations 3 and 4?

The partitioning and interpretation of equations 3 and 4 are unaffected by the s coefficients depending upon variable densities; the s coefficients have to be partitioned as they are from consideration of the single locus dynamics (see equation (1)); and this then forces epistasis to be what it is.

Yet, evolutionary dynamics can be more complicated in evolutionary epidemiology than in classic population genetic models as epidemiological feedbacks can change the sign and magnitude of the selective coefficients and epistasis. Although this creates richer dynamics, the fundamental population genetic principles still apply, and provide a useful interpretative framework.

2. If the authors disagree with Reviewer 1's questions about the interpretation that variation in susceptible density explains the effect in example 1 and in Lehtinen et al. 2019, more explanation is needed for the following points:

a) Why variation in clearance rate affects selection for resistance when there is a multiplicative cost on clearance rate, but no cost on transmission rate (for the re-interpretation of Lehtinen specifically)?

Consider the evolution of resistance to drug A in a metapopulation in which allele b is fixed. From (1), resistance is selected for in population x whenever the selection coefficient is positive, $s_{A}^{x} > 0$ . So resistance is selected for whenever the benefits $(τ_{A}^{x})$ exceed the costs $(c_{β_{A}}^{x} S^{x} + c_{μ_{A}}^{x})$ . This is true when the system has reached an endemic equilibrium but also during transient dynamics.

Now how is variation between the populations generated? This could occur a number of ways; but in all cases, the key is that the selection coefficient, $s_{A}^{x},$ differs between populations. For example, variation in treatment rate or in the parameters controlling the costs are sufficient. However, in keeping with Lehtinen et al., 2017, 2019, we assume that $τ_{A}^{x} = τ_{A}$ and that the cost parameters are population independent.

Then there are two ways in which variation can be produced, depending upon if the cost parameters are additive or are multiplicative (i.e., do they interact or not with the epidemiological parameters?)

i) Additive costs $(c_{β_{A}}^{x} = c_{β_{A}}, c_{μ_{A}}^{x} = c_{μ_{A}})$ : here costs are ‘independent’ of the population, except insofar as the transmission costs also depend upon density of susceptibles.

From $s_{A}^{x}$ , variation in resistance (due to selection) across populations is only possible if both c_βA > 0 and the susceptibles, S^x, show variation.

This was the example we focused upon in the manuscript; the prediction here is that any factor negatively correlated with S^x will lower the costs of resistance and so select for resistance. For example, populations with a longer duration of carriage (smaller µ^x or longer duration of carriage), or higher transmissibility (larger β^x), will more substantially deplete available susceptibles, S^x, lowering the costs of resistance, $c_{β_{A}} S^{x} + c_{μ_{A}},$ and so select for resistance.

ii) Multiplicative costs: here the cost parameters multiplicatively interact with the epidemiological parameters. This was true in the model of Lehtinen et al., 2017, 2019 (detailed above); for this model, by assumption (due to the multiplicative interaction), populations with a longer duration of carriage (smaller µ^x) pay disproportionately lower costs to acquire resistance than populations with shorter duration of carriage (larger µ^x), and so are more likely to become associated with resistance.

Another way of looking at this is if we suppose the populations correspond to serotype, then by specifying an interaction between $c_{μ_{A}}$ and µ^x, we are creating epistasis between population (serotype) x and resistance. Epistasis produces LD (in this case between serotype and resistance).

Therefore there are two distinct mechanisms generating variation: through the density of susceptibles, or by assuming multiplicative interactions between the cost parameters and the epidemiological parameters. We focused upon additive costs, because we wanted to emphasize the correlation between selection coefficients due to susceptible densities, but this does not mean that multiplicative costs are not valid.

For example, in the case of S. pneumoniae, there may be a physiological reason why the mechanism of resistance interacts with capsule, such that certain capsule types reduce the efficacy of antibiotic resistance over others, or make resistance more costly. Since different capsule types are associated with different durations of carriage, it is possible that such an interaction might mean that capsules with shorter duration of carriage amplify the costs of resistance.

In the Sup. Mat. on lines 716-738 we have clarified this issue; we have also modified the language in the main text (lines 247-250, lines 290-292) when we are discussing the Lehtinen et al., example to make clear that we are primarily focused upon correlations in variation mediated through susceptibles.

b) What is happening in example 2, where the effect is currently explained in terms of serotype-specific susceptible density, but the model does not include serotype-specific susceptibles?

The pool of available susceptibles depends upon serotype due to the stabilizing function intended to mimic adaptive immunity (this is the function from Lehtinen et al., 2017). So the ν(I,x) generates the fraction of susceptibles that are available to serotype x, i.e., S^x = ν(I,x)S. Thus there is variation in available susceptibles between serotypes and so there is serotype-specific susceptible density.

It would be possible to mechanistically model each of the populations of susceptibles, however, this would be computationally difficult and not add anything, as the stabilizing function ensures that there is: (1) the potential for serotype-specific differences in available susceptibles and (2) the variation in available susceptibles is linked to attributes of the serotypes (specifically, fitness).

In the Sup. Mat. on lines 756-763 we have added text which addresses this issue.

3. For examples 2 and 3, please either explain why additive transmission costs are reasonable (given the concerns highlighted below) or change the modeling of cost.

‘Additive’ transmission costs are always reasonable, given suitable constraints (e.g., ensuring transmission cannot be negative, as in the example provided by the reviewer); but the requirement for suitable constraints is also true for ‘multiplicative’ costs (e.g., the constraint 0 ≤ c_β ≤ 1). In the analysis we conducted in examples 2 and 3, appropriate constraints were attached. The key distinction is the costs have different evolutionary implications (as reviewer 1 correctly pointed out).

More generally, we have reframed the relevant sections so it does not seem as if we are advocating for either additive or multiplicative costs; our intent was to point out that they have different evolutionary interpretations. For example, one consequence of multiplicative costs is that they tend to produce epistasis.

We also have amended the text to stress that the inclusion of epistasis in a model is not ‘problematic’ (as our wording incorrectly suggested in the Sup. Mat.); in fact, in many circumstances it is necessary. Our point was that awareness of the presence/absence of epistasis can guide the interpretation.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed as outlined below in edited comments from Reviewer Sonja Lehtinen, as outlined in points 2-5 below. In addition, I have reviewed the discussion with the eLife editorial office about the article structure and have considered this question myself, ultimately coming to point 1.

1. Please do create a methods section for this paper. There is no need to remove anything from your current main text. Rather, in the new methods section please first overview the methods even if they are already described in more detail elsewhere in the paper (some repetition is ok, and I agree with the level of detail you have provided in the Results section, given the type of paper) in one or more subheadings, and then import all of the sub-headings and the full content of the appendix. Since we don't have page limitations and since this is important to the paper, I think it should be included in the methods section of the main text of the paper, rather than as an appendix.

We have now created a methods section that is the appendix. We have added a paragraph on lines 589600 providing an overview of the contents of the methods. In the main text, we now have changed all the references to refer to the relevant sections in the methods.

2. The point the authors are making about density of susceptibles and that this point is not limited to additive costs is understood. However, the section on equilibrium patterns of MDR still should be updated because the difference between additive and multiplicative costs is not explicitly flagged in the main text. (I think the discussion in the Sup Mat is very helpful, but most readers will not read the Supplement). Specifically, when considering equilibrium patterns of MDR, the paragraph starting line 283 is only true for additive costs. Without reading the supplement, this paragraph comes across as a general result. If costs are multiplicative, variation in duration of carriage is both necessary (i.e. variation in transmission rate does not produce LD at equilibrium) and sufficient (explicit transmission costs are not needed). Thus, overall, it is true to say that variation in duration of carriage is not necessary, but not that it is not sufficient. It would be helpful to explicitly make this difference clear somewhere in this section. Please also rephrase the sentence lines 270-273 to avoid implying that the costs in Lehtinen et al., are additive.

We have made a number of changes to address these points. First, we have added a box to the main text that discusses the costs of resistance (specifically, additive vs multiplicative) and their consequences for epistasis on pg. 5.

Next, we have made a number of changes to the wording to clarify the points the reviewer is concerned about. Specifically:

– On line 251, we identify that we are using additive costs

– On line 252 we state that this differs from Lehtinen et al., 2019 who use multiplicative costs

– On line 276-277 we identify again that our result here differs from Lehtinen et al., 2019 based upon the choice of costs

– On line 287 we state that our results are for additive costs

– In the caption for figure 3 we identify that we are dealing with additive costs (and no epistasis)

3. The main text would benefit from a discussion of additive costs. Specifically:

a. The authors make a good point that cost parameters are subject to constraints however they are modelled. The reviewer's point was that these constraints are qualitatively different for additive transmission costs: reasonable additive costs for single resistance can lead to negative transmission rate for dual resistance (as in the example in my review) in the absence of cost epistasis. This suggests that for cost parameters where this is the case, assuming no epistasis cannot be appropriate. Although the specific parameter values the authors use don't give rise to this problem, the authors also present more general results for a model with additive costs and no epistasis, so it would be helpful to highlight this constraint somewhere.

In box 1 we have highlighted the constraint, it is that $0 \le \β_ij^x \le \β_ab^x$. This constraint applies to any form of cost.

In the case of additive costs, there are a variety of ways this constraint could be implemented. We discuss this on lines 762-772 and provide two examples – one of which can lead to epistasis.

b. The authors' careful consideration of how the specification of costs affects interpretation is indeed very useful and interesting. In addition to this conceptual point, the authors say that they are aiming to highlight mechanisms – e.g. variation in transmission rate – that could plausibly give rise to patterns of MDR. In the case of equilibrium patterns, variation in transmission rate only gives rise to LD when transmission costs are additive. The plausibility of variation in transmission rate as an explanation for equilibrium patterns of MDR therefore depends on the plausibility of additive transmission costs. Therefore, if the authors want to suggest that variation in transmission rate is a plausible explanation for patterns of MDR, it is necessary to include a discussion of whether/how additive transmission costs might arise.

The magnitude of the costs and epistasis will be dictated by the biology of the system. For example, when discussing populations separated based upon capsular serotype, the structure of the capsule has a role in determining transmission/duration of carriage. If the resistance mechanism does not interact with these differences in capsule structure, then carriage of the resistance allele will come with additive costs. However, it is also possible that the differences between capsule serotype may interact with the resistance mechanism. For example, more transmissible capsular serotypes may pay lower or higher costs of resistance. If more transmissible serotypes pay higher costs, this is multiplicative costs (or at least how we have been using it); if more transmissible capsular serotypes pay lower costs, this would require a slightly different form of the costs – but this would also yield epistasis. Similar logic applies to duration of carriage. We have added a section in the Methods (lines 754-761) to address this.

More generally, without specific evidence of an interaction between differences between capsular serotypes and the resistance mechanism, then a priori it is reasonable to assume there is no interaction (additive costs).

4. Add "per carriage episode" to the summary of the biological interpretation of the duration of carriage effect in Lehtinen et al., 2017/2019 (lines 244-245), to make it clear that they were not suggesting that (host) populations with longer durations of pathogen carriage have greater antibiotic exposure. (This point was missed in the original review).

Whoops, yes we have added that. Thanks for pointing it out.

5. One final observation, which is not an essential revision and the authors should only address it if they think it will improve the paper. Points 2 and 3 are both related to the result that a model with a multiplicative transmission cost predicts, at equilibrium, LD between duration of carriage and resistance, but not transmission rate and resistance. Would the authors explain this in terms of variation in the density of susceptibles giving rise to the LD with duration of carriage, but this effect being offset by the epistatic interaction between transmission rate and cost of resistance leading to no LD between transmission rate and resistance? This might be interesting to explicitly discuss in the manuscript.

If we understand this scenario correctly, the reviewer is suggesting a model in which the selection coefficients are of the form :

$s_d^x = -c_{\β_d}^x\β^xS^x + \tau_d^x$ and there is variation in duration of carriage (but no costs of resistance to duration of carriage).

In this model, there is no way to isolate the susceptible density so that it is the only thing in the selection coefficient varying across populations without also setting $\β^x = \β$. This is because the parameter controlling the `costs’ in this case is not just $c_{\β_d}^x$, but is actually the composite parameter $c_{\β_d}^x\β^x$. This is related to our point about the interpretation of multiplicative costs: they assume a dependence upon other life-history characteristics.

At equilibrium, as the reviewer pointed out, there is cancellation between the $\β^x$ term in the costs and the $\β^x$ term in the denominator of the susceptible density. But this is really just a product of the mathematical assumptions; and it would be incorrect to refer to c_{\β_d} as the `costs’ because the cost `parameter’ actually depends upon $\β^x$.

Put another way, in the most general setting we could have simply stated we have transmission rates \β_ij^x without specifying any relationship whatsoever between `costs’ and ‘transmission’; in this case it is easy to see that only if we make specific mathematical assumptions will cancellation occur.

For example, if we were interested in a scenario in which populations with higher transmission rates pay disproportionately higher costs of resistance (which is the same biologic assumption as multiplicative costs), we would be equally justified in choosing the function

\β_{Ab}^x = \β^x/(1 + c_{\β_A} \β^x) (*)

as we would be in choosing multiplicative costs

\β_{Ab}^x = (1 – c_{\β_A})\β^x (#)

But using the function (*), at equilibrium \β^x won’t disappear from the invasion condition as it will with the function (#).

We have now highlighted in Box 1 that costs need not necessarily be ‘additive’ or ‘multiplicative’; they could just as easily take a more complex form as in our example (*) above. We simply want to emphasize that the choice of costs can produce epistasis, which is an important consideration for modelling multilocus dynamics.

https://doi.org/10.7554/eLife.65645.sa2

Article and author information

Author details

David V McLeod

Centre D'Ecologie Fonctionnelle & Evolutive, CNRS, Univ Montpellier, EPHE, IRD, Montpellier, France

Contribution
Conceptualization, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

For correspondence
david.mcleod@cefe.cnrs.fr

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2551-7877
Sylvain Gandon

Centre D'Ecologie Fonctionnelle & Evolutive, CNRS, Univ Montpellier, EPHE, IRD, Montpellier, France

Contribution
Conceptualization, Supervision, Investigation, Methodology, Writing - review and editing

For correspondence
sylvain.gandon@cefe.cnrs.fr

Competing interests
No competing interests declared

Funding

Natural Sciences and Engineering Research Council of Canada (Postdoctoral fellowship)

David V McLeod

Agence Nationale de la Recherche (ANR-17-CE35-0012)

Sylvain Gandon

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Sonja Lehtinen, Brian J Arnold and Stephen M Kissler for helpful comments. This work was supported by a NSERC-CRSNG postdoctoral fellowship to DVM. SG acknowledges support from the Agence Nationale de la Recherche, grant ANR-17-CE35-0012 ‘EVOMALWILD’.

Senior and Reviewing Editor

George H Perry, Pennsylvania State University, United States

Reviewers

Sonja Lehtinen, ETH Zurich, Switzerland
Stephen M Kissler, Harvard TH Chan School of Public Health, United States

Publication history

Received: December 10, 2020
Accepted: May 28, 2021
Accepted Manuscript published: June 1, 2021 (version 1)
Version of Record published: June 16, 2021 (version 2)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

1,088

Page views
155

Downloads
6

Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Article PDF

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

David V McLeod
Sylvain Gandon

(2021)

Understanding the evolution of multiple drug resistance in structured populations

eLife 10:e65645.

https://doi.org/10.7554/eLife.65645

Share this article

Cite this article

Notation used in main text.

Schematic of the dynamics of system (2).

Box 1.

The effect of migration upon LD at equilibrium depends upon the scale at which LD is measured.

Duration of carriage is one of many potential explanations for MDR over-representation at equilibrium.

Linkage disequilibrium for different drug pairs in Streptococcus pneumoniae.

Transient dynamics coupled with epistasis can explain patterns of serotype LD in Streptococcus pneumoniae.

Different antibiotic prescription strategies generate different patterns of LD at equilibrium.

Time between drug rotations affects the evolution of both single- and multi-drug resistance.

Author details

David V McLeod

Contribution

For correspondence

Competing interests

Sylvain Gandon

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading