Coupled environmental and demographic fluctuations shape the evolution of cooperative antimicrobial resistance

There is a pressing need to better understand how microbial populations respond to antimicrobial drugs, and to find mechanisms to possibly eradicate antimicrobial-resistant cells. The inactivation of antimicrobials by resistant microbes can often be viewed as a cooperative behaviour leading to the coexistence of resistant and sensitive cells in large populations and static environments. This picture is, however, greatly altered by the fluctuations arising in volatile environments, in which microbial communities commonly evolve. Here, we study the eco-evolutionary dynamics of a population consisting of an antimicrobial-resistant strain and microbes sensitive to antimicrobial drugs in a time-fluctuating environment, modelled by a carrying capacity randomly switching between states of abundance and scarcity. We assume that antimicrobial resistance (AMR) is a shared public good when the number of resistant cells exceeds a certain threshold. Eco-evolutionary dynamics is thus characterised by demographic noise (birth and death events) coupled to environmental fluctuations which can cause population bottlenecks. By combining analytical and computational means, we determine the environmental conditions for the long-lived coexistence and fixation of both strains, and characterise a fluctuation-driven AMR eradication mechanism, where resistant microbes experience bottlenecks leading to extinction. We also discuss the possible applications of our findings to laboratory-controlled experiments.


Introduction
Environmental conditions, such as temperature, pH, or available resources, endlessly change over time and shape the fate of natural populations.For instance, microorganisms often live in volatile environments where resource abundance fluctuates between mild and harsh conditions, and regimes of feast alternate with periods of famine [SK98,MK18,HM20].How environmental variability (EV), generally referring to changes not caused by the organisms themselves (e.g., supply of abiotic resources), affects species diversity is a subject of intense debate and research, see, e.g., [Gri73, MRS11, Fox13, CW81, Che94, Che00, ESAH19, MHRet al.21, VAME10, GPW10].Demographic noise (DN) arising from randomness in birth and death events in finite populations is another source of fluctuations.DN is negligible in large populations and strong in small ones, where it can lead to species fixation, when one species takes over the population, or to extinction, and hence can permanently set the makeup of a community [Ewe04,CK09,BM07,WSWP23].The dynamics of the population composition (evolutionary dynamics) is often coupled with that of its size (ecological dynamics) [Rou79], resulting in its eco-evolutionary dynamics [PGH09, KA14,WFM18].
When EV influences the size of a population, it also modulates the DN strength, leading to a coupling of DN and EV [WFM17, WFM18, WM20, TWAM20, SMM21, TWMA23].This interdependence is potentially of great relevance to understand eco-evolutionary dynamics of microbial communities.The coupling of DN and EV can lead to population bottlenecks, where new colonies consisting of few individuals are prone to fluctuations [WGSV02,PW09,BBG07,Bro07], and plays an important role in the eco-evolutionary dynamics of antimicrobial resistance (AMR) [CPL + 18, MB20].
The rise of AMR is a global threat responsible for millions of deaths [O'N16].Understanding how AMR evolves and what mechanisms can possibly eradicate the resistance to antimicrobials are therefore questions of great societal relevance and major scientific challenges.A common mechanism of antimicrobial resistance involves the production by resistant cells, at a metabolic cost, of an extra or intracellular enzyme inactivating antimicrobial drugs [Dav94, Wri05, YCD + 13].When the number of resistant cells exceeds a certain threshold, there are enough drug-inactivating enzymes, and the protection against antimicrobial drugs is shared with sensitive cells that can thus also resist antimicrobial drugs at no metabolic cost.However, below the resistant population threshold, only resistant microbes are protected against the drug (enzyme availability is limited and it can only inactivate the drug in the vicinity of resistant cells).AMR can hence be viewed as a thresholded cooperative behavior where widespread antimicrobial inactivation is a form of public good.This results in the spread of resistant microbes below the threshold, while sensitive cells thrive under high enzymatic concentration (above threshold).Hence, in static environments and large populations, both sensitive and resistant strains survive antimicrobial treatment and coexist in the long run [YCD + 13, VG14, MSL + 15, BWB16].In this work, we show that this picture can be greatly altered by the joint effect of demographic and environmental fluctuations, often overlooked, but ubiquitous in microbial communities that commonly evolve in volatile environments, where they can be subject to extreme and sudden changes [WGSV02, BBG07, Bro07, PW09, SPA + 12, SLKF12, CPL + 18].
Motivated by the problem of the evolution of AMR, here we study the coupled influence of EV and DN on the eco-evolutionary dynamics of a population of two species, one antimicrobial resistant strain and the other sensitive to antimicrobials.In our model, we assume that AMR is a cooperative behavior above a certain threshold for the number of resistant microbes, and the microbial community is subject to environmental fluctuations that can cause population bottlenecks.Here, EV involves random switches of the carrying capacity, causing the population size to fluctuate, while the antimicrobial input is kept constant.We thus study how the joint effect of EV and DN affects the fixation and coexistence properties of both strains, determining under which environmental conditions either of them prevail or if they both coexist for extended periods.This allows us to identify and fully characterise a fluctuation-driven antimicrobial resistance eradication mechanism, where environmental fluctuations generate transients that greatly reduce the resistant population and DN can then lead to the extinction of AMR.
In the next section, we introduce the model and discuss our methods.We present our results in section 3, where we first describe the main properties of the (in silico) model evolving under a fluctuating environment, and then study its properties analytically.In sections 3.1 to 3.3, we analyse the population dynamics in the large population limit, and then the model's fixation properties in static environments.In section 3.4, we characterise the fixation and coexistence of the strains in fluctuating environments, and discuss in detail the fluctuation-driven eradication of antimicrobial resistance arising in the regime of intermediate switching.Section 4 is dedicated to the discussion of the influence of environmental variability on the strains fraction and abundance (section 4.1), and to a review of our modeling assumptions (section 4.2).Our conclusions are presented in section 5. Technical and computational details are given at the end, in the annex supplemental material (SM).

Methods & Models
Microbial communities generally evolve in volatile environments: they are subject to suddenly changing conditions [SPA + 12, SLKF12], and fluctuations can play an important role in their evolution [Gri73, MRS11, Fox13, Che94, Che00, CPL + 18, ESAH19, MHRet al.21].For instance, fluctuating nutrients may be responsible for population bottlenecks leading to feedback loops and cooperative behavior [WGSV02, BBG07, Bro07, PW09], while sensitivity to antimicrobials depends on cell density and its fluctuations [Bro04, YCD + 13, VG14, MSL + 15, BWB16].Here, we study the eco-evolutionary dynamics of cooperative AMR by investigating how a well-mixed microbial community evolves under the continued application of a drug that hinders microbial growth when the community is subject to fluctuating environments.The evolutionary dynamics of the microbial community is modeled as a multivariate birth-and-death process [Gar02,vK92,BM07], whereas to model the fluctuating environment we assume that the population is subject to a time-varying binary carrying capac-ity [Ben06, HL06, TWAM20, TWMA23, HLGM16].

Microbial model
We consider well-mixed co-cultures composed of an antimicrobial resistant cooperative strain (denoted by R) and a defector type sensitive to antimicrobials (labeled S), under a constant input of antimicrobial drug, inspired by a chemostat laboratory set-up.The population, of total size N , hence consists of N R resistant and N S sensitive microbes, with N = N R +N S .Note that, since we later introduce EV as switches in the carrying capacity, the total population will fluctuate accordingly.A frequent mechanism of antimicrobial resistance relies on the production of an enzyme hydrolysing the antimicrobial drug in their surroundings [Dav94, Wri05, YCD + 13].Here we assume that each R cell produces the enzyme at a constant rate, regardless of the antimicrobial concentration, which is inspired by typical lab experiments, e.g., with resistance gene-bearing plasmids [YCD + 13, BWB16].When the number of R is high enough, the overall concentration of resistance enzyme in the medium suffices to inactivate the drug for the entire community: the enzyme hydrolyses the drug and sets it below the Minimum Inhibitory Concentration (MIC), therefore acting as a public good and protecting S as well.This mechanism can hence lead to antimicrobial resistance as a cooperative behavior [Dav94, Wri05, YCD + 13], for instance, by means of the β-lactamase resistance enzyme for the general β-lactam family of antibiotics [Bro09] (see section 4.2 for non-shared resistance mechanisms).
Here, we model this AMR mechanism by assuming that R acts as a cooperative strain when the number of R cells (proxy for resistance enzyme concentration) exceeds a fixed threshold N th , i.e., R cells are cooperators when N R ≥ N th , while they retain for themselves the benefit of producing the protecting enzyme when N R < N th [Bro04, YCD + 13, VG14, MSL + 15, BWB16].The effective regulation of public good production by means of a population threshold has been found in a number of microbial systems, see, e.g., [BJ01, CB06, BBG07, SG13, VG14], and is consistent with a slower microbial growth cycle with respect to the fast time scale of enzyme production and dispersion.In this work, we study the AMR evolution as a form of cooperative behavior under demographic and environmental fluctuations.Assuming fixed-volume fluctuating environments, the threshold for AMR cooperation is here set in terms of R abundance (rather than its concentration), see section 4.
In our model, R microbes have a constant birth rate independent of the biostatic drug hindering microbial growth [HA12, APNK07, SMM17]1 , with fitness f R = 1 − s, where 0 < s < 1 captures the extra metabolic cost of constantly generating the resistance enzyme.The birth rate of S depends on the public good abundance: when N R < N th , the enzyme concentration is low (below cooperation threshold) and the antimicrobial drug is above the MIC, the S fitness f S is thus lower than f R , with f S = 1−a, where 1 > a > s and a encodes growth rate reduction caused by the drug.When N R ≥ N th , the R abundance is above the cooperation threshold.This triggers the AMR cooperative mechanism: the drug is inactivated (below MIC), and the S birth rate, with f S = 1, is then higher than that of R, see figure 1a.Denoting by x ≡ N R /N the fraction of R in the population, here S fitness is where θ[z] is the Heaviside step function, defined as θ[z] = 1 if (z > 0) and θ[z] = 0 otherwise, and x th (N ) ≡ N th /N is the fraction of R at the cooperation threshold.The average population fitness is f = f R N R /N + f S N S /N .In this setting, this population evolves according to the multivariate birth-death process [Gar02,vK92,Ewe04] defined by the reactions occurring with transition rates [WFM17, WFM18, TWAM20, SMM21] with growth limited by the logistic death rate N/K (so that the total population N follows the standard logistic dynamics in the mean field limit, see equation ( 6)), where K is the carrying capacity, that is here assumed to be a time-fluctuating quantity, see below.Moreover, we have normalised f R/S by the average fitness f for mathematical convenience, without loss of generality (see section 4.2).This corresponds to the growth rate of each strain to be given by its fitness relative to the average population's fitness, a common assumption in the context of biological and evolutionary processes [Ewe04, CW81, TH09], which allows us to establish a neat relationship between our multivariate birth-death process and the classical Moran process.The latter is the reference birth-death-like process used to model the evolution of idealised populations of constant total size [Mor62, Ewe04, AS06, BM07, CMF11].The link with the Moran process enables us to take advantage of its well known properties, in particular the exact results for the fixation probability and mean fixation time [Mor62, Ewe04, AS06, BM07], to characterise analytically many features of our eco-evolutionary model (see sections 3.2, 3.3, and annex supplemental section D).
For simplicity, we consider that K(t) is driven by the colored dichotomous Markov noise (DMN) ξ(t) = {−1, 1} that randomly switches between K − and K + .The DMN is an important example of bounded noise, with finite correlation time, that is easy to simulate accurately (see annex supplemental section A) and amenable to analytical progress, and hence often employed in modelling evolutionary processes in fluctuating environments [Ben06, HL06, RDL11, HLGM16, WFM17, WFM18, WM20, SMM21, TWAM20].The dynamics of the DMN is defined by the simple reaction [Ben06, HL06, RDL11] endlessly occurring at rate (1 − δξ)ν, where −1 < δ < 1.Here, we always consider the DMN at stationarity where ξ = ±1 with probability (1 ± δ)/2.The stationary DMN ensemble average is thus , where ν is both half the inverse of the correlation time and average switching rate.We thus consider that the binary switching carrying capacity is [WFM17, WFM18, WM20, TWAM20] and K(t) thus switches from a state where resources are abundant (K + ) to another state where they are scarce (K − ) with rates ν + ≡ ν(1 − δ) and ν − ≡ ν(1 + δ) according to Environmental statistics can be characterised by the mean switching rate ν ≡ (ν − + ν + )/2 and by δ ≡ (ν − − ν + )/(ν − + ν + ) that encodes the environmental switching bias: when δ > 0, on average, more time is spent in the environmental state ξ = 1 than ξ = −1, and thus K = K + is more likely to occur than K = K − (symmetric switching arises when δ = 0).The time-fluctuating carrying capacity (4) modeling environmental fluctuations is responsible for the time-variation of the population size, and is coupled with the birth-and-death process (1) and (2).The master equation (ME) giving the probability P (N R , N S , ξ, t) for the population to consist of N R and N S cells in the environmental state ξ at time t is [Gar02]: where E ± R/S are shift operators such that E ± R/S f (N R/S , N S/R , t) = f (N R/S ± 1, N S/R , t), and the probabilities are set to P (N R , N S , ξ, t) = 0 whenever N R < 0 or N S < 0. The last line on the righthand-side of (5) accounts for the random environmental switching, see black line in figure 1b.Since T ± R/S = 0 whenever N R = 0 or N S = 0, this indicates that there is extinction of R (N R = 0) and fixation of S (N S = N ), or fixation of R (N R = N ) and extinction of S (N S = 0).When one strain fixates and replaces the other, the population composition no longer changes while its size continues to fluctuate2 .The multivariate ME (5) can be simulated exactly using standard stochastic methods (see annex supplemental section A), and encodes the eco-evolutionary dynamics of the model whose main distinctive feature is the coupling of the population size N and its composition x = N R /N , with DN coupled to EV, see (2) and below.

Results
In this section we analyse how the coupled demographic and environmental fluctuations shape the evolution of the fraction of R in cooperative AMR [CPL + 18].Our main goals are to establish the conditions under which EV and DN facilitate the eradication of R, and reduce the size of the remaining pathogenic microbial population (see also section 4).

Coupled environmental and demographic noise induces regimes of coexistence and dominance
The eco-evolutionary long-lived behavior of a microbial community is chiefly captured by: (I) the expected duration of the strains coexistence (Mean Coexistence Time, MCT, that here coincides with the unconditional mean fixation time [Ewe04,AS06]; see annex section B.2); and (II) by the fixation (or extinction) probability of each strain, i.e., the chance that a single strain eventually takes over the entire population (or that the strain is fully replaced by others).These properties have been extensively studied in populations of constant total size, e.g., in terms of the Moran process [Mor62, Gar02, Ewe04, AS06, BM07, TH09, PRS22], but are far less known in communities of fluctuating size when DN is coupled to EV.To gain some insight into the behavior of microbial co-cultures under coupled eco-evolutionary dynamics defined by equation ( 5), we compute in silico the R fixation probability, denoted by ϕ, and the strains coexistence probability, labeled by P coex , when the external conditions fluctuate between harsh (K − = 120, scarce resources) and mild (K + = 1000, abundant resources).Here, P coex is defined as the probability that both strains still coexist for a time exceeding twice the average stationary population size t > 2⟨N ⟩3 .In our simulations, we consider a wide range of the switching rate ν and bias δ, with ∼ 10 3 − 10 4 realizations for each dynamic environment, and different values of the cooperation thresholds, with N th ∼ 100.In our simulations, we respectively use s ∼ 0.1 − 0.2 and a ∼ 0.25 − 0.5 as plausible values for the resistance metabolic cost and the impact of the drug on S [vdHSS + 11, MWK15].Our choice of K ± ensures that the dynamics is not dominated mainly by DN or EV, but by the interplay of DN and EV, and the values of the cooperation threshold N th < K − guarantee that the fixation of either strain or their coexistence are all scenarios arising with finite probabilities in our simulations, see below and annex supplemental section A. Note that, as discussed in section 4.2 and annex section D.3, the behavior reported here can also be observed in big, realistic populations of N > 10 6 .Figure 2a-c shows the in silico ν − δ phase diagrams corresponding to the various fixation and coexistence scenarios arising for different cooperation thresholds.For small thresholds relative to EV (N th ≲ 10K + /K − , see annex supplemental section D.3), S displays a high fixation probability (red region) at intermediate ν and non-extreme δ, where R is most likely to be eradicated.Under high/low values of ν (when δ is not too low), the red region in figure 2a-c is surrounded by dark areas where the long-lived coexistence of the strains is most likely.When the threshold N th is closer to K − , R is most likely to prevail in the blue region of figure 2b-c, where the environment is predominantly in the harsh state (δ < 0).As N th increases, the blue region expands and gradually replaces the red and black areas: the fixation of R is likely to occur in most of the ν − δ diagram.In addition to the population makeup, the average population size is a decreasing function of ν at fixed δ, and increasing with δ at fixed ν; see figure 4d-e  In what follows, we analyse the different phases of figure 2a-c, focusing particularly on the characterization of the red area, and also determine how N R varies with the environmental parameters in the different phases.This allows us to determine the most favorable environmental conditions for the eradication of R and for the reduction of the population of pathogenic cells, which are issues of great biological and practical relevance.

Weak demographic noise promotes coexistence
To gain an intuitive understanding of the model's eco-evolutionary dynamics, it is useful to discuss the sample paths of figures 1b-c and 2d-f in terms of the population size N and the R fraction It is instructive to first consider the case of very large population arising with a constant and large carrying capacity K(t) = K 0 ≫ 1.In this setting, corresponding to a static environment, we ignore all forms of fluctuations and the system evolves according to the mean-field (deterministic where the dot indicates the time derivative.It is clear from equation (7) that the dynamics of the population composition, given by x, is coupled to that of its size N .According to the logistic equation ( 6), the population size reaches N = K 0 on a time scale t ∼ 1 independently of x, while the population composition is characterised by a stable equilibrium x = x th ≡ N th /N = N th /K 0 reached on a time scale of t ∼ 1/s or ∼ 1/(a − s) from x > N th /N or < N th /N , respectively.When s < a ≪ 1, there is a timescale separation, with N relaxing to its equilibrium much faster than x.We note that the coexistence equilibrium in terms of R and S is N eq R = N th and N eq S = K 0 − N th .Clearly, this suggests that S would unavoidably be wiped out if N th was greater than the carrying capacity, and hence we always consider that the latter exceeds the cooperation threshold (K > N th ).
When the population is large enough for demographic fluctuations to be negligible (1/ √ N → 0) and the sole source of randomness stems from the time-fluctuating environment (random switches of the carrying capacity), the dynamics becomes a so-called piecewise deterministic Markov process (PDMP) [Dav84].Between each environmental switch, the dynamics is deterministic and given by equations (6), with K 0 replaced by K ± in the environmental state ξ = ±1, and (7).Here, the PDMP is thus defined by 3 for the behavior at much larger populations and thresholds.Stronger blue (red) depicts a higher fixation probability of R (S).Darker color indicates higher coexistence probability, defined as the probability to not reach any fixation before t = 2⟨N ⟩, where we take the average total population in its stationary state.The area enclosed within the green solid line indicates the optimal regime for the eradication of R, see section 3.4.The white asterisks in (b) depict the environmental statistics for each of the bottom panels.(d-f ) Sample paths for the carrying capacity (K, black), number of R (N R , blue), number of S (N S , red), and fixed cooperation threshold N th = 80 (dashed blue) for the environmental parameters (ν, δ) depicted by the corresponding white asterisk in (b).
The high environmental switching frequency in (f) results in an effectively constant carrying capacity (K = K, dotted line); see section 3.2.
where the fluctuating carrying capacity K(t) is given by equation (4), coupled to (7).Sample paths of this PDMP are shown as solid lines in figures 1b-c and 2d-f.These realizations illustrate that N (t) tracks the switching carrying capacity K(t) independently of x, while x(t) evolves towards the coexistence equilibrium at the cooperation threshold x th (t) = N th /N (t), which changes in time as N varies.Hence, x increases when N R < N th , and it decreases when N R > N th .For extremely high environmental switching rate ν → ∞, the microbial community experiences a large number of switches, between any update of the population make-up.In this case, N is not able to track K(t), but experiences an effectively constant carrying capacity K = K ≡ 1/⟨1/K(t)⟩ obtained by self-averaging the environmental noise over its stationary distribution (see [WFM17, WFM18, WM20, SMM21, TWMA23]), leading to Hence, when ν → ∞, the community size is approximately N ≈ K and, provided that δ is not too close to −1 (K not too close to K − ), long-lived coexistence of both strains is likely (with abundances N R ≈ N th and N S ≈ K − N th ), as shown in figure 2f.

Antimicrobial resistance is robust to demographic noise in static environments
When EV causes a population bottleneck, DN about the coexistence equilibrium may cause the extinction of one strain and the fixation of the other (see figures 1b-c and 2d-e).To elucidate the fate of microbial communities under fluctuating environments, it is therefore necessary to first understand how a small community is able to fixate, or avoid extinction, in a static environment, when it is subject to a constant carrying capacity K 0 , with 1 ≪ K 0 ∼ K − ≪ K + .This condition ensures both fixation of one strain or long-lived coexistence are possible, i.e., demographic fluctuations, of order O(1/ √ K 0 ), matter but do not govern the dynamics.
Since the community composition tends to the coexistence equilibrium x → x th , see equation ( 7), the faster N dynamics reaches its steady state N → K 0 before any fixation/extinction events occur, see equation ( 6).Therefore, we assume a fixed N = K 0 .The evolutionary dynamics is thus modeled by the analytically tractable Moran process [Mor62, Ewe04, BM07, CMF11, WFM17, WFM18], where the population composition evolves stochastically by balancing each birth/death of R by the simultaneous death/birth of a S, according to the reactions with the effective transition rates ]. Due to DN, the R fraction fluctuates around x th until the eventual extinction of a strain.Therefore, from the classic Moran results (equation (S6) in the annex), we can derive a simplified, approximated expression for the R fixation probability by setting any initial composition directly at coexistence x 0 = x th , which yields where we now assumed (1 − a) N th , which is in line with our choices 0 < s < a < 1 and N th < K 0 .Here K * 0 is the microbial population size giving the same fixation probability 1/2 to R and S. In our examples, s = 0.1 and a = 0.25 (see section 3.1), and fixation equiprobability is reached at K * 0 ≈ 3N th , where the R and S abundance in the long-lived coexistence equilibrium are respectively N eq R ≈ K 0 /3 and N eq S ≈ 2K 0 /3. Figure 3a shows the excellent agreement between the approximation (9) (solid lines) and the exact R fixation probability of the underlying Moran process of annex supplemental equation (S6) (dotted lines), for different cooperation thresholds (N th = 20 − 100)4 .Equation (9) and figure 3a-b show that, in a static environment, the relative magnitude of the carrying-capacity-to-threshold ratio K 0 /N th with respect to K * 0 /N th ≈ ln (1 − a)/ ln (1 − s) clearly determines whether R fixates (for smaller K 0 /N th ), becomes extinct (larger K 0 /N th ), or coexists with S for a long time (larger K 0 /N th and large populations).To interpret these results we remember that the mean field behavior tends to N R = N th and N S = K 0 − N th .When K 0 /N th ∼ 1, N S ≈ K 0 − N th is small, and S is prone to extinction (R fixates).As the total population K 0 increases (at fixed cooperation threshold N th ), the equilibrium value N S ≈ K 0 − N th increases, making S less likely to go extinct.The fixation probability of S thus rises, and overcomes that of the strain R when K 0 > K * 0 .However, the expected time for the fixation of S increases exponentially with K 0 and, for large enough cooperation thresholds (typically for N th > 50), fixation takes too long and is unobservable in practice; see figure 3b.Note that the expected S fixation time (or R extinction time) saturates for K 0 > K * 0 because the equilibrium value N R ≈ N th is independent of the total population K 0 .For all examples in figure 2a-c we have K 0 = K + = 1000 > K * 0 when δ = 1 and K 0 = K − = 120 < K * 0 when δ = −1.This explains the dark areas (coexistence) in panels b-c where δ → 1, and the blue regions (R fixation) where δ → −1.In figure 2a, we observe dark regions (coexistence) for both δ = ±1 as the MCT is always larger than the coexistence threshold (t > 2K 0 ).Therefore, it appears that in static environments AMR always dominates or, at least, survives for extended periods.

