## Figures

## Abstract

Antibiotic-resistant infections are a growing threat to human health, but basic features of the eco-evolutionary dynamics remain unexplained. Most prominently, there is no clear mechanism for the long-term coexistence of both drug-sensitive and resistant strains at intermediate levels, a ubiquitous pattern seen in surveillance data. Here we show that accounting for structured or spatially-heterogeneous host populations and variability in antibiotic consumption can lead to persistent coexistence over a wide range of treatment coverages, drug efficacies, costs of resistance, and mixing patterns. Moreover, this mechanism can explain other puzzling spatiotemporal features of drug-resistance epidemiology that have received less attention, such as large differences in the prevalence of resistance between geographical regions with similar antibiotic consumption or that neighbor one another. We find that the same amount of antibiotic use can lead to very different levels of resistance depending on how treatment is distributed in a transmission network. We also identify parameter regimes in which population structure alone cannot support coexistence, suggesting the need for other mechanisms to explain the epidemiology of antibiotic resistance. Our analysis identifies key features of host population structure that can be used to assess resistance risk and highlights the need to include spatial or demographic heterogeneity in models to guide resistance management.

## Author summary

The burden of drug-resistant bacterial infections is rising, and the fear that we are nearing a “post-antibiotic” era has seeped into the public consciousness. Scientists and public health officials often rely on mathematical models to predict changes in resistance levels over time and the effects of hypothetical interventions. However, most models struggle to reproduce common trends seen in real-world data, limiting their practical use. Here we propose a simple model to account for variations in the likelihood of taking antibiotics if infected, which arise within and between regions due to factors like drug-prescribing practices, healthcare access or care-seeking behavior, or the co-occurrence of other diseases. This model extension robustly reproduces trends seen in data, such as sustained coexistence of both drug-resistant and drug-sensitive strains of bacteria, and differences in resistance levels between similar or adjacent regions.

**Citation: **Krieger MS, Denison CE, Anderson TL, Nowak MA, Hill AL (2020) Population structure across scales facilitates coexistence and spatial heterogeneity of antibiotic-resistant infections. PLoS Comput Biol 16(7):
e1008010.
https://doi.org/10.1371/journal.pcbi.1008010

**Editor: **Matthew (Matt) Ferrari,
The Pennsylvania State University, UNITED STATES

**Received: **February 18, 2020; **Accepted: **June 2, 2020; **Published: ** July 6, 2020

**Copyright: ** © 2020 Krieger et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability: **The data presented in Fig 1, S1 and S2 Figs were collected from publicly available sources that are cited in the figure captions. All equations and algorithms needed to generate the results of this paper are described in the Methods and S1 Text, and example code to simulate the model equations is available at https://github.com/alsnhll/SpatialAbResistance.

**Funding: **This work was supported by the National Institutes of Health (DP5OD019851, ALH) and the Bill & Melinda Gates Foundation (OPPOPP1148627, ALH). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests: ** The authors have declared that no competing interests exist.

## Introduction

Antibiotic resistance is a major threat to our ability to treat bacterial infections. Over the past century, resistance to each new class of drugs has appeared soon after clinical use began. Today, drug-resistant infections are estimated to cost perhaps $20 billion annually [1]. Individual bacteria that are resistant to multiple classes of antibiotics are now common in species such as *Streptococcus pneumoniae*, *Pseudomonas aeruginosa*, and *Clostridium difficile* [2], and nearly untreatable strains of *Neisseria gonorrheae* [3], *Klebsiella pneumoniae* [4], and *Acinetobacter baumannii* [5] have recently been identified. These trends have led to speculations about a “post-antibiotic future”, in which routine medical procedures such as surgeries, childbirth, and dental work might become as high-risk as they were in the pre-WWII era due to a lack of effective antibiotic prophylaxis and treatment [6, 7].

Beyond sensational headlines, there is widespread interest among experts in predicting the future morbidity, mortality, and economic impact of drug-resistant infections, so that appropriate investments to counter these trends can be encouraged from government and industry. Mathematical modeling has traditionally played a key role in predicting the dynamics of infectious diseases [8], and has a long history of application to antibiotic-resistant infections [9]. Despite this, there are several recurrent trends in the spatiotemporal dynamics of drug-resistant bacteria that are difficult to explain. Firstly, for many bug-drug pairs, resistant strains have not completely displaced drug-sensitive ones, but instead coexist stably at intermediate levels, in contrast to predictions of standard infection dynamics and the ecological principle of competitive exclusion (Fig 1i). Secondly, regions which prescribe similar levels of antibiotic can have very different levels of resistance, and these differences persist over time (Fig 1ii). Finally, even neighboring regions can have large and persistent differences in resistance frequencies, at the scale of neighborhoods all the way up to countries (Fig 1iii). Since models must, at minimum, be able to explain current trends before being trusted for forecasting, these fundamental disagreements with data have either discouraged efforts to make predictions or led to widespread suspicion of existing predictions (e.g. the Review on Antimicrobial Resistance [10, 11]).

**I)** Percentage of *S. pneumoniae* isolates in Spain which are resistant to erythromycin (a macrolide), over time. After an initial increase, drug-sensitive and drug-resistant strains have appeared to coexist stably for ∼ 15 years. [12]. **II)** The same amount of antibiotic use can lead to very different levels of resistance in different regions. Percentage of *S. pneumoniae* isolates from 2016 that were resistant to macrolides versus the number of daily doses of macrolides administered per capita, for each country participating in the European Centres for Disease Prevention and Control [13–16]. **III)** Neighboring regions can have vastly different rates of resistance which persist for long time periods. Time-averaged percent of *K. pneumoniae* isolates resistant to carbapenems in Austria, Belgium, Croatia, France, Germany, Italy, Luxembourg, and Slovenia from 2013-2016 and Switzerland from 2015-2016 [14, 17] (blue/yellow), along with *S. pneumoniae* isolates resistant to penicillin in six provinces of Spain from 1990-1998 (green/red) [18, 19]. The year-to-year deviation from this average is less than 3% (see S1 Fig). Each country is labeled with the resistance level (%). Note that for *S. pneumoniae* and macrolides or penicillin, “% resistant” includes isolates classified as “non-susceptible”.

In this paper we propose a mechanism to explain these perplexing spatiotemporal trends in antibiotic resistance levels that have eluded standard infection models. We consider that host populations may be structured, with heterogeneities in mixing patterns within and between regions, as well as in the distribution of antibiotic use. Spatial differences in treatment rates could arise, for example, from local differences in prescribing guidelines or norms, availability of care, or care-seeking behavior; from the presence of hospitals or other facilities with higher rates of antibiotic usage, or, an increased use of antibiotics following spatially-localized viral epidemics. We find that in contrast to well-mixed populations, structured populations often support long-term coexistence of drug-resistant strains at intermediate levels. We observe that the predicted prevalence of resistance depends not only on the frequency of antibiotic use, but also on how drug use is distributed and details of the transmission network. This mechanism can also reproduce the observation of sharp gradients in resistance levels between neighboring regions, even when infections can spread directly between them. We discuss how these ideas can be used to better understand factors promoting or hindering antibiotic resistance and to construct more realistic models to evaluate antibiotic stewardship policies.

## Results

### Observed spatiotemporal trends in antibiotic resistance

#### I) Frequencies of resistance over time seem to defy competitive exclusion.

For a very large number of bacteria-antibiotic combinations, the prevalence of resistant strains has not reached 100% in the population, but has rather appeared to saturate at an intermediate frequency (Fig 1i). Prominent examples include *Escherichia coli* and aminopenicillins in Europe over the past decade [14, 20], methicillin-resistant *Staphylococcus aureus* in the United States [2], and penicillin resistance in *S. pneumoniae* in Spain from the 1980’s to 1990’s [18, 19] as well as in other locations/timeframes [21, 22]. This long-term coexistence of drug-sensitive and drug-resistant strains is difficult to explain. Resistant strains are clearly selected for in treated individuals, but generally carry a fitness cost such that sensitive strains do better in untreated individuals [23–25]. Mathematical models of infection dynamics under treatment predict that in most conditions the population will move towards an equilibrium where only a single strain persists, with drug-resistant strains reaching fixation when treatment is common enough and costs of resistance are low enough, and drug-sensitive strains dominating otherwise [26]. More generally, long-term coexistence of two strains would seem to defy the ecological principle of competitive exclusion [27–29], which dictates that when two species compete for a single resource (here, the susceptible population), only one may survive.

Dozens of previous mathematical models have been constructed in an attempt to explain this coexistence. While many studies have focused in particular on *S. pneumoniae*, these and others have attempted to be as general as possible so the conclusions can be extended to other bacterial species. Studies have repeatedly found that well-mixed populations in which individuals can only be infected with one strain (either sensitive or resistant) at a time do not generally support coexistence [26, 30–32]. Models that include more details of within-host infection dynamics have suggested some potential mechanisms of coexistence. Drug resistant and sensitive strains can coexist if dual infection of individuals with both is possible, but only for a very narrow range of values for both the treatment coverage level and the cost of resistance [26, 30, 31]. If *de novo* generation of one strain in an individual infected with the other can occur (presumably via point mutation or horizontal gene transfer from non-pathogenic co-colonizers), then coexistence is again possible [26], but such a process is likely to lead to low-level mutation-selection balance, not coexistence at near intermediate levels. Recently, Davies et al. [32] have suggested that co-colonization of hosts with subsequent within-host competition results in a type of frequency-dependent selection that helps maintain coexistence. This occurs because resistant strains have an advantage when they are rare, because they will mostly co-colonize drug-sensitive hosts who, when treated, will be rid of competitors. This mechanism is quite robust to parameter values, but only relevant to commensal bacteria for which co-colonization is common.

Another set of mechanisms shown to promote coexistence is host population structure. If there are two completely isolated sub-populations [26], or if treated and untreated individuals interact extremely rarely [30], then sensitive and resistant strains can both persist, but this is unlikely to apply to any realistic scenario. Blanquart *et al*. [33] recently extended this idea to populations consisting of an arbitrary number of completely interconnected sub-populations, and similarly found that coexistence could occur in some regions of parameter space but only when transmission between the sub-populations was weak. A study simulating transmission in a population in which both contact patterns and the probability of being treated are age-dependent and informed by data for *S. pneumoniae* also reveals more opportunity for coexistence [31]. However, this explanation is still far from general, since coexistence only occurs for costs of resistance <8%, does not display the full range of resistance prevalence levels—and their dependence on antibiotic use—observed in data, and is somewhat reliant on the division of the bacterial population into discrete serotypes which each support different resistance levels.

There are also ecological effects known to promote coexistence of species more generally [34], but their relevance to bacterial resistance to antibiotics remains unclear. For example, coexistence can occur when species compete more strongly among themselves than with other species [30]. The obvious candidate for such a mechanism in the context of infectious diseases is strain-specific immunity [35], with the effect that hosts are less susceptible to re-infection by a strain with which they have previously been infected. Strain-specific immunity leads to balancing selection, since low frequency strains have an advantage. However, there is generally not expected to be a connection between the resistance status of a strain and its immunogenicity (e.g. serotype). Lehtinen et al. [36] have recently shown that such a connection may not be necessary, as linkage between a locus under balancing selection and a polymorphic locus that influences the relative fitness of resistant and sensitive strains (such as duration of carriage for *S. pneumoniae*) can promote coexistence. The relevance of these mechanisms to most antibiotic resistant bacteria remains uncertain.

Overall, despite some progress in explaining the long-term coexistence of drug-sensitive and drug-resistant strains, many proposed mechanisms only reproduce coexistence for small regions of parameter space or species-specific scenarios, and the ultimate set of causes for real-world coexistence is far from fully understood.

#### II) Different frequencies of resistance are seen in regions with similar rates of drug prescription.

Another puzzling trend in spatiotemporal antibiotic resistance data is the appearance of very different levels of resistance in regions that seem to have the same levels of drug usage. While overall, there is a strong correlation between the amount of antibiotics consumed in a region and the prevalence of drug resistant strains [37, 38], there are many cases that diverge from this trend. For instance, comparing countries participating in the European Centre for Disease Prevention and Control (ECDC), the rate of resistance can vary by as much as 30% even between countries that have equivalent rates of outpatient prescription of penicillin, such as Spain and Poland [13–16] (Fig 1ii). Other examples can be found at smaller spatial scales, for instance between cities in the American South [39] and provinces in Spain [19]. While surveillance practices can differ in across regions and nations, the fact that these trends are seen in both ECDC data (where there is significant but not complete standardization) as well as in the smaller-scale studies conducted by single research teams with unified protocols over long time periods, solidifies that they are not artifacts of recording practices.

