Coexistence of Competing Microbial Strains under Twofold Environmental Variability and Demographic Fluctuations

Microbial populations generally evolve in volatile environments, under conditions fluctuating between harsh and mild, e.g. as the result of sudden changes in toxin concentration or nutrient abundance. Environmental variability (EV) thus shapes the long-time population dynamics, notably by influencing the ability of different strains of microorganisms to coexist. Inspired by the evolution of antimicrobial resistance, we study the dynamics of a community consisting of two competing strains subject to twofold EV. The level of toxin varies in time, favouring the growth of one strain under low drug concentration and the other strain when the toxin level is high. We also model time-changing resource abundance by a randomly switching carrying capacity that drives the fluctuating size of the community. While one strain dominates in a static environment, we show that species coexistence is possible in the presence of EV. By computational and analytical means, we determine the environmental conditions under which long-lived coexistence is possible and when it is almost certain. Notably, we study the circumstances under which environmental and demographic fluctuations promote, or hinder, the strains coexistence. We also determine how the make-up of the coexistence phase and the average abundance of each strain depend on the EV.

In the absence of detailed knowledge about the timevariation of external factors, EV is generally modelled by means of noise terms affecting the species growth and/or death rates [26,.Demographic noise (DN) is another important source of fluctuations: it can lead to fixation, which is the phenomenon arising when one strain takes over the entire community.The effect of DN is significant in communities of small size, and becomes negligible in large populations [41,42,[67][68][69].Significantly, the time development of the size and composition of populations are often interdependent [70][71][72][73][74][75][76][77][78], with fluctuations of the population size modulating the strength of DN [57, 61-63, 65, 66, 79].The interplay between EV and DN is crucial in shaping microbial communities, but the quantitative effects of their coupling are as yet still mostly unknown.
Environmental and demographic fluctuations play a crucial role in the evolution of AMR [13,38,40,60,64].When treatments reduce a microbial community to a very small size, but fail to eradicate the microorganisms resistant to the drugs, resistant cells may replicate and restore infection, hence possibly leading to the spread of antibiotic resistance.Moreover, within a small population, demographic fluctuations may also lead to the extinction of resistant cells.Understanding how the coexistence of cells resistant and sensitive to antibiotics is affected by the joint effect of environmental variability and demographic fluctuations, and how the fraction of resistant cells varies with environmental conditions with the possibility of eradication, are thus central questions in the effort to understand the evolution of AMR [13,36,38,40,79].
It is worth noting that considerable efforts have recently been dedicated to study the mechanisms underpinning the coexistence of competing species under various scenarios, see e.g.[31,32,34,65].The influence of different kinds of variability (e.g.quenched disorder, heterogeneous rates, "spillover protection") on species coexistence in ecosystems exhibiting cyclic dominance [59,80] has notably been investigated in [47,[81][82][83][84][85].Here, inspired by the AMR evolution in a chemostat setup [19,65], we study the eco-evolutionary dynamics of an idealised microbial community consisting of two competing strains subject to a time-varying level of toxin, with the growth of one strain favoured under low toxin level and a selective advantage to the other strain under high toxin level.We also assume that the resource abundance varies according to a time-switching carrying capacity that drives the fluctuating size of the community.In most of previous works, EV is either encoded in fluctuating growth rates, with the size or carrying capacity of the population kept constant [26, 41-45, 47-51, 53-56, 58, 60, 64, 86? ], or EV is modelled by a time-varying carrying capacity that affects the species death rates and drives the population size [57,[61][62][63]66] (see also [65]).The distinctive feature of this study is therefore the twofold environmental variability accounting for fluctuations stemming from the variation of the toxin level and the switches of the carrying capacity resulting in the coupling of DN and EV; see Fig. 1.As main results, we obtain the fixation-coexistence diagrams of the system, and these allow us to determine the environmental conditions under which long-lived coexistence of the strains is possible or certain, and when one strain dominates the other.We also analyse the make-up of the population when the strains coexist, and their average abundance.
The organisation of the paper is as follows: the model is introduced in Sec.II.Sec.III is dedicated to the study of the case with a constant carrying capacity (subject to a static or varying toxin level) by means of a mean-field analysis and a mapping onto a suitable Moran process.The twofold influence of time-varying fitness and carrying capacity on the coexistence and fixation of the species is analysed in Sec.IV.Sec.V is dedicated to the influence of the EV on the make-up of the coexistence phase and strains abundance.We present our conclusions in Sec.VI.Additional technical details are given in the appendix.

II. MODEL
We consider a well-mixed population of fluctuating size N (t) = N R (t) + N S (t) which, at time t, consists of N R bacteria of strain R and N S of type S, which compete for the same resources.The former refers to a strain that pays a metabolic cost to be resistant to a certain toxin, and the latter to microorganisms sensitive to that toxin.Based on mounting evidence showing that microbial communities generally evolve in volatile environments [5,6,9,12,25], we study the ecoevolutionary dynamics of this population under twofold environmental variability: external conditions fluctuate between harsh and mild, and affect the level of toxin and resources that are available in the population; see Fig. 1.
For concreteness, we assume that the toxin is biostatic and reduces the growth rate of the sensitive strain, but does not affect the resistant bacteria [87][88][89] 1 .In this setting, resistant R bacteria have a constant fitness f R , FIG. 1. Cartoon of the model characterised by twofold EV.A microbial community consisting of two strains, denoted by R (resistant, blue (charcoal), fitness fR) and S (sensitive, red (grey), fitness fS), evolves in a time-varying environment, as illustrated by the arrows: the level of toxin, ξT , stochastically switches with rates ν ± T (vertical arrows) between low (ξT = +1) and high (ξT = −1), and the amount of available resources, modelled by the carrying capacity K(t), stochastically switches with rates ν ± K (horizontal arrows) between scarce and abundant.The strain R fares better than S under high toxin level (fR > fS), while S grows faster under low toxin level (fR < fS).The carrying capacity K(t) = K+ when there are abundant resources, while K(t) = K− < K+ when available nutrients are scarce.The four environmental states characterising the twofold EV are indicated by the coloured background: striped / solid green (light grey) refers to low toxin level and scarce / abundant resources, while striped / solid red (grey) indicates high toxin level with scarce / abundant resources.See the text for the greyscale coding, and the table in Sec.SM1 of the appendix for detailed notations and definitions.
whereas the sensitive S bacteria have an environmentdependent fitness f S (ξ T ), where ξ T (t) is a time-varying environmental random variable encoding the toxin level: ξ T > 0 for the low toxin level and ξ T < 0 for the high toxin level.As in previous studies [29,30,[91][92][93], we here consider where s > 0 denotes the selection bias favouring the strain S when ξ T > 0, and strain R when ξ T < 0. The parameter s therefore encodes both the selection and strength of the environmental variability associated with the changes in toxin level (T -EV).As in many recent theoretical studies [56, 57, 61-63, 65, 91, 94], T -EV is here modelled by coloured dichotomous Markov noise here.This is not particularly limiting since the same drug can often act as a biostatic or biocidal toxin at low/high concentration [90].
In this context, DMN, that is easy to simulate accurately and whose relationships with more complex forms of noise have been studied (see Sec. SM4 in the appedix and [66,95,96]), provides us with mathematicallyamenable models that can be regarded as the theoretical counterparts of commonly-used laboratory experimental chemostat set-ups [8,11,13,19].
The environmental effect on the level of nutrients (K-EV), fluctuating between scarcity and abundance, is modelled by a binary switching carrying capacity K(t) ∈ {K − , K + } that is driven by the binary random variable (also following a DMN process) ξ K (t) ∈ {−1, 1}.The state ξ K = −1 thus corresponds to a harsh state with scare resources, where the carrying capacity is K − , whereas nutrients are abundant in the mild state ξ K = +1 where the carrying capacity is K + > K − ≫ 1.As in [57, 61-63, 65, 66], this is encoded in the time-switching carrying capacity which, with K 0 ≡ K++K− 2 and γ ≡ K+−K− 2K0 , can conveniently be written as This randomly switching carrying capacity drives the population size N (t) and is hence responsible for its fluctuations, with K(t) ≫ 1 ensuring that the population dynamics is never solely dominated by demographic fluctuations.It is worth noting the choice of ξ K as a DMN ensures that the carrying capacity is always bounded and physical, i.e.K(t) ∈ {K − , K + } and K(t) > 0, which is not the case if EV is modelled by unbounded (e.g.Gaussian) noise, see Sec.SM4 in the appendix.
The population thus evolves under twofold EV encoded in the environmental states ξ T (t), ξ K (t), see Fig. 1, subject to time switching according to the reactions and where ν ± α are the switching rates of the α-DMN, with α ∈ {T, K} indicating the relevant environmental noise.It is also useful to define the average switching rate ν α and switching bias δ α for each α-DMN as . This means that δ T > 0 corresponds to a bias towards low toxin level (mild T state, ξ T = +1) favouring the S strain, whereas δ T < 0 indicates a bias towards high toxin level (harsh T state, ξ T = −1) where the growth of S is hampered and the spread of R is favoured, see Fig. 1.Similarly, δ K > 0 corresponds to bias towards the environmental state rich in nutrients (where K = K + ), while δ K < 0 is associated with a bias towards an environment where nutrients are scarce (K = K − ).In all cases, we consider α-DMN at stationarity, where ξ α = ±1 with probability (1 ± δ α )/2, yielding the average ⟨ξ α ⟩ = δ α and autocovariance (autocorrelation up to a constant) , where ⟨•⟩ denotes the α-DMN ensemble average.We notice that the correlation time of the α-DMN, 1/(2ν α ), is half of the inverse of the average switching rate of ξ α [95][96][97].From Eq. ( 1), we find that the average carrying capacity is ⟨K⟩ = K 0 (1 + γδ K ) and the variance of , with the amplitude of K-EV thus scaling as K 0 γ, while the variance and amplitude of the T -EV increase with s; see Sec.SM2 the appendix.
For concreteness, we here assume that ξ T and ξ K are totally uncorrelated.In our motivating example, this corresponds to the reasonable assumption that nutrient and antibiotic uptake are independent processes.The case where ξ T and ξ K are fully correlated or fully anticorrelated, with ξ T = ξ K = ξ or ξ T = −ξ K = ξ, where ξ is a single DMN process, is briefly discussed in Sec.SM8 of the appendix.
The system considered here best translates to a chemostat setup whereby toxin and nutrient levels can be maintained at a constant level through time and switched via changing concentrations of medium coming into the system [19,65].The switch ξ T → −ξ T with ξ T = −1 can thus be envisioned as corresponding to switching the concentration of an antibiotic drug from above the minimum inhibitory concentration (MIC), where the growth of the sensitive strain is hampered, to a concentration below the MIC where the S strain can spread at the expense of R [36,38,64].
At time t the fraction of R-types in the system is x(t) ≡ N R (t) N (t) and the average population fitness is f (x, ξ T ) ≡ x + (1 − x) exp(sξ T ), which depends on the population composition x and the toxin state ξ T , .We assume that mutation rates between strains are negligible, and seek to characterise the population dynamics by the evolution of its size and composition according to the multivariate birth-death process [98-100] where the time-dependent birth and death transition rates are respectively The per-capita birth rates f R/S /f (where we normalise with f in line with the standard Moran process [41,[67][68][69]101]) thus vary with the toxin level and population composition, while the logistic-like per-capita death rate N/K varies with nutrient level and population size.
With N ≡ (N R , N S ), the master equation giving the probability P (N, ξ T , ξ K , t) for the population to consist of N R and N S bacteria of type R and S, respectively, in the environmental state (ξ T , ξ K ) at time t is where are shift operators such that We note that P (N, ξ T , ξ K , t) = 0 whenever N R < 0 or N S < 0, and the last two lines on the right-hand-side of Eq. ( 4) account for the random environmental switching of toxin (ξ T → −ξ T ) and carrying capacity (ξ K → −ξ K ).The dynamics encoded by the multivariate master equation Eq. ( 4) can be simulated exactly by a stochastic algorithm as explained in [102,103], see Sec.SM4 in the appendix.Since T ± R/S = 0 whenever N R/S = 0, there is extinction of R and fixation of S (N R = 0, N = N S ), or fixation of R and extinction of S (N S = 0, N = N R ).When one strain fixates and replaces the other, the population composition no longer changes while its size continues to fluctuate 2 .Fixation of one strain and extinction of the other is expected when strains compete for the same resources (competitive exclusion principle), and always occur in a finite population even when its size fluctuates [57,61,63].In stark contrast, here we show that environmental fluctuations can lead to the long-lived coexistence of competing species and nontrivially shape the abundance distribution of both strains.
The notation and definitions of the model parameters and physical quantities of interest are conveniently summarised in the table of Sec.SM1 in the appendix. 2Finally, the population will settle in the absorbing state N R = N S = 0 corresponding to the eventual extinction of the entire population.This occurs after a time that grows exponentially with the system size [57,61,63,104].This phenomenon, irrelevant for our purposes (since we always have K(t) ≫ 1), is not considered here.

III. CONSTANT CARRYING CAPACITY: MEAN-FIELD ANALYSIS AND MORAN PROCESS
Since ξ T and ξ K are independent, it is useful to first consider the case of a constant carrying capacity, with environmental variability stemming only from the fluctuations of the toxin level in the birth rates of Eq. (3).
In this section, we thus assume that the carrying capacity is constant and large: K(t) = K 0 ≫ 1.After a short transient the population size fluctuates about K 0 , with N ≈ K 0 .When K 0 ≫ 1, we can approximate the population size by N = K 0 and make analytical progress by using the well-known results of the Moran process [41,[67][68][69]101].In this approximation, the population is kept constant, which requires the simultaneous birth and death of individuals of either species, and the population evolves according to a fitness-dependent Moran process [57,[61][62][63], defined in terms of Eqs. ( 2) and (3) by the reactions corresponding, respectively, to the simultaneous birth of an R and death of an S with rate T + R , and death of an R and birth of an S with rate T − R , where A. Mean-field analysis We now consider the case where N = K 0 → ∞, and thus ignore demographic fluctuations.In this case, the population composition evolves according to the meanfield equation [98]: where the dot denotes the time derivative.It is important to notice that, owing to environmental noise ξ T , Eq. ( 7) is a mean-field stochastic differential equation that defines a so-called "piecewise deterministic Markov process" (PDMP) [56,57,61,63,105].According to this PDMP, after a switch to an environmental state ξ T , x evolves deterministically with Eq. ( 7) and a fixed value of ξ T , until a switch ξ T → −ξ T occurs, see Sec.V A. We consider Eq. ( 7) in the regimes of (i) low, (ii) high, and (iii) intermediate switching rate ν T : (i) Under low switching rate, ν T → 0, the population settles in its final state without experiencing any Tswitches.In this regime, the population reaches its final state in its initial toxin level ξ T (0), i.e. ξ T (0) = ξ T (∞) = ±1 with probability (1 ± δ T )/2.In this regime, Eq. ( 7) thus boils down to ẋ = − x(1−x)(e s −1) with probability 1−δ T 2 .
Since s > 0, with a probability (1 + δ T )/2 we have ξ T (0) = ξ T (∞) = +1 and x → 0 (R vanishes), while with a probability (1 − δ T )/2 we have ξ T (0) = ξ T (∞) = −1 and x → 1 (S vanishes).In either case, the mean-field dynamics are characterised by the dominance of one of the strains.Therefore, in the absence of demographic fluctuations, there is never long-lived coexistence of the strains R and S under low switching rate of the toxin level.
(ii) Under high switching rate, ν T ≫ 1, the population experiences a large number of T -switches before relaxing into its final state; see below.In this case the T -DN selfaverages [29,53,57,61,63,66] and we are left with a Moran process defined by the effective rates obtained by averaging ξ T over its stationary distribution, yielding (8) When N → ∞, the mean-field (MF) rate equation associated with this effective Moran process thus reads [68,98,99]: where the right-hand-side (RHS) can be interpreted as the RHS of Eq. ( 7) averaged over ξ T .In addition to the trivial fixed points x = 0, 1, Eq. ( 9) admits a coexistence equilibrium when − tanh s 2 < δ T < tanh s 2 .This equilibrium stems from the T -DMN and thus is a fluctuation-induced coexistence point.In the case of large s we have that coth s 2 → 1, and x * exists (0 < x * < 1) for all values of δ T .Since ) < 0 is negative, linear stability analysis reveals that x * is the sole asymptotically stable equilibrium of Eq. ( 9) when it exists (x = 0, 1 are thus unstable).When s ≪ 1, coth s 2 → 2 s , and x * exists only for This means that for s ≪ 1, coexistence is essentially possible only under symmetric switching (δ T = 0), see Sec.SM2 the appendix.In what follows, we focus on the less restrictive case s = O(1), for which coexistence is possible for a broad range of parameters (ν T , δ T ).
(iii) In the regime of intermediate switching rate, where ν T ∼ 1, the population experiences a finite number of Tswitches prior to settling in its final state.Depending on this number, as well as the selection strength s and the population size, the dynamics may be closer to either the low or high ν T regime, with dominance or coexistence possible but, in general, not certain; see Fig. 3 below.

B. Finite populations -fixation and long-lived coexistence
From the MF analysis, we have found that when N → ∞ species coexistence is feasible under fast T -EV switching, whereas only R or S dominance occurs under slow switching.We now study how these results nontrivially morph when the population is fixed and finite.
Since the model is defined as a finite Markov chain with absorbing boundaries, see Eqs. ( 5) and ( 6), its final state unavoidably corresponds to the fixation of one strain and the extinction of the other, i.e. the population eventually ends up in either the state (N R , N S ) = (N, 0) or (N R , N S ) = (0, N ) [41,68,100,101].This means that, strictly, the finite population does not admit stable coexistence: when it exists, x * is metastable [106][107][108].
In fact, while it is guaranteed that eventually only one of the strains will finally survive, fixation can occur after a very long time and can follow a long-term coexistence of the strains, as suggested by the MF analysis of the regime with ν T ≫ 1.It is thus relevant to study under which circumstances there is long-lived coexistence of the strains.
The evolutionary dynamics is characterised by the fixation probability of the strain R, here denoted by ϕ.This is the probability that a population, consisting initially of a fraction x 0 of R bacteria, is eventually taken over by the strain R. A related quantity is the unconditional mean fixation time (MFT), here denoted by τ , which is the average time for the fixation of either species to occur.
In what follows, we study how the R fixation probability ϕ(ν T ) and the MFT τ (ν T ) vary with the average switching rate of the T -EV for different values of K 0 , δ T , and s (treated as parameters), and determine when there is long-lived coexistence of the strains.

appendix.
For a given initial resistant fraction3 , the fixation probability in the slow-switching regime, ϕ(ν T → 0) is obtained by averaging ϕ MA (N ) ξ T , denoting the R fixation probability in the realm of the MA in static environment ξ T , over the stationary distribution of ξ T [57,61,63,66]: When N ≫ 1 and ξ T = +1 the strain S is always favoured and ϕ MA (N ) ξ T =+1 ≈ 0, whereas R is favoured when ξ T = −1 and in this case ϕ MA (N ) The probability that S fixates when ν In the fast-switching regime the fixation probability is that of the Moran process defined by the effective rates in Eq. ( 8).Using Eq. ( 8) with x = n/N , we thus find [41,101]: see Sec.SM3.A the appendix.A similar analysis can also be carried out for τ , see Sec.SM3.B the appendix.
Results reported in Fig. 2 show that Eqs. ( 11) and ( 12) accurately capture the behaviour of ϕ in the limiting regimes ν T → 0, ∞, see Fig. 2(a).Fig. 2(b) shows that the predictions for τ when ν T → 0, ∞, are also in good agreement with simulation results, with a much larger MFT under high ν T than under low switching rate (at fixed δ T ).In Fig. 2(b), the MFT when ν T ≫ 1 for δ T = 0 is significantly larger than under δ T ̸ = 0.This stems from x * = x 0 = 1/2 being the attractor of Eq. ( 9) when δ T = 0, but not being an equilibrium when δ T = 0.3 or δ T = 0.5.Fig. 2(a,b) also illustrate the excellent agreement between the predictions of the MA with N = K 0 and those obtained from stochastic simulations with K = K 0 constant.Note that typical error bars are shown for δ T = 0 in Fig. 2(b).These are found to be small and almost coincide with the markers.
For the sake of readability, we have thus omitted similar error bars from the other panels and figures.
The MF analysis and results of Fig. 2 suggest that under sufficiently high switching rate ν T there is longlived coexistence of the strains.We can rationalise this picture by noting that in the regime of dominance of one strain the MFT scales sublinearly with the population size N , while the MFT grows superlinearly (exponentially when N = K 0 ≫ 1, see Fig. 2(c)) in the regime of long-lived coexistence [82,101,109,110].The dominance and long-lived coexistence scenarios are separated by a regime where the MFT scales with the population size, i.e. τ ∼ N , where the dynamics is governed by random fluctuations.This leads us to consider that long-lived coexistence of the R and S strains arises whenever the MFT exceeds 2⟨N ⟩, i.e. when τ > 2⟨N ⟩, where ⟨N ⟩ is the mean population size at (quasi-)stationarity; see below 4 .This is illustrated in the provided videos of [111] commented in Sec.SM7 of the appendix.When, as in this section, N = K 0 or N fluctuates about the constant carrying capacity K 0 (N ≈ K 0 ), we simply have ⟨N ⟩ = K 0 .The criterion τ > 2⟨N ⟩ = 2K 0 thus prescribes that long-lived coexistence occurs when the MFT scales superlinearly with K 0 and hence exceeds the double of the average population size, 2K 0 ≫ 1.
The conditions under which the long-lived coexistence criterion, τ > 2⟨N ⟩, is satisfied can be estimated by noting that, from the MF analysis, we expect coexistence to occur when ξ T self-averages under sufficiently high switching rate ν T .Since the average number of Tswitches by t = 2⟨N ⟩ scales as ν T ⟨N ⟩, self-averaging occurs when ν T ⟨N ⟩ ≫ 1.We thus consider that there is fast T -EV switching when ν T ≫ 1/⟨N ⟩, while ν T ≪ 1/⟨N ⟩ corresponds to the slow T -EV regime.To ensure long-lived coexistence, the necessary condition ν T ≫ 1/⟨N ⟩ is supplemented by the requirement that s ∼ 1.This ensures enough environmental variability and a regime of coexistence where the MFT is generally τ ∼ e c⟨N ⟩ (where c is some positive constant) when s = O(1) [41,42,101,[107][108][109]112] guaranteeing τ > 2⟨N ⟩.
Hence, the expected conditions for long-lived coexistence are ν T ≫ 1/⟨N ⟩ (fast T -switching) and s = O(1) (enough EV), which are satisfied in the examples considered here when ν T ∼ 1, s ∼ 1 or greater, and ⟨N ⟩ ≫ 1.
We have studied the influence of T -EV on the fixation and coexistence properties of the model with constant carrying capacity K = K 0 and selection strength s by running a large number of computer simulations up to a time t = 2K 0 across the ν T − δ T parameter space.When just after t = 2K 0 both species are present, the run for (ν T , δ T , s) is characterised by longlived coexistence which is RGB coded (0, 1, 0).There is no long-lived coexistence for the run (ν T , δ T ) if one of the species fixates by t ≤ 2K 0 : either the strain R, which is RGB coded (1, 0, 0), or the strain S, which is coded by (0, 0, 1).This procedure is repeated for 10 3 realisations for different (ν T , δ T ) and, after sample averaging, yields the RGB-diagram of Fig. 3(a-c); see Sec.SM4 the appendix.In greyscale, the RGB coding translates into red → grey, blue → charcoal, and green → light grey.In what follows, the crossover regimes are coded in magenta and faint green with dark grey and faint light grey as their respective greyscale counterparts, see below.
It is also useful to study the effect of the T -EV in the realm of the MA by means of numerically exact results.For this, with the transition rates Eq. ( 6), we notice that when N = K 0 is constant and there are initially n cells of type R, the R-strain fixation probability, ϕ ξ T n , in the environmental state ξ T satisfies the first-step analysis equation [53,98,100] ( subject to the boundary conditions ϕ ξ T 0 = 0 and ϕ ξ T N = 1.The mean fixation time in the environmental state ξ T , τ ξ T n , satisfies a similar equation: with boundary conditions τ ξ T 0 = τ ξ T N = 0. Eqs. ( 13) and ( 14) are thus solved numerically, and the fixation probability and MFT are obtained after averaging over the stationary distribution of ξ T , yielding In our examples, we always consider x 0 = 1/2, and henceforth write ϕ MA (N ) ≡ ϕ N/2 for the R-fixation probability and τ MA (N ) ≡ τ N/2 for the MFT in the realm of the MA.For each triple (ν T , δ T , s), we numerically solved Eq. ( 14) and, in the region of the the parameter space where τ MA (N ) > 2K 0 , there is longlived coexistence, which is coded by (0, 1, 0) in the RGBdiagram of Fig. 3(d-f).When τ MA (N ) ≤ 2K 0 , there is dominance of one of the species, characterised by the fixation probabilities ϕ MA (N ) and 1 − ϕ MA (N ) of R and S, respectively, obtained from Eq. ( 13) and coded by (ϕ MA (N ), 0, 1 − ϕ MA (N )) in Fig. 3(d-f).Exact numerical results for the MA with N = K 0 in Fig. 3(d-f) are in excellent agreement with those of simulations obtained for K = K 0 in Fig. 3(a-c).In line with the MF analysis, we find that long-lived coexistence, occurs for T -EV of sufficiently large magnitude, i.e. s ∼ 1 or higher, and under high enough switching rate, i.e. ν T ∼ 1 or higher, shown as green (light grey) areas in Fig. 3.The region of coexistence separates regimes dominated by either species, especially at high ν T when ϕ ≈ 0 when δ T 0 while ϕ ≈ 1 when δ T < 0. In Fig. 3, the boundaries between the regimes of R/S dominance, coded in blue (charcoal) / red (grey), and coexistence, areas in green (light grey)), are interspersed by crossover regimes where both species are likely to fixate (magenta (dark grey) in Fig. 3), or coexist with probability between 0 and 1 (faint green (faint light grey) in Fig. 3), as coded in Fig. 3(g).

IV. TWOFOLD ENVIRONMENTAL VARIABILITY: COEXISTENCE AND FIXATION UNDER TIME-VARYING FITNESS AND SWITCHING CARRYING CAPACITY
We have seen that under a constant carrying capacity, long-lived coexistence of the strains is feasible when s and ν T are of order 1 or higher (enough T -EV and fast Tswitching).We now study how this picture morphs when, in addition to the time-variation of f S and f , the carrying capacity K(t) switches according to Eq. ( 1) and drives the fluctuating population size N .EV is thus twofold, and the population evolves under the joint effect of T -EV and K-EV.
We consider where N is independent of s and affected only by K-EV via ξ K in Eq. ( 1), while the evolution of x in Eq. ( 16b) is impacted by ξ K , ξ T , and The population composition is hence coupled to the evolution of the population size (eco-evolutionary dynamics), while the statistics of N , like its average, denoted by ⟨N ⟩, are obtained by ensemble averaging over ξ K .The stochastic logistic differential equation Eq. (16a) defines an N -PDMP whose properties allow us to characterise the distribution of N [56,57,61,105].

A. N -PDMP approximation
The PDF p(N ) captures well the main effects of the K-EV on the QPSD, which is bimodal under low ν K and becomes unimodal when ν K ≫ 1, see Fig. 4. In the realm of the N -PDMP approximation, p(N ) aptly reproduces the location of the QPSD peaks and the transition from a bimodal to unimodal distribution as ν K increases: The distribution is sharply peaked around N ≈ K ± when ν K → 0 (with probability (1∓δ K )/2), flattens when ν K ∼ 1, and then sharpens about Since p(N ) ignores DN, it cannot capture the width of the QPSD about the peaks, but it provides an accurate description of the mean population size, see Fig. 4, that is well approximated by with ⟨N ⟩ ≈ K 0 (1 + γδ K ) when ν K ≪ 1 and ⟨N ⟩ ≈ K when ν K ≫ 1 [57,[61][62][63], as shown in Fig. 4.
The PDF p(N ) is particularly useful to obtain the fixation/coexistence diagrams under the effects of both T -EV and K-EV.Theoretical/numerical predictions of the fixation and coexistence probabilities can indeed be derived in the vein of [57,[61][62][63]66] by focusing on situations where coexistence occurs when x and N relax on similar timescales.Long-lived coexistence of the strains thus arises when s ∼ 1, with fixation typically occurring after N has settled in the QPSD, see Sec.SM7 the appendix and videos in [111].Hence, for the analytial description of the fixation/coexistence diagrams the R fixation probability (with x 0 = 1/2) can be suitably approximated by averaging ϕ MA (N ) over p(N ) as follows: where ϕ MA (N ) is obtained from solving the corresponding Eq. ( 13) for the R fixation probability of the associated Moran process, as seen in Sec.III B. We can use the PDF p(N ) and the results for the mean fixation time τ MA (N ), obtained from solving Eq. ( 14), to determine the probability of coexistence in the realm of the N -PDMP approximation.For this, we first solve Eq.( 14) for τ MA (N * ) = 2⟨N ⟩, where ⟨N ⟩ is given by Eq. ( 18).Since τ MA is an increasing function of N , see Fig. 2(c), we have τ MA (N ) > 2⟨N ⟩ for all N > N * , which is the long-lived coexistence condition.Within the N -PDMP approximation, the lowest possible value of ).We then determine the probability η that this condition is satisfied by integrating p(N ) over [max(N * , K − ), K + ]: (20) where N * depends on both T -EV and K-EV (via ⟨N ⟩), while the integrand depends only on K-EV.Clearly, η → 1 when N * → K − .Hence, long-lived coexistence is almost certain when N * ≈ K − , i.e. whenever the meanfixation time of the population of fixed size N = K − exceeds 2 ⟨N ⟩.Based on the results of Sec.III B, η increases with ν T and s, and thus for sufficiently large ν T and s (but not too large δ T ), we expect N * → K − and η → 1.

B. Fixation/coexistence diagrams under T -EV and K-EV
The fixation/coexistence diagrams under joint effect of T -EV and K-EV are obtained as in Sec.III B, with the difference that long-lived coexistence arises when t > 2 ⟨N ⟩, a condition that depends on (ν K , δ K ), see Fig. 2(c).
In our simulations, we considered different values of ν K (letting δ K = 0 for simplicity), and ran simulations until t = 2 ⟨N ⟩.Each run in which both species still coexist just after t = 2 ⟨N ⟩  19) and (20), see text for details.All diagrams are colour-coded (greyscale-coded) as in Fig. 3.
The comparison of the top and bottom rows of Fig. 5 shows that the theoretical RGB (greyscale) diagrams quantitatively reproduce the features of those obtained from simulations.In general, we find that coexistence regions are brighter in diagrams obtained from N -PDMP based approximation than in those stemming from simulations.This difference stems from former ignoring demographic fluctuations which slightly broaden the crossover (magenta (dark grey) and faint green (faint light grey)) regimes in the latter.
The regions of Fig. 5 where |δ T | → 1 are characterised by dominance of one of the strains, and essentially reduces to the model studied in [57,61,63], and we can therefore focus on characterising the coexistence phase.
When K 0 is large, under sufficient environmental variability (s = 0.5, γ = 2/3 in Fig. 5), the joint effect of T -EV and K-EV on the phase of long-lived coexistence in the RGB (greyscale) diagrams of Fig. 5 can be summarised as follows: (i) when ν K → 0, a (bright green) region where η ≈ 1 and coexistence is almost certain is surrounded by a faint green (faint light grey)"outer shell" where coexistence is possible but not certain (0 < η < 1), see Fig. 5(a,d); (ii) at low, but nonvanishingly small, values of ν K , the outer-shell where 0 < η < 1 fades, and there is essentially only a bright green (bright light grey) region of coexistence where η ≈ 1, see Fig. 5(b,e); (iii) when ν K ≫ 1, the coexistence region corresponds essentially to η ≈ 1 (bright green / bright light grey) and is broader than under low ν K , Fig. 5(c,f).In all scenarios (i)-(iii), η increases with ν T ≳ 1 (for not too large δ T ) and hence all the green (bright and faint light grey) coexistence phases in Fig. 5) become brighter as ν T is raised and η → 1.
These different scenarios can be explained by the dependence of the QPSD on ν K , well captured by the PDF Eq. (17).In regime (i) where ν K ≪ 1/K 0 , the QPSD and p(N ) are bimodal, N ≈ K ± with probability (1 ± δ K )/2, and any K-switches by t = 2 ⟨N ⟩ ∼ K 0 are unlikely, yielding the faint green (faint light grey) outer shell of Fig. 5(a,d) corresponding to longlived coexistence arising only when N ≈ K + , with a probability η ≈ (1 + δ K )/2.In regime (ii), where 1/K 0 ≪ ν K ≪ 1, the QPSD and p(N ) are still bimodal but some K-switches occur by t ∼ K 0 , resulting in longlived coexistence arising almost only when ν T is high enough to ensure η ≈ 1 when N ≈ K − .In regime (iii), where ν K ≫ 1 the QPSD and p(N ) are unimodal with average ⟨N ⟩ ≈ K ≥ K − , which results in a long-lived coexistence region where η ≈ 1 that is broader than in (i) and (ii), Fig. 5(c,f).The size of the coexistence region in regime (iii) actually depends nontrivially on ν K , as revealed by the modal value of the PDF Eq. ( 17) when with lim ν K →∞ N = ⟨N ⟩ = K.We notice that N is an increasing function of ν K when γ > δ K , and it decreases if γ < δ K (remaining constant when γ = δ K ).As a consequence, the long-lived coexistence region under high K switching rate grows with ν K when γ > δ K , as in Fig. 5(c,f), and, when γ < δ K , shrinks as ν K is increased, see Sec.SM6 Fig. (S3) the appendix.

C. Influence of the K-EV amplitude on coexistence
We have seen that increasing the selection bias s, raises the amplitude of the T -EV and facilitates the emergence of long-lived coexistence.Here, by keeping K 0 constant, we investigate the influence of the parameter γ, which controls the amplitude of K-EV, on the fixation/coexistence diagrams.When γ → 1 and K 0 ≫ 1, there is K-EV of large amplitude, with the population subject to a harsh population bottleneck (K − → 0) accompanied by strong demographic fluctuations.The latter facilitate fixation of either strain and counter the effect of T -EV that drives the community to coexistence.Results of Fig. 6 illustrate the influence of γ under low and high K-switching rate (δ K = 0): -Under low ν K , the probability of long-lived coexistence η decreases together with the value of K − = K 0 (1 − γ) when γ is increased (all other parameters being kept fixed).As a consequence, the bright green (bright light grey) region in Fig. 6(a) where long-lived coexistence is almost certain (η ≈ 1) shrinks with γ and is gradually replaced by a faint green (faint light grey) area where coexistence occurs with a lower probability (η = (1 + δ K )/2 < 1), see Fig. 6(b,c).
-Under high ν K , we have N ≈ K and the effect of γ is encoded in the expression of K = K 0 (1 − γ 2 )/(1 − γδ K ).When δ K ≤ 0, K and η decrease with γ, and as a result the bright green (bright light grey) region in Fig. 6(d) shrinks and is eventually replaced by a smaller faint green (faint light grey) region where coexistence is possible but not certain (0 < η < 1), see Fig. 6(e,f).When δ K > 0, there is a bias towards K = K + and K increases with γ until γ = γ ≡ (1 − 1 − δ 2 K )/δ K and then decreases, with K < K 0 , when γ > δ K .This results in a non-monotonic dependence of the coexistence region where η ≈ 1: under ν K ≫ 1 and δ K > 0, the long-lived (bright-green) coexistence region grows with γ up to γ and shrinks when γ > γ.
We have thus found that the environmental fluctuations have opposite effects on the species coexistence: increasing the amplitude of T -EV (by raising s) prolongs the coexistence of the strains and expands the coexistence region, but raising the amplitude of K-EV (by raising γ) can significantly reduce the probability of long-lived coexistence for all values of ν K .

V. MAKE-UP OF THE COEXISTENCE PHASE AND STRAINS AVERAGE ABUNDANCE
Having characterised in detail the conditions under which long-lived coexistence and fixation occur, we now study the make-up of the coexistence phase and then use this result to determine the stationary average abundance of each strain.

A. Coexistence phase make-up
We are interested in the characteristic fraction of the resistant strain R in the coexistence phase, here defined as x * .This is the fraction of R expected, given that we have coexistence at t = 2 ⟨N ⟩.According to the mean-field theory, the fraction of the strain R in the coexistence phase is given by the expression Eq. ( 10) of x * .It turns out that deep into the coexistence region whereby η ≈ 1 and ν T is sufficiently high, there is good agreement between theory and simulations, see Fig. 7(a,b).In addition, even when η < 1, the mean-field prediction x * remains remarkably close to the value of fraction of R measured in the coexistence state obtained in simulations, with small deviations arising as η approaches 0. We notice that the characteristic fraction of R, for given δ T , is almost independent of ν T .
We can also predict the fraction of R regardless of coexistence or fixation, here denoted by ⟨x⟩.The quantity ⟨x⟩ thus characterises the fraction of R in the coexistence, fixation, and crossover regime where both coexistence and fixation are possible, with respective probabilities η and ϕ, but neither is certain.Making use of Eqs. ( 10), (19), and (20) we thus define ⟨x⟩ as This captures well the dependence of ⟨x⟩ on ν T and reduces to the fraction of R in the coexistence phase, ⟨x⟩ = x * , when η ≈ 1 and long-lived coexistence is almost certain (see SM5, Fig. S2).As shown in Sec.SM5 of the appendix, a closed-form alternative to ⟨x⟩ is provided by the modal value of the stationary PDF of the PDMP defined by Eq. ( 7), which, while less accurate than ⟨x⟩, matches qualitatively well to simulations.

B. Strain average abundance
In this section, we study the (quasi-)stationary average abundance of the strains R and S, respectively denoted by ⟨N R ⟩ and ⟨N S ⟩.Since ⟨N S ⟩ = ⟨N ⟩ − ⟨N R ⟩, and ⟨N ⟩ is well described by Eq. ( 18), see Fig. 4, we only need to focus on studying ⟨N R ⟩.
In fact, while the evolution of N is governed by K-EV and is well-captured by the stochastic logistic equation Eq. (16a) and the corresponding N -PDMP, the dynamics of the abundance of the R strain depends on both T -EV and K-EV.In the mean-field limit, where demographic fluctuations are neglected, we indeed have [98] which is a stochastic differential equation depending on both ξ K and ξ T , and coupled to the N -and x-PDMPs defined respectively by Eqs.(16a) and (16b).
In the dominance regimes, ⟨N R ⟩ ≈ 0 (S dominance) or ⟨N R ⟩ ≈ ⟨N ⟩ (R dominance), which can be obtained from Eq. ( 18).However, finding ⟨N R ⟩ in the coexistence phase is a nontrivial task.Progress can be made by noticing that, ξ K and ξ T being independent, we can write where ⟨N ⟩ ηx * is the contribution to ⟨N R ⟩ when there is coexistence (with probability η), and ⟨N ⟩ (1 − η)ϕ is the contribution arising when there is fixation of the strain R (with probability (1 − η)ϕ).In our theoretical analysis, x * , ⟨N ⟩ , ϕ and η are obtained from Eqs. ( 10) and ( 18)- (20).Eq. ( 23) thus captures the behaviour of ⟨N R ⟩ in each regime: the dominance regime where η ≈ 0 and we have ⟨N R ⟩ ≈ ⟨N ⟩ ϕ, deep in the coexistence phase where we have η ≈ 1 and ⟨N R ⟩ ≈ ⟨N ⟩ x * , and where 0 < η < 1 and coexistence is possible but not certain where we have ⟨N R ⟩ ≈ ⟨N ⟩ ⟨x⟩.
In Fig. 8, we find that the theoretical predictions based on Eq. ( 23) agree well with simulation results over a broad range of ν K and ν T , and for different values of δ K and δ T .The dependence of ⟨N R ⟩ on ν K reflects that of ⟨N ⟩ shown in Fig. 4: ⟨N R ⟩ decreases with ν K at fixed δ K , see Fig. 8(a), and we have ⟨N ⟩ ≈ K when ν K → ∞ yielding ⟨N R ⟩ ≈ Kx * deep in the coexistence phase where ν T ≫ 1, and similarly ⟨N ⟩ ≈ K 0 (1 + γδ K ) when ν K → 0 yields ⟨N R ⟩ ≈ K 0 (1 + γδ K )x * .Not shown in Fig. 8(a) is the case of ν T ≪ 1, whereby we have only dominance such that ⟨N R ⟩ ≈ ⟨N ⟩ (1 − δ T )/2.Fig. 8(b) shows that the dependence of ⟨N R ⟩ on ν T can be non-monotonic and exhibit an extreme dip (δ T > 0) or peak (δ T < 0) at intermediate T -switching rate, ν T ∼ 1.This behaviour can be understood by referring to the diagrams of Fig. 6: as ν T is raised from ν T = 0 with δ T < 0 kept fixed, the R fixation probability first slowly increases across the slightly R-dominant phase where coexistence is unlikely (η ≈ 0) and ⟨N R ⟩ ≈ ⟨N ⟩ ϕ.When ν T is increased further and R is the strongly dominant species (blue (charcoal) phases in Fig. 6), with ϕ ≈ 1 and ⟨N R ⟩ ≈ ⟨N ⟩ is maximal; coexistence then becomes first possible (0 < η < 1, faint green (faint light gray) in Fig. 6) and then almost certain (η ≈ 1, bright green (bright light gray) in Fig. 6) when ν T is increased further, which results in a reduction of the R abundance to ⟨N R ⟩ ≈ ⟨N ⟩ x * < ⟨N ⟩.A similar reasoning holds for the S strain when δ T > 0 and results in a maximal value ⟨N S ⟩ ≈ ⟨N ⟩ and therefore a dip of the R abundance, with a minimal value ⟨N R ⟩ ≈ 0, when ν T ∼ 1.
The results of this section hence show that the twofold EV has nontrivial effects on the make-up of the coexistence phase, and on the average number of cells of each strain, as shown by Fig. 7 and the nonmonotonic dependence of ⟨N R ⟩ on ν T in Fig. 8.

VI. CONCLUSION
Microorganisms live in environments that unavoidably fluctuate between mild and harsh conditions.Environmental variability can cause endless changes in the concentration of toxins and amount of available nutrients, and thus shapes the eco-evolutionary properties of microbial communities including the ability of species to coexist.Understanding under which circumstances various microbial species can coexist, and how their coexistence and abundance vary with environmental factors, is crucial to shed further light on the mechanisms promoting biodiversity in ecosystems and to elucidate the evolution of antimicrobial resistance (AMR).
Motivated by these considerations, and inspired by the antimicrobial resistance evolution in a chemostat setup, we have studied the eco-evolutionary dynamics of an idealised microbial community of fluctuating size consisting of two strains competing for the same resources under twofold environmental variability (T -EV and K-EV): the level of toxin and the abundance of nutrients in the community both vary in time.One of the strains is resistant while the other is sensitive to the drug present in the community, and both compete for the same resources.
Environmental variability is thus assumed to affect the strains growth and death rates, and is modelled by means of binary randomly time-switching fitness (T -EV) and carrying capacity (K-EV).Under harsh conditions, the level of toxin is high and resources are scarce, while environmental conditions are mild when the level of toxin is low and resources are abundant.In this setting, the strain resistant to the drug has a selective advantage under high toxin-level, whereas it is outgrown by the sensitive strain when the level of toxin is low.Moreover, the time-switching carrying capacity drives the fluctuating size of the microbial community, which in turn modulates the amplitude of the demographic fluctuations, resulting in their coupling with the variation of the available resources.When the environment is static, there is no lasting coexistence since one species dominates and rapidly fixates the entire population.Here, we have shown that this picture changes radically in fluctuating environments: we have indeed found that long-lived species coexistence is possible in the presence of environmental fluctuations.Using stochastic simulations and the properties of suitable piecewise-deterministic and Moran processes, we have computationally and analytically obtained the fixationcoexistence phase diagrams of the system.These have allowed us to identify the detailed environmental conditions under which species coexist almost certainly for extended period of times, and the phases where one species dominates, as well as the crossover regimes where both coexistence and fixation are possible but not guaranteed.
We have found that long-lived coexistence requires sufficient variation of the toxin level, while resource variability can oppose coexistence when strong K-EV leads to population bottlenecks responsible for large demographic fluctuations facilitating fixation.More generally, our analysis has allowed us to assess the influence of the population size distribution, whose shape changes greatly with the rate of K-EV, on the fixation-coexistence phase diagram.We have also determined how the make-up of the coexistence phase and average abundance of each strain depend on the rates of environmental change.
Environmental variability generally comes about in many forms in a variety of settings throughout biology and ecology, and the conundrum of coexistence within a system is impacted by it, alongside demographic fluctuations.This leads to complex eco-evolutionary dynamics.In particular, how microbial communities evolve subject to environmental variability is vital when considering the issue of AMR, so that the effectiveness of treatments can be maximised, while minimising the harmful effects.In considering twofold environmental variations, we have shown that these can have qualitative effects on the population evolution as they can either promote or jeopardise lasting species coexistence.
In summary, our analysis allows us to understand under which circumstances environmental variability, together with demographic fluctuations, favours or hinders the long-lived coexistence of competing species, and how it affects the fraction and abundance of each strain in the community.This work hence contributes to further elucidate the role of fluctuations on the maintenance of biodiversity in complex ecosystems.
In particular, our findings demonstrate the influence of environmental fluctuations on biodiversity in microbial communities, and may thus have potential impacts on numerous applications.For instance, the model studied here is well suited to describe the in vitro evolution of antimicrobial resistance in a chemostat setup where the level of antibiotics would fluctuate below and above the minimum inhibitory concentration.In this context, the model is able to predict, under a broad range of external constraints, the best conditions to avoid the fixation of the strain resistant to the drug and when both strains coexist.A more realistic model of AMR evolution would take into account that the drug resistance is often mediated by a form of public goods [36,79], and that there may exist more than two competing species and various toxins.The eco-evolutionary dynamics of communities consisting of multiple species, resources and toxins can generally not be simply inferred from those of two-species eco-systems, even though in some cases simple models can be illuminating [65].Another potential application of the model considered here, with varying drug levels, concerns the so-called adaptive therapy used in cancer treatment to prevent or delay the cancer from becoming completely drug resistant [113].

APPENDIX
In this appendix, we provide some further technical details and supplementary information in support of the results discussed in the main text.We also provide additional information concerning the mean-field, Moran and PDMP approximations used in the main text, the simulation methods, we illustrate our main findings by discussing typical sample paths, and briefly discuss the generalisation of the model to correlated/anticorrelated EV.

SM1. TABLE OF PARAMETERS AND PHYSICAL QUANTITIES OF INTEREST
The various parameters of the model and physical quantities of interest are summarised in the following table.Scaled difference in carrying capacity values / strength of nutrient environmental variability γ ≡ Average switching rate and half of the inverse correlation time of ξα, with α ∈ {T, K} να Switching rate from + into − (+) or from − into + (−) of the environmental noise ξα, with α ∈ {T, K} ν ± α Switching bias of the environmental noise ξα, with α ∈ {T, K} δα ≡ Fraction of resistant bacteria out of total population x Average fitness of current population f Birth (+) and death (−) rates for bacteria of type i ∈ {R, S} T ± i Birth (+) and death (−) rates for the Moran approximation The selection strength s does not only shape the dynamics of the composition, see Eq. ( 7), but it also determines the amplitude of the T -EV fluctuations, i.e. its variance, which is linked to stronger coexistence in the fast-switching regime (see Sec. III).To show this, we first consider the normalised fitness of the resistant subpopulation .
Since we here focus on how s may shape coexistence through the toxin level fluctuation size, and coexistence dominates in the fast toxin-switching regime, we assume ν T ≫ 1.As discussed in the main text, in the fast-switching regime, the per capita growth rate (normalised fitness) of strain R averaged over the stationary distribution of ξ T reads: .
To derive its variance, we also have to compute the average of its square as Combining both, we obtain the variance of the normalised resistant fitness due to the environmental fluctuations in the toxin level A similar analysis for the sensitive strain, which has a normalised fitness exp(ξ T s) times that of the resistant strain, provides In both cases, we conclude that the variance (arising from the T -EV) indeed increases with the selection strength s.Therefore, s does shape the strength of coexistence.
Note that the amplitude of the fluctuations of the carrying capacity K-EV (at no bias δ K = 0, for simplicity) increases with γ.However, in this case, the larger the K-EV the weaker the coexistence, as the harsh state K − → 0 further promotes extinction.In conclusion, and as discussed in Sec.IV B and Fig. 6, the environmental variability can either promote coexistence (T -EV) or jeopardise it (K-EV), and both parameters s and γ determine long-lived coexistence.
To discuss the coexistence timescale, we consider the small s regime by expanding the numerator and denominator of the right-hand-side of Eq. ( 9) to order O(s 2 ), which yields ẋ ≈ x(1 − x) 2 where x * = 1 2 − δ T s is the coexistence equilibrium under s ≪ 1.The equilibrium x * is physical (and stable) only when − s 2 < δ T < s 2 ; i.e.only in the special case of almost symmetric switching (δ T = O(s) or smaller).From Eq. (S1), the coexistence equilibrium is approached slowly, with a relaxation time on the order of ∼ 1/s 2 ; thus, taking into account demographic noise, the expected fixation time is τ ∼ e ⟨N ⟩s 2 when ⟨N ⟩ s 2 ≫ 1.

A. Moran fixation probability
The exact Moran fixation probability ϕ N 0 R , N, s, ξ T for the resistant strain to take over the entire population of constant size N in the static toxin environment ξ T , starting with N 0 R resistant individuals, is [41,98,99,101] and N 0 R = 1, 2, ..., N. (S2) , and we take into account all the possible reactions that can take place.In the case of the full model this means: (1) the four possible birth or death reactions with rates {T + R (t), T − R (t), T + S (t), T − S (t)} that depend on the variables {N R (t), N S (t), K(t), ξ T (t)} and the constant parameter s; (2) the nutrient level switches stochastically with constant rate ν ± K for the state K(t) = K ± ; and (3) the toxin level switches stochastically with constant rate ν ± T for the state ξ T (t) = ±1.We perform efficient stochastic simulations by implementing the Next Reaction Method [102,103].
We have determined the standard errors for the data sets shown in Figs. 2, 4, and 8, in each case on the sample mean over 10 3 realisations.In the main text, we show the standard error on the mean for the case δ T = 0.0 in Fig. 2(b) where we expect the highest variance in simulation result.This illustrates that typical error bars in simulation results reported in Figs. 2, 4, and 8 do not ever exceed the size of the markers.Hence to maintain readability of the figures, and without loss of information, we have omitted error bars in the figures others than Fig. 2(b) for δ T = 0.0.
Regarding Fig. 5, the direct numerical implementation of Eqs.(19) and (20) to predict fixation and coexistence probabilities is not feasible, because calculating ϕ(N ) and the coexistence probability for each integer N ∈ [K − = 400, K + = 2000] is computationally too expensive.Therefore, we capitalise on the exact Moran results for static environments (see Eqs. (S3) and (S4)), and the N -PDMP PDF p(N ) (see Eq. ( 17) and Fig. 4), to provide the theoretical prediction for the fixation and coexistence probabilities.For the first regime at very low ν K (Fig. 5(a,d)), the starting environment is unlikely to switch, and the distribution of the total population is bimodal; see left inset in Fig. 4. The fixation and coexistence probabilities of Fig. 5(d) are then computed as the weighted average of the Moran probabilities for the static environment cases N = K ± , with weights (1 ± δ K ) /2.For the regime of Fig. 5(b,e), irrespective of the starting environment, the population size will, on average, spend significant time on both N = K − and K + .Since relative demographic fluctuations are larger at K − , the simulated fixation and coexistence probabilities of Fig. 5(b) are captured in Fig. 5(e) by a Moran process under total population N = K − .Finally, for high ν K , the environmental noise averages out and the behaviour of the system corresponds to that of a Moran process at N = K; see right inset in Fig. 4.
An additional comment on environmental noise is in order.As explained in the main text, we have chosen to model environmental variability using dichotomous Markov noise (DMN).This choice allows us to make mathematical progress while keeping the theoretical modelling close to laboratory experimental conditions.In fact, the properties of DMN and their relations with other forms of noises have been extensively studied [66,[95][96][97].The advantage of modelling environmental variability with DMN is well illustrated by considering the time-varying carrying capacity where γ ≡ (K + − K − )/(2K 0 ) and 0 < K 0 ≡ (K + + K − )/2 (with 0 < γ < 1 and 0 < K 0 < ∞).When ξ K (t) ∈ {−1, 1} is a dichotomous noise, K(t) ∈ [K − , K + ] is bounded and always physical, and the evolution of the community is fully characterised by the master equation of Eq. ( 6) that can be simulated exactly by the methods described above.However, if ξ K (t) was an unbounded noise of continuous range (−∞, ∞), e.g. a Gaussian white noise (zero correlation time) or an Ornstein-Uhlenbeck noise (finite correlation time), we would face a number of challenging problems: (i) The carrying capacity would no longer be bounded, and could take unphysical (negative) values (when ξ K < −1/γ).(ii) The evolutionary process resulting from the coupling of the birth-death process defined by Eq. ( 4) with the dynamics of K(t) driven by an (unbounded or bounded) environmental noise with continuous range would generally result in a non-Markovian process [66].There would therefore be no corresponding master equation governing the population dynamics.The mathematical analysis would thus be difficult, and even simulating the process would be a challenging task since the above methods could not be used directly.Circumventing these issues generally requires some approximations (e.g.truncating a Gaussian noise, finding a way to restore the Markov property) whose accuracy, validity and implications may be difficult to control, see e.g.[51,55,66] and in particular the discussion in the Supplemental Material of [66].

SM6. COEXISTENCE UNDER K-EV AND FINAL FIXATION
In Fig. 5(a-c) we find that the region of coexistence grows when ν K is increased.This behaviour is expected for an effective population size that would increase with ν K , as suggested by the MFT that increase with the population size (Fig. 2(c)).However, this seems at odds the average population size ⟨N ⟩ decreasing with ν K as shown by Fig. 4. A more suitable characterisation of the influence of ν K on coexistence phase is thus provided by the modal value N of the N -PDMP PDF, given by Eq. (21).Indeed, Fig. S3 illustrates that N , unlike ⟨N ⟩, increases ν K for δ K < γ in line with the results reported in Fig. 5(a-c).
For a complete characterisation of the coexistence regime we now briefly discuss the final state that is attained after after a long transient ensuring that a large fluctuation drive the system from the long-lived metastable coexistence to one of the two absorbing states where a single strain takes over the entire population [107,108].Fig. S4 illustrates the full fixation outcome diagram (under no K-EV for simplicity), even in the phase of long-lived coexistence.Fixation in this regime is determined by the sign of δ T , with a sharp transition at δ T = 0.This is because fixation occurs FIG.S3.Modal value N of the PDF Eq. ( 17) against νK using Eq. ( 21) (solid line) and points from simulations averaged over 10 2 realisations (×) obtained by tracking the modal value of the QPSD.Parameters used are K0 = 1000, γ = 0.5, δK = 0.0, νT = 10.0, δT = 0.0.As discussed in Sec.IV B, we find that for δK < γ there is an increase of N with νK up to K = 750 (fast-switching limit), thus allowing for a greater range of δT values that give coexistence.most likely in the absorbing state closest to the coexistence equilibrium x * given by Eq. ( 10): here, the toxin bias eventually imposes the fixation of S (δ T > 0) or R (δ T < 0).

SM7. VIDEOS OF SAMPLE PATHS
In this section we provide example trajectories for representative realisations of the full model under both fluctuations in the toxin level and in the carrying capacity.The corresponding movies are provided online in [111].
Fig. S5 captures the dynamics of the system in six example realisations under different EV parameters (ν T , δ T , ν K ) for a selection strength s = 0.5, which sets the magnitude of T -EV; and mean carrying capacity K 0 = 1200 with γ = 2/3 (K − = 400 and K + = 2000), setting the magnitude of K-EV.For simplicity, we keep δ K = 0.
-The example of Fig. S5(a) illustrates an unbiased fluctuating environment with intermediate T -switching rate (white background: low toxin level, grey: high toxin level; ν T = 0.5, δ T = 0) and K-switching rate (solid black line; ν K = 0.2).We notice that total population N rapidly attains quasi-stationarity well before a fixation event occurs (S fixation in this example).
-In Fig. S5(b), we show faster environmental T and K switching, with ν T = 10 and ν K = 5.There is also a bias towards the harsh T state favouring the strain R (δ T = −0.15,high toxin level) that is responsible for a robust offset in strain abundance eventually leading to extinction of S and fixation of R.
In this example, the K-EV switching rate is sufficiently high to ensure its self-averaging: in this fast K switching regime, the population size is distributed about the effective carrying capacity K = K 0 (1 − γ 2 )/(1 − γδ K ) (here N ≈ K = 667), see Fig. 4 and Sec.IV A.
-Fig.S5(c) illustrates the case of fast T -EV (ν T = 10) with a bias towards the mild/low T state favouring the strain S (δ T = 0.15).In this example, the K-EV switching rate is low (ν K = 0.05) and the population size N follows K(t) and fluctuates about K − and K + .The T -EV bias (δ T > 0) is here responsible for a systematic offset in the strain abundance, with N S > N R , resulting in the fixation of S and extinction of R, see Fig. 2(a).In this example, K-EV is responsible for larger demographic fluctuations in the environmental state ξ K = −1, where N ≈ K − and S fixation is more likely than when N ≈ K + (ξ K = +1).
- -Fig.S5(e,f) illustrate the dynamics under fast T -EV (ν T = 10) and extremely slow K-switching rate (ν K = 5 × 10 −6 ), with an initial carrying capacity K(0) = K − in (e) and K(0) = K + (f).There is also a small bias towards the mild T state favouring S (δ T = 0.1).This choice of parameters corresponds to a point in the faint green region in the diagram of Fig. 5(a), where 0 < η < 1 and long-lived coexistence is possible but not certain.In panel (e), the population size fluctuates about K(0) = K − (N ≈ K − ) and is not able to sustain long-lived coexistence: demographic fluctuations yield S fixation in a time t ≲ 2K − = 800.In panel (f), we have N ≈ K + and, as demographic fluctuations are unable to cause the fixation/extinction of either strain in a time less than 2 ⟨N ⟩ = 2K + = 4000 (for clarity, the time series in Fig. S5 have been truncated at t = 2400), we have long-lived coexistence.
In all examples in the panels of Fig. S5, the population size N reaches quasi-stationarity (settles in the QPSD) a significant time before the fixation of one strain and extinction of the other, in line with the considerations underpinning Eqs.(19) and (20).

SM8. FULLY CORRELATED AND ANTI-CORRELATED ENVIRONMENTAL VARIABILITY
In the case of fully correlated/anti-correlated T -EV and K-EV, environmental variability is no longer twofold since we have ξ ≡ ξ T and ξ K = ξ (fully correlated EV) or ξ K = −ξ (fully anti-correlated EV).The switching carrying capacity can thus be written as where γ = γ in the fully correlated case, and γ = −γ when T -EV and K-EV are fully anti-correlated.For instance, this implies that, in the correlated case, the environmental state ξ = +1 corresponds to f s = e s > 1 and K = K + (low toxin level, abundant resources), while ξ = −1 is associated with f S = e −s < 1 and K = K − (high toxin level, scarce resources).As said, under fully correlated/anti-correlated T /K-EV, environmental variability is no longer twofold: ξ simultaneously drives the level of toxin and the abundance of resources.Hence, we can characterise the effect of fully correlated/anti-correlated EV in terms of ν ≡ ν T and δ ≡ δ T , and the fully anti-correlated case is related to completely correlated EV via γ → − γ.An example of fully-anticorrelated EV modelled in terms of a dichotomous process driving the level of toxins and resources in the context of competitive exclusion is considered in [65].
Fig. S6 shows the comparison between the uncorrelated T -EV and K-EV studied in the main text (top row), and the fully correlated (middle row) and fully anti-correlated (bottom row) cases reported here, all under K 0 = 1000 and γ = 0.9, and for different selection strengths s ∈ {0.2, 2, 20} (left to right columns).For the uncorrelated case, the parameters of K-EV are independent from those of T -EV parameters, and are here chosen to be (ν K , δ K ) = (10 −4 , 0), i.e. unbiased slow-switching K-EV (similar to the example of Fig. 5(a)).
Since the fully correlated and anti-correlated cases (middle and bottom rows) are mirror images through the horizontal axis (symmetric respect to δ T = 0) and under a red-blue colour change, we focus on the correlated case only (middle row).For this case, γ = γ = 0.9, with ξ K = ξ T ≡ ξ, ν K = ν T , and δ K = δ T .In the fully correlated case, we observe that both the blue (resistant fixation) and the bright green (coexistence) regions shift upwards (to higher toxin level biases δ T ) as the selection strength increases (from left to right column).We can understand this phenomenon in light of the T -correlated K-EV fluctuations.This is, since δ K = δ T , a lower value of δ T in the diagrams implies longer cumulative periods in the harsh toxin level, but also in the low carrying capacity environment.Therefore, lower δ T provides the selective advantage to resistant strain at the same time that it shrinks the total population.Demographic fluctuations being stronger in smaller populations, correlated T -EV and K-EV thus provide higher R fixation probability (blue region shifted up), as well as lower coexistence probability (green region shifted up).Moreover, since total population increases with the bias towards positive carrying (here δ K = δ T ), see Sec.IV A, the MFT thus increases with δ T (see Fig. 2(c)), and the coexistence probability thus shifts upwards.The magnitude of the upwards shift in the correlated case, is small but increases with selection strength s that increases the amplitude of the T -EV fluctuations.
In summary, we obtain the same qualitative results for fully correlated/anti-correlated T /K-EV as when ξ T /K are independent (uncorrelated environmental noise, twofold EV), with some minor quantitative differences, as shown in Fig. S6 and discussed above.We conclude that the similar behaviour observed for uncorrelated and (anti-)correlated T /K-EV indicates that our findings are robust against the detailed model specifications: the results are expected to be valid for the general case of twofold environmental variability where T /K-EV are neither completely independent nor fully correlated/anti-correlated.

FIG. 3 .
FIG. 3. Fixation/coexistence diagrams under T -EV in the (νT , δT ) parameter space for a small system of constant carrying capacity K0 = 50, with selection bias s = 0.1 (a,d), s = 1 (b,e), and s = 10 (c,f), after a time t = 2K0.(a)-(c): phase diagrams obtained from stochastic simulations of the model with a constant K0 (N fluctuates about K0) over 10 3 realisations and coded according to the RGB colourmap of panel (g): red / blue (grey / charcoal) corresponds to the likely fixation of S/R (red (grey): S dominance; blue (charcoal): R dominance), regions indicated as "S/R fixation" in panel (a); magenta (dark grey) indicates where fixation of R or S is likely, area between "S/R fixation" regions in panel (a); green (light grey) indicates where long-lived coexistence is most likely, area highlighted as "S + R coex." in panel (c).(d)-(f): same as in (a)-(c) but from numerically exact solutions of Eq. (15) for the corresponding Moran process with a constant population size N = K0, see text.
N ) Starting fraction of resistant bacteria x0 Fixation probability in state ξT with Moran approximation (initially n bacteria of type R) ϕ ξ T n Mean fixation time in state ξT with Moran approximation (initially n bacteria of type R) τ ξ T n Fixation probability with the Moran approximation (x0 = 0.5) ϕMA Mean fixation time with the Moran approximation (x0 = 0.5) τMA Fixation probability of strain R ϕ Coexistence probability η Mean population size ⟨N ⟩ Modal population size N Equilibrium fraction of R in coexistence state x * Fraction of resistant strain R (regardless of coexistence/fixation) ⟨x⟩ (⟨x⟩ ≈ x * when η ≈ 1) Population density function of the N -PDMP p(N ) Mean population size of type R ⟨NR⟩ SM2.T -SWITCHING COEXISTENCE: FLUCTUATION AMPLITUDE AND TIME SCALE FIG. S1.Predictions of x-PDMP stationary density ρ(x) against x for δT ∈ {−0.5, 0.0, 0.5} from left to right.Theoretical predictions (solid lines) from Eq. (S7) are in excellent agreement with simulation results (×).Blue vertical line shows predicted mean x * from Eq. (10) and green vertical line shows predicted mode x from Eq. (S8).Parameters used are s = 10, νT = 5, and fixed K0 = 100.Simulation results have been averaged over 10 3 realisations.
Fig. S5(d) shows typical trajectories in the long-lived coexistence phase in the case of unbiased fast switching T -EV and K-EV, with ν T = 10, δ T = 0 and ν K = 50.This set of parameters essentially corresponds (here ν K = 50 instead of ν K = 5) to a point in the bright green region in the diagram of Fig. 5(c), where η ≈ 1 and long-lived coexistence is almost certain.As in panel (b), the fast K switching (solid black line not shown in Fig. S5(d)) leads to fluctuations of the total population size about K, i.e.N ≈ K = 667.In this example x * = 1/2, see Fig. 7, and the number of R and S cells fluctuates about their averages: ⟨N R ⟩ ≈ ⟨N S ⟩ ≈ K/2 ≈ 333.5, see Sec.V B.

TABLE I .
Table of parameters and physical quantities of interest.