Demographic noise can eradicate antimicrobial resistance in fluctuating environments
Under low and high environmental switching rates, the community behaves as in static total populations of size K ± and K, respectively; see supplemental sections D.1-D.2.Richer and novel dynamical behavior arises at intermediate switching rate (in figure 2a-c, see red areas around ν = 10 −2 − 10 0 ), when there are several environmental switches prior to fixation, and the quantities ϕ and P coex cannot be simply expressed in terms of their counterparts in a population of constant effective size.This switching regime is characterised by the full interplay of the ecological and evolutionary dynamics: as shown in figures 1b and 2e, environmental switches can thus lead to transient "bumps" and "dips" in N R (after the carrying capacity increases K − → K + or decreases K + → K − , respectively).The transient N R dips, together with demographic fluctuations caused by the population bottleneck (K + → K − ), can thus lead to the rapid eradication of R with the fixation of S (red areas in figure 2a-c).Each dip has a small but non-negligible probability to eradicate R and hence reduces the expected coexistence time.Therefore, the ingredients for this fluctuation-driven AMR eradication mechanism are: (I) intermediate environmental switching, so that the total population N fluctuates by tracking K(t) without lagging behind, see green lines in figure 1b; (II) a slower population composition x coupled to the faster N , so that the R population N R = xN experiences transient bumps and dips about its equilibrium N R ≈ N th , see blue lines in figure 1b-c; and (III) a small number of R at the bottom of transient dips N R ∼ 1, so that DN can drive R to extinction, see blue noisy line in figure 1b.
Here, we are interested in characterising the transient N R dips as the main fluctuation-driven mechanism leading to the possible eradication of R. To study their properties, it is useful to consider the PDMP description of the transient R behavior in large populations: , where K(t) and N (t) are respectively given by equations ( 4) and ( 8), and we assume N R < N th .We note that, after a switch from the mild to harsh environment (K + → K − ) in the absence of DN, R always survives the ensuing transient dip, and N R rises towards the coexistence equilibrium N R = N th ; see thick solid line in figure 1b.However, when K − ≪ K + and the microbial community experiences a population bottleneck, a transient dip to a small value of N R can form.When this occurs, R is prone to extinction caused by non-negligible demographic fluctuations (stronger when N R is small).
To characterise the region of the ν−δ phase diagram where transient N R dips cause eradication of R, we need to estimate t dip , defined as the time from the onset of the dip to when N R reaches its minimal value according to equation (10), see figure 1b.To determine t dip from (10) we require ṄR (t dip ) = 0 which, assuming From the solution of equation ( 6) with the initial condition N (t = 0) ≈ K + , we find: Ignoring DN, we can thus estimate the R population at the bottom of the transient dip N dip R , reached at t = t dip .This is, we find the R fraction at the bottom of the dip x (t dip ) in the small x limit of equation ( 7) and combine it with the above N (t dip ) to obtain (see annex supplemental section D.3): where we assumed that R started from N R (t = 0) = N th .Demographic fluctuations at the bottom of a dip are of the order N dip R .For DN to possibly drive R to extinction, and the fluctuation-driven eradication scenario to hold, it is necessary that N dip R ∼ 10 or lower.This condition is certainly satisfied when K − and N th are of comparable size (with K − > N th ), and each of order K + , which can also hold for realistically large populations of N > 10 6 , see section 4.2 and annex section D.3.
Under these sufficient requirements, the optimal environmental conditions to rapidly eradicate R in large but fluctuating populations can be estimated from equations ( 6)-( 7), (10), and (11).First, in the mild environment (K = K + ), R needs to be able to evolve to the coexistence equilibrium N R = N th , requiring a longer average duration of the mild environment ν −1 + as compared to the evolutionary time scale s −1 , i.e., ν −1 + ≳ s −1 .Second, after the switch from mild to harsh environment (K + → K − ), R needs to reach the bottom of the transient dip and experience demographic fluctuations, which imposes an average duration of the harsh environment ν −1 − longer than the average time to reach the bottom of the dip t dip , that is, ν −1 − ≳ t dip .Third, if R survives the dip, the environment should go back to the mild ξ = 1 state to rule out the extinction of S when the environment stays in the harsh state ξ = −1.For this, we require the harsh environment to be short, while ensuring that the dip is not interrupted by a switch; see figure 1b-c.This enforces ν −1 − ≲ 2 ln ( K+ K− )(a − s) −1 , where the right-hand-side, derived from the small x limit of equation ( 7), is twice the expected time to reach the equilibrium in the harsh state N R = N th and N S = K − − N th .As a fourth condition, we demand that this cycle should be as fast as possible to maximise the number of transient dips (while still allowing the population to evolve back to N R = N th after a bump), yielding ν −1 + ≲ 2 ln ( K+ K− )s −1 , which, similarly as in the previous condition, is twice the average time needed to return to equilibrium in the mild state.Using the environmental parameters ν and δ, the above lead to s The green contour lines in figure 2a-c enclose the predicted optimal region for the fast eradication of R under s = 0.1, a = 0.25, K − = 120 and K − = 1000, and fall in the red areas observed in silico.The borders of these regions depend on N th .This stems from the dependence of ϕ and MCT on N th (see figure 3b) and the criterion for long-lived coexistence (t > 2⟨N ⟩).The conservative prediction (13) ignores any dependence on N th .
In summary, DN can eradicate antimicrobial resistance in fluctuating environments when the population make-up x evolves on a much slower timescale than the population size N , which requires relatively small values of s and a.Moreover, the variability in the carrying capacity K + /K − needs to be of the order of the cooperation threshold N th or larger; the threshold has to fall below the lowest value of the carrying capacity K − ; and the switching rate ν has to be of order s and hence comparable to the rate of relaxation of the population composition.Note that all conditions above can be met in biologically relevant systems of any size; see section 4.2 and annex supplemental section D.3.