While this finding has not yet been specifically examined in the context of theoretical models, all of the models mentioned in the previous section [26, 30–33, 36] admit only one stable equilibrium for a given set of parameters, including antibiotic consumption level. Therefore, if each region is considered to be an isolated and well-mixed population, extra mechanisms are needed to explain different resistance levels for the same overall antibiotic use. It could be that one or both of the regions in question had not yet reached a steady-state prevalence of resistance by the time of sampling, though when longitudinal surveys of these same resistant strains are available, they generally suggest that levels have approximately equilibrated (e.g. ECDC data [14]). Another possibility is that the cost of resistance varies between regions, for example due to different mechanisms of resistance or different genetic backgrounds (e.g. spa types in *S. aureus* [40] or serotypes in *S. pneumoniae* [41]), but experimental evidence to support or refute this idea is lacking so far. Finally, the discrete regions considered during resistance surveillance are in reality neither completely isolated nor well-mixed. Movement of individuals and microbes between regions leads to interdependence of resistance dynamics. Within a region, there could be heterogeneous distributions of treatment and infections which impacts the overall resistance level. Consequently, models that attempt to explain this trend must consider the connected nature of regions across spatial scales.

#### III) Neighboring regions can show very different frequencies of resistance.

Another confounding element of antibiotic resistance epidemiology is the high degree of heterogeneity in resistance levels between neighboring populations at many differing scales. This is slightly different from the preceding point (II): not only do we find two regions with the same levels of antibiotic use but with different levels of resistance, we find that these regions can be bordering each other. Neighboring regions are likely to exchange infected individuals frequently, which would be expected to ameliorate differences in resistance levels over time, just like the predictions of well-mixed models. In contrast, they can actually sustain sharp gradients in the frequency of resistance seemingly indefinitely.

The scales at which this phenomenon is observed range from the size of nations down to neighborhoods of a city (Fig 1iii). For instance, in Europe, the frequency of carbapenems resistance in *K. pneumoniae* has been much much higher in Italy than in neighboring states for many years [14]. Similar trends have led to differences in the frequency of certain resistant strains across the United States [42, 43]. At a smaller scale, in the long-term study of *S. pneumoniae* strains in Spain revealed a large difference in frequency between penicillin resistance in Aragon and its neighboring regions [12, 18, 19, 44]. Among municipalities around Copenhagen, Denmark, the prevalence of resistance to trimethoprim and other antibiotics used to treat urinary tract infections in *E. coli* isolates varied by 3-fold [45]. At an even smaller scale, recent (2015) data on the United States revealed that in three neighborhoods in Greater Miami, rates of ciprofloxacin resistance in *Proteus mirabilis* differed by a staggering 62% overall, with rates of 13% (Fort Lauderdale, 172 cases), 41% (Miami, 264 cases), and 75% (Pembroke Pines, 8,979 cases) [39]. The fact that these large differences in resistance levels appear frequently in data at different spatial scales and between neighbors that disburse similar amounts of drug therefore require a unifying explanation.

### A general model for the evolution of drug resistance in a structured population

To better understand the mechanisms that could be responsible for the spatiotemporal motifs seen in drug resistance data, we developed a simple model for competition between strains of an infection in a structured population. We assume that the total population is divided up into multiple subpopulations (also known as ‘demes’ [46, 47]) (Fig 2). These demes could represent any subdivision of a human population of interest, such as different countries, regions within a nation, neighborhoods within a city, demographic groups, households, and so forth. Within this population we consider the concurrent spread of drug-sensitive and drug-resistant strains of an infection. Infected individuals spread the disease to uninfected individuals at rate *κ* if they are in the same deme (‘within-deme’ transmission rate) and rate *β* if they are in different but connected demes (‘between-deme’ transmission rate), with *κ* ≥ *β*. We do not allow for any super-infection or co-infection: only susceptible individuals can be infected. Infected individuals recover and become susceptible again at rate *g*.

**A)** An example population, divided into six equal subpopulations (“demes”, black squares) of five individuals (circles). Infection can spread within demes, and also between demes that are connected (black lines). Individuals are categorized based on their infection status (uninfected: open circle, infected with drug-sensitive strain: green circle, infected with drug-resistant strain: red circle). The deme where an individual is located may determine whether or not they will receive drug treatment (blue shading), or more generally, their probability of receiving treatment. **B)** Untreated deme: The wild-type strain (green) is transmitted at rate *κ* within a deme (bottom) and rate *β* between demes (top). Individuals infected with any strain recover with a rate *g*. The resistant strain (red) pays a cost *c* in its transmission rate with or without drug. **C)** Drug-treated deme: Transmission of the resistant strain is unaffected by whether the source individual is receiving treatment, but transmission of the wild-type strain (green) both within (bottom) and between (top) demes is reduced by a factor (1 − *ϵ*) if the source individual is treated, where *ϵ* is the drug efficacy.

Each individual has some probability of receiving drug treatment dependent on the deme they belong to. For individuals infected with the drug-sensitive strain, treatment reduces the rate at which they transmit the infection to others (in the same or connected demes) by a factor (1 − *ϵ*), where *ϵ* is the efficacy of the drug, 0 < *ϵ* < 1. We assume that the resistant strain is perfectly resistant to the drug, so that transmission rates of individuals infected with it are unaffected by the presence of treatment. However, the drug-resistant strain incurs a fitness cost *c*, which results in lower transmission rates relative to the wild-type strain with or without treatment: *κ* → *κ*(1 − *c*), *β* → *β*(1 − *c*). Another possibility for the effect of treatment is that it accelerates disease recovery, and we consider this possibility in S1 Text.

The “infected” status in our model could represent a symptomatic disease state caused by a pathogenic microbe, or an asymptomatic “colonized” state due to a commensal organism. Like others [26, 30–33, 36], we do not explicitly model virulence or mortality. Our model is also agnostic as to whether treatment is received only in response to infection (e.g. if a symptomatic infection triggers administration of antibiotics) or independently of infection status (e.g a commensal colonizer is subject to bystander selection from antibiotics administered for another condition [48, 49]), since we assume treatment does not impact susceptibility to catching an infection.

Our model is an extension of the classic Susceptible-Infected Susceptible (SIS) model for disease transmission [50] to a case with two disease strains (sometimes called an *SI*_{1}*I*_{2}*S* model [51]). We assume that each deme is large enough that the dynamics can be considered deterministically. This model is structurally neutral [52], meaning that if two identical strains are simulated (e.g. if *ϵ* = 0 and *c* = 0), they will continue to persist *ad infinitum* at the same level at which they were initiated.

### Absence of coexistence in an unstructured population

We first consider the dynamics of competition between drug-resistant and drug-sensitive strains in a single large, well-mixed deme (an unstructured population) (Fig 3A). A fraction *f* of all individuals infected by either strain will receive drug treatment for the duration of infection. The dynamics are described by a set of four differential equations, tracking the proportion of total individuals in each of four infected states (drug-sensitive untreated, drug-sensitive treated, drug-resistant untreated, drug-resistant treated), while the remaining individuals are uninfected (see Methods).

**A)** A population of individuals in a single well-mixed deme, in which a fraction *f* receive drug treatment (blue haloes). Individuals may be uninfected (hollow blue circles), or infected with either the wild-type (green circles) or drug-resistant (red circles) strain. **B)** The total prevalence of infection (wild-type + drug-resistant) as a function of the fraction of treated individuals (*f*) for different costs of resistance (*c*). **C)** The % of infections that are drug-resistant as a function of the fraction of treated individuals (*f*) for different parameters. Infection switches between 0% and 100% resistant when *f* = *c*/*ϵ*. Coexistence never occurs. **D)** Schematic of a two-deme population (left-untreated, right-treated) and the two strains (green-wild type, red-resistant) considered in the model. **E-G)** Each panel shows the infection level (shading) as a function of the relative connectivity between demes (*β*/*κ*) and the cost of resistance (*c*). **E)** The % of all infections that are drug-resistant strain across the entire population. **F)** The % of individuals in each deme who are infected with the wild type strain. **G)** The % of individuals in each deme who are infected with the resistant strain.
For all results, the transmission rate is *κ* = 0.25/day, the recovery rate is *g* = 0.1/day, and the treatment efficacy is *ϵ* = 0.9.

The basic reproductive ratio (*R*_{0}), defined as the average number of secondary infections caused by a single infected individual in an otherwise uninfected population, can be defined for each disease strain in this model and completely determines the equilibrium behavior (details in S1 Text). If both infections have *R*_{0} < 1, then neither strain of infection will persist and the equilibrium consists only of uninfected individuals. Otherwise, only the strain with the larger *R*_{0} will persist (at a prevalence 1 − 1/*R*_{0}, Fig 3B) while the other will go extinct. Based on the formulas for *R*_{0}, the resistant strain persists if and only if the coverage and efficacy of treatment are enough to offset the cost of resistance, *ϵf* > *c* (Fig 3C). Therefore, when the population is well-mixed, the heterogeneity introduced by partial treatment coverage (0 < *f* < 1) is insufficient to allow coexistence of wild-type and drug-resistant infections, and competitive exclusion always occurs, in agreement with other studies [26, 30–32].

### Limited coexistence in two connected subpopulations with heterogeneous treatment distribution

We next consider the case of two separate but connected subpopulations of equal size (Fig 3D). As an extreme example, we assume that individuals in one of the demes always receive treatment, while individuals in the other deme are never treated. The system is described by a set of equations tracking the proportion of individuals of each infection status (uninfected, infected with drug sensitive, infected with drug resistant) in each deme (treated vs untreated). The basic reproductive ratio can be derived for each strain ( and ) (see Methods).

This population supports qualitatively different infection dynamics than the single-deme case (Fig 3E–3G). While most parameter regimes still lead to persistence of only the drug-sensitive or the drug-resistant strain (even when *R*_{0} > 1 for both), there is now a stable equilibria where there is a mixture of both strains coexisting (Fig 3D). This region of coexistence occurs when the between-deme infection rate is relatively low compared to the within-deme rate (*β*/*κ* < 0.7 when *ϵ* = 0.9) and for intermediate values of the cost of resistance *c*. The mathematical derivations of these boundary regions and their physical intuition is included in S1 Text. Interestingly, in this mixed equilibrium, both demes contain a mixture of sensitive and resistant strains, a stricter type of strain coexistence than a situation in which each strain dominates a different deme.

This simple example, consisting of only two preferentially-mixing subgroups of equal size, demonstrates that co-clustering of transmission and propensity for antibiotic use can lead to stable coexistence of sensitive and resistant strains at not only the population level but also the *sub*population level, suggesting a resolution for the ubiquity of coexistence in empirical data (observation **(I)**). This host population structure produces a situation in which even though a strain is disfavored in a global sense by having a smaller *R*_{0} value than its competitors, it can be locally favored and thus avoid extinction. However, in this simple example, the region in parameter space where coexistence is seen is relatively small, even for this extreme example of treatment clustering, motivating the question of whether more complex structures will expand the stability of coexistence. Moreover, this two-deme example cannot explain the observation that regions receiving the same amount of treatment or regions neighboring one another often support very different levels of resistance (observations **II** and **III**). Replicating real-world data requires finding situations in which neighbors with *the same* level of drug treatment still have very different levels of resistance. In the subsequent sections, we examine whether more complex population structures can recreate these observed patterns.

### Robust coexistence in more complex population structures

Our results for simple one- and two-deme scenarios suggest that host population structure can promote coexistence of drug-sensitive and resistant strains of an infection. To investigate this relationship further, we extended our model to include an arbitrary number of demes and connectivity patterns (Fig 2A). For simplicity, we first assumed that each deme in the population was the same size and was connected randomly to a fixed number of other “neighbor” demes, therefore creating a collection of random regular graphs [53] (Fig 4A). Although in reality much more complex human population structures are possible, these simplified networks provide a good base to examine the impact of specific properties on infection dynamics. Each deme was independently assigned to be “treated” (meaning individuals in that deme receive treatment, while individuals in “untreated” demes do not). We relax these assumptions in later sections.

**A)** Randomly generated population structures on which infection was simulated. Each node represents a deme (a well-mixed sub-population of individuals), and each edge indicates that infection can spread in either direction between those two demes. Ten example populations were selected out of 1000 total simulated, each with twenty demes randomly connected to three neighbors each, to represent a broad range of outcomes. **B)** Fraction of infections that are resistant in the entire population (*y*-axis) versus fraction of demes treated, *ρ* (*x*-axis). Each color represents a different parameter set (blue background—baseline, red background—lower cost of resistance, teal background—more between-deme connectivity). Numbers show data points for the ten example populations. The colored envelope is created by shading between sigmoidal curves that encompass all the data. **C)** For each population structure shown (*y*-axis) and each treatment level (*x*-axis), the proportion of simulations that resulted in *robust* coexistence between drug-sensitive and drug-resistant strains is shown (by the colored area of the box). Robust coexistence was defined as at least 80% of demes supporting both strains at frequencies above 10%. **D)** Differences in resistance levels (% of all infections that are with the drug resistant-strain) are measured between all pairs of directly-connected untreated demes. **E)** Histograms showing the distribution of pairwise differences in resistance for a given population structure. Lighter shaded histograms combine results from all population graphs. All simulations used kinetic parameters *κ* = 0.25/day, *g* = 0.1/day, and *ϵ* = 0.9, and pooled results from 100 simulations with different random allocation of treatment across demes. Pairwise differences were calculated with 30% treatment.

