Population structure across scales facilitates coexistence and spatial heterogeneity of antibiotic-resistant infections

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.

The transmission rate k represents the rate at which each infected individual contacts and infects each susceptible individual. Although we use the terms S t and S u to track susceptible individuals who are "treated" or "untreated", our model is agnostic as to whether these uninfected individuals actually receive drug, or if these markers simply indicate if these individuals will (or will not) receive treatment when infected. While here we assume a pre-determined fraction f of the total population will receive treatment if infected, similar models could be constructed assuming the decision to receive treatment is made instantaneously upon infection with probability f , or that there is some rate of receiving treatment, such that a fraction f will receive treatment before recovery. The analysis of these models gives identical steady state results. Because in our model the total population size N , as well as the fraction of individuals who will (or will not) receive treatment if infected, are both constant, we can remove the variables for susceptible individuals using S u = N (1 − f ) − W u − R u and S t = N f − W t − R t . To further reduce the variables in this system and remove any dependence on total population size, we let κ = kN and x i = X i /N (the fraction of individuals instead of the total number, for X = W, R and i = u, t), which then leads to Eq. 1. By keeping κ = kN as the constant parameter, we are implicitly assuming that the total rate at which a single individual can contact and infect others when the population is completely susceptible (kN ) is independent of the population size. This is a reasonable assumption for many infections which require some type of active contact, and here allows us to reduce the number of parameters required. None of our results depend on this scaling, and the effect of changing population size is equivalent to changing κ in this deterministic model.
There are three possible equilibria of Eq. 1. In each of these equilibria, there is no possible coexistence of infections. These equilibria are as follows: This uninfected equilibrium is only stable when 2. Only wild-type (drug-sensitive) strain survives This infected equilibrium is stable when where the later condition implies c > f .

Only drug-resistant strain survives
This infected equilibrium is stable when where the former condition implies c < f .

Analytic results for two connected subpopulations
In an effort to be more general than in the main text, we consider the case of a population of total size N divided into two demes, with a fraction f of the total population is contained in the drug-treated deme. The system describing the number of individuals in each state can then be written as: Since the size of each deme is fixed, we can remove the variables for susceptible individuals using To create a nondimensional version of the system, we instead track the fraction of individuals in each deme who have a particular infection status, so , and κ = kN and β = bN , which gives the altered system:ẇ Using the next-generation technique [?], we arrive at the following values for the basic reproductive ratios, In general, this system has nine equilibria. However, only one of these is stable and physically realizable (e.g. having between 0% and 100% of individuals infected) for any given parameter set. These are the values that are plotted in the main text. These outcomes fall into one of two qualitative categories. Either one strain drives the other to competitive extinction (true in Fig 3 for c > 0.4 and c < 0.06, for instance) or both strains exist in both demes simultaneously. The boundaries separating this region from the competitiveexclusion region are R r 0 = R w 0 , the lower bound when viewed as a function of the cost of resistance c, and the curve which solves the equilibrium values. This second curve, in words, is where the drug-treated population is infected at the same rate by both drug-resistant and drug-sensitive strains.
For the parameters chosen here and used throughout the paper, it is impossible to have R w 0 < 1. For sufficiently high c, we can have R r 0 < 1. However, this line in parameter space lies deep within the region where the population is infected with the wild-type strain only.

Derivation of a multi-deme model
Here we derive the general case for the spread of wild-type and resistant strains of infection in a population consisting of multiple connected demes (with connectivity described by adjacency matrix ∆), where demes may have different sizes (D i ) and different treatment levels (f i ). The dynamics of infection are described by the following system of equations where X it represents the total number of individuals in deme i infected with strain X (or susceptible, S) who will be treated if infected, whereas X iu is those who won't be treated if infected.
Since the population size D i of each deme is constant and we assume it is predetermined who will get treatment if infected, we can remove the equations for the susceptible individuals using To derive a dimensionless form the the equations, we let x i = X i /D i , and define σ i = D i /D whereD is the average deme size. We define κ = kD and β = bD, to derive the simplified and scaled equationṡ In the limit of either 0% or 100% treatment in a deme, this system reduces tȯ

(S.13)
In the limit of both constant deme size and 0% or 100% treatment in a deme, this system reduces to Eq. 3.