Discussion
The results of the previous section characterise the long-term microbial population makeup under random switches between mild and harsh environmental conditions (high and low carrying capacity, K = K + and K − , respectively), for a broad range of the exogenous parameters (mean switching frequency ν and switching bias δ).Another important aspect of the time evolution of microbial population concerns the nontrivial impact of the environmental variability on the fraction and abundance of drug-resistant (R) and drug-sensitive (S) microbes in the different regimes, and especially in their phase of coexistence.It is also important to review to what extent our modeling assumptions are amenable to experimental probes.

Impact of environmental variability on the strains fraction and abundance
It was recently shown that in the fluctuating environment considered here, the average size of the microbial community ⟨N ⟩ is a decreasing function of the random switching rate ν (with δ kept fixed), that ⟨N ⟩ decreases with lower δ (keeping ν fixed), and that ⟨N ⟩ → K ± as δ → ±1 [WFM17, WFM18, TWAM20, SMM21]; see figure 4d-e.As a consequence, in the blue and red areas of the phase diagrams of figure 2, where only one strain survives (figure 2a-c), the surviving pathogenic population can be reduced by increasing the environmental switching frequency ν and/or the time spent in harsh state (by enforcing δ → −1).Moreover, since the R fraction x is directly coupled to N through the cooperation threshold N th , see equation ( 7), environmental variability (EV) non-trivially shapes the R fraction in the coexistence regime (colored areas in figure 4f).Under low switching frequency relative to the rate of evolutionary dynamics (see ν ≪ resistance extra metabolic cost s ∼ 10 −1 in all figures), R cells starting in the mild environment (K + ) are able to reach the coexistence equilibrium N R = N th before experiencing a switch; see figure 1c.However, if the starting environment is harsh (K − ), demographic noise (DN) can rapidly eradicate S and destroy coexistence (see figure 2b-c).The distributions of N R , N S and N in the regime ν → 0 are thus approximately bimodal because they combine both mild and harsh (effectively constant) environments; with N R ≈ N th , N S ≈ K + − N th , and N ≈ K + for the former; and N R ≈ K − , N S ≈ 0, and N ≈ K − for the latter; see figure 4a.As the switching rate is increased to ν ≲ s, fixation dominates, and the N R and N S bimodal distributions become approximately trimodal, that is, N R/S ≈ 0, K − or K + ; see figure 4b.The relative weight of the peaks at N R/S ≈ 0 is set by S/R fixation probability, which is modulated by the environmental bias δ, see supplemental equation (S9) in the annex.The total population distribution is still bimodal about N = K ± since its relaxation dynamics (of timescale ∼ 1, see equation ( 6)), is faster than the evolutionary timescale ∼ 1/s.Finally, when ν is increased further (ν ≫ s), we enter the coexistence regime characterised by an effective carrying capacity K = K [WFM17, WFM18, WM20, TWAM20, SMM21, TWMA23], and all distributions become unimodal about the coexistence equilibrium N R ≈ N th , N S ≈ K − N th , and N ≈ K; see figure 4c.As a consequence, if R is eradicated, imposing high EV (ν ≫ s) and harsh conditions δ → −1 would considerably reduce the abundance of the surviving community of pathogenic S cells; see figure 4d, green solid line, and figure 4e.However, if R survives, imposing ν ≫ 1 and δ < 0 would not only decrease the abundance of both strains but it would also increase the R fraction, and risk further antimicrobial resistance (AMR) spreading, see figure 4f (magenta /bluish areas).