Using these “meta-populations”, we found that coexistence is possible for a broad range of costs of resistance, treatment levels, and population structures. As observed in data, the overall prevalence of resistance in a population is roughly correlated with the treatment coverage, increasing gradually as treatment levels increase (Fig 4B). Moreover, when resistance levels were intermediate, coexistence also persisted at the level of individual demes. We defined ‘robust’ coexistence to mean that at least 80% of demes had at least 10% of infections caused by each strain. For the baseline parameter values we used a 20% cost of resistance, treatment levels between 20 and 40% supported robust coexistence in at least some population structures and treatment allocation schemes. For a lower cost of resistance (*c* = 0.05), robust coexistence occurred at a lower range of treatment levels (Fig 4C). When the amount of infection spread (*β*) arising from contact between “neighboring” demes increases relative to the within-deme spread (*κ*), resistance increases more sharply with treatment level, but regions of robust coexistence remain (teal color, and more parameter values in S14 and S15 Figs).

These results show that even relatively simple host population structures can reconcile mathematical models with empirical observations of long-term coexistence between drug-sensitive and drug-resistant strains. These results also provide a connection between our model and the data for another spatio-temporal observation, **(II)**: the same overall rate of drug consumption can lead to different prevalence of drug resistance in different populations (even if treatment efficacy and cost of resistance are the same). For example, in Fig 4B, a particular treatment level (e.g. 20% for the red parameter set) can be associated with very different frequencies of resistance in different population structures (up to 30% difference between graph 6 and 9, similar to the data in Fig 1). A more detailed study of the role of higher-order network properties is presented in a later section.

### Large differences in resistance levels are possible even between connected regions

With these complex multi-deme populations, we examined whether neighboring regions with equal levels of antibiotic consumption could sustain vastly different amounts of resistance (observation **(III**)). For each population structure, random treatment allocation, and parameter set, we chose all nearest-neighbor pairs which were both untreated, and generated a distribution of the pairwise differences in resistance levels between these neighbors (Fig 4D and 4E). We found that large differences in resistance levels between neighbors were common, with 34% of pairs receiving no treatment differing by more than 10% resistance prevalence, for our baseline parameter values when averaging across all population structures (Fig 4E). In the same simulations, 2% of pairs overall—and up to 15% in some structures—differed by greater than 30% resistance frequency, a value near the upper limit of observations from data shown in Fig 1. Differences of up to 60% between neighbors were observed in some individual simulations.