The apparent absence of multiple equilibria
Although the generic system of equations modeled by Eq. 3 admits multiple equilibria, we actually never observed any dynamic phenomena other than a steady convergence to a single stable fixed point, even when initial conditions were wildly different from that considered in the main text. It is not uncommon for infection models to admit equilibria which are physically unrealizable, i.e., taking values outside the unit simplexit may therefore be the case that any other stable equilibria are simply unphysical.
However, to comprehensively demonstrate this fact -either theoretically or numerically -is outside the scope of the current work. To prove such a result theoretically would necessarily draw on ideas of nonlinear algebra or algebraic geometry, such as intersection theory, which are outside the present scope; to demonstrate numerically that for all parameter sets and graph structures there is only one stable equilibrium is computationally infeasible, even on graphs of a given small size.
To illustrate the idea that there seems to be only one stable equilibrium per system, we fixed an exact distribution of treatment with ρ = 0.35 to be used on every random regular graph used in the main text. This means that, since all nodes in all graphs are labelled (and the random regular graphs in the main text have twenty demes each), seven labelled demes (here 1, 3, 8, 11, 15, 19, 20) were drug-treated irrespective of graph structure. We then performed, for all 30 graphs, 100 trials with uniformly random initial conditions over every deme, i.e., every deme can start with any possible arrangement of wild-type and resistant infectives (so long as the sum of the two does not exceed 1). We then compared the equilibria seen with the equilibria found for the initial condition used in the main body of the text.
There are many possible metrics which would illustrate the relative closeness of all the equilibria seen. In S17 Fig we have elected to plot the maximum matrix norm (the 2-norm) of the difference between the equilibrium seen in any of the 100 trials and the equilibrium from the main text. That is to say, all of the trial equilibria differ in norm from the main-text equilibrium by an amount not more than that plotted in the figure for a given graph structure. The parameters used were our default: κ=0.25, β=0.05/day, g=0.1/day, =0.9, c=0.2.
In all of these cases, the maximum norm deviation is extremely small. This could be due to two possible reasons. The first is due simply to aggregated numerical error from the many steps involved of a) calculating the equilibrium (which involves a threshold of stopping when the time-derivatives in all demes fall below a certain threshold, which can allow for some very small freedom in what an "equilibrium" is in a given trial), b) calculating the matrix norm (which amplifies these small errors). The second explanation is that while multiple stable equilibria may exist, they are extremely close together, with almost clinically-unmeasurable differences in the rates of wild-type and resistant infections within a given deme in a given community.

Model with effects of treatment on recovery rates
To demonstrate the robustness of the results described in the main text, we explored a second model in which the drug, rather than inhibiting the spread of the infection, instead speeds up recovery. The effect of treatment is to increase the rate of recovery from g to g + τ . In this model formulation a perfectly efficacious drug would correspond to τ → ∞ (instantaneous clearance of infection). A finite value of τ takes into account the fact that it takes some time for an infected individual to seek a clinician, obtain a prescription for drug, and clear the infection as a result of treatment.
This model matches qualitatively with that discussed in the main text, in that in reproduces the spatial motifs in the same way. The two models are also quantitatively quite similar for reasonably-akin parameter sets. Here we reproduce the main results for this alternative model.

A single well-mixed deme
Here, the system can be written aṡ (S.14) The basic reproductive ratio R 0 for each strain alone is For the purpose of comparing this model, the value of τ for which R w 0 in the main model (with our baseline = 0.9) equals R w 0 in this alternative model is τ = 0.9. This system admits only three stable, physical (having fractions of treated individuals between 0 and 1) solutions: one in which all infections die out (R w 0 < 1, R r 0 < 1), one in which only wild-type infections persist (R w 0 > 1, R r 0 < R w 0 ) and one in which only resistant infection persists (R r 0 > 1, R w 0 < R r 0 ). There is no possibility for coexistence. The condition for resistance to persist (instead of wild type)is c < f τ g+τ .

(S.16)
Using the next-generation technique, we were able to derive the values of R 0 in a two-deme system for the alternative model, as an analog to Eq. S.10: (S.17) The expressions are slightly more complicated, but the curves generated are similar to the main text. S7 Fig shows the analog of Fig 3. As in the main text, here the boundaries of the coexistence region are the curves R w 0 = R r 0 (lower boundary) and the implicit curve where more wild-type than resistant infected individuals are produced per unit time in the drug-treated deme (upper boundary).

Multi-deme population
The scaled version of the multi-deme model presented in Eq. 3 for the case of treatment acting on the recovery rates isẇ