Review of the modeling assumptions
Since we study an idealised microbial model, it is important to review our modeling assumptions in light of realistic laboratory experimental conditions.A key assumption to consider is the effectively sharp cooperation threshold N th , which is based on a number of experimental observations of microbial cooperation; see [Dav94, Wri05, YCD + 13, SG13].Accordingly, we have assumed that EV changes chemical concentrations (e.g., nutrient density) while the volume of the microbial ecosystem is kept constant [SG13].The cooperation threshold is then fixed at a constant number of R microbes N R = N th because, at constant volume, the resistance enzyme concentration is proportional to the number of public good producers R.This crucial ingredient fixes the stable number of R at N th across fluctuating environments, and is responsible for the transient dips which are at the origin of the novel eco-evolutionary mechanism for the eradication of AMR reported here.The complementary scenario, where the cooperation threshold is set by a fixed R fraction x th is also relevant (for a different set of microbial ecosystems), and is a topic for future research.Furthermore, in some microbial cases, R could regulate the production of resistance enzyme by quorum sensing [PTY12], but its impact on cooperative AMR remains an open problem.We also note that some resistance mechanisms can show anti-cooperative behaviour, such as efflux-pumps, which could result in an enhanced exposure of sensitive cells to the drug [Poo07,Sot13].In the case of non-shared resistance mechanisms, our model reduces to the eco-evolutionary processes studied in [WFM17, WFM18, TWAM20, SMM21].Further analytical results for the non-shared resistance models are discussed in [UH11], as well as in [Lam06, PQ07, PW08] in the case of a static environment.
A second assumption to review concerns the simulation results obtained here, for systems with K ± ∼ 10 2 − 10 3 and N th ∼ 100, that we are able to computationally probe (see annex supplemental section A) but that correspond to populations of relatively small size.In the supplemental section D.3 we provide a detailed discussion on how the rich microbial behavior and novel eco-evolutionary AMR eradication mechanism reported here can be translated to larger, more realistic, microbial communities of size of order N ≳ 10 6 [SG13] to N ≳ 10 8 [CPL + 18, SWJ56, Can56, Fel76, PDH + 07], or higher.In our discussion we argue that, as long as N th K − /K + ≲ 10 and 0 < s < a ≲ 10 −1 − 10 −2 , regardless of the magnitude of K ± or N th , the transient dips studied here will drag R close to extinction, where demographic fluctuations are instrumental for the likely and rapid eradication of AMR.Note that, for very fast/slow fluctuating environments, where transient dips are hindered (see section 3.4), R and S populations will always coexist unless K − − N th ≲ 10.
A specificity of our study is its focus on biostatic antimicrobial drugs.However, since most antimicrobials gradually change from acting as biostatic to biocidal as their concentration in the medium grows [HA12, APNK07, SMM17], our approach is consistent with a low antimicrobial concentration scenario.Conveniently, the combined biostatic effect of the drug and the normalization of strain fitness in equation (2) [Ewe04,BM07] decouples the total population N from its composition x.If any of the above conditions would not hold, N would then directly depend on x, a case already studied for a simpler model in [WFM17,WFM18].It is worth noting that we have confirmed that the main findings reported here are robust, as they do not depend crucially on the detailed choice of the transition rates in equation (2), and in particular they are found to be essentially independent of the normalization by the average fitness; see sections 2.1 and 3.3.We also note that the values used in our examples for the extra metabolic cost to generate the resistance enzyme (s ∼ 10% to 25%), and for the impact of the antimicrobial drug on S growth (a ∼ 25% to 50%), while only indicative, are plausible figures [vdHSS + 11, MWK15].
For the sake of simplicity, we have focused on modelling EV through binary switches of the carrying capacity.These switches capture sudden changes in the available resources (as in feast and famine cycles [SK98, BSO18, MK18, HM20]) that can also occur in presence of antimicrobial drugs, e.g., in polluted environments or during drug treatment.In the context of evolutionary processes, environments that fluctuate via random switches are commonly modelled in terms of dichotomous Markov noise (DMN), also known as telegraph noise [TvO04, Ben06, HL06, HLGM16, HSM17, WFM17, MB20, WM20, TWAM20, SMM21, TWMA23].Moreover, binary switching is the standard way to implement EV in laboratory-controlled experiments, where the concentration of nutrients can be regulated in a chemostat set-up [AMvO08, LK14, ARTG21, NLGS21].Although laboratory experiments are often carried out with periodically switching environments (e.g., [AMvO08,LK14]), and natu-ral environmental conditions often vary continuously in time and magnitude (e.g., [NLGS21]), the relationship between DMN and other commonly used forms of EV has already been extensively studied [Ben06,HL06,TWMA23].Therefore, our choice of modelling EV with DMN is natural, convenient and non-limiting: it allows us to make mathematical progress while keeping the theoretical modelling close to laboratory experimental conditions.The literature suggests that the essence of our findings are expected to hold for general fluctuating environments with a time-varying carrying capacity, but the extent to which other and more complex forms of EV than binary random switching may alter our results for microbial communities exhibiting cooperative AMR remains a problem to be studied.
Finally, we note that the novel eco-evolutionary mechanisms reported in this study to eradicate cooperative AMR, and to reduce the total pathogenic microbial community, or minimise the coexistence fraction of R, all take place at a biologically and clinically relevant range of environmental switching rates.Indeed, although our theoretical study does not set a specific timescale of microbial growth, a plausible rough estimate for a single replication cycle of a microbe could be of the order of ∼ 1 hour.The novel AMR eradication mechanism at ν ∼ s then comes into play when a single environmental phase lasts, on average, 1/s ∼ 10 hours.This could be consistent with the periodic administration of a treatment that enforces microbial population bottlenecks, and is a feasible time scale for laboratory experiments.Our idealised model however assumes a homeostatic influx of antimicrobial drug in all environments.Thus, an interesting approach for future work would involve the joint application of antimicrobial drug and population bottlenecks (in the harsh environment), with no drug administered in the mild environmental state.