While populations with higher levels of between-deme mixing (higher *β*/*κ*, teal color) often supported persistence of both drug-sensitive and drug-resistant strains at equilibrium (Fig 4B), they supported smaller differences in resistance levels between demes receiving the same treatment (Fig 4E). In this parameter regime, sensitive and resistant strains are more segregated between untreated and treated demes, respectively, explaining why simulations on these structures rarely produced “robust” coexistence as we defined it (over 80% of demes supporting at least 10% of each strain) (Fig 4C). We also found that populations with weak connections between all demes were less likely to have large pairwise differences in prevalence of resistance (S11 Fig). Overall, these findings corroborate the intuition that adding edges between demes or raising the between-deme infection rate of connected demes brings the system closer to well-mixed population and hinders coexistence.

### Properties contributing to higher resistance

In order to understand how specific properties of host population structure contribute to the frequency of resistance and the likelihood of coexistence, we simulated infection dynamics on a large set of transmission networks and then statistically analyzed the results. To first generate a large ensemble of population structures which varied in many graph-theoretic properties, we used the Watts-Strogatz [54] algorithm to create 1000 unique networks of 50 demes each, which had different average degree, variance in degree, clustering, centrality, efficiency, and mean path length (see Methods for definitions of these properties and S3 Fig for their values). For each graph, we generated 50 different allocations of treatment, keeping the overall population treatment level the same. Then, we collected the results and employed LASSO (least absolute shrinkage and selection operator) regression to identify the most important properties contributing to the frequency of resistant infections (Table 1). The regression was conducted both at the level of individual demes (i.e. to identify local properties of demes that influenced the level of resistance within a deme) and at the level of the whole population (i.e. to identify properties of the whole population that influenced the overall level of resistance). The simulations and the regression were repeated for a “low resistance” case, in which overall 24% of individuals in the population were treated and as a result less than half of infections were resistant at equilibrium, and a “high resistance” case, where 40% were treated and over half of infections were resistant (S4 Fig). Properties were classified based on whether they were always associated with an increase or decrease in resistance, or whether they pushed resistance towards intermediate levels (i.e. facilitated coexistence) or towards extreme levels (i.e. hindered coexistence).

We first examined predictors of the resistance level within a deme (Table 1, left side). The most important contributor to resistance was whether or not the deme itself was treated, and the distribution of treatment among the deme’s first, second, and third nearest neighbors were also important– albeit less so—predictors. The next strongest predictor was the total number of first degree neighbors of a deme (the “degree”). Higher degree is associated with lower resistance when resistance is rare, but higher resistance when it is common, and hence pushes resistance to more extreme levels, acting against coexistence of both strains within a deme (− /+ in Table 1). The same is true for the eigencentrality of the deme, a measure of its overall level of connectivity within the network. These results imply that highly-interconnected, hub-like demes are more likely to harbor the more common strain, while demes that are on the periphery of a network are more likely to support coexistence of the strain that is rarer in the population overall. Hence, when resistance is rare, it may be preferentially found in less connected demes, whereas when it is more common, it may cluster in more central demes.

We next evaluated predictors of resistance at the level of the entire population (Table 1, right side). Surprisingly, the simplest statistics of a graph—including the mean and variance in the number of neighbors of each deme—were the least important contributors to resistance. Rather, the best predictors involved a mix of treatment-clustering properties, like the proportion of sets of three connected demes where two of the demes were treated (“U-T-T”), and node-clustering properties, such as measures of “efficiency”. Treatment promotes consistently higher resistance levels when it is distributed more evenly throughout the population (e.g. treated and untreated demes are interspersed in the network, “U-T-T”, “U-T”), but more extreme clustering of communities receiving the most treatment (“T-T-T”) instead facilitates coexistence. This suggests that the effect of increasing antibiotic consumption on resistance levels will depend strongly on which communities experience the increases. The clustering coefficient and local efficiency are both measures of how interconnected the neighbors of a given deme are, on average. Global efficiency and mean path length are both measures of the ease of moving between any two random demes in the whole population. Coexistence is promoted by higher mean path length, which leads to graphs with more segregated transmission clusters, but hindered by global and local efficiency. The more interconnected a network is, the harder it is for multiple strains to coexist. The reason why the clustering coefficient always acts against resistance is not entirely clear, but its interpretation is complicated by the fact that the treatment-clustering properties are partially determined by the level of clustering in the graph.

### Robustness of results to mechanism of treatment action

Throughout this paper, we have modeled the effect of antibiotic treatment as reducing the ability of an infected individual to spread infection to others. Hence, treated individuals are still infected, and therefore immune to infection with the other strain, but do not contribute to transmission (or contribute less, if treatment is imperfect, *ϵ* < 1). Alternatively, treatment could instead accelerate the rate at which an individual recovers from infection, a mechanism that has been used in many previous models of antibiotic resistance [30, 31]. To determine if the mechanism of treatment action had any influence on our findings about the spatiotemporal patterns of resistance and to facilitate comparison of our results with those of other models, we considered a variant of the model where in treated individuals, the recovery rate increases from *g* → *g* + *τ* (see S1 Text. for details). We found that all of the results reported in the main text are recapitulated when treatment increases recovery rate, including the absence of coexistence in partially-treated well-mixed populations (Fig 1C), a limited though larger parameter regime of coexistence in two-deme populations (S7E Fig), common coexistence in multi-deme populations (S8B and S8C Fig) and divergent resistance levels between neighboring demes receiving the same treatment level (S8E Fig).

### Robustness of results to variations in treatment distribution and community structure

The results presented so far consider communities of demes where each deme is connected to the same number of other demes (e.g. “uniform networks”), where each deme is of the same size, and where there is either a 0% or 100% chance of receiving treatment if infected. While these simplifications allowed us to evaluate whether population structure alone could recreate observed spatiotemporal resistance trends in a minimal model, it is obviously not capturing many sources of variability that occur in the real world. Therefore, we also considered how relaxing each of these assumptions impacted our results. Populations described by graphs with heterogeneous connectivity (mean degree = 3, variance up to 10, S9 Fig) supported similar levels of coexistence and between-neighbor differences in resistance levels as homogeneous graphs. Variation in deme size also did not seem to change the overall results (coefficient of variation up to 30%, S10 Fig). We also created networks in which every deme was connected to every other deme at a fractional weight relative to the main network connections (S11 Fig). Only as this level approached 1 did the likelihood of coexistence shrink. Intermediate levels of this background mixing (e.g. 0.01, 0.1) actually promoted robust coexistence (S11C Fig), by allowing wild type and resistant strains to distribute throughout the population while still maintaining high and low treatment niches that selected for one or the other. Finally, we considered continuous distributions of treatment levels across demes, ranging from bimodal distributions where low and high levels of treatment were more common than intermediate levels, and unimodal distributions where the average treatment level was most likely and there was only a small variation between demes (S12 Fig). In general, extreme differences in treatment levels support coexistence over a wider range of parameters (S13 and S14 Figs), but are not strictly necessary for it to occur and be distributed across most demes in the population (S15 Fig). When differences in antibiotic consumption across groups is small, coexistence only occurs when the average population-level drug efficacy is roughly equal to the cost of resistance (*ρϵ* ∼ *c*, assuming groups are otherwise equal). In this case, other mechanisms may be needed to explain the pervasiveness of coexistence in the real world.

## Discussion

The prevalence of antibiotic-resistant bacterial infections displays several surprising spatial and temporal trends that defy canonical predictions of infection dynamics. The goal of this paper was to examine a simple and flexible model for competition between drug-sensitive and drug-resistant strains of an infection in a structured host population to see whether it could capture several of these trends. The first perplexing pattern observed in the data is the long-term coexistence of both sensitive and resistant strains at intermediate levels, which many other studies have attempted to explain [26, 30–33]. Our simulations reproduce coexistence for a wide range of parameter values and suggest that it is a natural feature to be expected in spatially (or otherwise) heterogeneous populations with competing infection. In addition, we also captured more nuanced trends observed in surveillance data that other models have ignored. Regions that administer similar amounts of a particular antibiotic may experience persistent differences in drug resistance levels for many years, and these “frozen gradients” are observed even if regions are neighboring one another. These phenomena were reproduced even in our simple model (Fig 5A), but also occurred for more realistically variable population structures and treatment distributions (S9–S13 Figs). We found that in a collection of partially-mixing subpopulations, subtle details about the pattern of connectivity and how treatment is distributed between them can lead to large differences in resistance prevalence between regions that otherwise seem similar (Fig 5A). We were able to identify specific covariates that contribute to resistance levels and to promoting coexistence with drug-sensitive strains. Similar to the use of “risk-mapping” to predict the distribution of vector-borne diseases based on environmental [55–58] and human [58–62] factors, these findings suggest that spatial risk-based assessment may be useful for the study of drug resistance. Surveillance on a finer scale will be needed to test these ideas on empirical data.

**A)** Our analysis suggests two possible explanations for large differences in resistance levels between two regions (orange and blue) receiving similar amounts of antibiotics. **Left:** The regions may differ in the distribution of antibiotic consumption within the population and its connection with the underlying transmission network. **Right:** The regions may differ in their connectivity to other regions which consume different amounts of antibiotics. **B)** Predictions about the impact of interventions on future resistance levels can be incorrect if population structure isn’t accounted for. **Left:** A population consisting of 100 sub-populations, with 20% treatment coverage at baseline. **Right:** Simulations of resistance levels in the “true” structured population (blue) with and without a hypothetical intervention (reduction to 15% treatment) applied at year 5. For comparison, predictions of a well-mixed model that with the same parameters (black), which approximates the dynamics before year 5 but diverges afterwards. The well-mixed model would predict large reductions in resistance in the short term but eventual fixation of resistance nonetheless (black dashed line), whereas the structured model predicts modest but sustained reductions which depend on the particular sub-populations targeted (blue dotted lines), and long-term coexistence. All simulations used kinetic parameters *κ* = 0.25/day, *β* = 0.05/day, *g* = 0.1/day, and *c* = 0.2. For the left of (A) and for (B), *ϵ* = 0.9, and for the right of (A), *ϵ* = 0.5. The population in (B) is a Watts-Strogatz network with degree 4 and re-wiring probability 0.1.

Increasing levels of antibiotic resistance are widely considered to be a major public health threat, and attempts are continually underway to forecast resistance levels with or without additional interventions using at least some type of mathematical modeling. Our findings highlight the fact that infection models which rely on assumptions of well-mixed populations are unlikely to be a useful tool for this task. As these models cannot reproduce ubiquitous trends observed in existing data on drug resistance epidemiology, they are unlikely to make predictions that are even qualitatively trustworthy (Fig 5B). For example, it is very important to know whether resistance for specific bug-drug combinations will tend towards 100% or settle at an intermediate level, and so any model must be able to explain the phenomena of coexistence for most strains. While clearly any degree of resistance complicates clinical management, the economic impact of partial vs complete resistance in a population could differ dramatically. In addition, by ignoring the details of spatial or demographic heterogeneity, we may miss an opportunity to more effectively target intervention or surveillance to particular subpopulations. The obvious difficulty is that such models may have to be highly-tailored towards a particular disease, and currently, there is a stark lack of data required to construct or calibrate realistic spatial models. While good-quality data is available on antibiotic consumption and drug resistance at the level of US states [63] and European countries [14, 16], and could be combined with other datasets tracking inter-region travel, this is likely too large of a spatial scale to capture the relevant heterogeneity for this mechanism. More routine high-resolution mapping of resistance levels over space and time [39, 64], combined with spatial or demographic patterns of antibiotic consumption [45, 65–68] and data to estimate population connectivity patterns over which disease spreads, are sorely needed.

Throughout this paper we have considered a relatively simple type of population structure, while in reality it could be much more complicated. There may be multiple levels of structure—from a household to countries—with much more complicated patterns of connectivity and a continuum of transmission rates. Structures may be dynamic, as movement and interaction patterns of humans, animals, or disease-carrying material change over time. The mechanism leading to population structure varies depending on the infection of interest. For hospital-acquired infections, the movement of medical staff or the spatial arrangement of patients may be most important, while for community-acquired infections, the nature of close household, school, or workplace contacts may be of primary interest. For sexually-transmitted infections, networks of sexual contacts determine disease spread, while for diseases of livestock, relatively well-mixed contact within crowded farms along with patterns of transfer between farms may create an overall “meta-population” structure. Host population structures are notoriously difficult to measure in practice, but some approximation of them is likely needed to create realistic models for the spread of drug resistant infections. Technological advancements that allow better tracking of human mobility, such as wearable proximity censors [69–71] or fitness trackers, cell phone location data [72–75], or large-scale transit networks [76, 77], provide one source of population structure data. Advances in genetic epidemiology allow us to retrospectively trace the spread of infection between individual hosts or communities using pathogen phylogenies [78–82], which also provides information on population structure. Only with better empirical characterization of transmission-relevant population structures and variations in antibiotic use can we test whether the mechanism described here can explain coexistence for specific bug-drug pairs.

Our model of competition between drug-sensitive and resistant strains is deterministic; we do not consider the stochastic early-stage invasion of an initial resistant mutant or the possibility of extinction of either strain. This approach is justified as we are interested in the dynamics of resistant strains that are already well-established and at high prevalence in the population. Since most antibiotics are closely related to natural products, resistance genes also naturally occur in environmental microbes, making “extinction” of resistant strains unlikely. Previous work by our group and others has examined how spatial heterogeneity or transmission network structure can accelerate the initial probability of emergence of drug resistance mutations [83–89] and also on the loss of resistant strains when drug is withdrawn [90]. Models for the emergence of methicillin-resistant Staphylococcus aureus (MRSA) in which demes represent hospital vs community settings have found that stochastic effects can create a new regime of dynamics, whereby sporadic subcritical outbreaks can sometimes overwhelm control resources and push the system to endemic infection [91]. Even in our context of established resistance, some important differences could occur in a stochastic version of our model. If deme sizes are relatively small, the balanced coexistence of sensitive and resistant strains we observe at equilibrium could actually be a dynamic process of repeated colonization and extinction at the deme level (as described by metapopulation theory [92]). Extinction times for competing epidemics in homogeneous populations have only recently been estimated for homogeneous populations [93], and the added complexity of our model suggests we are far from an analytic theory of invasion and extinction for infections in structured populations.

Our paper focuses on explaining trends observed for antibiotic resistance, but drug resistance is also a problem for antimicrobial therapy more broadly, including antivirals and antiparasitic drugs. To our knowledge, long-term coexistence of resistant and sensitive strains despite constant drug exposure has not been observed for HIV or malaria, the viral and parasitic infections for which targeted therapy is most widespread. In the case of HIV, rapid within-host evolution, including both *de novo* generation of resistant strains and reversion to susceptible strain in the absence of treatment would be expected to promote coexistence even in the absence of population structure. However, the prevalence of resistance appears to be changing rapidly along with dramatic changes in antiretroviral treatment availability and quality [94]. For example, the prevalence of resistance at the time of diagnosis (before treatment administration) has been increasing rapidly in South Africa and other low and middle income countries [95, 96], where widespread access to treatment is recent but follow-up care and resistance testing are still limited. In contrast, this primary resistance has been steadily decreasing in Europe, where treatment has been common since the 1990s albeit originally with suboptimal drugs (e.g. [97]). For malaria, it appears that resistance levels in a population continually increase under drug exposure, as extremely high levels of chloroquine or sulfadoxine/pyrimethamine resistance have been observed (e.g. [98]). However, antimalarial treatment is not generally individualized and international guidelines suggest that country or region-wide treatment recommendations switch to new drug classes once suspected resistance is above a threshold of around 15-25%. Consequently, the removal of selection for those resistance strains and their subsequent decline(e.g. [99]) likely precludes the occurrence of stable coexistence. Interestingly, a recent modeling study examining the role of co-infection with resistant and sensitive strains of malaria showed how this can suppress selection for resistance [100], very similar to recent ideas suggested for promoting coexistence in commensal bacterial infections [32].

Spatial heterogeneity has been understood by population geneticists to be an important factor in the spread of genes since the days of Fisher [101] and later Kimura [102]. Other classic work uncovered mechanisms by which the coexistence—as opposed to competitive exclusion—of species could be promoted by spatial heterogeneity [103, 104]. More recently, Gavrilets & Gibson [105] examined a two-deme system with two competing types, each of which has higher fitness than the other in one of the two demes, and observed a phase diagram very similar to that for our two-deme infection model (Fig 3), in which polymorphic equilibria were possible. This idea was later extended to multiple connected demes [106], where it was demonstrated that spatial heterogeneity can preserve a species in spite of fitness differences which in a well-mixed model would drive it to extinction. Other work involving more complex population structures [107] or more complex fitness distributions across space [108, 109] has examined how the rate of invasion of rare mutants depends on spatial heterogeneity. This previous work all involved traditional population-genetic models, such as Moran or Wright-Fisher processes, and before this study it was not clear which results would generalize to infection dynamics models. For example, the role of network structure in modulating the fixation probability of new strains differs dramatically between Moran processes [107] and infection models [88].

Several other recent studies have examined the role of population structure in facilitating coexistence of drug resistant infections. Kouyos et al. [110] considered populations divided into community members and hospitalized patients, and examined if such structure could support the coexistence of two separate strains of methicillin-resistant *Staphylococcus aureus* (MRSA). Similar to our model one strain was more resistant yet more costly, and was favored in more highly-treated hospitalized individuals, while the other strain was less resistant but more transmittable. They found a relatively narrow parameter regime supporting coexistence, which could be increased only somewhat if they additionally structured the population into different age groups. Interestingly, they found that beyond the parameter region that supported coexistence at equilibrium, there was a much wider regime leading to decades-long transient coexistence, suggesting that our estimates of the extent of coexistence are likely very conservative underestimates. Cobey et al. [31] examined the role of age structure in supporting coexistence of drug sensitive and resistant strains of *S. pneumoniae*, and found that this mechanism could lead to coexistence but only for a relatively narrow range of costs of resistance. In addition, their model included a mixture of different serotypes, and they found that apparent population-level coexistence often occurred due to summing many unique serotypes which were each nearly all resistant or all sensitive. Similar to us, Blanquart et al. [33] considered a population consisting of a set of interconnected “demes” with different treatment levels, but in contrast, for mathematical tractability, they restricted their analysis to structures in which either all demes are connected to all others with equal weight, or are completely isolated. Their model was parameterized such that the maximal force of infection from all outside demes equaled to the force of infection from within the same deme, whereas in our model the force of infection from outside can be much larger than from inside. They also concluded that population structure can support coexistence, but in contrast to our findings, under their model assumptions coexistence only occurs with very weak inter-deme mixing. By considering larger and more variable transmission network structures and treatment distributions, our work thus suggests the potential for a more prominent role of population structure in supporting coexistence. Compared to previous work, our criteria for the occurrence of coexistence in models is more stringent, since we require at least 5% prevalence of each strain (compared to >0% [33, 110] or >2% [31]), and additionally require coexistence within each individual subpopulation (“robust coexistence”). Our paper additionally differs from these previous works by going beyond coexistence to examine other common spatial trends in resistance data (i.e. trends **(II)** and **(III)**), and by comparing the role of a range of spatial metrics describing population structure and treatment allocation patterns in promoting coexistence (Table 1).

While our results show that population structure can explain observed spatiotemporal patterns of antibiotic resistance, it is unlikely to be the only mechanism contributing to either long-term coexistence or differences between regions. Persistence of both drug-sensitive and drug-resistant strains is also promoted by co-infection and superinfection [30–32], by overlaps between resistance status and serotype [30, 31] or other traits under balancing selection [36], and by heterogeneities in transmission or recovery rates as opposed to treatment coverage [30, 31]. Regions receiving similar overall levels of drug but experiencing different frequencies of resistant infections could alternatively be explained by local differences in i) prescribed vs consumed doses, ii) who receives drugs (e.g. age group, hospital vs outpatient) and for what type of condition, iii) genetic background of the pathogen and the cost of resistance, iv) transmission or recovery rates, or v) strain-specific control measures. For example, take two regions with identical transmission networks, drug consumption, prescribing behavior, and infection burden. If one region institutes a control policy whereby individuals infected with the drug-resistant strain are isolated, then in that region the cost of resistance (*c*) is effectively higher—even for a genetically identical strain—and the resistance level will be lower. It is difficult to judge precisely how realistic our results are, since our models of infection dynamics and population structure are both dramatically simplified, but it likely cannot explain coexistence over the entire range of parameters for which it is observed empirically. In addition, in our results it is relatively rare to see large differences in resistance prevalence between regions with similar treatment rates unless they have very different connectivity patterns within or between them, and unless there are large spatial or demographic heterogeneities in treatment use that overlap with clustering of disease transmission. In reality, multiple mechanisms probably act in concert to explain the observed patterns. Moreover, there is unlikely to be a single explanation for each different bug-drug pair for which these trends are observed. Other disease-specific factors that may be relevant to the ecology of antibiotic resistance include whether the bacteria is primarily pathogenic or also a commensal colonizer, whether infection is more common in the community or in healthcare facilities, whether resistance is carried on plasmids or chromosomally, and whether there are environmental reservoirs and if they are exposed to antibiotics.

## Methods

### Model assumptions

Our two-strain model of infectious disease dynamics makes a number of simplifying assumptions. We assume that individuals can transmit the infection immediately upon becoming infected, hence ignoring any sort of “latent” phase or “exposed” class. Birth, death, and migration into and out of the population are not modeled. The infection is assumed to be non-lethal, and may even represent a non-pathogenic “colonized” state with a commensal organism. Infected individuals all recover to be susceptible again, which presumes that there is no long-lasting immunity to re-infection. We do not allow for any co-infection (simultaneous infection with both strains) or super-infection (replacement of one strain within an individual by another), which implies that there is complete cross-immunity during infection. We model resistance as a binary trait: there are only two strains, either completely sensitive to the drug, or completely resistant. In reality, there may be strains with intermediate levels of resistance.

Individuals who are assigned a status of “treated” will receive treatment immediately upon being infected, every time they are infected. This ignores the possibility that there may be a delay to receiving treatment or that treatment behavior may vary within an individual over time. Treatment is assumed to reduce transmission of the infection, but may not be perfectly effective at doing this (which could approximate the effect of treatment delay or imperfect drug efficacy). In S1 Text. we consider an alternative model in which treatment acts to increase the rate of recovery from infection instead. In reality, treatment likely has both effects, since it generally acts to reduce pathogen loads within a host, which leads to reduced contagiousness and reduced duration of infection. Although we track uninfected (susceptible) individuals as either “untreated” or “treated”, our results are independent of whether uninfected individuals actually ever receive drug.

Our population structure consists of *M* subpopulations (“demes”) of equal size *D* for a total population of size *N*. The structure is assumed to be static. For each model, the transmission rates *κ* and *β* were scaled so that the model behavior did not depend on *N*, *D*, or *M*. In the one- and two-deme case, *κ* was multiplied by *N*, and in the multi-deme cases, both *κ* and *β* were multiplied by *D*. A more detailed derivation of this scaling is given in S1 Text. This scaling is equivalent to saying that individuals are limited in how many contacts they can have with others, and that this limit is independent of the total population size.

By using differential equations to model infection, we assume that the population sizes of individual demes are large enough so that variation from the average behavior is not important. We additionally assume that resistance is pre-existing and that stochastic extinction of either strain from the population is not possible. For the most part, we examine only the equilibrium behavior of the models.

### Model equations

**Well-mixed population (single deme):** Infection spreads between all individuals with the within-deme transmission rate *κ*. A fraction *f* of all individuals are treated with drug immediately upon infection by either strain. We track the proportion of the total population who are infected with the wild type and untreated (*w*_{u}), infected with the wild type and treated (*w*_{t}), infected with the resistant strain and untreated (*r*_{u}), and infected with the resistant strain and treated (*r*_{t}). The number of uninfected (susceptible) individuals is given by *s* = 1 − *w*_{u} − *w*_{t} − *r*_{u} − *r*_{t}. The dynamics are then described by the following set of differential equations (see S1 Text. for a derivation):
(1)

**Two deme population:** Infection spreads between all individuals in the same deme with rate *κ*, and between any two individuals in different demes at rate *β*. The demes are of equal size *D*. Treatment is assigned to all individuals within only one of the two demes. We track the proportion of individuals in the treated deme that are infected with the wild-type/drug-sensitive strain (*w*_{u}) and infected with the resistant strain (*r*_{u}), and the same proportions for the treated deme (*w*_{t} and *r*_{t}). The number of uninfected (susceptible) individuals in the untreated deme is *s*_{u} = 1 − *w*_{u} − *r*_{u} and in the treated deme is *s*_{t} = 1 − *w*_{t} − *r*_{t}. The system is then described by the following set of differential equations (see S1 Text. for a derivation and generalization to unequal deme sizes):
(2)

**Multi-deme population:** Infection spreads between all individuals in the same deme with rate *κ*, and between any two individuals in different but connected demes at rate *β*. There are *M* demes each of size (*D*). Each deme may only be able to spread infection to a sub-set of other demes, and the connectivity of the population is described by the adjacency matrix Δ_{ij} (Δ_{ij} = 1 if an individual in deme *j* can be infected by an individual in deme *i*). The *degree* of deme *i*, *d*_{i}, is the number of neighbors it is connected to, *d*_{i} = *Σ*_{i}Δ_{ij}. Treatment is assigned at the level of the deme, described by indicator variable *T*_{i} (where *T*_{i} = 1 if deme *i* is drug-treated and *T*_{i} = 0 otherwise). The fraction of demes that are treated is *ρ*. The system of equations describing the fraction of individuals in each deme who are infected, with either the wild-type (*w*_{i}) or drug-resistant strain (*r*_{i}) is:
(3)

A derivation that is generalized to the case of demes of unequal size and with arbitrary treatment levels in each deme is provided in S1 Text. Note that in constructing the adjacency matrices Δ_{ij} we ensure that the graph contains a single giant component, and not multiple isolated parts.

Derivations for the equilibria, stability conditions, and basic reproductive ratio (using the next-generation technique [111]) are given in S1 Text.

### Numerical results

The differential equations were numerically integrated in MATLAB, using the Runge-Kutta solver provided in the function ode45. The initial condition for each parameter-population structure combination was each deme having a level of infection with the wild type that would occur if there were no resistant strain (a fraction 1 − 1/*R*_{0} of individuals are infected with the wild type), and a very small level of infection with the resistant strain (10^{−3}). The rest of the population was uninfected.

The system was integrated until an equilibrium was reached, which was defined as the point at which the sum of the time derivatives for the fraction of each strains in all demes became less than 10^{−4}/*M* per day, where *M* is the number of demes. To ensure that the equilibrium reached by the solver is stable, once an equilibrium was reached by the above standard, we applied a small random fluctuation with magnitude not exceeding 1% to the value of all strains in all demes. If the same equilibrium was arrived at after repeating this process three times, the equilibrium was taken to be sufficiently stable and used as the result of that simulation. We also explored using more stringent thresholds, but found that our results were nearly always within 0.01% of the value that would be achieved and reduced simulation times (S16 Fig). To determine whether there was only a single positive and stable equilibrium for the multi-deme system for a given structure and parameter set, we simulated the same systems many times with uniformly-at-random initial conditions. For each deme, the initial fractions [*w*_{i}, *r*_{i}] were both drawn from the uniform distribution on [0, 0.7]. If their sum exceeded 1, both were re-drawn. Previous work has proven mathematically that stable equilibria exist for these systems [112, 113]. We found that the same equilibrium was always reached (S17 Fig). More details are given in S1 Text.

### Network generation

**Random regular graphs:** Random regular graphs are randomly-generated graphs in which every node has the same degree (number of incident edges). There are multiple existing algorithms to create random regular graphs. Our random regular graphs were generated according to a pairwise-construction mechanism [114] implemented in MATLAB [115]. For most results presented in the paper we used networks with *M* = 20 total nodes and degree *d* = 3.

**Watts-Strogatz graphs:** The Watts-Strogatz algorithm is a method of constructing networks with *M* total nodes and *M***d* total edges, but with heterogeneous properties. In the original construction, a ring lattice is formed with each of the *M* nodes connected with *d* edges to neighbors in a ring formation. Then, with probability *p*, the target of each edge is rewired to a uniformly random node. For our 1000 graphs, we selected *M* = 50, but allowed degree *d* to be selected uniformly-at-random for each graph from {3, 4, 5, or 6}, and allowed *p* to be selected uniformly-at-random for each graph from [0.1, 0.6]. The resulting graphs have a broad distribution all graph properties of interest (S3 Fig) and of the frequency of resistance (S4 Fig).

### Network statistics

A set of graph-theoretical properties were tested for their relationship to resistance levels within demes and within populations. A description of these properties is given below.

**Eigencentrality**: For a node *i*, its eigencentrality is a measure of not only its centrality (having a high degree, for example), but is also weighted according to how high the centrality scores of its neighbors are; as an example, an academic paper with many citations may have high centrality, but if the papers citing it have themselves low centrality then the eigencentrality of the original paper may be low. Formally, the eigencentrality is the leading right eigenvalue of the adjacency matrix, and so satisfies the implicity equation *E*_{i}: *E*_{i} = (1/λ)∑_{j} Δ_{ij} *E*_{j}, where λ is a normalization factor.

**Global efficiency**: In general, the average efficiency of a graph is defined as , where *M* is the number of demes/nodes and *d*_{ij} is the smallest number of edges between nodes *i* and *j*. Therefore, graphs that have higher efficiency have on average shorter path lengths between nodes. Since *M* = 50 is fixed in all of our structures examined in this section, the only varying factor is the distances between nodes; the smaller the distances between all given nodes, the higher the efficiency. Global efficiency involves comparing the graph to the most efficient possible graph with as many nodes. Since this simply involves another prefactor which is equal for all of our graphs, our global efficiency is proportional to the average efficiency.

**Local efficiency**: Local efficiency is defined as Eff_{l}(*G*) = *M*^{−1}∑_{i} Eff(*G*_{i}), where *G*_{i} is the subgraph in which we keep the neighbors of node *i* but remove *i* itself, and the efficiency is defined as above.

**Clustering coefficient**: The network average clustering coefficient was calculated by averaging the local clustering coefficient over all nodes. The local clustering coefficient is calculated for a node *i* as the fraction of all possible edges that could exist between the neighbors of *i* which actually do exist in the network. Said another way, the local clustering coefficient is determined by counting all triplets (sets of three connected nodes) centered on node *i* and then calculating the fraction of these triplets that are triangles (all three nodes connected to each other).

**Mean path length**: Given any two nodes, we can calculate the path length as the smallest number of edges that must be traversed to move from one to the other. The average path length is the average over the distribution that arises from calculating this smallest number for all possible pairs of nodes.

### LASSO regression

LASSO (least absolute shrinkage and selection operator) [116] is an advanced regression technique that favors sparse models. This is accomplished by forcing the sum of all the regression coefficients to be less than a certain number; in optimal fitting, some coefficients are then set to zero (i.e., deemed unimportant) rather than assigned a non-zero value. This is usually accomplished by adding the sum of all the regression coefficients, multiplied by a regularization parameter λ, to the quantity being minimized (such as sum of the squared error). When λ is large, having few covariates in a model becomes more important than the model fitting the data accurately. Therefore, the error-minimizing model usually occurs for some intermediate value of λ.

We employed Matlab’s lasso (and lassoPlot) function, which imposes the LASSO constraint on the **L**_{1} norm of all fit coefficients onto a linear regression. In the “deme-level” regression, the outcome variable was the proportion of infections that were with the resistant strain within that deme, and the possible predictors were local properties of that deme. In the “population-level” regression, the outcome variable was the proportion of resistant infections across the entire population, and the possible predictors were properties of the entire population graph. At each level, the regression was conducted both under low resistance level conditions (24% demes treated, leads to < 50% resistant) and high resistance level conditions (40% demes treated, leads to > 50% resistant).

Although a rigorous and robust significance test for the LASSO is still lacking, we can compute the mean-squared error of the predictions of the model to the actual data for cross-validation. Here we choose five-fold validation, meaning that the data is randomly split into five chunks of equal size. Four of these chunks are used to fit the regression, and then the mean-squared error between the actual fifth chunk and its prediction by the regression is measured. If this mean-squared error is sufficiently small, the model is a good fit. The results of this cross-validation procedure for LASSO regression fits to both deme and population-level properties are shown in S5 Fig.

The rank of each property displayed in Table 1 is determined from the magnitude of the regression coefficient at the value of the regularization constraint λ that minimizes the mean-squared error after five-fold cross-validation (S5 and S6 Figs), averaged over both the low and high treatment case.

## Supporting information

### S1 Text. Supplementary methods.

Additional information on model derivation and analysis, equations for models with alternative drug action, and computational methods.

https://doi.org/10.1371/journal.pcbi.1008010.s001

(PDF)

### S1 Fig. Stable spatial heterogeneity in resistance levels in Europe.

Percent of *Klebsiella pneumoniae* isolates resistant to carbapenems in Austria, Belgium, Croatia, France, Germany, Italy, Luxembourg, and Slovenia from 2013-2016 as reported in the ECDC’s European Antibiotic Resistance Surveillance Network [13, 14] and Switzerland as reported by the Swiss Centre for Antibiotic Resistance from 2015-2016. Each country is labeled with the resistance level (%). The year-to-year deviation in Italy from the average value shown in Fig 1 is less than 1% for all years. For all other countries, the frequency of resistance never reaches even 1/10 of the average value seen in Italy.

https://doi.org/10.1371/journal.pcbi.1008010.s002

(TIFF)

### S2 Fig. Spatial heterogeneity in resistance levels in Europe for more bug-drug pairs.

Percent of isolates resistant to particular antibiotics as reported in the ECDC’s European Antibiotic Resistance Surveillance Network [13, 14] in 2017. Each country is labeled with the resistance level (%). Four example bug-drug pairs are shown: *Staphylococcus aureaus* + methicillin, *Pseudomonas aeruginosa* + ceftazadime, *Eschericia coli* + cephalosporins, *Enterococcus faecium* + vancomycin.

https://doi.org/10.1371/journal.pcbi.1008010.s003

(TIF)

### S3 Fig. The distribution of graph-theoretic properties of simulated populations used in LASSO regression analysis.

Populations were created as random networks of 50 demes using a variant of the Watts-Strogatz algorithm. Treatment was randomly assigned to a portion *ρ* of demes so that an overall desired fraction treated was achieved (here, *ρ* = 0.24). Distributions show properties for 1000 such networks. The first six properties are intrinsic to the network structure, while the latter six describe how treatment is allocated over the network. Quantities of the form X-Y give the proportion of all pairs of connected demes in which one deme has treatment status X and one has treatment status Y (U = untreated, T = treated). Quantities of the form X-Y-Z give the same information for triples of connected demes (order ignored). A description of the methods for creating the networks and for calculating each property is given in the **Methods**.

https://doi.org/10.1371/journal.pcbi.1008010.s004

(TIF)

### S4 Fig. The distribution of the fraction of all infections which are drug-resistant in the simulated populations used in the LASSO regression analysis.

Populations were created as random networks of 50 demes using a variant of the Watts-Strogatz algorithm. Treatment was randomly assigned to a portion *ρ* of demes so that an overall desired fraction treated was achieved (on the left *ρ* = 0.24, and on the right *ρ* = 0.4). Infection dynamics were simulated (Eq (3)) until an equilibrium was reached, at which the proportion of all infections with the drug-resistant strain was recorded. Distributions show results for 1000 such simulations each with a unique random network and treatment allocation. Parameters used were *κ* = 0.25/day, *β* = 0.05/day, *g* = 0.1/day, *ϵ* = 0.9, *c* = 0.2.

https://doi.org/10.1371/journal.pcbi.1008010.s005

(TIF)

### S5 Fig. Five-fold cross-validation of LASSO regression models.

Mean-squared error of model predictions vs sub-sample of 1/5 of data after model fitting to 4/5 of data, as a function of the internal parameter penalizing the *L*_{1} norm for fitted coefficients (λ). Red dots show the mean-squared error of the model, with grey bars showing the standard error of the mean. The green circle/line indicates the value of λ that gives the minimum cross-validation error, whereas the blue circle/line locates the λ value where the mean-squared error is one standard error above that of the λ value with minimum cross-validation error.

https://doi.org/10.1371/journal.pcbi.1008010.s006

(TIF)

### S6 Fig. Results of LASSO regression to identify predictors of resistance.

Plots show the regression coefficients obtained as a function of the internal constraint on *L*_{1} norm of coefficients (“regularization parameter”, λ). Each curve is a different predictor variable (described in Table 1). Analysis were separately run to predict deme-level resistance from deme-level properties (left side) or population-level resistance from graph-level properties (right side). Each analysis was conducted for two different levels of drug coverage in the population (either a fraction *ρ* = 0.24 or *ρ* = 0.4 of demes treated). The labels on the curves match the names of the predictors in Table 1.

https://doi.org/10.1371/journal.pcbi.1008010.s007

(TIF)

### S7 Fig. Competition between drug-sensitive and drug-resistant strains in one and two deme populations, when treatment increases recovery rate.

**A)** A population of individuals in a single well-mixed deme, in which a fraction *f* will receive drug treatment when infected (blue haloes). Individuals may be uninfected (hollow blue circles), or infected with either the wild-type (green circles) or drug-resistant (red circles) strain. **B)** The total prevalence of infection (wild-type + drug-resistant) as a function of the fraction of treated individuals (*f*) for different costs of resistance (*c*). **C)** The % of infections that are drug-resistant as a function of the fraction of treated individuals (*f*) for different parameters. Infection switches between 0% and 100% resistant when . Coexistence never occurs. **D)** Schematic of a two-deme population (left-untreated, right-treated) and the two strains (green-wild type, red-resistant) considered in the model. **E-G)** Each panel shows the infection level (shading) as a function of the relative connectivity between demes (*β*/*κ*) and the cost of resistance (*c*). **E)** The % of all infections that are drug-resistant strain across the entire population. **F)** The % of individuals in each deme who are infected with the wild type strain. **G)** The % of individuals in each deme who are infected with the resistant strain. For all results, the transmission rate is *κ* = 0.25/day, the recovery rate is *g* = 0.1/day, and the increase in recovery rate due to treatment is treatment efficacy is *τ* = 0.9/day.

https://doi.org/10.1371/journal.pcbi.1008010.s008

(TIF)

### S8 Fig. Dynamics of drug-resistant infections in populations consisting of networks of inter-connected demes, when treatment increases recovery rate.

**A)** Randomly generated population structures on which infection was simulated. Each node represents a deme (a well-mixed sub-population of individuals), and each edge indicates that infection can spread in either direction between those two demes. Ten example populations were selected out of 1000 total simulated, each with twenty demes randomly connected to three neighbors each, to represent a broad range of outcomes. **B)** Fraction of infections that are resistant in the entire population (*y*-axis) versus fraction of demes treated, *ρ* (*x*-axis). Each color represents a different parameter set (blue background—baseline, red background—lower cost of resistance, teal background—more between-deme connectivity). Numbers show data points for the ten example populations. The colored envelope is created by shading between sigmoidal curves that encompass all the data. **C)** For each population structure shown (*y*-axis) and each treatment level (*x*-axis), the proportion of simulations that resulted in *robust* coexistence between drug-sensitive and drug-resistant strains is shown (by the colored area of the box). Robust coexistence was defined as at least 80% of demes supporting both strains at frequencies above 10%. **D)** Differences in resistance levels (% of all infections that are with the drug resistant-strain) are measured between all pairs of directly-connected untreated demes. **E)** Histograms showing the distribution of pairwise differences in resistance for a given population structure. Lighter shaded histograms combine results from all population graphs. All simulations used kinetic parameters *κ* = 0.25, *g* = 0.1, and *τ* = 0.9, and pooled results from 100 simulations with different random allocation of treatment across demes. Pairwise differences were calculated with 30% treatment.