Conclusion
Understanding how environmental variability affects the demographic and ecological evolution of microbes is central to tackle the threat of AMR, an issue of pressing societal concern [O'N16].Central questions in studying AMR involve how the fraction of resistant microbes changes in time, and by what mechanisms these can possibly be eradicated.
It is well established that AMR is an emergent property of microbial communities, shaped by complex interactions.In particular, certain resistant cells able to inactivate antimicrobials can, under certain conditions, protect the entire microbial community.This mechanism can hence be viewed as an AMR cooperative behavior.Moreover, microbial populations are subject to changing conditions.For instance, the size of a microbial population can vary greatly with the variation of the nutrients or toxins, and can, e.g., experience bottlenecks.As a result of evolving in volatile environments, microbial communities are prone to be shaped by fluctuations.In general, these stem from environmental variability (EV, exogenous noise) and, chiefly in small populations, from demographic noise (DN).The underlying eco-evolutionary dynamics, characterised by the coupling of DN and EV, is ubiquitous in microbial ecosystems and plays a key role to understand the AMR evolution, but is still rather poorly understood.
In this work, we have studied an idealised model of cooperative AMR where a well-mixed, microbial population consisting of sensitive and resistant cells is treated with an antimicrobial (biostatic) drug, hindering microbial growth, in a fluctuating environment.The latter is modeled by a binary switching carrying capacity that fluctuates between two values corresponding to mild and harsh conditions (high/low values, respectively).Based on a body of experimental work [YCD + 13, VG14, MSL + 15, BWB16, Dav94, Wri05], we assume that resistant cells produce, at a metabolic cost, an enzyme that inactivates the antimicrobial drug.Importantly, the abundance of resistant microbes is thus a proxy for the concentration of the drug-inactivating enzyme, which above a certain abundance threshold, becomes a public good by providing drug protection, at no metabolic cost, to the sensitive strain.Above the cooperative threshold, the latter hence have a fitness advantage over the resistant strain, whereas, below the threshold, the drug is responsible for a reduced fitness of the sensitive cells.In this setting, the evolution of AMR can be viewed as a public good problem in a varying environment, whose outcome is shaped by the coupling of environmental and demographic fluctuations.
We have identified three regimes characterising the eco-evolutionary dynamics of the model, associated with the fixation of the resistant or sensitive microbes, or with the long-lived coexistence of both strains.Our analysis shows that, while AMR generally survives, and often prevails, in static environments, a very different scenario can emerge under environmental variability.In fact, we demonstrate that fluctuations between mild and harsh conditions, coupled to DN, can lead to "transient dips" in the abundance of resistant microbes, which can then be driven to extinction by demographic fluctuations.Here, we determine that this fluctuation-driven AMR eradication mechanism occurs when the rate of environmental change is comparable to that of the relaxation of the evolutionary dynamics (ν ∼ s).By computational means, we show that this fluctuation-driven mechanism speeds up the eradication of resistant cells, and argue that it holds also for large microbial communities, comparable to those used in laboratory experiments (N > 10 6 ).We have also studied how EV non-trivially affects the strain abundance in the various regimes of the model, and in particular have determined the complex long-lived distribution of the fraction of resistant cells when both strains coexist and the environment fluctuates.
In conclusion, we have shown the existence of a biophysically plausible novel mechanism, driven by the coupling of EV and DN, to eradicate resistant microbes, and have demonstrated how EV shapes the long-lived microbial population in the possible scenarios of strains coexistence or fixation.Our work thus paves the way for numerous possible applications, for instance, in microbial experiments with controlled environmental fluctuations (it is currently possible to track even individual microbes, e.g., see [BLB + 21, MSCD + 21]), which might shed light on new possible treatments against AMR in real-world clinical infections.

B.2 Exact General Mean Coexistence Time (MCT)
It is also possible to exactly compute the general mean duration of coexistence regardless of the final state (either fixation or extinction of R), i.e., the Mean Coexistence Time (MCT) t N 0 R , N , when the transition rates T ± R are time-independent.It is worth noting that the MCT here coincides with the unconditional mean fixation time (and mean extinction), since t N 0 R , N gives the mean time after which coexistence is lost due to the fixation of one strain and the extinction of the other.
The exact formula for the MCT reads [Gar02, vK92, Ewe04, AS06, TH09]: where γ N 0 R , N and T ± R N 0 R , N are defined as in equations (S1) and (S2).For clarity, we split in two the sum over k in the first numerator (from k = 1 to N 0 R − 1, and from k = N 0 R to N − 1), and then rearrange the equation as where ϕ ≡ ϕ N 0 R , N is the R fixation probability starting from N R (t = 0) = N 0 R in a population of overall size N , and is given by (S1).

B.3 Coexistence probability
To derive an expression for the coexistence probability P coex at a fixed total population N , cooperation threshold N th , and starting with the number of R cells at equilibrium N 0 R = N th , we first compute the exact Moran MCT t(N th , K 0 ) from equation (S8); see main figure 3b and supplemental figure S1b for N = K 0 , dotted lines.Since the microbial community tends to a coexistence equilibrium, the fixation of a strain occurs on a slow time scale, driven by fluctuations.We thus assume that the full density of coexistence times roughly approximates an exponential distribution of mean t(N th , K 0 ), a known property of systems exhibiting metastability [AM17].P coex is thus the exponential cumulative probability remaining after an elapsed time 2K 0 , i.e., P (t > 2 ⟨N ⟩ = 2K 0 ), as previously computed in silico (but for a wide range of ν); see main section 3.1 and figure 2a-c.

C Full model in a static environment with constant carrying capacity
In this section we relax the Moran approximation of a fixed total population size, and allow N to fluctuate around a constant capacity, here denoted by K 0 .In this case, the behavior of the microbial community does not directly correspond to the Moran process, and follows a bivariate process in terms of the number N R/S of R/S individuals (but it does not depend on ξ since the environment is here static).The probability P (N R , N S , t) that the population consists of N R and N S at time t, now satisfies the ME where the transition rates T ± R/S are given by main manuscript, equation (2), and the carrying capacity is now constant as K → K 0 .
For direct comparison, in figure S1 we display exact simulations of the above behavior on top of the exact Moran results of main figure 3 (dotted lines).We observe that the theoretical predictions of the Moran process quantitatively capture the results in silico with some minor systematic deviations.These small discrepancies arise from the fact that, below K * 0 /N th = ln (1 − a)/ ln (1 − s) ≈ 3 (see main manuscript, section 3.3), the MCT depends exponentially on K 0 /N th .Therefore, it is exponentially more probable to observe faster fixation under smaller populations, i.e., when random fluctuations drive N below its expected value ⟨N ⟩ = K 0 .
The Moran approximation, based on assuming N = K 0 , misses demographic fluctuations of order √ K 0 about K 0 , which results in underestimating ϕ and overestimating the MCT.The relative amplitude of these deviations scale with those of the standard deviations of N/K 0 which are of order O( √ K 0 /K 0 ) and hence become vanishingly small when K 0 → ∞.Hence, the analytical predictions from the Moran process at N = K 0 are relevant for static environments because they quantitatively capture both the fixation probability and MCT in silico within a small relative error.Moreover, this error decreases the bigger the total population and cooperation threshold as the fluctuations about fixed K 0 /N th become negligible, which is the case for more realistic, biologically plausible microbial populations.D Dynamic environments: additional analytical derivations D.1 Infrequent environmental switching limit ν → 0 When ν → 0, the average number of environmental switches prior to fixation of one strain (extinction of the other) is very low, and the fixation probability can be obtained by averaging its constant-N counterpart over the stationary distribution of ξ.We can indeed assume that the community evolves subject to a static carrying capacity set by the starting environment K(t = 0) ≡ K 0 , where K 0 = K ± with probability (1 ± δ)/2.For these infrequently switching environments, the R fixation probability at an arbitrary environmental bias δ is the average To this end, we first notice that, as discussed in section 4 of the main text, the environmental parameters must fulfill 1 < K − /N th to avoid that S dives into extinction following a switch to the harsh environment.For this, the stable number of S, N S = K − −N th (see section 3.2 in the main text), must be high to resist demographic fluctuations, for which a reasonable estimate is K − −N th > 10.The second point to notice is that the model tends to two possible coexistence equilibria 0 < x = N th /K ± < 1 as long as K − > N th (see section 3.2).Hence, the expected behavior for large populations in the two possible static environments K ± is coexistence, which has a duration that scales exponentially with the cooperation threshold N th ; see main figure 3b.Since typically the resistant cooperation threshold is of the order N th ≲ K − , the duration of microbial coexistence is thus expected to grow exponentially with K − .Therefore, fixation in static environments is never observed in realistically big communities where K − ≳ 10 6 .
However, as discussed in the main section 3.4, the coupled eco-evolutionary dynamics in fluctuating environments generates significant transient N R dips when switching from mild to harsh environments (K + → K − ) at an intermediate switching rate ν ∼ s.The frequency and depth of these transient dips are maximised in a certain range of environmental parameters ν and δ, derived in section 3.4.The rapid eradication of R in this optimal dynamic environment regime is shown in the green-enclosing areas of figure 2a-c.We now ask whether realistically big values of {N th , K − , K + } actually enhance the extinction of cooperative AMR in the fluctuation-driven optimal AMR eradication regime {ν, δ}.
In such an optimal AMR-eradication regime, the minimum possible expected number of R, reached in the transient dips, is given by the main equation (12).To derive this equation, on the one hand we take the low R fraction limit x → 0 in the main equation ( 7) and thus, assuming x(t = 0) = x eq + ≡ N th K+ , where t = 0 is the time at the K + → K − environmental switch, we obtain the microbial composition dynamics at short times On the other hand, we can exactly solve the total population logistic dynamics of main equation (6) in the K − ≫ 1 environment after the switch, with N (0) = K + , as N (t) = K + K − e t K + (e t − 1) + K − .
Taking K − /K + ≪ 1 in t dip given by main equation (11); evaluating x (t = t dip ) above; noting that N (t = t dip ) ≡ K − (1 − s)/(1 − a) from α R ≃ 0 in main equation (10), see main section 3.4; and multiplying both resulting expressions, provides the final estimate of the R number N dip R at the bottom of the transient dip: Since we focus on the case where DN eradicates R, for the transient-dip fluctuation-driven eradication mechanism to possibly work, we then need large demographic fluctuations (of order N dip R ) relative to N dip R .This typically suggests to consider N dip R ∼ 10 or lower.Big cooperative AMR microbial communities in ecosystems with antimicrobial drugs could present populations of, for instance, N ≈ K + ∼ 10 12 in nutrient abundance conditions.Moreover, biophysically plausible values for the remaining parameters could be s = 0.1, a = 0.25, and a resistant cooperation threshold N th = 2 • 10 6 for an example fixed environmental volume.Sudden and drastic ecological bottlenecks in events of nutrient scarcity (or additional toxins) could then kill most of the community and reduce the total population by a factor of, e.g., few in a million; where the total population would decrease from N ∼ 10 12 to N ≈ K − ∼ 5 • 10 6 .These realistic parameters would fulfill the condition 1 < K − /N th = 2.5 with K − − N th ∼ 10 6 ≫ 10.Crucially, these plausible values would fulfill N dip R ≈ 1.7 • 10 ≲ 10 in equation (S12), so that there would be a significant chance that R becomes extinct during each transient N R dip in dynamic environments.
Finally we note that, as shown above, the relative magnitude of the population bottleneck K − /K + is critical to enhance the extinction of R during transient dips.Any increase in the population drop between abundance and scarcity environments, i.e., K − /K + → 0, boosts the eradication of AMR; whereas a reduction in the drop size, i.e., K − /K + → 1, hinders R extinction and promotes strain coexistence.Tuning the nutrient abundance or scarcity levels in each environment modulates the relative magnitude of the bottleneck K − /K + .But, additionally, introducing an intermediate environmental step between harsh and mild regimes could reduce the population bottleneck and boost coexistence [SG13].

Figure 1 :
Figure 1: Microbial community model.(a) Top: When the abundance of R (blue microbes) is below the cooperation threshold N th , antimicrobial drug hinders the growth rate of S (red microbes) and R cells have a growth advantage.Bottom: AMR becomes cooperative when the number of R exceeds N th and these generate enough resistance enzyme (public good in green shade) to hydrolyse the antimicrobial drug below the MIC for the whole medium, so that protection against the drug is shared with S (with green shields).(b) Temporal eco-evolution dynamics of the microbial community for example parameters s = 0.2, a = 0.5, K − = 50, K + = 250, ν = 0.2, and δ = 0.6; thick black line shows the sample path of the time-switching carrying capacity K(t), with a cooperation threshold N th = 30 (dashed blue line); thick solid lines depict the N → ∞ piecewise deterministic (deterministic between two switches of K) process defined by equations (7) and (8) for the total microbial population (N , green), number of R (N R = N x, blue), and number of S (N S = N (1 − x), red); noisy lines show an example stochastic realization of the full model under the joint effect of demographic and environmental fluctuations.In the absence of DN, R can experience bumps and dips (thick blue line), and t dip indicates the mean time to reach the bottom of a dip from its inception; see section 3.4.In the presence of DN, fluctuations about the dip can lead to the extinction of R (blue arrow).(c) R fraction x = N R /N for the same sample path of varying environment as in (b); line styles as in panel (b); the dashed black line shows the stable R fraction in each environment as K(t), driven by ξ(t), switches in time.

Figure 2 :
Figure2: Eco-evolutionary dynamics in the phase diagram of the joint fixation and coexistence probability.(a-c) Fixation and coexistence joint probability in silico at a given environmental bias δ and mean switching frequency ν for s = 0.1, a = 0.25, K − = 120, and K + = 1000 at resistant cooperation thresholds N th = 60, 80, and 100; see the discussions in sections 3.3, 4.2, and annex supplemental section D.3 for the behavior at much larger populations and thresholds.Stronger blue (red) depicts a higher fixation probability of R (S).Darker color indicates higher coexistence probability, defined as the probability to not reach any fixation before t = 2⟨N ⟩, where we take the average total population in its stationary state.The area enclosed within the green solid line indicates the optimal regime for the eradication of R, see section 3.4.The white asterisks in (b) depict the environmental statistics for each of the bottom panels.(d-f ) Sample paths for the carrying capacity (K, black), number of R (N R , blue), number of S (N S , red), and fixed cooperation threshold N th = 80 (dashed blue) for the environmental parameters (ν, δ) depicted by the corresponding white asterisk in (b).The high environmental switching frequency in (f) results in an effectively constant carrying capacity (K = K, dotted line); see section 3.2.

Figure 3 :
Figure3: Moran theory for R fixation probability and Mean Coexistence Time (MCT) in static environments.(a) R fixation probability ϕ in terms of the total microbial population normalised by the resistant cooperation threshold K 0 /N th for five example thresholds, from N th = 20 (dark green) to 100 (yellow green); the starting microbial composition is set at the coexistence equilibrium x th = N th /K 0 ; solid lines show the approximated prediction of equation (9); dotted lines depict the exact Moran behavior of annex supplemental equation (S6), only distinguishable for the smallest threshold; open diamonds illustrate the predicted K * 0 /N th that provides fixation equiprobability for each resistant cooperation threshold, see section 3.3.(b) Mean Coexistence Time vs K 0 /N th in log-linear scale; solid lines show the exact Moran MCT, computed from supplemental equation (S8); legend and symbols as in panel (a).

Figure 4 :
Figure 4: Total population, strain abundance, and coexistence composition in fluctuating environments.(a-c) In silico probability distributions of the total population (N , green), number of R (N R , blue), and number of S (N S , red), with parameters s = 0.1, a = 0.25, K − = 120, K + = 1000, and N th = 80, under no environmental bias (δ = 0) and for mean switching rates in slow ν = 10 −4 (a), intermediate ν = 10 −1 (b), and fast ν = 10 2 (c) conditions.Histograms are smoothed by a Gaussian filter of width σ = 10 in cell number.(d) Average overall population (number of individuals on the vertical axis) and strain abundances under no bias, i.e., δ = 0, as a function of switching rate ν; colors as in (a-c).Lines are smoothed by a log-scale Gaussian filter of width σ = 10, i.e., one frequency decade.(e) Average overall population size in dynamic environments.(f ) Coexistence composition and fixation probability (of any strain) in dynamic environments.Stronger blue (red) depicts a higher coexistence fraction of R (S).Lighter color indicates lower coexistence probability, defined as the probability for no fixation event to occur before t = 2⟨N ⟩.The white and black asterisks in (e-f) depict the environmental statistics for each of the top panels.All panels are computed at quasistationarity reached after a time t = 2⟨K⟩, ensuring that N reaches its (quasi-)stationary state, where ⟨N ⟩ ≤ ⟨K⟩; see text.

Figure S1 :
Figure S1: Full model simulations in static environments against exact Moran theory for R fixation probability and Mean Coexistence Time (MCT).(a) R fixation probability ϕ in terms of the total microbial population normalised by the resistant cooperation threshold K 0 /N th for three example thresholds, N th = 20 (dark green), 40 (green), and 60 (yellow green); the starting microbial composition is set at the coexistence equilibrium x 0 = x th = N th /K 0 ; dotted lines depict the exact Moran behavior of equation (S6), solid lines show simulation data of the full model averaged over 10 3 realizations.(b) Mean Coexistence Time vs K 0 /N th in log-linear scale; dotted lines show the exact Moran MCT, computed from equation (S8), solid lines show averaged simulation data of the full model over 10 3 runs; legend and symbols as in panel (a).