https://doi.org/10.1371/journal.pcbi.1008010.s009

(TIF)

### S9 Fig. Dynamics of drug-resistant infections in populations consisting of networks of inter-connected demes with heterogeneous connectivity.

**A)** Randomly generated population structures on which infection was simulated. Each node represents a deme (a well-mixed sub-population of individuals), and each edge indicates that infection can spread in either direction between those two demes. Each network had twenty demes and the number of neighbors of each deme was drawn from a gamma distribution with mean three and a different variance. Ten example populations were selected out of 1000 total simulated to represent a broad range of variances in connectivity: from graph 1:4.4, 2:5.7, 3:6.2, 4:6.6, 5:6.9, 6:7.2, 7:7.5, 8:8.0, 9:8.6, 10:10.5. **B)** Fraction of infections that are resistant in the entire population (*y*-axis) versus fraction of demes treated, *ρ* (*x*-axis). Each color represents a different parameter set (blue background—baseline, red background—lower cost of resistance, teal background—more between-deme connectivity). Numbers show data points for the ten example populations. The colored envelope is created by shading between sigmoidal curves that encompass all the data. **C)** For each population structure shown (*y*-axis) and each treatment level (*x*-axis), the proportion of simulations that resulted in *robust* coexistence between drug-sensitive and drug-resistant strains is shown (by the colored area of the box). Robust coexistence was defined as at least 80% of demes supporting both strains at frequencies above 10%. **D)** Differences in resistance levels (% of all infections that are with the drug resistant-strain) are measured between all pairs of directly-connected untreated demes. **E)** Histograms showing the distribution of pairwise differences in resistance for a given population structure. Lighter shaded histograms combine results from all population graphs. All simulations used kinetic parameters *κ* = 0.25/day, *g* = 0.1/day, and *ϵ* = 0.9, and pooled results from 100 simulations with different random allocation of treatment across demes. Pairwise differences were calculated with 30% treatment.

https://doi.org/10.1371/journal.pcbi.1008010.s010

(TIF)

### S10 Fig. Dynamics of drug-resistant infections in populations consisting of networks of inter-connected demes, with varying deme size.

**A)** Randomly generated population structures on which infection was simulated. Each node represents a deme (a well-mixed sub-population of individuals), and each edge indicates that infection can spread in either direction between those two demes. Ten example populations were selected out of 1000 total simulated, each with twenty demes randomly connected to three neighbors each, to represent a broad range of outcomes. Deme sizes were chosen from a normal distribution with coefficient of variation of *σ* = 0, 0.1, of 0.3. **B)** Fraction of infections that are resistant in the entire population (*y*-axis) versus fraction of demes treated, *ρ* (*x*-axis). Each color represents a different level of heterogeneity in deme size (n navy blue—uniform sizes (*σ* = 0), medium blue—*σ* = 0.1, light blue—*σ* = 0.3). Numbers show data points for the ten example populations. The colored envelope is created by shading between sigmoidal curves that encompass all the data. **C)** For each population structure shown (*y*-axis) and each treatment level (*x*-axis), the proportion of simulations that resulted in *robust* coexistence between drug-sensitive and drug-resistant strains is shown (by the colored area of the box). Robust coexistence was defined as at least 80% of demes supporting both strains at frequencies above 10%. **D)** Differences in resistance levels (% of all infections that are with the drug resistant-strain) are measured between all pairs of directly-connected untreated demes. **E)** Histograms showing the distribution of pairwise differences in resistance for a given population structure. Lighter shaded histograms combine results from all population graphs. All simulations used kinetic parameters *κ* = 0.25/day, *β* = 0.05/day, *g* = 0.1/day, *c* = 0.2, and *ϵ* = 0.9, and pooled results from 100 simulations with different random allocation of deme sizes and treatment across demes. Pairwise differences were calculated with 30% treatment.

https://doi.org/10.1371/journal.pcbi.1008010.s011

(TIF)

### S11 Fig. Dynamics of drug-resistant infections in populations consisting of networks of inter-connected demes, with varying levels of background mixing between demes.

**A)** Randomly generated population structures on which infection was simulated. Each node represents a deme (a well-mixed sub-population of individuals), and each edge indicates that infection can spread in either direction between those two demes. Each population consisted of twenty demes randomly connected to three neighbors each, and additionally connected with relative weight *α* to all other demes. Ten example populations were selected out of 1000 total simulated to represent a broad range of outcomes. **B)** Fraction of infections that are resistant in the entire population (*y*-axis) versus fraction of demes treated, *ρ* (*x*-axis). Each color represents a different level of background connectivity between all demes (navy blue—no background mixing (*α* = 0), dark blue—*α* = 0.01, medium blue—*α* = 0.1, light blue—*α* = 1). Numbers show data points for the ten example populations. The colored envelope is created by shading between sigmoidal curves that encompass all the data. **C)** For each population structure shown (*y*-axis) and each treatment level (*x*-axis), the proportion of simulations that resulted in *robust* coexistence between drug-sensitive and drug-resistant strains is shown (by the colored area of the box). Robust coexistence was defined as at least 80% of demes supporting both strains at frequencies above 10%. **D)** Differences in resistance levels (% of all infections that are with the drug resistant-strain) are measured between all pairs of directly-connected untreated demes. **E)** Histograms showing the distribution of pairwise differences in resistance for a given population structure. Lighter shaded histograms combine results from all population graphs. All simulations used kinetic parameters *κ* = 0.25/day, *β* = 0.05/day, *g* = 0.1/day, *c* = 0.2, and *ϵ* = 0.9, and pooled results from 100 simulations with different random allocation of treatment across demes. Pairwise differences were calculated with 30% treatment.

https://doi.org/10.1371/journal.pcbi.1008010.s012

(TIF)

### S12 Fig. Alternate distributions of treatment across demes.

Each panel shows a distribution across demes of the fraction of individuals within the deme receiving treatment if infected. Treatment levels were drawn from Beta distributions with varying mean (*ρ*, rows) and shape parameter (*α*, columns). The left most column corresponds to the results presented in the rest of the paper, where individuals in each deme either have a 0 or 100% chance of being treated if infected.

https://doi.org/10.1371/journal.pcbi.1008010.s013

(TIF)

### S13 Fig. Frequency of drug resistant infections for alternate distributions of treatment across demes.

The percent of infections over the whole population that are caused by the the drug resistant strain as a function of the average treatment level across demes (*ρ*, y-axis) and the Beta-distribution parameter describing the spread of treatment levels (*α*, x-axis). Dotted lines border the region where resistance levels fall between 5% and 95%. Each panel shows a distribution across demes of the fraction of individuals within the deme receiving treatment if infected. Treatment levels were drawn from Beta distributions with varying mean (*ρ*, rows) and shape parameter (*α*, columns). Example treatment distributions are shown in S12 Fig. The left most edge of the graph (*α* = 0) corresponds to the results presented in the rest of the paper, where individuals in each deme either have a 0 or 100% chance of being treated if infected. Higher *α* values correspond to more continuous, unimodal treatment distributions with lower variance. The population consisted of twenty demes randomly connected to three neighbors each. Results were averaged over 1000 graphs with with 1000 random treatment allocations for each. For all results, the transmission rates are *κ* = 0.25/day (intra-deme) and *β* = 0.05/day (inter-deme), the recovery rate is *g* = 0.1/day, the cost of resistance is *c* = 0.1, and the treatment efficacy is *ϵ* = 0.9.

https://doi.org/10.1371/journal.pcbi.1008010.s014

(TIF)

### S14 Fig. Frequency of drug resistant infections for alternate distributions of treatment across demes.

Each panel shows the average fraction of infections that are resistant in the entire population as a function of the relative connectivity between demes (*β*/*κ*) and the cost of resistance (*c*). Dotted lines border the region where resistance levels fall between 5% and 95%. The mean and variance in the amount of treatment per deme varied between panels and is given by the distributions in S12 Fig. The left most column corresponds to the case where individuals in each deme either have a 0 or 100% chance of being treated if infected. Higher *α* values correspond to more continuous, unimodal treatment distributions with lower variance. The population consisted of twenty demes randomly connected to three neighbors each. Results were averaged over 1000 graphs with with 1000 random treatment allocations for each. For all results, the intra-deme transmission rate is *κ* = 0.25/day, the recovery rate is *g* = 0.1/day, and the treatment efficacy is *ϵ* = 0.9.

https://doi.org/10.1371/journal.pcbi.1008010.s015

(TIF)

### S15 Fig. Distributed coexistence for alternate distributions of treatment across demes.

Each panel shows the average fraction of demes that supported both strains at frequencies of at least 10%, as a function of the relative connectivity between demes (*β*/*κ*) and the cost of resistance (*c*). The mean and variance in the amount of treatment per deme varied between panels and is given by the distributions in S12 Fig. The left most column corresponds to the case where individuals in each deme either have a 0 or 100% chance of being treated if infected. The population consisted of twenty demes randomly connected to three neighbors each. Results were averaged over 1000 graphs with with 1000 random treatment allocations for each. For all results, the intra-deme transmission rate is *κ* = 0.25/day, the recovery rate is *g* = 0.1/day, and the treatment efficacy is *ϵ* = 0.9.

https://doi.org/10.1371/journal.pcbi.1008010.s016

(TIF)

### S16 Fig. Sensitivity of results to the simulation stopping parameter.

We tested the sensitivity of our results to the value of the derivative which we considered ``at equilibrium'', defined as the point where the sum of the derivatives of all strains in each deme was less than the stopping parameter divided by the number of demes. **(Top row)** We considered 10^{−7}, the smallest value tested, as the ``true'' solution and measured this solution against other solutions with all parameters held equal, but the stopping parameter increased. Box plots show two representatives of the error from this solution as the stopping parameter is increased up to our parameter 10^{−4} used in the main-text results and beyond: **(left)**, the matrix 2-norm (also known as **L**_{2} norm, Frobenius norm, or Euclidean distance). **(Right)**: the maximum demewise error (**L**_{∞} norm). **(Bottom row:** For these same simulations, we show: **(left)**, the time in internal units (simulated days) until equilibrium was reached; **(right)**, the physical computation time of the simulation. The measures of error for our chosen parameter are small given the advancement in computation time afforded by that choice.

https://doi.org/10.1371/journal.pcbi.1008010.s017

(TIF)

### S17 Fig. Difference between equilibrium values obtained for infection dynamics in multi-deme systems as initial conditions vary.

For 10 different random regular graphs with 20 demes and degree (*d*) 3, 4, or 5, treatment was allocated randomly across demes with an overall proportion treated *ρ* = 0.35. For each graph-treatment allocation combination, 100 different sets of initial conditions (levels of resistant and sensitive infections in each deme) were chosen uniformly at randomly, and infection dynamics were run until an equilibrium was reached as described in the Methods. The maximum difference (in matrix 2-norm, also known as **L**_{2} norm, Frobenius norm, or Euclidean distance) was calculated between the equilibria seen in any of 100 trials and the equilibrium values used for results reported in the main text. Very low maximum norm values suggest there is a single stable equilibrium value for each parameter set.

https://doi.org/10.1371/journal.pcbi.1008010.s018

(TIF)

## Acknowledgments

We thank Jeff Gerold, Sam Sinai, Kamran Kaveh, Anjalika Nande, Gabriel Leventhal, and participants in the 2018 Gordon Research Conference on Drug Resistance for helpful discussions. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.

## References

- 1. Roberts RR, Hota B, Ahmad I, Scott RD II, Foster SD, Abbasi F, et al. Hospital and Societal Costs of Antimicrobial-Resistant Infections in a Chicago Teaching Hospital: Implications for Antibiotic Stewardship. Clinical Infectious Diseases. 2009;49(8):1175–1184. pmid:19739972
- 2. McDonald LC. Trends in Antimicrobial Resistance in Health Care: Associated Pathogens and Effect on Treatment. Clinical Infectious Diseases. 2006;42(Supplement):S65–S71.
- 3. Unemo M, Nicholas RA. Emergence of multidrug-resistant, extensively drug-resistant and untreatable gonorrhea. Future Microbiology. 2012;7(12):1401–1422.
- 4.
Neuner EA, Yeh JY, Hall GS, Sekeres J, Endimiani A, Bonomo RA, et al. Treatment and outcomes in carbapenem-resistant
*Klebsiella pneumoniae*bloodstream infections. Diagnostic Microbiology and Infectious Disease. 2011;69(4):357–362. https://doi.org/10.1016/j.diagmicrobio.2010.10.013 pmid:21396529 - 5. Perez F, Hujer AM, Hujer KM, Decker BK, Rather PN, Bonomo RA. Global Challenge of Multidrug-Resistant Acinetobacter baumannii. Antimicrobial Agents and Chemotherapy. 2007;51(10):3471–3484.
- 6. Reardon S. WHO warns against ‘post-antibiotic’ era. Nature News. 2014.
- 7.
Molteni M. The Post-Antibiotic Era Is Here. Now What? Wired. 2017;.
- 8.
Anderson R, May R. Infectious Diseases of Humans—Dynamics and Control. Oxford, UK: OUP; 1992.
- 9. Opatowski L, Guillemot D, Boelle PY, Temime L. Contribution of mathematical modeling to the fight against bacterial antibiotic resistance. Current Opinion in Infectious Diseases. 2011;24(3):279.
- 10.
O’Neill J. The Review on Antimicrobial Resistance: Final Report and Recomendations; 2016. Available from: https://amr-review.org/sites/default/files/160525_Final%20paper_with%20cover.pdf.
- 11. Kraker MEAd, Stewardson AJ, Harbarth S. Will 10 Million People Die a Year due to Antimicrobial Resistance by 2050? PLOS Medicine. 2016;13(11):e1002184.
- 12.
Liñares J, Ardanuy C, Pallares R, Fenoll A. Changes in antimicrobial resistance, serotypes and genotypes in
*Streptococcus pneumoniae*over a 30-year period. Clin Microbiol Infect. 2010;16:402–410. - 13.
European Centre for Disease Prevention and Control. Surveillance of antimicrobial resistance in Europe. Stockholm: ECDC; 2018. Available from: https://ecdc.europa.eu/sites/portal/files/documents/EARS-Net-report-2017-update-jan-2019.pdf.
- 14.
European Centre for Disease Prevention and Control. European Centre for Disease Prevention and Control;. Available from: https://ecdc.europa.eu/en/antimicrobial-resistance/surveillance-and-disease-data/data-ecdc.
- 15.
European Centre for Disease Prevention and Control. Antimicrobial consumption. Stockholm: ECDC; 2018. Available from: https://www.ecdc.europa.eu/sites/portal/files/documents/AER_for_2017-antimicrobial-consumption.pdf.
- 16.
European Centre for Disease Prevention and Control. European Centre for Disease Prevention and Control;. Available from: https://ecdc.europa.eu/en/antimicrobial-consumption/database/trend-country.
- 17.
Swiss Centre for Antibiotic resistance. Interactive Database—Antibiotic resistance;. Available from: http://www.anresis.ch/index.php/Interactive-database.html.
- 18.
Fenoll A, Bourgon CM, Munoz R, Vicioso D, Casal J. Serotype Distribution and Antimicrobial Resistance of
*Streptococcus pneumoniae*Isolates Causing Systemic Infections in Spain, 1979-1989. Reviews of Infectious Diseases. 1991;13(1):56–60. - 19.
Fenoll A, Jado I, Vicioso D, Perez A, Casal J. Evolution of
*Streptococcus pneumoniae*Serotypes and Antibiotic Resistance in Spain: Update (1990 to 1996). Journal of Clinical Microbiology. 1998;36(12):3447–3454. - 20.
Bergman M, Nyberg ST, Huovinen P, Paakkari P, Hakanen AJ. Association between Antimicrobial Consumption and Resistance in
*Escherichia coli*. Antimicrobial Agents and Chemotherapy. 2009;53(3):912–917. - 21.
Song JH, Jung SI, Ko KS, Kim NY, Son JS, Chang HH, et al. High Prevalence of Antimicrobial Resistance among Clinical
*Streptococcus pneumoniae*Isolates in Asia (an ANSORP Study). Antimicrobial Agents and Chemotherapy. 2004;48(6):2101–2107. pmid:15155207 - 22.
Camargos P, Fischer GB, Mocelin H, Dias C, Ruvinsky R. Penicillin resistance and serotyping of
*Streptococcus pneumoniae*in Latin America. Paediatric Respiratory Reviews. 2006;7(3):209—214. - 23. Andersson DI. The biological cost of mutational antibiotic resistance: any practical conclusions? Curr Opin Microbiol. 2006;9:461–465.
- 24. Andersson DI, Hughes D. Antibiotic resistance and its cost: is it possible to reverse resistance? Nat Rev Microbiol. 2010;8:260–271.
- 25. Melnyk AH, Wong A, Kassen R. The fitness costs of antibiotic resistance mutations. Evolutionary Applications. 2014;8(3):273–283.
- 26. Spicknall IH, Foxman B, Marrs CF, Eisenberg JNS. A Modeling Framework for the Evolution and Spread of Antibiotic Resistance: Literature Review and Model Categorization. American Journal of Epidemiology. 2013;178(4):508–520.
- 27. Gauze GF. Experimental Studies on the Struggle for Existence. Journal of Experimental Biology. 1932;9(4):389–402.
- 28. Hardin G. The competitive exclusion principle. Science. 1960;131:1292–1297.
- 29. Hutchinson GE. The paradox of the plankton. Amer Nat. 1961;95:137–145.
- 30. Colijn C, Cohen T, Fraser C, Hanage W, Goldstein E, Givon-Lavi N, et al. What is the mechanism for persistent coexistence of drug-susceptible and drug-resistant strains of Streptococcus pneumoniae? Journal of The Royal Society Interface. 2009; p. rsif20090400.
- 31. Cobey S, Baskerville EB, Colijn C, Hanage W, Fraser C, Lipsitch M. Host population structure and treatment frequency maintain balancing selection on drug resistance. Journal of The Royal Society Interface. 2017;14(133):20170295.
- 32. Davies NG, Flasche S, Jit M, Atkins KE. Within-host dynamics shape antibiotic resistance in commensal bacteria. Nature Ecology & Evolution. 2019;3(3):440.
- 33. Blanquart F, Lehtinen S, Lipsitch M, Fraser C. The evolution of antibiotic resistance in a structured host population. Journal of The Royal Society Interface. 2018;15(143):20180040.
- 34. Chesson P. Mechanisms of Maintenance of Species Diversity. Annu Rev Ecol Syst. 2000;31:343–366.
- 35. Cobey S, Lipsitch M. Niche and neutral effects of acquired immunity permit coexistence of pneumococcal serotypes. Science. 2012;335:1376–1380.
- 36. Lehtinen S, Blanquart F, Croucher NJ, Turner P, Lipsitch M, Fraser C. Evolution of antibiotic resistance is linked to any genetic mechanism affecting bacterial duration of carriage. Proceedings of the National Academy of Sciences. 2017;114(5):1075–1080.
- 37. Goossens H, Ferech M, Vander Stichele R, Elseviers M, ESAC Project Group. Outpatient antibiotic use in Europe and association with resistance: a cross-national database study. Lancet. 2005;365(9459):579–587.
- 38. Bell BG, Schellevis F, Stobberingh E, Goossens H, Pringle M. A systematic review and meta-analysis of the effects of antibiotic consumption on antibiotic resistance. BMC Infectious Diseases. 2014;14(1):13.
- 39. MacFadden DR, Fisman D, Andre J, Ara Y, Majumder MS, Bogoch II, et al. A Platform for Monitoring Regional Antimicrobial Resistance, Using Online Data Sources: ResistanceOpen. The Journal of Infectious Diseases. 2016;214(Supplement 4):S393–S398. pmid:28830108
- 40. Grundmann H, Aanensen DM, Wijngaard CCvd, Spratt BG, Harmsen D, Friedrich AW, et al. Geographic Distribution of Staphylococcus aureus Causing Invasive Infections in Europe: A Molecular-Epidemiological Analysis. PLOS Medicine. 2010;7(1):e1000215. pmid:20084094
- 41. Yahiaoui RY, Bootsma HJ, den Heijer CDJ, Pluister GN, John Paget W, Spreeuwenberg P, et al. Distribution of serotypes and patterns of antimicrobial resistance among commensal Streptococcus pneumoniae in nine European countries. BMC Infectious Diseases. 2018;18(1):440. pmid:30157780
- 42. Whitney CG, Farley MM, Hadler J, Harrison LH, Lexau C, Reingold A, et al. Increasing Prevalence of Multidrug-Resistant Streptococcus pneumoniae in the United States. New England Journal of Medicine. 2000;343(26):1917–1924. pmid:11136262
- 43.
McCormick A, et al. Geographic diversity and temporal trends of antimicrobial resistance in
*Streptococcus pneumoniae*in the United States. Nature Medicine. 2003;9(4):424–431. pmid:12627227 - 44. Fenoll A, Asensio G, Jado I, Berrón S, Camacho MT, Ortega M, et al. Antimicrobial susceptibility and pneumococcal serotypes. Journal of Antimicrobial Chemotherapy. 2002;50(Suppl. S2):1–8.
- 45. Nygaard Jensen J, Melander E, Hedin K, Bjerrum L, Kornfalt Isberg H, Holm A, et al. Comparison of antibiotic prescribing and antimicrobial resistance in urinary tract infections at the municipal level among women in two Nordic regions. Journal of Antimicrobial Chemotherapy. 2018;73(8):2207–2214. pmid:29757408
- 46. Kimura M, Weiss GH. The stepping stone model of population structure and the decrease of genetic correlation with distance. Genetics. 1964;49:561–576.
- 47. Ball F. Stochastic and deterministic models for SIS epidemics among a population partitioned into households. Math Biosci. 1999;156:41–67.
- 48. Tedijanto C, Olesen SW, Grad YH, Lipsitch M. Estimating the proportion of bystander selection for antibiotic resistance among potentially pathogenic bacterial flora. Proceedings of the National Academy of Sciences of the United States of America. 2018;115(51):E11988–E11995.
- 49. Morley VJ, Woods RJ, Read AF. Bystander Selection for Antimicrobial Resistance: Implications for Patient Health. Trends in Microbiology. 2019;27(10):864–877.
- 50.
Keeling MJ. Modeling infectious diseases in humans and animals. Princeton: Princeton University Press; 2008.
- 51. Humplik J, Hill AL, Nowak MA. Evolutionary dynamics of infectious diseases in finite populations. Journal of Theoretical Biology. 2014;360:149–162.
- 52. Lipsitch M, Colijn C, Cohen T, Hanage WP, Fraser C. No coexistence for free: neutral null models for multistrain pathogens. Epidemics. 2009;1:2–13.
- 53.
Bollobás B. Random Graphs. Cambridge, UK: Cambridge University Press; 2001.
- 54. Watts DJ, Strogatz SH. Collective dynamics of small-world networks. Nature. 1998;393:440–442.
- 55. Hay SI, Snow RW. The Malaria Atlas Project: Developing Global Maps of Malaria Risk. PLOS Medicine. 2006;3(12):e473.
- 56. Bhatt S, Gething PW, Brady OJ, Messina JP, Farlow AW, Moyes CL, et al. The global distribution and burden of dengue. Nature. 2013;496(7446):504–507. pmid:23563266
- 57. Messina JP, Kraemer MU, Brady OJ, Pigott DM, Shearer FM, Weiss DJ, et al. Mapping global environmental suitability for Zika virus. eLife. 2016;5:e15272. pmid:27090089
- 58. Kraemer MUG, Hay SI, Pigott DM, Smith DL, Wint GRW, Golding N. Progress and Challenges in Infectious Disease Cartography. Trends in Parasitology. 2016;32(1):19–29.
- 59.
Tatem AJ, Smith DL. International population movements and regional
*Plasmodium falciparum*malaria elimination strategies. Proceedings of the National Academy of Sciences. 2010;107(27):12222–12227. - 60. Wesolowski A, Eagle N, Tatem AJ, Smith DL, Noor AM, Snow RW, et al. Quantifying the impact of human mobility on malaria. Science. 2012;338(6104):267–270. pmid:23066082
- 61.
Carrasco-Escobar G, Gamboa D, Castro MC, Bangdiwala SI, Rodriguez H, Contreras-Mancilla J, et al. Micro-epidemiology and spatial heterogeneity of
*P. vivax*parasitaemia in riverine communities of the Peruvian Amazon: A multilevel analysis. Scientific Reports. 2017;7(1):8082. pmid:28808240 - 62.
Kang SY, Battle KE, Gibson HS, Cooper LV, Maxwell K, Kamya M, et al. Heterogeneous exposure and hotspots for malaria vectors at three study sites in Uganda. bioRxiv. 2018;.
- 63.
CDC. Centers for Disease Control and Prevention;. Available from: https://arpsp.cdc.gov/.
- 64. Livermore DM, Pearson A. Antibiotic resistance: location, location, location. Clinical Microbiology and Infection. 2007;13:7–16.
- 65. Goyal MK, Johnson TJ, Chamberlain JM, Casper TC, Simmons T, Alessandrini EA, et al. Racial and Ethnic Differences in Antibiotic Use for Viral Illness in Emergency Departments. Pediatrics. 2017;140(4).
- 66. Olesen SW, Grad YH. Racial/Ethnic Disparities in Antimicrobial Drug Use, United States, 2014-2015. Emerging Infectious Diseases. 2018;24(11). pmid:30334733
- 67. Olesen SW, Barnett ML, MacFadden DR, Brownstein JS, Hernandez-Diaz S, Lipsitch M, et al. The distribution of antibiotic use and its association with antibiotic resistance. eLife. 2018;7:e39435. pmid:30560781
- 68. Volpi C, Shehadeh F, Mylonakis E. Correlation of antimicrobial prescription rate and county income in medicare part D. Medicine. 2019;98(22):e15914.
- 69. Salathe M, Kazandjieva M, Lee JW, Levis P, Feldman MW, Jones JH. A high-resolution human contact network for infectious disease transmission. Proceedings of the National Academy of Sciences. 2010; p. 201009094.
- 70. Vanhems P, Barrat A, Cattuto C, Pinton JF, Khanafer N, Regis C, et al. Estimating Potential Infection Transmission Routes in Hospital Wards Using Wearable Proximity Sensors. PLoS ONE. 2013;8(9):e73970. pmid:24040129
- 71. Voirin N, Payet C, Barrat A, Cattuto C, Khanafer N, Regis C, et al. Combining High-Resolution Contact Data with Virological Data to Investigate Influenza Transmission in a Tertiary Care Hospital. Infection Control & Hospital Epidemiology. 2015;36(03):254–260.
- 72. Wesolowski A, Buckee CO, Pindolia DK, Eagle N, Smith DL, Garcia AJ, et al. The Use of Census Migration Data to Approximate Human Movement Patterns across Temporal Scales. PLOS ONE. 2013;8(1):e52971. pmid:23326367
- 73. Tizzoni M, Bajardi P, Decuyper A, King GKK, Schneider CM, Blondel V, et al. On the Use of Human Mobility Proxies for Modeling Epidemics. PLOS Computational Biology. 2014;10(7):e1003716. pmid:25010676
- 74. Wesolowski A, O’Meara WP, Eagle N, Tatem AJ, Buckee CO. Evaluating Spatial Interaction Models for Regional Mobility in Sub-Saharan Africa. PLOS Computational Biology. 2015;11(7):e1004267.
- 75. Wesolowski A, Buckee CO, Engo-Monsen K, Metcalf CJE. Connecting Mobility to Infectious Diseases: The Promise and Limits of Mobile Phone Data. The Journal of Infectious Diseases. 2016;214(suppl_4):S414–S420.
- 76. Brockmann D, Helbing D. The Hidden Geometry of Complex, Network-Driven Contagion Phenomena. Science. 2013;342(6164):1337–1342.
- 77. Nunes MRT, Palacios G, Faria NR, CS E Jr, Pantoja JA, Rodrigues SG, et al. Air Travel Is Associated with Intracontinental Spread of Dengue Virus Serotypes 1–3 in Brazil. PLOS Neglected Tropical Diseases. 2014;8(4):e2769. pmid:24743730
- 78. Pybus OG, Suchard MA, Lemey P, Bernardin FJ, Rambaut A, Crawford FW, et al. Unifying the spatial epidemiology and molecular evolution of emerging epidemics. Proceedings of the National Academy of Sciences. 2012;109(37):15066–15071.
- 79. Snitkin ES, Zelazny AM, Thomas PJ, Stock F, Henderson DK, Palmore TN, et al. Tracking a Hospital Outbreak of Carbapenem-Resistant Klebsiella pneumoniae with Whole-Genome Sequencing. Science Translational Medicine. 2012;4(148):148ra116–148ra116. pmid:22914622
- 80. Frost SDW, Pybus OG, Gog JR, Viboud C, Bonhoeffer S, Bedford T. Eight challenges in phylodynamic inference. Epidemics. 2015;10:88–92.
- 81.
Taylor AR, Schaffner SF, Cerqueira GC, Nkhoma SC, Anderson TJC, Sriprawat K, et al. Quantifying connectivity between local
*Plasmodium falciparum*malaria parasite populations using identity by descent. PLOS Genetics. 2017;13(10):e1007065. pmid:29077712 - 82. Dudas G, Carvalho LM, Bedford T, Tatem AJ, Baele G, Faria NR, et al. Virus genomes reveal factors that spread and sustained the Ebola epidemic. Nature. 2017;544(7650):309–315. pmid:28405027
- 83. Kepler TB, Perelson AS. Drug concentration heterogeneity facilitates the evolution of drug resistance. Proceedings of the National Academy of Sciences. 1998;95(20):11514–11519.
- 84. Greulich P, Waclaw B, Allen RJ. Mutational Pathway Determines Whether Drug Gradients Accelerate Evolution of Drug-Resistant Cells. Physical Review Letters. 2012;109(8):088101.
- 85. Hermsen R, Deris JB, Hwa T. On the rapidity of antibiotic resistance evolution facilitated by a concentration gradient. Proceedings of the National Academy of Sciences. 2012;109(27):10775–10780.
- 86. Moreno-Gamez S, Hill AL, Rosenbloom DIS, Petrov DA, Nowak MA, Pennings PS. Imperfect drug penetration leads to spatial monotherapy and rapid evolution of multidrug resistance. Proceedings of the National Academy of Sciences. 2015;112(22):E2874–E2883.
- 87. Fu F, Nowak MA, Bonhoeffer S. Spatial Heterogeneity in Drug Concentrations Can Facilitate the Emergence of Resistance to Cancer Therapy. PLOS Computational Biology. 2015;11(3):1–22.
- 88. Leventhal GE, Hill AL, Nowak MA, Bonhoeffer S. Evolution and emergence of infectious diseases in theoretical and real-world networks. Nature Communications. 2015;6.
- 89. De Jong MG, Wood KB. Tuning Spatial Profiles of Selection Pressure to Modulate the Evolution of Drug Resistance. Physical Review Letters. 2018;120(23):238102.
- 90. Liechti JI, Leventhal GE, Bonhoeffer S. Host population structure impedes reversion to drug sensitivity after discontinuation of treatment. PLOS Computational Biology. 2017;13(8):e1005704.
- 91. Cooper BS, Medley GF, Stone SP, Kibbler CC, Cookson BD, Roberts JA, et al. Methicillin-resistant Staphylococcus aureus in hospitals and the community: stealth dynamics and control catastrophes. Proceedings of the National Academy of Sciences of the United States of America. 2004;101(27):10223–10228. pmid:15220470
- 92. Hanski I, Ovaskainen O. Metapopulation theory for fragmented landscapes. Theoretical Population Biology. 2003;64(1):119–127.
- 93.
Lopes F, Luczak M. Extinction time for the weaker of two competing SIS epidemics. arXiv:180204037 [math]. 2018;.
- 94.
WHO. HIV Drug Resistance Report 2017. Geneva, Switzerland: World Health Organization; 2017. Available from: http://www.who.int/hiv/pub/drugresistance/hivdr-report-2017/en/.
- 95. Chimukangara B, Lessells RJ, Rhee SY, Giandhari J, Kharsany ABM, Naidoo K, et al. Trends in Pretreatment HIV-1 Drug Resistance in Antiretroviral Therapy-naive Adults in South Africa, 2000–2016: A Pooled Sequence Analysis. EClinicalMedicine. 2019;9:26–34. pmid:31143879
- 96. Gupta RK, Gregson J, Parkin N, Haile-Selassie H, Tanuri A, Forero LA, et al. HIV-1 drug resistance before initiation or re-initiation of first-line antiretroviral therapy in low-income and middle-income countries: a systematic review and meta-regression analysis. The Lancet Infectious Diseases. 2018;18(3):346–355. pmid:29198909
- 97. Olson A, Bannert N, Sonnerborg A, de Mendoza C, Price M, Zangerle R, et al. Temporal trends of transmitted HIV drug resistance in a multinational seroconversion cohort. AIDS. 2018;32(2):161. pmid:29112061
- 98. Talisuna AO, Bloland P, D’Alessandro U. History, Dynamics, and Public Health Importance of Malaria Parasite Resistance. Clinical Microbiology Reviews. 2004;17(1):235–254.
- 99.
Okombo J, Kamau AW, Marsh K, Sutherland CJ, Ochola-Oyier LI. Temporal trends in prevalence of
*Plasmodium falciparum*drug resistance alleles over two decades of changing antimalarial policy in coastal Kenya. International Journal for Parasitology: Drugs and Drug Resistance. 2014;4(3):152–163. - 100. Bushman M, Antia R, Udhayakumar V, Roode JCd. Within-host competition can delay evolution of drug resistance in malaria. PLOS Biology. 2018;16(8):e2005712.
- 101. Fisher RA. The Wave of Advance of Advantageous Genes. Annals of Eugenics. 1937;7(4):355–369.
- 102. Kimura M, Weiss GH. The Stepping Stone Model of Population STructure and the Decrease of Genetic Correlation with Distance. Genetics. 1964;49(4):561–576.
- 103. Levins R. Theory of Fitness in a Heterogeneous Environment. I. The Fitness Set and Adaptive Function. The American Naturalist. 1962;96(891):361–373.
- 104. Holt RD. Spatial Heterogeneity, Indirect Interactions, and the Coexistence of Prey Species. The American Naturalist. 1984;124(3):377–406.
- 105. Gavrilets S, Gibson N. Fixation probabilities in a spatially heterogeneous environment. Popul Ecol. 2002;44:51–58.
- 106. Whitlock MC, Gomulkiewicz R. Probability of Fixation in a Heterogeneous Environment. Genetics. 2005;171:1407–1417.
- 107. Lieberman E, Hauert C, Nowak MA. Evolutionary dynamics on graphs. Nature. 2005;433(7023):312–316.
- 108. Manem VSK, Kaveh K, Kohandel M, Sivaloganathan S. Modeling Invasion Dynamics with Spatial Random-Fitness due to Micro-Environment. PLoS One. 2015;10:e0140234.
- 109.
Kaveh K, McAvoy A, Nowak MA. The Effect of Spatial Fitness Heterogeneity on Fixation Probability. arXiv. 2017; p. 1709.03031.
- 110. Kouyos R, Klein E, Grenfell B. Hospital-Community Interactions Foster Coexistence between Methicillin-Resistant Strains of Staphylococcus aureus. PLOS Pathogens. 2013;9(2):1–14.
- 111. Heffernan JM, Smith RJ, Wahl LM. Perspectives on the basic reproductive ratio. J R Soc Interface. 2005;2:281–293.
- 112. Li MY, Shuai Z. Global-stability problem for coupled systems of differential equations on networks. Journal of Differential Equations. 2010;248(1):1–20.
- 113. Shuai Z, van den Driessche P. Global Stability of Infectious Disease Models Using Lyapunov Functions. SIAM Journal on Applied Mathematics. 2013;73(4):1513–1532.
- 114. Kim JH, Vu VH. Generating Random Regular Graphs. Combinatorica. 2006;26(6):683–708.
- 115.
Pundak G. Random Regular Graph Generator. MATLAB Central File Exchange. 2010;.
- 116. Tibshirani R. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society Series B (Methodological). 1996;58(1):267–288.