Population Dynamics in a Changing Environment: Random versus Periodic Switching
PPopulation Dynamics in a Changing Environment: Random versus Periodic Switching
Ami Taitelbaum § , Robert West § , Michael Assaf, ∗ and Mauro Mobilia † Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds LS2 9JT, U.K.
Environmental changes greatly influence the evolution of populations. Here, we study the dynam-ics of a population of two strains, one growing slightly faster than the other, competing for resourcesin a time-varying binary environment modeled by a carrying capacity switching either randomly or periodically between states of abundance and scarcity. The population dynamics is characterizedby demographic noise (birth and death events) coupled to a varying environment. We elucidatethe similarities and differences of the evolution subject to a stochastically- and periodically-varyingenvironment. Importantly, the population size distribution is generally found to be broader underintermediate and fast random switching than under periodic variations, which results in markedlydifferent asymptotic behaviors between the fixation probability of random and periodic switching.We also determine the detailed conditions under which the fixation probability of the slow strain ismaximal. PACS numbers:
The evolution of natural populations is influenced byvarying environmental conditions: the abundance of nu-trients, toxins, or external factors like temperature aresubject to random and seasonal variations, and have animportant impact on population dynamics [1–3].Several models of a population response to a changingenvironment assume that external conditions vary eitherperiodically or stochastically in time [4–26]. These ex-ternal variations are often modeled by taking a binaryenvironment that switches between two states [26–46].In finite populations, demographic noise (DN) is anotherform of randomness that can lead to fixation (one speciestakes over the population [47, 48]). DN is strong in smallpopulations and negligible in large ones. Importantly,the evolution of a population composition is often cou-pled with the dynamics of its size [49–55]. This can leadto coupling between DN and environmental variability(EV), with external factors affecting the population size,which in turn modulates the DN strength. The interplaybetween EV and DN plays a key role in microbial com-munities [56–63]: the variations of their composition andsize are vital to understand the mechanisms of antimi-crobial resistance [26, 58], and may lead to populationbottlenecks , where new colonies consisting of few indi-viduals are prone to fluctuations [56, 59, 61–63]. Inter-actions between microbial communities and environmenthave also been found to influence cooperative behavior in
Pseudomonas fluorescens biofilms [60–62]. EV and DNare also important in ecology, e.g. , in modeling tropicalforests [14–16], and in gene regulatory networks [17, 21].In most studies, there is no interdependence betweenthe fluctuations stemming from DN and EV, with growthrates often assumed to vary independently of the popula-tion size [6, 7, 11–13, 17–22, 24, 27–29, 31, 36–38, 41, 43]. § Equally contributed to this work. ∗ Electronic address: [email protected] † Electronic address: [email protected]
Hence, there is as yet no systematic comparison of thedynamics under random and periodic switching: someworks report that they lead to similar evolutionary pro-cesses while others find differences, see e.g. , Refs. [29, 40].Here, we systematically study the coupled influence of EVand DN on the dynamics of a population, where slow- andfast-growing strains compete for resources subject to arandomly- and periodically-switching carrying capacity.A distinctive feature of this model is that it accountsfor the stochastic or periodic depletion and recovery ofresources via a binary environment, varying with a fi-nite correlation time or period, and the
DN and EV cou-pling , see Fig. 1. This setting is simple enough to enableus to scrutinize whether environmental perturbations ofdifferent nature lead to the same dynamics, and includesmany features (switching environment, varying popula-tion size) that can be tested in controlled microbial ex-periments [28, 52, 55, 57, 63].To address the fundamental question of evolution un-der stochastic and deterministic variations, we considerrandom and periodic environmental switching. This al-lows us to elucidate the influence of EV on the populationsize distribution (PSD) and the fixation properties. Weanalytically show that the PSD is generally broader un-der intermediate and fast random switching than underperiodic variations, leading to markedly different fixationprobabilities. We also determine the switching conditionsfor which the slow strain’s fixation probability is maxi-mized.We consider a well-mixed population of time-fluctuating size N ( t ) = N S ( t ) + N F ( t ) consisting of twostrains. At time t , N S ( t ) individuals are of a slow-growing strain S , corresponding to a fraction x = N S /N of the population, and N F are of a fast-growing species F . The respective per-capita growth rates of S and F are(1 − s ) / ¯ f and 1 / ¯ f , which sets the model’s time scale [64].Here, ¯ f = (1 − s ) x + 1 − x = 1 − sx is the popula-tion average fitness and 0 < s (cid:28) F over S [44, 45, 50, 51]. a r X i v : . [ q - b i o . P E ] J u l FIG. 1: (a) K vs. time t : asymmetric random (pink/lightgray) and periodic (black dashed) switching between K + and K − yield fluctuating population composition and size(large/small circles), see text. (b, c) Typical realizationsof N (black), N S (red/gray) and K (black dashed) vs. t under random (b) and periodic (c) switching: composi-tion changes until fixation occurs. Here ( s, K , ν, γ, δ, x ) =(0 . , , . , . , . , . Growth is limited by a logistic death rate
N/K , where K (cid:29) K ) yielding a constant or logistically-varying N [48, 65–68]. Here, we instead consider a popula-tion of fluctuating size subject to a time-varying envi-ronment, and obeying the birth-death process [45, 64]: N S/F T + S/F −−−→ N S/F + 1 and N S/F T − S/F −−−→ N S/F −
1, withtransition rates T + S = (1 − s ) N S / ¯ f , T + F = N F / ¯ f and T − S/F = (
N/K ( t )) N S/F . We model EV via a switchingcarrying capacity K ( t ) = K [1 + γξ α ( t )] , ξ α ( t ) ∈ {− , +1 } , (1)where K ≡ ( K + + K − ) / γ ≡ ( K + − K − ) / (2 K ),while α ∈ { r, p } and γ = O (1). Here, resources vary ei-ther randomly ( α = r ) or periodically ( α = p ), betweenstates of scarcity, K = K − ( ξ α = − K = K + ( ξ α = +1), where K + > K − (cid:29)
1, causing fluc-tuations of population size and composition, see Fig. 1.This specific choice of birth-death process coupled to atime-varying binary environment is arguably the simplestbiologically-relevant model to study population dynam-ics under the joint influence of EV and DN, see Sec. S1.1in [64].When K ( t ) switches randomly , ξ r is a colored asym-metric dichotomous (telegraph) Markov noise (ADN) [69,70], with the transition ξ r → − ξ r occurring at rate ν ± when ξ r = ±
1. The (average) switching rate is ν = ( ν + + ν − ) / δ = ( ν − − ν + ) / (2 ν ) measures theswitching asymmetry ( | δ | <
1, with δ = 0 for symmet-ric switching). In this model, the ADN is a stationarynoise of mean (cid:104) ξ r ( t ) (cid:105) = δ and autocorrelation function (cid:104) ξ r ( t ) ξ r ( t (cid:48) ) (cid:105) − (cid:104) ξ r ( t ) (cid:105)(cid:104) ξ r ( t (cid:48) ) (cid:105) = (1 − δ ) e − ν | t − t (cid:48) | ( (cid:104)·(cid:105) de-notes ensemble averaging). When K ( t ) switches periodi-cally , ξ p is a rectangular wave defined by the rectangular function, rect( · ) [71], of period T = (1 /ν + ) + (1 /ν − ) =2 / [(1 − δ ) ν ]: ξ p ( t ) = ∞ (cid:88) j = −∞ (cid:34) rect (cid:32) t + ν + + jT /ν + (cid:33) − rect (cid:32) t − ν − + jT /ν − (cid:33)(cid:35) , which becomes the square wave ξ p ( t ) = − sign { sin ( πνt ) } when δ = 0. In our simulations, ξ p ( t ) averaged over aperiod T has the same mean and variance as ξ r ( t ). Hence,the mean and variance of K ( t ) are the same for α ∈{ r, p } : (cid:104) K ( t ) (cid:105) = K (1 + γδ ) and var( K ) = ( γK ) (1 − δ ) [72].The model considered here gives rise to a long-lived population size distribution (PSD) followed by an even-tual extinction of the entire population which occursafter a very long time (practically unobservable when K (cid:29) t = O ( s − ), a timescale on which one species is likely tohave gone extinct and the other fixated the populationthat is in its long-lived PSD [64]. We show that the fixa-tion probabilities strongly depend on the PSD which is en-coded in the underlying master equation [45, 46, 74, 75],see [64] for details.Insight into the dynamics is gained by ignoring fluctu-ations and considering the mean-field picture of a verylarge population with constant K = K . Here, N and x evolve according to dN/dt ≡ ˙ N = N (1 − N/K ) and˙ x = − sx (1 − x ) / (1 − sx ) [50, 51, 64], with x decayingon a timescale t ∼ s − (cid:29) N ( t ) = O ( K ) after t = O (1) [76]. Thus, a timescale separation occurs: therelaxation of x is much slower than that of N .However, when dealing with a finite population, DN(random birth/death events) must be taken into account,yielding the fixation of one of the species. The S fixationprobability, given a fixed population size N , and an initialfraction x = N S (0) /N (0) of S individuals, is [48, 67, 75] φ ( x ) | N = (cid:2) e − Nx ln(1 − s ) − (cid:3) / (cid:2) e − N ln(1 − s ) − (cid:3) , (2)which exponentially decreases with N . For s (cid:28) N − / (cid:28) φ ( x ) | N (cid:39) ( e − Ns (1 − x ) − e − Ns ) / (1 − e − Ns ) [44, 45, 66].While Eq. (2) provides a good approximation for the fix-ation probability also when N fluctuates about constant K = K , this picture changes drastically when, in addi-tion to DN, the population is subject to a time-varying K ( t ), see Fig. 1. Below we study the joint influence ofEV and DN on the PSD and fixation properties. Population size distribution.
Simulations show thatthe marginal quasi-stationary PSD, P ( α ) ν ( N ) (uncondi-tioned of ξ α ), is characterized by different regimes de-pending on the switching rate ν , with markedly differentfeatures in the case of random and periodic variationswhen ν = O (1) and ν (cid:29)
1, see Fig. 2.The case of random switching can be treated as in[44, 45] for δ = 0. Upon ignoring DN, N ( t ) is there-fore subject only to ADN according to the piecewise-deterministic Markov process (PDMP) [64, 77, 78] FIG. 2: P ( r ) ν ( N ) (blue/dark gray) and P ( p ) ν ( N ) (red/gray)for different ν : (a) ν = 0 .
05, (b) ν = 17 .
5, (c) ν = 1 .
4, (d) ν = 1. Symbols are from simulations; solid black lines in (a)-(d) are from P PDMP ν , those in cyan/light gray are from P ( N )in (a), P Kap ν in (b), and P PPP ν in (c,d); vertical lines show N = K ± (dashed) in (a,c) and N = N min / max (cyan/lightgray) in (c,d), see text and Sec. S2.3 in [64]; horizontal dashedlines are eyeguides. Here ( s, K , γ, x ) = (0 . , , . , . δ = 0 . δ = − . defined by the stochastic differential equation ˙ N = N [1 − ( N/ K )(1 − γξ r ) / (1 − γδ )], where K ≡ K (1 − γ ) / (1 − γδ ). When ν → ∞ , the ADN self-averages, ξ ν →∞ −−−−→ (cid:104) ξ (cid:105) = δ , and N ν →∞ −−−−→ K . The marginal PSD ofthis PDMP has support [ K − , K + ] and can be computedexplicitly [45, 69]: its expression P PDMP ν ( N ) is given byEq. (S22) of [64]. Although P PDMP ν only accounts for EV,when K (cid:29) γ = O (1), it captures the peaks of P ( r ) ν and the average population size, see Figs. 2 and S3(b)in [64]. However, P PDMP ν ignores DN and cannot cap-ture the width of P ( r ) ν about its peaks, see Fig. 2(a,c,d).Yet, this can be remedied, by a linear noise approxima-tion, see [45] and Sec. S3.2 in [64]. We can also obtain aPDMP-like approximation (ignoring DN) [70, 79] of theperiodic PSD by solving the mean-field equation for N ( t )with periodic K ( t ). By inverting N ( t ) we then obtain thepiecewise periodic process (PPP) approximation P PPP ν of P ( p ) ν , given by (S19) in Sec. S2.3 of [64], which is validover a broad range of switching rates, see Fig. 2(c,d) andbelow.Furthermore, for periodic switching, the full P ( p ) ν canbe found analytically in the limits of very slow ( ν → ν (cid:29)
1) variations. For ν → i.e. K ( t ) (cid:39) K (0). The PSD is thus the samefor periodic and random switching: P ( p )0 = P ( r )0 ≡ P ,and can be computed from the master equation. Assum-ing K (cid:29) γ = O (1), the PSD is bimodal withpeaks about N = K ± , whose intensity depends on δ [64]: P ( N ) (cid:39) [(1 + δ ) K N +1+ e − K + + (1 − δ ) K N +1 − e − K − ] / [2 N · N !]. This result excellently agrees with simulations, seeFig. 2(a). Under fast periodic switching, P ( p ) ν differsmarkedly from its random counterpart, see Fig. 2(b).An approximate expression of P ( p ) ν to leading order in1 /ν , here denoted by P Kap ν , and peaked at N = K when ν → ∞ is given by Eq. (S15) in [64]. P Kap ν is obtainedfrom the master equation by using the WKB approx-imation [80] and the Kapitza method [12, 24, 81], i.e. separating the dynamics into fast and slow variables,see Sec. 2.2 of [64]. In Fig. 2(b), we notice that both P ( p ) ν (cid:39) P Kap ν and P ( r ) ν (cid:39) P PDMP ν are unimodal andpeaked about N ≈ K when ν (cid:29)
1, but P Kap ν is muchsharper and narrower than P PDMP ν . In fact, the varianceof P PDMP ν scales as K /ν when 1 (cid:28) ν (cid:28) K , and ismuch larger than that of P Kap ν , see Sec. S4.3 in [64].Note that while P and P Kap ν account for DN and EV, P PDMP ν and P PPP ν only account for EV. Yet, DNis negligible compared to EV when 1 (cid:46) ν (cid:28) K and1 (cid:46) ν (cid:28) √ K in the random and periodic cases, re-spectively [64]. P PDMP ν and P PPP ν are therefore suitableapproximations of P ( α ) ν in those regimes.In particular, P PDMP ν and P PPP ν allow us to char-acterize interesting phenomena arising in the interme-diate asymmetric switching regime where ν (cid:38) ν − > ν + <
1, or ν − < ν + > i.e. when1 / (1 + | δ | ) < ν < / (1 − | δ | ). In the former case ( δ > P ( r ) ν has a peak at N ≈ K + and, under sufficiently strongEV, exhibits also a peak N ∗ between K − and K + ( i.e. K − < N ∗ < K + ), whose position is aptly captured by P PDMP ν , see Fig. 2(c) and Sec. S3.1 in [64]. In Fig. 2(c), P ( p ) ν is less broad than P ( r ) ν and has also two peaks wellreproduced by P PPP ν whose support is narrower than thatof P PDMP ν [64]. When ν (cid:38)
1, with ν − < ν + > δ < P ( r ) ν and P ( p ) ν exhibit a single peak at N ≈ K − ,well predicted by P PDMP ν and P PPP ν , with the latter beingnarrower than the former in Fig. 2(d). In fact, Figs. 2(b)and 2(c) show that the transition from bimodal to uni-modal PSD (slow to fast switching) is generally moreabrupt under periodic than under random switching. Fixation probability.
We denote by φ α the slow ( S )species fixation probability subject to α -switching ( α ∈{ r, p } ). As aforementioned, when s (cid:28) t (cid:38) O (1),the system has settled in its long-lived PSD. Thus, given x , φ α can be approximated by averaging φ ( x ) | N over P ( α ) ν/s ( N ), upon rescaling ν → ν/s [44, 45] φ α ( ν ) (cid:39) (cid:90) ∞ P ( α ) ν/s ( N ) φ ( x ) | N dN, α ∈ { r, p } . (3)This result is valid under weak selection, 1 /K (cid:28) s (cid:28)
1, when there are O ( ν/s ) switches prior to fixation [44,45, 64]. The difference between φ r and φ p stems fromthe different ν -dependence of P ( r ) ν and P ( p ) ν , see Fig. 2.Approximations of φ r and φ p are obtained by respectivelysubstituting P ( α ) ν/s by P PDMP ν/s and P PPP ν/s into Eq. (3). Thisyields expressions (S38) and (S39) of [64] which are validover a broad range of ν [44, 64], see Fig. 3 and S2(c,d) of[64]. Notably, when ν/s (cid:29) φ p is better approximatedby substituting P ( p ) ν/s by P Kap ν/s in Eq. (3), see below and[64].When ν → P ( α ) ν/s is peakedat N = K ± . Hence, with Eq. (3), lim ν → φ α ( ν ) (cid:39) φ (0) =[(1 − δ ) φ ( x ) | K − + (1 + δ ) φ ( x ) | K + ] /
2. Fig. 3(d) confirmsthat φ r and φ p approach φ (0) when ν/s (cid:28) FIG. 3: (a)-(d) fixation probability for random/periodicswitching (circles/squares): symbols are from simulations.In (a) solid lines are from (S38) (blue/dark gray, α = r ),(4) (red/gray, α = p ) and (S39) (black, α = p ) of [64].(a) φ α versus ν with δ = 0 .
2; dashed line shows φ ( ∞ ) .Inset: φ r /φ p versus ν with δ = 0 .
2. (b,c) ln ( φ α /φ ( ∞ ) )versus s/ν for random (b) and periodic (c) switching with δ = 0 . δ = 0 (blue/dark gray). Dashedgray lines are eyeguides ∝ s/ν in (b) and ∝ ( s/ν ) in (c).(d) Non-monotonic φ α ( ν ) with δ = 0 . δ = 0 . φ (0 , ∞ ) . φ r ( ν ) and φ p ( ν ) are maximal at ν = ν ∗ r ≈ . ν = ν ∗ p ≈ .
07, see text. (e) Heatmapof ν ∗ r (see Sec. S5.1 in [64] for details and heatmap of ν ∗ p ): ν ∗ r → , ∞ in the black and white areas, respectively; φ r ( ν )is non-monotonic in the red-yellow/gray area, with ν ∗ r ≈ . ν ∗ r ≈ . δ = 0 . δ = 0 . s, K , γ, x ) = (0 . , , . , .
5) in (a)-(c) and (0 . , , . , .
6) in (d,e).
When ν/s (cid:29) P ( α ) ν/s is sharplypeaked at N (cid:39) K , see Fig. 2(b), and to leading orderlim ν →∞ φ α ( ν ) (cid:39) φ ( ∞ ) = φ ( x ) | K [44, 45]. Simulation re-sults of Fig. 3 confirm that at ν (cid:29) s , φ r ( ν ) and φ p ( ν ) con-verge to φ ( ∞ ) . Thus, the fixation probability under fastrandom/periodic switching is the same to lowest order in1 /ν . Yet, the rate of convergence differs, see Fig. 3(a).This is explained by computing the next-to-leading or-der of φ α in ν/s (cid:29)
1. For this, we use Eq. (3) withEq. (2) and P PDMP ν/s and P Kap ν/s for random and periodicswitching, respectively. A saddle-point calculation, with1 /K (cid:28) s (cid:28)
1, yields (see Sec. S4 in [64])ln (cid:18) φ α ( ν ) φ ( ∞ ) (cid:19) (cid:39) (cid:40) A r ( s/ν ) ( α = r ) A p ( s/ν ) ( α = p ) . (4)Here φ ( ∞ ) = e m/ , m ≡ K (1 − x ) ln (1 − s ), and A r = m (4 + m )(1 − δ )( γ/ (1 − γδ )) /
16 while A p = K (1 − (1+ m/ K ) )( γ/ (1 − γδ )) /
72. Thus, when K s (cid:29) φ α ( ν ) converges to φ ( ∞ ) much faster in the periodic than in the random case, see Fig. 3(a)-(c). The differentasymptotic behavior can be understood by noting that P ( r ) ν is generally broader than P ( p ) ν , with respective vari-ances scaling as ν − and ν − . N can thus attain smallervalues under random than periodic switching, which en-hances φ r with respect to φ p [82]. When ν/s (cid:29) φ r,p isdetermined by the mean (cid:104) N (cid:105) (cid:39) K of P ( α ) ν , and the differ-ent rate of convergence to φ ( ∞ ) stems from the deviationsof (cid:104) N (cid:105) from K , which decrease as ν − when α = r and ν − when α = p , see Sec. S4.3 in [64]. Another signatureof the different asymptotic behavior is the sharp peak ofthe ratio φ r /φ p at a nontrivial ν , see Fig. 3(a, inset).Under intermediate (rescaled) switching, φ α exhibits arich behavior, see Fig. 3(d). When the switching asym-metry is sufficiently large, φ α is a non-monotonic func-tion of ν in a nontrivial region γ > γ c ( s ), δ > δ c ( γ, s ) ofthe parameter space that can be found from Eq. (3), seeFig. 3(d,e) and Sec. S5.1 in [64]. The PDMP- and PPP-based approximations [Eqs. (S38) and (S39) in [64]] ad-equately capture the ν -dependence of φ α in this regime,and its maximum at ν ∗ α ∼ s . This optimal switch-ing rate, which maximizes the S species fixation prob-ability at given ( γ , δ , s ), corresponds to O (1) switchesprior to fixation. The relative increase in φ α ( ν ), givenby φ α ( ν ∗ α ) / max( φ (0) , φ ( ∞ ) ) − ν ∗ p (cid:46) ν ∗ r , and φ p ( ν ∗ p )is narrower around the peak than φ r ( ν ∗ r ), see Figs. 3(d,e)and S2(e) of [64]. When the asymmetry is not toolarge ( | δ | < δ c ), φ α ( ν ) is a monotonic function: it in-creases/decreases with ν below/above a critical selectionintensity s c (with γ, δ fixed), see Sec. S5.2 and Fig. S2(d)in [64]. Remarkably, transitions between monotonic andnon-monotonic behavior of φ α ( ν ) are also found when S produces public goods benefiting the entire population,see Sec. S7 in [64].Inspired by the evolution of microbial communities influctuating environments, we have studied the dynam-ics of a population of two strains competing for resourcessubject to a binary carrying capacity, switching randomly or periodically in time. We have analyzed how the cou-pling of demographic noise and environmental variabilityaffects the population size and fixation properties. Wehave shown that the population size distribution is gener-ally broader under random variations than under periodicchanges in the intermediate/fast switching regime, whichlead to markedly different asymptotic behaviors of thefixation probabilities. We have also determined the con-ditions under which the probability that the slow speciesprevails is maximal. Our work sheds light on the similar-ities and differences of evolution in stochastically- versusdeterministically-varying environments, and is thus rele-vant to microbial communities, often subject to frequentand extreme environmental changes.We are grateful to E. Frey, A. M. Rucklidge, andK. Wienand for useful discussions. AT and MA ac-knowledge support from the Israel Science Foundationgrant No. 300/14 and the United States-Israel Bina-tional Science Foundation grant No. 2016-655. The sup-port of an EPSRC Ph.D. studentship to RW (Grant No. EP/N509681/1) is also gratefully acknowledged. [1] C. R. Morley, J. A. Trofymow, D. C. Coleman, andC. Cambardella, Microbiol. Ecol. , 329 (1983).[2] C. A. Fux, J. W. Costerton, P. S. Stewart, and P. Stood-ley, Trends Microbiol. , 34 (2005).[3] Caporaso et al., ,Genome Biology :R50 (2011).[4] H. Beaumont, J. Gallie, C. Kost, G. Ferguson, and P.Rainey, Nature , 90 (2009).[5] P. Visco, R. J. Allen, S. N. Majumdar, and M. R. Evans,Biophys. J. , 1099 (2010).[6] R. M. May, Stability and complexity in model ecosystems (Princeton University Press, Princeton, USA, 1973).[7] S. Karlin and B. Levikson, T. Pop. Biol. , 383 (1974).[8] P. L. Chesson and R. R. Warner, American Naturalist , 923 (1981).[9] M. Loreau and C. de Mazancourt, American Naturalist , E49 (2008).[10] B. K. Xue and S. Leibler, Phys. Rev. Lett. , 108103(2017).[11] E. Kussell and S. Leibler, Science , 2075 (2005).[12] M. Assaf, A. Kamenev and B. Meerson, Phys. Rev. E.78, 041123 (2008).[13] M. Assaf, A. Kamenev and B. Meerson, Phys. Rev. E.79, 011127 (2009).[14] R. A. Chisholm et al. Ecology Letters : 855 (2014).[15] D. A. Kessler and N. Shnerb, J. Theor. Biol. , 1(2014).[16] M. Kalyuzhny, R. Kadmon and N. M. Shnerb, EcologyLetters : 572 (2015).[17] M. Assaf, E. Roberts, Z. Luthey-Schulten, and N. Gold-enfeld, Phys. Rev. Lett. , 058102 (2013).[18] Q. He, M.Mobilia, and U. C. T¨auber, Phys. Rev. E ,051909 (2010).[19] U. Dobramysl, and U. C. T¨auber, Phys. Rev. Lett. ,048105 (2013).[20] M. Assaf, M. Mobilia, and E. Roberts, Phys. Rev. Lett. , 238101 (2013).[21] E. Roberts, S. Be’er, C. Bohrer, R. Sharma and M. Assaf,Phys. Rev. E. 92, 062717 (2015).[22] A. Melbinger and M. Vergassola, Scientific Reports ,15211 (2015).[23] M. Assaf and B. Meerson, J. Phys. A: Math and Theo.50, 263001 (2017).[24] O. Vilk and M. Assaf, Phys. Rev. E. 97, 062114 (2018).[25] U. Dobramysl, M. Mobilia, M. Pleimling, and U. C.T¨auber, J. Phys. A: Math. Theor. , 063001 (2018).[26] L. Marrec and A.-F. Bitbol, PLoS Comput. Biol. :e1007798 (2020).[27] E. Kussell, R. Kishony, N. Q. Balaban, and S. Leibler,Genetics , 1807 (2005).[28] M. Acar, J. Mettetal, and A. van Oudenaarden, NatureGenetics , 471 (2008).[29] M. Thattai and A. Van Oudenaarden, Genetics , 523(2004).[30] S. P. Otto and M. C. Whitlock, Genetics , 723 (1997).[31] B. Ga´al, J. W. Pitchford, and A. J. Wood, Genetics ,1113 (2010).[32] K. Wienand, MSc Thesis (Ludwig-Maximilians- Universit¨at M¨unchen, 2011).[33] P. Patra and S. Klumpp, PLoS One (5): e62814 (2013).[34] P. Patra and S. Klumpp, Phys. Biol. , 683 (2013).[36] P. Ashcroft, P. M. Altrock, and T. Galla, J. R. Soc. In-terface , 20140663 (2014).[37] M. Danino and N. M. Shnerb, J. Theor. Biol. , 84(2018).[38] P. G. Hufton, Y. T. Lin, and T. Galla, J. Stat. Mech.Theory Exp. (2018).[39] Q. Su, A. McAvoy, L. Wang, and M. A. Nowak, Proc.Natl. Acad. Sci. USA , 25398 (2019).[40] I. Meyer and N. M. Shnerb, e-print: arXiv1912.06386.[41] P. G. Hufton, Y. T. Lin, T. Galla, and A. J. McKane,Phys. Rev. E , 052119 (2016).[42] J. Hidalgo, S. Suweis, and A. Maritan, J. Theor. Biol. , 1 (2017).[43] R. West, M. Mobilia, and A. M. Rucklidge, Phys. Rev. E , 022406 (2018).[44] K. Wienand, E. Frey, and M. Mobilia, Phys. Rev. Lett , 158301 (2017).[45] K. Wienand, E. Frey, and M. Mobilia, J. R. Soc. Interface , 20180343 (2018).[46] R. West and M. Mobilia, J. Theor. Biol. , 110135(2020).[47] J. F. Crow and M. Kimura, An Introduction to Popu-lation Genetics Theory (Blackburn Press, New Jersey,2009).[48] W. J. Ewens,
Mathematical Population Genetics (Springer, New York, 2004).[49] J. Roughgarden,
Theory of Population Genetics and Evo-lutionary Ecology: An Introduction (Macmillan, NewYork, 1979).[50] A. Melbinger, J. Cremer, and E. Frey, Phys. Rev. Lett. , 178101 (2010).[51] J. Cremer, A. Melbinger, and E. Frey, Phys. Rev. E ,051921 (2011).[52] J. Cremer, A. Melbinger, and E. Frey, Sci. Rep. , 281(2012).[53] A. Melbinger, J. Cremer, and E. Frey, J. R. Soc. Interface , 20150171 (2015).[54] C. S. Gokhale and C. Hauert, Th. Pop. Biol. , 28(2016).[55] J. S. Chuang, O. Rivoire, and S. Leibler, Science ,272 (2009).[56] L. M. Wahl, P. J. Gerrish, and I. Saika-Voivod, Genetics , 961 (2002).[57] K. Wienand, M. Lechner, F. Becker, H. Jung, and E.Frey, PloS one, 10(8), e0134300 (2015).[58] J. Coates, B. R. Park, D. Le, E. S¸im¸sek, W. Chaudhry,and M. Kim, eLife :e32976 (2018).[59] Z. Patwas and L. M. Wahl, Evolution , 1166 (2009).[60] P. B. Rainey and K. Rainey, Nature 425, 72 (2003).[61] M. A. Brockhurst, A. Buckling, and A. Gardner, Curr.Biol. , 761 (2007).[62] M. A. Brockhurst, PLoS One , e634 (2007). [63] J. Cremer, A. Melbinger, K. Wienand, T. Henriquez, H.Jung, and E. Frey, e-print: arXiv:1909.11338.[64] See the Supplementary Material provided as an appendixto this e-print (pages 7-22); also available on Figshare atthe URL https://doi.org/10.6084/m9.figshare.12613370,where additional supporting resources, including videosillustrating Figures 2, 3(a), S1 and S4(b)-(e), are pro-vided[65] P. A. P. Moran, The statistical processes of evolutionarytheory (Clarendon, Oxford, 1962).[66] R. A. Blythe and A. J. McKane, J. Stat. Mech.
P07018 (2007).[67] T. Antal and I. Scheuring, Bull. Math. Biol. , 1923(2006).[68] R. M. Nowak, Evolutionary Dynamics (Belknap Press,Cambridge, USA, 2006).[69] W. Horsthemke and R. Lefever,
Noise-Induced Transi-tions (Springer, Berlin, 2006).[70] I. Bena, Int. J. Mod. Phys. B , 2825 (2006).[71] A rectangular function is defined as follows: rect( x ) = 1 if | x | < /
2, rect( x ) = 0 if | x | > /
2, while rect( ± /
2) = 0.[72] Here, (cid:104)·(cid:105) is the ensemble average in the case of randomswitching ( α = r ), and the average obtained by inte-grating over T , i.e. (cid:104) K ( t ) (cid:105) = (1 /T ) (cid:82) t + Tt K ( τ ) dτ underperiodic switching.[73] After the exinction of one species, the other has a logistic-like dynamics, see text, with a mean time to extinctionscaling exponentially with its carrying capacity [23, 74].[74] M. Assaf and B. Meerson, Phys. Rev. E. 81, 021116(2010).[75] S. Redner, A Guide to First-Passage Processes , (Cam-bridge, New York, 2001).[76] When N and x are not coupled and s (cid:28)
1, the initialcondition N (0) is here irrelevant. We set N (0) = K or N (0) = K (1 − γ ) / (1 − γδ ) = K and confirmed that ourresults were independent of N (0).[77] K. Kitahara, W. Horsthemke, and R. Lefever, Phys. Lett. , 377 (1979).[78] M. H. A. Davis, J. R. Stat. Soc. B , 353 (1984).[79] C. R. Doering and W. Horsthemke, J. Stat. Phys. ,763 (1985).[80] V. Elgart and A. Kamenev, Phys. Rev. E 70, 041106(2004).[81] L. D. Landau and E. M. Lifshitz, Mechanics (Pergamon,Oxford, 1976).[82] This occurs for s > s c (generic case), when φ α ( ν ) is adecreasing monotonic function. On the other hand, φ r (cid:46) φ p when s < s c and γ < γ c and δ < δ c , see below andSec. S5.2 in [64].[83] In Sec. S2.2, we first consider a general periodic function ξ p ( t ) = ξ p ( t + T ). In all the other sections of the maintext and this Supplemental Material, we have specificallyfocused on the case of a rectangular wave of period T andduty cycle (1 + δ ) /
2, see Sec. S5.4.[84] D. T. Gillespie, J. Comput. Phys. , 403 (1976).[85] D. F. Anderson, J. Chem. Phys. , 214107 (2007).[86] M. Dykman, E. Mori, J. Ross, and P. Hunt, J. Chem.Phys. , 5735 (1994).[87] C. M. Bender and S. A. Orszag, Advanced Mathemat-ical Methods for Scientists and Engineers (New York:Springer, 1999).[88] When b >
0, the average number of switches prior tofixation is O ( ν/s ) as in the cases δ = 0 [45] and b = q = 0, see Sec. S6 and Fig. S3(c). Appendix: Supplementary Material to
Population Dynamics in a Changing Environment:Random versus Periodic Switching
In this Supplemental Material, we provide some further technical details and supplementary information in supportof the results discussed in the main text. We also provide additional information concerning the population’s meanfixation time (MFT), and the generalization of the model in a scenario where the slow strain is a public goods producer.In what follows, unless stated otherwise, the notation is the same as in the main text and the equations andfigures refer to those therein. This document and additional supporting resources are available at the following URL:https://doi.org/10.6084/m9.figshare.12613370.
S1. MODEL DESCRIPTION, MASTER EQUATION AND SIMULATION METHODS
In this section, we describe in detail the model and discuss our modelling choices. We then give the master equation(ME) of the birth-death process according to which the population evolves, and describe the methods used to simulatethe population dynamics in the case of random and periodic switching.
S1.1. Model description
As explained in the main text, the population evolves according to a multivariate birth-death process where repro-duction of
S/F individuals, N S/F → N S/F + 1, occurs at a transition rate T + S/F , and death N S/F → N S/F −
1, occursat a transition rate T − S/F , with [44, 45] T + S = f S ¯ f N S , T + F = f F ¯ f N F and T − S = NK ( t ) N S , T − F = NK ( t ) N F . (S1)In the main text we explicitly consider f S = 1 − s and f F = 1, with 0 < s (cid:28)
1, yielding the population’s average(relative birth) fitness ¯ f = ( N S f S + N F f F ) /N = 1 − sx , where x ≡ N S /N is the fraction of S individuals (slowgrowers). In the transition rates (S1), the carrying capacity K ( t ) varies in time either randomly or periodicallyaccording to Eq. (1) of the main text, and switches with rates ν ± , see also below. It is worth noting that ourchoice of f i , i ∈ { S, F } sets the typical time scale of the dynamics. In a more general setting, the biological factorsdetermining the per capita growth and death rates can be written as the product of a global and relative terms: T + i = g ( x, N ) f i ( x ) N i / ¯ f and T − i = d ( x, N ) w i ( x ) N i / ¯ w , where ¯ w = ( N S w S + N F w F ) /N . In this formulation, g ( x, N )and d ( x, N ) are respectively referred to as the global birth fitness and global weakness and are species independent(acting similarly on both strains), whereas f i ( x ) and w i ( x ) are the species-dependent relative birth fitness and relativeweakness, respectively [50, 51]. In this general setting, g and f i affect the strains’ birth rates, while d and w i determinetheir survival or viability. Within this framework, various evolutionary scenarios can be investigated, see below andRefs. [44–46, 50–53].In this work, as in many applications, see e.g. Refs. [50–53], we have assumed that S and F (slow and fast growers)have equal survival chances and are subject to a logistic growth, and hence we set w S = w F = 1, and d ( x, N ) = N/K for the global weakness. For the sake of simplicity, we have assumed that the relative birth fitness (referred to as“fitness” for brevity) is constant for each species, with f S = 1 − s and f F = 1, while the global birth fitness is g = 1 in the main text, where we focus on the “pure resource competition scenario” of Refs. [44, 45]. In Sec. 7 of thisSupplemental Material (SM), we also consider a “public good scenario” in which the slow growers (strain S ) are publicgood (PG) producers, and the global growth birth fitness (global growth rate) is g ( x ) = 1 + bx (with b > x of S individuals inthe population. This choice corresponds to the “balanced growth scenario” considered in Refs. [50–53] with a constantcarrying capacity. In such a scenario, birth and death events balance each other, and the population size fluctuatesabout its carrying capacity after a short transient. Interestingly, the “balanced growth scenario” (with PG production, b >
0) has been used in Ref. [52] to explain the Simpson’s paradox found in the microbial experiments of Ref. [55]. Thisframework also allows us to model the effect of bacteriostatic (biostatic) and bactericidal (biocidal) antimicrobials onthe time evolution of sensible microorganisms in communities of sensible and resistant cells: bacteriostatic suppressessensible cells growth, and hence affects f i (but neither d nor w i ), while bactericidal induces sensible cells death andthus affects w i (but neither g nor f i ), see, e.g. , Refs. [26, 58].While different other model formulations are of course possible, studying the birth-death process defined by Eqs. (S1)and (1) of the main text, is arguably the simplest way to investigate analytically, in a biologically simple and relevantsetting, the effect of demographic noise (random birth/death events) coupled to environmental variability. Namely,this coupling is achieved via the switching carrying capacity that drives the dynamics of the population size. At thispoint, it is useful to summarize the main properties of the birth-death process defined by Eqs. (S1) and (1):- As reported in the main text, at mean-field level (constant K = K (cid:29)
1, large population), the population sizeobeys the logistic equation dN/dt ≡ ˙ N = (cid:80) i ( T + i − T − i ) = N [1 − ( N/K )], while the population compositionevolves according to the replicator-like equation [68] ˙ x = ( T + S − T − S ) /N − x ( ˙ N /N ) = − x (1 − x )[ f S − ¯ f ] / ¯ f = − sx (1 − x ) / (1 − sx ) [44, 45]. This model, and its generalization (see Sec. 7 of this SM), therefore have a soundeco-evolutionary dynamics.- When the population size is constant ( N = K ) and there is no environmental variability (only demographicnoise), the dynamics can be mapped onto that of the well-known fitness-dependent Moran model [48, 66, 68]defined by the reactions SF → SS and SF → F F , respectively occurring at rates (cid:101) T + S = T + S T − F /N = (1 − s ) x (1 − x ) N/ (1 − sx ) and (cid:101) T − S = T − S T + F /N = x (1 − x ) N/ (1 − sx ), see Ref. [45]. This allows us to obtain Eq. (2)in the main text, used in Eq. (3) to compute the fixation probability when K varies in time.- The model studied here is conceptual, but many of its features are biologically relevant . With modern bio-engineering techniques, it is in fact possible to perform controlled microbial experiments in settings allowingto test the theoretical predictions of models featuring switching environment, time-varying population size, PGproduction, cooperation dilemma, see, e.g. , Refs. [28, 55, 57, 63].- The birth-death process underpinning this model can generalized in different ways. In addition to the scenariowith PG production, see above and Sec. 7 of this SM, a possible generalization is the “dormancy scenario” ofRef. [51] where g ( x, N ) = 1 + x − ( N/K ) , d ( x, N ) = 0 and same f i , w i as here. The above general frameworkcan also accommodate more realistic and complex processes in which g ( x, N ) and f i ( x ), and/or w i ( x ), dependon ξ α , with α ∈ { r, p } and hence, also vary with the environment along with d = N/K . S1.2. Master equation of the underlying birth-death process
Using (cid:126)N = ( N S , N F ) and ± as a shorthand notation for ξ r = ±
1, the ME for the birth-death process definedby (S1), where the carrying capacity K ( t ) varies randomly by switching according to K + → K − with rate ν + and K − → K + with rate ν − [see Eq. (1) in the main text, with α = r ], reads dP ( r ) ν ( (cid:126)N , + , t ) dt = ( E − S − T + S P ( r ) ν ( (cid:126)N , + , t )] + ( E − F − T + F P ( r ) ν ( (cid:126)N , + , t )] (S2a)+ ( E + S − T − S P ( r ) ν ( (cid:126)N , + , t )] + ( E + F − T − F P ( r ) ν ( (cid:126)N , + , t )] + ν − P ( r ) ν ( (cid:126)N , − , t ) − ν + P ( r ) ν ( (cid:126)N , + , t ) ,dP ( r ) ν ( (cid:126)N , − , t ) dt = ( E − S − T + S P ( r ) ν ( (cid:126)N , − , t )] + ( E − F − T + F P ( r ) ν ( (cid:126)N , − , t )] (S2b)+ ( E + S − T − S P ( r ) ν ( (cid:126)N , − , t )] + ( E + F − T − F P ( r ) ν ( (cid:126)N , − , t )] + ν + P ( r ) ν ( (cid:126)N , + , t ) − ν − P ( r ) ν ( (cid:126)N , − , t ) , where E ± S/F are shift operators such that E ± S f ( N S , N F , ξ, t ) = f ( N S ± , N F , ξ, t ) and similarly for E ± F . Clearly,Eqs. (S2a) and (S2b) are coupled and the terms on the 2nd lines’ right-hand-side account for environmental switching.For periodic switching, the carrying capacity K ( t ) = K [1 + γξ p ( t )] , varies deterministically with ξ p ( t ) ≡ ξ p ( t + T ),where the shape of ξ p ( t ) is taken to be a rectangular wave of period T = (1 /ν + ) + (1 /ν − ) [83]. In this case, the MEof the birth-death processs (S1) with periodically switching K ( t ) reads dP ( p ) ν ( (cid:126)N , t ) dt = ( E − S − T + S P ( p ) ν ( (cid:126)N , t )] + ( E − F − T + F P ( p ) ν ( (cid:126)N , t )]+ ( E + S − T − S ( ξ p ) P ( p ) ν ( (cid:126)N , t )]] + ( E + F − T − F ( ξ p ) P ( p ) ν ( (cid:126)N , t )] , (S3)where T − S/F ( ξ p ) are now the time-dependent transition rates given by (S1) that vary periodically with ξ p . Note, thatin both MEs (S2)-(S3), P ( α ) ν ( (cid:126)N , t ) = 0 whenever N S < N F < S1.3. Simulation methods
While the MEs (S2) and (S3) fully describe the population dynamics in the case of random and periodic switching,respectively, in general, they cannot be solved analytically. However, to gain insight into to the stochastic dynamics,one can employ efficient numerical simulations. In the case of random switching, process (S2) defined by the birth-death (S1) and switching ξ r → − ξ r reactions, can be exactly simulated using the standard Gillespie algorithm [84]. Inthe case of periodic switching, it is convenient to simulate the birth-death process (S3) with time-dependent (periodic)transition rates (S1) using the simulation method outlined below. S1.3.1. Simulation of the periodic switching case with the modified next reaction method
In the periodic case we used the modified next reaction method [85], which is a suitable algorithm for systems withexplicit time dependent rates. Unlike the classic Gillespie Algorithm, this version considers all possible birth/deathprocesses as independent reactions. We can calculate the time step ∆ t i in which the next reaction occurs by generatinga random number from a uniform distribution r i ∈ U (0 ,
1) for the probability that reaction i did not occur after timeinterval ∆ t i . Here, we have four stochastic reactions i ∈ { , . . . , } (birth/death of S and F ) each with a propensityfunction a i ∈ { T + S , T − S ( ξ p ( t )) , T + F , T − F ( ξ p ( t )) } , and thus we have r i = exp (cid:104) − (cid:82) t +∆ t i t a i ( t (cid:48) ) dt (cid:48) (cid:105) .We start the simulation at time t = 0, and for each reaction we set the “internal time” T i = 0 and the quantity P i = ln (1 /r i ). We also set the initial number of each species, the environmental state (with probability determinedby the duty cycle), and the initial time to the next switch ∆ t switch . Here, the time step ∆ t i is found by computing (cid:82) t +∆ t i t a i ( t (cid:48) ) dt (cid:48) = P i − T i , which can be easily solved, since K is discrete and thus in each iteration it is constant.At this point we find the reaction that has the minimal time step ∆ t µ = min i { ∆ t i } , propagate time t → t + ∆ t µ ,and update the population size, the internal times T i → T i + (cid:82) t +∆ t µ t a i ( t (cid:48) ) dt (cid:48) , and P i → P i + δ i,µ ln (1 /r i ). Then werecalculate the rates a i , generate another random number r i ∈ U (0 , ξ → − ξ , that occurred during a period of1 /ν ± , as follows: if ∆ t switch < ∆ t µ , we switch ξ → − ξ and propagate the time t → t + ∆ t switch . S2. APPROXIMATIONS OF THE QUASI-STATIONARY POPULATION DENSITY: PERIODICSWITCHING
In this section we compute the quasi-stationary population size distribution (PSD), P N , in the slow switchingregime, as well as under fast and intermediate periodic switching. This is done by first computing the PSD in the caseof constant carrying capacity, assuming a static environment ξ α ( t ) = ξ and carrying capacity K ( t ) = K . To do so,we start with the ME for P ( N ) | K – the probability that the total population size is N given a carrying capacity KdP ( N ) | K dt = ( N − P ( N − | K + ( N + 1) K P ( N + 1) | K − (cid:18) N + N K (cid:19) P ( N ) | K . (S4)The PSD can be found by putting ˙ P ( N ) | K = 0 and demanding a reflecting boundary condition at N = 1. The latterassumes that the probability flux to the extinction state P ( N = 0) is negligibly small, which is legitimate since themean time to extinction is assumed to be much larger than the time scales we are interested in here, see main text.The normalized solution of the resulting recursion equation reads P ( N ) | K = 1Ei ( K ) − γ − ln ( K ) K N N N ! (cid:39) KN K N e − K N ! , (S5)where Ei ( x ) = − (cid:82) ∞− x dte − t /t is the exponential integral function, γ = 0 . ... is the Euler–Mascheroni constant, andthe last approximation holds when K (cid:29) S2.1. Quasi-stationary PSD under slow switching
When ν →
0, on average there are no switches prior to fixation, and the population evolves in a static environment ξ α = ±
1, with ξ α that is distributed with a probability p ( ξ α ) = (1 ± δ ) /
2. Namely, if ξ α = ±
1, the population is0subject to a constant carrying capacity K = K ± . Hence, using Eq. (S5), the PSD under slow switching reads P ( N ) = (cid:88) ξ α = ± P ( N | ξ α ) p ( ξ α ) (cid:39) (cid:18) δ (cid:19) K + N K N + e − K + N ! + (cid:18) − δ (cid:19) K − N K N − e − K − N ! . (S6)As explained in the main text, this result is valid both for periodic and random switching. S2.2. Quasi-stationary PSD under fast periodic switching: Kapitza method
In the opposite limit ν (cid:29)
1, the carrying capacity K rapidly oscillates around K . To find the PSD in the case offast periodic switching, we employ the Kapitza method [12], valid for a general periodic ξ p ( t ), which involves separatingthe dynamics into fast and slow variables, and averaging the fast variables over the period of variation.Our starting point is ME (S4), but now with K = K [1 + γξ p ( t )], i.e., explicitly time-dependent rates. To treatEq. (S4) semi-classically, we define the probability generating function G ( p, t ) = (cid:80) ∞ m =0 P ( m, t ) p m , where p is anauxiliary variable. Conservation of probability yields G (1 , t ) = 1. The definition of G is useful since P ( N, t ) = 1 N ! ∂ N G ( p, t ) ∂p N | p =0 . (S7)Multiplying Eq. (S4) by p N and summing over all N ’s, we obtain a second-order partial differential equation for G ∂G∂t = (1 − p ) (cid:18) − p ∂G∂p + 1 K ( t ) ∂G∂p + pK ( t ) ∂ G∂p (cid:19) . (S8)This equation cannot be solved in general. An approximate solution can be found by using the fact that the typicalcarrying capacity is large, K (cid:29)
1, and employing the WKB ansatz G = G exp [ − K S ( p, t )] in Eq. (S8) [23]. Keepingleading- and subleading-order terms with respect to O ( K ), we arrive at the following Hamilton-Jacobi equation − ∂S∂t = q (1 − p ) (cid:18) − p + 1 K ( t ) + K K ( t ) pq (cid:19) ≡ H ( p, q ) , (S9)where H is the Hamiltonian, S is the action associated with the Hamiltonian, and we have defined q = − ∂S∂p as thecoordinate conjugate to the variable p , see [86].Let us separate the fast and slow time scales by denoting q ( t ) = X ( t ) + ζ ( t ) and p ( t ) = Y ( t ) + η ( t ). Here X and Y are slow variables, while ζ and η are small corrections (to be verified a-posteriori) that rapidly oscillate around0 [12]. Expanding the Hamiltonian (S9) up to second order around q = X and p = Y we find H ( q, p, t ) (cid:39) H ( X, Y, t ) + ζ ∂H∂X + η ∂H∂Y + ζ ∂ H∂X + η ∂ H∂Y + ζη ∂ H∂X∂Y ≡ ˜ H ( X, Y, t ) . Using the Hamilton equations ˙ q = ˙ X + ˙ ζ (cid:39) ∂ Y ˜ H ( X, Y, t ) and ˙ p = ˙ Y + ˙ η (cid:39) − ∂ X ˜ H ( X, Y, t ), and equating the rapidlyoscillating terms yields in the leading order in K (cid:29) ζ (cid:39) (cid:0) X − X Y (cid:1) ( B/ν ) and η (cid:39) − (cid:0) XY − XY (cid:1) ( B/ν ).Here B ( t ) = O (1) is defined in Eq. (S12), and in the calculation X and Y were considered as constants during theperiod of rapid oscillations. In addition, we have neglected terms of order ζ and η , but kept their time derivatives(proportional to ν (cid:29) q, p ) to the new ( X, Y ) variable q (cid:39) X + X (1 − Y ) Bν + 2 X (1 − Y ) B ν , p (cid:39) Y − (cid:0) Y − Y (cid:1) X Bν − X (cid:0) Y − Y + 2 Y (cid:1) B ν , (S10)which can be obtained using the generating function F ( q, Y, t ) = qY − q (cid:0) Y − Y (cid:1) ( B/ν ). This transformation iscanonical up to second order in the small parameter 1 /ν , as the Poisson brackets satisfy { q, p } ( X,Y ) = 1 + O (cid:0) ν (cid:1) . Using Eqs. (S9) and (S10) and defining H (cid:48) = H + ∂F ∂t , by averaging over a period of a rapid oscillation, we find H ( X, Y ) = XY (1 − Y ) (cid:20) AX − X C (cid:0) − Y + 4 Y − XA (cid:1) ν (cid:21) , (S11)1where we have defined the following O (1) variables A ( γ, δ ) = K K ( t ) = 1 T (cid:90) t + Tt K K ( t ) dt ; C ( γ, δ ) = B = 1 T (cid:90) t + Tt B ( t ) dt ; B ( t ) = K ν (cid:90) dt (cid:20) K ( t ) − K ( t ) (cid:21) , (S12)and used the fact that B ( t ) is periodic. It can be shown that the O ( ν ) terms in Hamiltonian (S11) yield the PSDin the constant environment case [Eq. (S5)].Having found the time- independent Hamiltonian (S11), which effectively takes into account the rapid environmentaloscillations, we can compute the PSD by finding the nontrivial zero-energy trajectory of H ( X, Y ). Up to second orderin ν , this trajectory is given by X ( Y ) = 1 /A − C/ ( A ν ) (2 Y − + O (cid:0) ν − (cid:1) . Thus, recalling that q = − ∂S∂p , and usingthe fact that the transformation ( q, p ) → ( X, Y ) is canonical, we find S ( Y ) = − (cid:82) XdY = − Y /A + C/ (6 A ν )(2 Y − .As a result, the generating function becomes G ( Y ) (cid:39) G exp [ − K S ( Y )] = G exp (cid:20) K YA − K C A ν (2 Y − (cid:21) , where G is a constant, see below. Therefore, the PSD can be found by employing the Cauchy theorem to Eq. (S7): P ( N ) = 12 πi (cid:122) G ( Y ) Y N +1 dY = G πi (cid:122) Y exp [ K g ( Y )] dY, where the integration has to be performed over a closed contour in the complex Y plane around the singular point Y =0, and we have defined g ( Y ) = − S ( Y ) − NK ln Y . This integral can be calculated using the saddle point approximation[87]. The saddle point, up to second order in 1 /ν , is found at Y ∗ = AN/K + [ CN/ ( AK )] (2 AN/K − ν − .Furthermore, since g (cid:48)(cid:48) ( Y ∗ ) > P ( N ) ∼ = (cid:104) G / (cid:16) Y ∗ (cid:112) πK | g (cid:48)(cid:48) ( Y ∗ ) | (cid:17)(cid:105) e K g ( Y ∗ ) . Note,however, that only the leading-order result can be taken into account here; accounting for the prefactor would be anexcess of accuracy since we have ignored the p -dependent prefactors in both G and in P . Putting it all together, wefinally obtain P ( N ) (cid:39) C exp (cid:20) N − N ln ( AN/K ) − K C A ν (2 AN/K − (cid:21) , (S13)where C is a normalization constant which can be found by demanding (cid:82) P ( N ) dN = 1. S2.2.1. Rectangular wave
Our derivation above has been carried out for a general periodic function ξ p ( t ). We now compute the PSD in theparticular case of a rectangular wave. Using the expression of ξ p ( t ) given in the main text, we find B ( t ) = γ − γ × (cid:40) ( δ − νt − / − ν + ≤ t ≤ δ + 1) νt − / ≤ t ≤ ν − , (S14)where the constant of integration was determined by the demand that B = 0. Plugging this into Eq. (S12) yields K ( t ) − = A/K = (1 − γδ ) / [ K (cid:0) − γ (cid:1) ] and C = (1 / γ / (cid:0) − γ (cid:1) . Using these results, Eq. (S13) becomes P Kap ν (cid:39) P ( N )exp (cid:34) − K ν (cid:18) γ − γ (cid:19) (cid:18) N −K K (cid:19) (cid:35) , (S15)which is the expression of P Kap ν used in the main text, with lim ν →∞ P Kap ν = P ( N ) ∝ exp[ N (1 − ln ( N/ K ))], peakedat K = K (1 − γ ) / (1 − γδ ). Hence, P Kap ν ( N ) is unimodal and peaked about N ≈ K when ν (cid:29)
1, see Fig. 2(b).
S2.3. Quasi-stationary PSD under intermediate periodic switching
We now consider the quasi-stationary PSD in the regime of intermediate periodic switching where ν = O (1). Inthis regime, progress can be made upon neglecting demographic noise, and by considering only the environmental2periodic modulation for N ( t ). This leads to an approximation of P ( p ) ν , here referred to as “piecewise periodic process”and denoted by P PPP ν , that is the periodic counterpart of the PDMP approximation, see Eq. (S22) and below. Thisapproach is similar in nature to that of Refs. [70, 79] (whose focus was on symmetric switching).Our starting point is the mean-field rate equation for the total population size, in the case of periodic switching,upon ignoring demographic noise. Using the definition of ξ p ( t ) from the main text, the equation reads˙ N = N (cid:20) − NK ( t ) (cid:21) = N (cid:26) − NK [1 + γξ p ( t )] (cid:27) , (S16)At t → ∞ , after the transient has decayed, the periodic solution reads N = f ( t ) = K (cid:0) − γ (cid:1) (cid:104) − γ + 2 γe − ˜ t − e − /ν − − e − T (cid:105) − < ˜ t < ν + ( ξ = 1) K (cid:0) − γ (cid:1) (cid:104) γ − γe − ˜ t e /ν + − − e − T (cid:105) − ν + < ˜ t < T ( ξ = − , (S17)where ˜ t = t − /ν − − (cid:98) t − /ν − T (cid:99) T , such that 0 ≤ ˜ t ≤ T = (1 /ν + ) + (1 /ν − ). As a result, for each segment of the solution,one can express ˜ t as function of N : ˜ t = g ± ( N ) = ln ( B ± ) + ln ( N ) − ln (1 − N/K ± ). Here, B ± is a cumbersomeexpression independent on N and hence irrelevant for our purposes, whereas the subscripts + and − stand for thefirst and second segment in each period, respectively. Therefore, we can approximate the PSD, P PPP ν , as P PPP ν ( N ) = (cid:90) T d ˜ t (cid:48) P PPP ν (cid:0) N, ˜ t (cid:48) (cid:1) , (S18)where P PPP ν (cid:0) N, ˜ t (cid:48) (cid:1) ∼ δ (cid:2) N − f (cid:0) ˜ t (cid:1)(cid:3) is the probability that the population size at time ˜ t is N , and we have omittedthe normalization constant. Here we have neglected demographic noise by assuming that the instantaneous totalpopulation size N is sharply peaked around its deterministic solution. Performing the integral in Eq. (S18) , we find P PPP ν ( N ) = C (cid:20)(cid:12)(cid:12)(cid:12)(cid:12) dg + dN (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) dg − dN (cid:12)(cid:12)(cid:12)(cid:12)(cid:21) = C (cid:20) K + − N + 1 N − K − (cid:21) , (S19)where C is a normalization constant. This expression is valid for N min ≤ N ≤ N max , where the boundaries N min = N (˜ t = 0) and N max = N (˜ t = 1 /ν + ) satisfy N min = K (cid:0) − γ (cid:1) (cid:20) − γ +2 γ − e − /ν − − e − T (cid:21) − , N max = K (cid:0) − γ (cid:1) (cid:20) − γ +2 γ e − /ν + − e − T − e − T (cid:21) − , (S20)while the normalization constant is given by C − = ln (cid:20) ( K + − N min ) ( N max − K − )( K + − N max ) ( N min − K − ) (cid:21) . (S21)The PPP approximation P PPP ν of the periodic PSD is shown in Fig. 2(c,d) of the main text, where it is found to agreewell with the simulation results and to reproduce the main features of P ( p ) ν in the intermediate switching regime. Italso accurately captures the average population size, as shown in Fig. S3(b). S3. QUASI-STATIONARY PSD FOR RANDOM SWITCHING: THE PDMP APPROXIMATION
When demographic noise is neglected, by assuming that the fluctuating population size is always large, and theonly source of noise stems from the randomly switching carrying capacity, we have seen that the PSD, P ( r ) ν ( N ), canbe described in terms of the marginal stationary probability density of the underlying piecewise-deterministic Markovprocess (PDMP). Upon omitting the normalization constant, this PSD reads [45, 69] P PDMP ν ( N ) ∝ N (cid:34)(cid:18) K + N − (cid:19) ν + − (cid:18) − K − N (cid:19) ν − − (cid:35) , (S22)where the dependence on γ, δ and ν is given by K ± = (1 ± γ ) K and ν ± = (1 ∓ δ ) ν . Clearly, P PDMP ν has support[ K − , K + ] and accounts for environmental noise, but ignores all demographic fluctuations. The expression of P PDMP ν gives a suitable description of P ( r ) ν in the intermediate switching regime where interesting phenomena arise (seeSec. S4.3 below for a detailed discussion of the validity of the PDMP-like approximations).3 FIG. S1: Phase diagram for the PSD, P ( r ) ν ( N ), and its approximations P PDMP ν and P LNA ν (insets), see Eq. (S26). We distinguishfour regions described in the text: In addition to a peak about K + , the PSD has always a local maximum K − < N ∗ < K + in the intermediate switching regime in I; in regime II and III, the PSD and P PDMP ν have a peak about K + and, dependingon ν , possibly another peak at some values K − < N ∗ < K + , see insets; the PSD and P PDMP ν have one single peak about K + in IV. Insets illustrate the form of P ( r ) ν ( N ), P PDMP ν and P LNA ν in regions I-III. In the insets, solid lines are from the P PDMP ν , given by Eq. (S22), dashed lines are from P LNA ν , given by Eq. (S26), solid areas are from computer simulations, andthe vertical dashed lines are eyeguides showing N = K ± . Parameters are: ( K , γ, s, x ) = (250 , . , . , .
5) and (inset I) δ = 0 . ν = (0 . , . , .
5) (pink, orange, blue); (inset II) δ = 0 . ν = (1 , , .
5) (purple, blue, green); (inset III) δ = 0 . ν = (1 , ,
12) (purple, blue, green). In inset I, N ∗ is in the intermediate regime for ν = 1 . N ∗ is in theintermediate regime for ν = 1 (purple) and ν = 6 . N ∗ is in the intermediate regime for ν = 1 (purple).We notice that the LNA excellently agrees with simulation results for the PSD: P ( r ) ν ( N ) and P LNA ν are almost indistinguishablein each inset. S3.1. PSD dependence on γ and δ in the intermediate switching regime The PSD, P ( r ) ν ( N ), and its PDMP approximation, P PDMP ν , are bimodal, with peaks about K ± , when ν <
1, andunimodal when ν > N ∗ that is the smaller solution to N − ( ν (1 − γδ ) + 1) K N + (1 − γ ) K ν = 0 , (S23)with N ∗ → K as ν → ∞ [44–46]. In addition, two other regimes can arise under asymmetric switching at intermediate rate when 1 / (1 + | δ | ) < ν < / (1 − | δ | ). Here, the PSD has a different form not found when δ = 0: When δ < / (1 − δ ) < ν < / (1 + δ ), P PDMP ν and P ( r ) ν have a peak at N (cid:39) K − . When δ > / (1 + δ ) < ν < / (1 − δ ), P PDMP ν and P ( r ) ν have a peak at N (cid:39) K + and, depending on δ, γ and ν , also a peak at N ∗ . The condition for theexistence of such a peak at K − < N ∗ < K + can be inferred from the PDMP approximation (S22) by noting that(S23) has real roots when (1 − γδ ) ν − γ ( δ − γ )) ν + 1 > . (S24)We thus distinguish four regions, I-IV, in the ( δ, γ ) - space, see Fig. S1:I: δ < γ , where N ∗ exists for all intermediate ν .II: γ < δ < γ γ , where N ∗ exists for all intermediate ν that lie outside the interval between the two solutions of(S24), here denoted by ν , (with ν ≥ ν ).4III: γ γ < δ < γ − γ , where N ∗ only exists if δ < ν < ν .IV: δ > γ − γ , where N ∗ does not exist.Simulation results of Figs. 2 and S1 confirm that the above analysis correctly reflects the properties of P ( r ) ν , see thevideos of the Figshare resources [64].As shown by Fig. 2 of the main text, the PSD under intermediate periodic switching is qualitatively characterizedby the same features as P ( r ) ν , with some generic quantitative differences: P ( p ) ν is generally narrower and has sharperpeak than P ( r ) ν . All these features are well captured by the PPP approximation (S19)-(S21) of P ( p ) ν . In particular, P PPP ν has a narrower support [ N min , N max ] than the support [ K − , K + ] of P PDMP ν , since N min > K − and N max < K + ,see Fig. 2 (c,d) and the Figshare resources of [64]. S3.2. Linear noise approximation about the PDMP solution
While P PDMP ν (S22) captures well the position of the peaks of the PSD and some of its main features, the PDMPapproximation fails to capture the width of P ( r ) ν ( N ). In order to account for the demographic noise responsible forthe shape of P ( r ) ν ( N ) near its peaks, we can perform a linear noise approximation (LNA) about the PDMP [45] ddt N = N (cid:20) − N K (cid:18) − γξ r − γδ (cid:19)(cid:21) , (S25)whose probability density P PDMP ν ( N, ξ r ) in the environmental state ξ r = ±
1, is given by [69] P PDMP ν ( N, ξ r ) ∝ δN (cid:104) K + N − (cid:105) ν + − (cid:104) − K − N (cid:105) ν − , ( ξ r = +1) − δN (cid:104) K + N − (cid:105) ν + (cid:104) − K − N (cid:105) ν − − , ( ξ r = − , where K ± = (1 ± γ ) K and ν ± = (1 ∓ δ ) ν . As in Refs. [41, 45], we also make the simplifying assumption thatdemographic noise is approximately the same in each environmental state, yielding the Gaussian distribution ∝ exp (cid:16) − ( N − (cid:101) N ) / (2 (cid:101) N ) (cid:17) / (cid:112) (cid:101) N for the demographic fluctuations N − (cid:101) N about the PDMP (S25). Proceeding as in thecase of symmetric switching ( δ = 0), see Ref. [45] where full details are provided, and omitting the normalizationconstant, we obtain the LNA of the marginal stationary probability density about the PDMP (S25) P LNA ν ( N ) ∝ (cid:90) K + K − e − ( N − (cid:102) N ) (cid:102) N (cid:101) N / (cid:40) (1+ δ ) (cid:20) K + (cid:101) N − (cid:21) ν + − (cid:20) − K − (cid:101) N (cid:21) ν − +(1 − δ ) (cid:20) K + (cid:101) N − (cid:21) ν + (cid:20) − K − (cid:101) N (cid:21) ν − − (cid:41) d (cid:101) N. (S26)The results shown in the insets of Fig. S1 illustrate that P LNA ν ( N ) is an excellent approximation of the PSD: itaccurately predicts all the details of the PSD P ( r ) ν ( N ) obtained from stochastic simulations. However, while P LNA ν significantly improves over P PDMP ν to describe the PSD, we have verified that computing φ r in the realm of thePDMP-based approximation [i.e. with Eq. (S38)] or by averaging φ ( x ) | N over P LNA ν , as an approximation of P ( r ) ν ,according to Eq. (3), yields essentially the same results: As shown in Fig. S2(c), the fixation probability calculatedusing P LNA ν gives only a minute improvement over the results obtained with P PDMP ν . The LNA approximation (S26)is thus useful to describe the PSD, but the PDMP approximation is sufficient to compute the fixation probability.Note, that while we have not carried it out explicitly, a similar LNA treatment can be done in the periodic case.This would allow us to accurately reproduce the PSD in the low and intermediate periodic switching regime. S4. FIXATION PROBABILITY UNDER FAST SWITCHING: SADDLE-POINT CALCULATIONS
In this section we perform a saddle-point approximation to find the fixation probability, φ α , in the fast switchingregime ν/s (cid:29)
1, and then discuss the validity of the PDMP-like (PDMP and PPP) approximations.To perform a saddle-point calculation of φ α under fast switching, we rewrite Eq. (3) of the main text in terms ofthe total population density y = N/K . Accounting for the normalization of the probability distribution, the fixation5probability can be written as φ α ( ν ) = (cid:82) ∞ P ( α ) ν/s ( y ) exp [ K (1 − x ) ln (1 − s ) y ] dy (cid:82) ∞ P ( α ) ν/s ( y ) dy ≡ (cid:82) ∞ exp (cid:104) f ( α )num ( y ) (cid:105) dy (cid:82) ∞ exp (cid:104) f ( α )den ( y ) (cid:105) dy , (S27)where we have defined f ( α )den ( y ) = ln P ( α ) ν/s ( y ), and f ( α )num ( y ) = f ( α )den ( y ) + K (1 − x ) ln (1 − s ) y , and α denotes either r (random) or p (periodic). Evaluating both integrals separately via the saddle point approximation, we obtain φ α ( ν ) (cid:39) (cid:113) κ ( α )1 /κ ( α )2 e f ( α )num (cid:16) y ( α )2 (cid:17) − f ( α )den (cid:16) y ( α )1 (cid:17) . (S28)Here y ( α )1 and y ( α )2 are the positions of the saddle points of the denominator and numerator, respectively, andsatisfy ( d/dy ) f ( α )den (cid:16) y ( α )1 (cid:17) = 0 and ( d/dy ) f ( α )num (cid:16) y ( α )2 (cid:17) = 0. In addition, κ ( α )1 = ( d /dy ) f ( α )den (cid:16) y ( α )1 (cid:17) and κ ( α )2 =( d /dy ) f ( α )num (cid:16) y ( α )2 (cid:17) represent the curvatures at the saddle point of the denominator and numerator, respectively. S4.1. Fast random switching
Here we compute Eq. (S28) in the case of randomly switching environment in the realm of the PDMP approximation,with P ( r ) ν/s (cid:39) P PDMP ν/s . To compute the denominator of Eq. (S27), with Eq. (S22), we define f ( r )den ( y ) = ln P PDMP ν/s ( y ) = − νs ln y + (cid:104) (1 − δ ) νs − (cid:105) ln (1 + γ − y ) + (cid:104) (1 + δ ) νs − (cid:105) ln ( y − γ ) . (S29)Thus, the saddle point is found at y ( r )1 (cid:39) (cid:0) − γ (cid:1) (1 − δγ ) (cid:34) γ ( δ − γ )(1 − δγ ) ν/s (cid:32) (cid:0) − γ + δγ (cid:1) (1 − δγ ) ν/s (cid:33)(cid:35) . As a result, we find f ( r )den (cid:16) y ( r )1 (cid:17) (cid:39) ( ν/s ) (cid:26) (1 + δ ) ln (cid:20) γ (1 + δ ) (1 − γ )(1 − δγ ) (cid:21) + (1 − δ ) ln (cid:20) γ (1 − δ ) (1 + γ )(1 − δγ ) (cid:21) − (cid:20) − γ − δγ (cid:21)(cid:27) + ln (cid:34) (1 − δγ ) γ (1 − δ ) (1 − γ ) (cid:35) + ( δ − γ ) (1 − δ ) (1 − δγ ) ν/s ,κ ( r )1 = d dy f ( r )den (cid:16) y ( r )1 (cid:17) (cid:39) − − δγ ) ν/s (1 − δ ) γ (1 − γ ) + 2 (1 − δγ ) (cid:0) δγ − δ γ − γ − δ (1 − γ ) (cid:1) (1 − δ ) (1 − γ ) γ . (S30)To compute the numerator of (S28) we define f ( r )num ( y ) = f ( r )den ( y ) + K ln(1 − s ) (1 − x ) y , and find the saddle point at y ( r )2 (cid:39) (cid:0) − γ (cid:1) (1 − δγ ) (cid:40) γ (cid:2) − δγ ) ( δ − γ ) + bγ (cid:0) − γ (cid:1) (cid:0) − δ (cid:1)(cid:3) − δγ ) ν/s (cid:20) − γ (cid:0) δ − δγ (cid:1) − bγ (1 − γ ) (cid:0) δ − γ + δ γ (cid:1) − δγ ) ν/s (cid:21)(cid:41) , b = K (1 − x ) ln (1 − s ). As a result, we find f ( r )num (cid:16) y ( r )2 (cid:17) (cid:39) ( ν/s ) (cid:26) (1 + δ ) ln (cid:20) γ (1 + δ ) (1 − γ )(1 − δγ ) (cid:21) + (1 − δ ) ln (cid:20) γ (1 − δ ) (1 + γ )(1 − δγ ) (cid:21) − (cid:20) − γ − δγ (cid:21)(cid:27) + ln (cid:34) (1 − δγ ) γ (1 − δ ) (1 − γ ) (cid:35) + b (cid:0) − γ (cid:1) − δγ + (cid:2) δ − γ ) (1 − δγ ) + bγ (cid:0) − γ (cid:1) (cid:0) − δ (cid:1)(cid:3) − δ ) (1 − δγ ) ν/s ,κ ( r )2 = d dy f ( r )num (cid:16) y ( r )2 (cid:17) (cid:39) − − − δγ ) ν/s (1 − δ ) γ (1 − γ ) − − δγ )(1 − δ ) (1 − γ ) γ (cid:20) (5 − b ) γ + 3 bγ + δ (cid:0) b ) γ − bγ (cid:1) − δ γ (cid:0) b (cid:0) − γ (cid:1)(cid:1) + δ γ (cid:0) − γ − b (cid:0) − γ (cid:1)(cid:1) + δγ (cid:0) b (cid:0) − γ (cid:1) − (cid:0) γ (cid:1)(cid:1) − (cid:21) . (S31) S4.2. Fast periodic switching
Here we compute Eq. (S28) in the case of periodically switching environment using P ( p ) ν/s (cid:39) P Kap ν/s with Eq. (S15).To compute the denominator of Eq. (S28) we define f ( p )den ( y ) = ln P Kap ν/s ( y ) = K (cid:34) y − y ln (cid:18) K K y (cid:19) − ν /s (cid:18) γ − γ (cid:19) (cid:18) y − K K (cid:19) (cid:35) . (S32)Thus, using Eq. (S12) the saddle point is found at y ( p )1 (cid:39) A − (cid:2) − C/ ( A ν ) (cid:3) , where A and C are given in Sec. 2.2.1.As a result, we find f ( p )den (cid:16) y ( p )1 (cid:17) (cid:39) A (cid:18) − C A ν (cid:19) , κ ( p )1 = d dy f ( p )den (cid:16) y ( p )1 (cid:17) (cid:39) − A − CAν . (S33)To compute the numerator of (S28) we define f ( p )num ( y ) = f ( p )den ( y ) + K ln(1 − s ) (1 − x ) y , and find the saddle point at y ( p )2 (cid:39) A − (1 − s ) − x (cid:26) − C/ ( A ν ) (cid:104) − s ) − x − (cid:105) (cid:27) . As a result, we find f ( p )num ( y ( p )2 ) (cid:39) (1 − s ) − x A − C A ν (cid:16) − s ) − x − (cid:17) ,κ ( p )2 = d dy f ( p )num (cid:16) y ( p )2 (cid:17) (cid:39) − A (1 − s ) − x (cid:26) CA ν (cid:104) − s ) − x − (cid:105) (cid:104) − s ) − x − (cid:105)(cid:27) . (S34)Thus, for both random and periodic switching (S28) predicts the same fixation probability φ r = φ p (cid:39) φ ( ∞ ) = φ ( x ) | K , for ν → ∞ . Yet, the asymptotic convergence to φ ( ∞ ) is markedly different [see Eq. (4) in the main text]:ln (cid:18) φ α φ ( ∞ ) (cid:19) = (cid:40) A r ( ν/s ) − for randomly switching environment A p ( ν/s ) − for periodically switching environment, with A r = (1 − x ) ln (1 − s ) K (cid:0) − δ (cid:1) γ − δγ ) (cid:18) − x ) ln (1 − s ) K (cid:19) , A p = K (cid:110) − [1 + 2 (1 − x ) ln (1 − s )] (cid:111) (cid:18) γ − γδ (cid:19) . (S35)These show that φ p ( ν ) approaches φ ( ∞ ) much faster than φ r ( ν ) as ν increases: the convergence towards the fastswitching limit is attained much quicker with periodic than random switching, see Figs. 3(a) and S2(c).7 S4.3. Validity of the PPP and PDMP approximations in the intermediate/fast switching regime
Simulation results show that P PPP ν and P PDMP ν are generally good approximations of P ( p ) ν and P ( r ) ν for a broadrange of ν , from slow to fast switching. We now combine the results of Sections S2 S2.3, S3 S3.2 and S4 S4.1 of thisSM to assess the theoretical validity of the PPP and PDMP approximations, P PPP ν ( N ) and P PDMP ν ( N ), given byEqs. (S19) and (S22), under intermediate/fast switching. This can be done by computing the variance of P PPP ν and P PDMP ν , i.e., σ and σ , and by comparing these results with K , which is the variance of the PSD when, inthe limit ν → ∞ , it is solely governed by demographic noise. Indeed, when ν → ∞ , the PSD [both the Kapitzaapproximation given by Eq. (S15) as well as the LNA given by Eq. (S26)] reduces to a Gaussian of mean and variance K , i.e., P ν →∞ ( N ) ∝ e − ( N −K ) / (2 K ) / √K . To compute the variances in the limit of ν (cid:29)
1, we perform a saddle-pointcalculation as in the previous section, and find to leading order in 1 /ν that σ = (cid:90) N max N min ( N − (cid:104) N (cid:105) PPP ) P PPP ν ( N ) dN = 112 (cid:18) γ − γδ (cid:19) K ν ,σ = (cid:90) K + K − ( N − (cid:104) N (cid:105) PDMP ) P PDMP ν ( N ) dN = 12 (cid:18) γ − γδ (cid:19) (cid:0) − δ (cid:1) K ν . (S36)Here we have used Eqs. (S19) and (S22), while (cid:104) N (cid:105) PPP = (cid:90) N max N min N P
PPP ν ( N ) dN = K (cid:20) γ − γδ ) ν (cid:21) , (cid:104) N (cid:105) PDMP = (cid:90) K + K − N P
PDMP ν ( N ) dN = K (cid:20) γ (1 − δ )2(1 − δγ ) ν (cid:21) , (S37)to leading order in 1 /ν . Notably, from Eqs. (S36) one can see that when ν (cid:29) σ ∼ ν − while σ ∼ ν − .This indicates that the PSD’s width under random switching is significantly larger than in the periodic case, whichallows the total population size to probe smaller values of N in the random than periodic case. This ultimately leadsto a larger fixation probability φ r than φ p (when s > s c , see Sec. S5.2). Importantly, since it is the PSD’s mean thatdetermines the fixation probability at high switching rates, the fact the PSD’s mean here converges at a different rateto K for periodic and random switching gives rise to the different asymptotic behavior of φ r and φ p when ν → ∞ ,yielding Eq. (4) in the main text, see also Eq. (S35). Notably, the convergence details φ α ν →∞ −−−−→ φ ( ∞ ) are expected togenerally depend on the underlying periodic/random processes ξ α ( t ).What is the regime of applicability of these PDMP-like approximations? According to Eqs. (S36), σ (cid:29) K for1 (cid:28) ν (cid:28) √ K , while σ (cid:29) K for 1 (cid:28) ν (cid:28) K ; in these regimes the variance stemming from periodic/randomswitching is much larger than the variance caused by demographic fluctuations. Hence, P PPP ν and P PDMP ν are accurateapproximations of P ( p ) ν and P ( r ) ν in the fast switching regime respectively when ν (cid:28) √ K and ν (cid:28) K . Remarkably,environmental noise also dominates over demographic fluctuations for slow/intermediate switching regime when ν (cid:46) P ν ( N )over a broad range of ν . It is worth noting that for ν (cid:29)
1, the variance of P Kap ν to leading order in 1 /ν satisfies σ = K (cid:2) O (cid:0) K /ν (cid:1)(cid:3) . Thus, as we have checked, the variances of P PPP ν and P Kap ν coincide in the leading order,for 1 (cid:28) ν (cid:28) √ K . However, σ (cid:29) σ when ν (cid:38) √ K , which indicates that the Kapitza-based approximationis superior to that of the PPP in this regime. This also reiterates that, at very high switching rates, one must takedemographic noise into account as is done using the Kapitza method (see previous section). Since the Kapitza-basedapproximation works well for any arbitrary switching rate ν (cid:29)
1, the calculation of Sec. S4.2 leading to the fastswitching asymptotic behavior of φ p ( ν ) has been carried out using P Kap ν/s instead of P PPP ν/s .However, while the PPP and PDMP approximations characterize well the PSD in the random and periodic cases,respectively when ν (cid:28) K and ν (cid:28) K / , they can still be aptly used for the purpose of calculating the fixationprobability, φ α , at arbitrary high switching rates according to (S39) and (S38), see Figs. 3(a,d) and S2 (d). This isbecause only the vicinity of the PSD’s maximum contributes to the leading-order calculation of φ α at high ν , whichis well captured, for any high switching rate, by these PDMP-like approximations. S5. FURTHER DETAILS ABOUT FIGURE 3(D,E) IN THE MAIN TEXT
Here we elaborate on our findings, see main text, that under certain conditions the fixation probability of the S species, φ α , at given s, γ, δ, K , x is optimal for a nontrivial switching rate ν ∗ α , see Fig. 3(d,e). We also discuss8the critical selection intensity below/above which φ α ( ν ) is an increasing/decreasing function under weak switchingasymmetry. S5.1. Region of the parameter space in which the fixation probability is nonmonotonic
Our starting point is Eq. (3) of the main text which, when substituting P ( r ) ν/s by its PDMP approximation, reads φ r ( ν ) (cid:39) (cid:90) K + K − P PDMP ν/s ( N ) φ ( x ) | N dN. (S38)We now give further details on how to determine from this equation the region of the parameter space of Fig. 3(e) in FIG. S2: (a) Triangular-like region in the parameter space in which φ r ( ν ) has a nontrivial maximum at ν = ν ∗ r for s =0 . , . , . , . , . , .
09 (red to blue, left to right) obtained from Eq. (S38). This region, defined by γ > γ c , δ > δ c ( γ, s ) isdelimited by the solid and dashed lines, see text and compare with Fig. 3(e) in the main text. (b) Critical selection intensity s c as a function of δ for γ = 0 . , . , . , . K = 250 and x = 0 .
6, see text. (c) φ r (circles) and φ p (squares) versus ν with ( s, K , γ, x ) = (0 . , , . , . δ = 0 . δ = 0 . P LNA ν to approximate the PSD in Eq. (3), seetext. (d) Colored symbols show φ r versus ν for ( s, K , x ) = (0 . , , .
6) and ( γ, δ ) = (0 . , − .
5) (orange), ( γ, δ ) = (0 . , . γ, δ ) = (0 . , .
6) (blue). Solid colored lines are from (S38) and dashed lines show φ (0 , ∞ ) for ν → , ∞ . Results for φ p and same values of ( γ, δ ) are shown as black squares, with solid black lines from Eq. (S39). (e) Heatmap of ν ∗ p : ν ∗ p → , ∞ in the black and white areas, respectively; φ p ( ν ) is non-monotonic in the red-yellow area, with the values of ν ∗ p given in thevertical bar; parameters are ( s, K , x ) = (0 . , , . ν ∗ p for γ = 0 . δ = 0 . δ = 0 . φ p shown in (c). To be compared with Fig. 3(e) of the main text showing the heatmap of ν ∗ r , see text. which φ α ( ν ) is non-monotonic, and how this region changes when s is increased. Using the diffusion approximation φ ( x ) | N (cid:39) ( e − Ns (1 − x ) − e − Ns ) / (1 − e − Ns ), and Eq. (S22), we compare the PDMP-based approximation of φ r ( ν )[Eq. (S38)] for different switching rates (slow, intermediate and fast switching) for a given set ( K , s, γ, δ, x ), anddetermine for which of these values φ r is maximal.If φ r ( ν (cid:28) s ) > φ r ( ν (cid:29) s ) , φ r ( ν ∼ s ), we say that the optimal fixation probability is φ r (0) = φ (0) at slow switching,i.e., ν ∗ r = 0. Similarly, if φ r ( ν (cid:29) s ) > φ r ( ν (cid:28) s ) , φ r ( ν ∼ s ), the optimal fixation probability is φ r ( ∞ ) = φ ( ∞ ) ,i.e., ν ∗ r = ∞ , see Fig. 3(a). Otherwise, the S fixation probability is maximal at a non-trivial switching rate ν ∗ r ∼ s that Eq. (S38) captures reasonably well. In this case, φ r varies non-monotonically with ν , see Fig. 3(d,e). Wehave performed extensive stochastic simulations of the model’s dynamics and found that this behavior arises in atriangular-like region in the subset of the parameter space where γ and δ exceed some critical values γ c ( s ) and δ c ( γ, s ), see Fig. 3(e). In order to determine the boundary ( γ c ( s ) , δ c ( γ, s )) of this triangular-like region at fixed s and9 K , we have systematically calculated the fixation probability as ν varies from 10 − (proxy for ν (cid:28) s ) to 1 (proxy for ν (cid:29) s ) for fixed ( γ, δ ), with γ (cid:38) .
7, and δ (cid:38) .
2. For each pair ( γ, δ ) we have then found the value of ν for whichit attains its maximum and store it in a matrix. In practice, the diagonal part of the boundary is then found bykeeping γ fixed and increasing δ until we find the first entry of the matrix for which 10 − < ν ∗ r <
1. This determines δ c ( γ, s ). The left hand side of the boundary, γ c ( s ) is found by finding the largest value of γ such that ν ∗ r = 1 or ν ∗ r = 10 − for all δ . Predictions of Eq. (S38) are in good agreement with simulation results which confirm that φ α ( ν )has a nontrivial maximum at ν ∗ α when γ > γ c ( s ) and δ > δ c ( γ, s ), see Figs. 3(d) and S2(c). As shown in Fig. S2 (a), γ c ( s ) is an increasing function of s while δ c ( γ, s ) changes little when γ and s are increased. As a result, when theselection intensity s is increased (at K fixed), the triangular-like region of the parameter space in which φ α ( ν ) hasa nontrivial maximum is “squeezed out”, as shown in Fig. S2(a), whereas we have verified that the optimal fixationprobability remains unaltered when s changes but K s is kept fixed.A similar analysis can be carried out for the periodic switching. In this case, the fixation probability is approximatedby substituting the PPP approximation (S19) with rescaled switching rate into the PSD in Eq. (3), yielding φ p ( ν ) (cid:39) (cid:90) N max N min P PPP ν/s ( N ) φ ( x ) | N dN. (S39)This expression gives a sound approximation of φ p , and correctly predicts the same qualitative behavior as underrandom switching, as shown in Figs. 3(d) and S2(d,e). We have used Eq. (S39) to obtain the heatmap of ν ∗ p ofFig. S2 (e) giving the switching rate ν ∗ p for which φ p ( ν ) is maximal. The comparison with the heatmap of ν ∗ r ofFig. 3(e) in the main text for the same parameters s, K , x , shows that both φ p ( ν ) and φ r ( ν ) are non-monotonicin qualitatively similar triangular-like regions of the γ − δ parameter space (triangular-like region of Fig. S2 (e) andFig. 3(e) are of similar size). We notice that ν ∗ p (cid:46) ν ∗ r for given γ and δ , which translates in the triangular-like regionof Fig. S2 (e) to be overall more “reddish” than the corresponding region of Fig. 3(e). S5.2. Critical selection intensity
As explained in the main text, it is useful to determine the critical selection intensity s c such that φ (0) = φ ( ∞ ) .Here, using the diffusion approximation we have φ (0) = [(1 − δ ) φ ( x ) | K − + (1 + δ ) φ ( x ) | K + ] / φ ( ∞ ) = φ ( x ) | K .By introducing z = exp( − sK − ), a = (1 + γ ) / (1 − γ ) and a = (1 + γ ) / (1 − δγ ), which yields K + = a K − and K = a K − , s c is obtained by solving φ (0) = φ ( ∞ ) , i.e., it is the solution of the transcendental equation(1 − δ ) (cid:18) z − x − − z (cid:19) + (1 + δ ) z a − (cid:18) z − a x − − z a (cid:19) − z a − (cid:18) z − a x − − z a (cid:19) = 0 . The numerical solutions of this equation, for s c = s c ( γ, δ ), are reported in Fig. S2(b), where we find that s c decreaseswith δ , and increases with γ . When s < s c , φ (0) < φ ( ∞ ) and φ (0) > φ ( ∞ ) if s > s c . This allows us to determine themonotonic behavior of φ α ( ν ) under weak switching asymmetry ( | δ | < δ c ): When s < s c , φ α is an increasing functionof ν , while it decreases with ν if s > s c (at given s, γ, δ ). In the examples of Fig. S2(d) we find that s c ≈ .
06 for( γ, δ ) = (0 . , . s c ≈ .
03 for ( γ, δ ) = (0 . , . s c ≈ .
095 for ( γ, δ ) = (0 . , − . s = 0 .
05, thiscorresponds to φ α ( ν ) being an increasing function of ν for ( γ, δ ) = (0 . , .
5) and ( γ, δ ) = (0 . , − . ν in the case ( γ, δ ) = (0 . , . < s (cid:28) K (cid:29)
1, as considered in this work, the generic case is s > s c and φ α ( ν )therefore generally decreases with ν as in Fig. 3(a).In the regime γ > γ c ( s ) and δ > δ c ( γ, s ) where φ α ( ν ) is non-monotonic, we can determine that φ α ( ν ) increasessteeper at slow/intermediate switching if s < s c while it is the opposite when s > s c . As confirmed by the resultsreported in Fig. S2(d), where φ r and φ p exhibit the same ν -dependence, these results hold for both random and periodicswitching, since in both cases φ α ( ν ) essentially coincides with φ (0) and φ ( ∞ ) for slow/fast switching, respectively. S5.3. Effective selection intensity under fast switching
As discussed in the main text, see also Sec. S4, under fast random and periodic switching φ α ( ν ) (cid:39) φ ( ∞ ) = φ ( ∞ ) | K , with φ ( ∞ ) = e m/ and m ≡ K (1 − x ) ln (1 − s ) = 2 (cid:2) (1 − γ ) / (1 − γδ ) (cid:3) K (1 − x ) ln (1 − s ), to leadingorder in 1 /ν . When 1 /K (cid:28) s (cid:28) / √ K and γ = O (1), the above expression simplifies: φ ( ∞ ) (cid:39) e −K s (1 − x ) = e − s [(1 − γ ) / (1 − γδ )] K (1 − x ) . Hence, in this regime, under fast random and periodic switching, the S fixation probabilityis the same as in a population subject to a constant carrying capacity K under a rescaled selection intensity s → FIG. S3: (a) Unconditional MFT under random (circles) and periodic (squares) switching versus ν for ( K , s, x ) =(250 , . , . γ, δ ) = (0 . , .
5) (purple), (0 . , − .
2) (orange), (0 . , .
7) (green) and (0 . , .
8) (teal). The MFT scales as1 /s ; the effect of random/periodic switching is to reduce the subleading corrections due to the decreasing average populationsize (cid:104) N (cid:105) , see text. (b) (cid:104) N (cid:105) versus ν for random (circles) and periodic (squares) switching, for the same parameters as in (a).Solid colored lines are given by (cid:104) N (cid:105) PDMP and solid black lines show (cid:104) N (cid:105) PPP , see Eq. (S37). (c) Average number of switchesdivided by ν prior to fixation versus s , for δ = 0 . , − . ν = 0 . , ,
10 (squares, circles, triangles) areshown to be O (1 /s ) with data for different values of ν collapsing together. Other parameters are: ( K , γ, x ) = (250 , . , . s (cid:48) = s (1 − γ ) / (1 − γδ ). This result yields the following remarkably simple and enlightening interpretation: in theabove regime, the effect of environmental variability, when δ < γ , is to effectively reduce the selection intensity withrespect to the static environment, yielding φ α > φ ( x ) | K under a selection intensity s . Similarly, there is an effectiveincrease of selection intensity ( s (cid:48) > s ) when δ > γ resulting in φ α < φ ( x ) | K . S5.4. Duty cycle and general effect of δ on φ α The parameter δ measures the asymmetry in the switching rate, and can be used to define the “duty cycle” as(1 + δ ) / K + and K − with period T = (1 /ν − ) + (1 /ν + ). The duty cyclegives the fraction (1 /ν + ) /T of one period spent in the environmental state ξ = 1. Clearly, when δ >
0, the populationspends more time in the environmental state ξ = +1 (with K = K + ) than in the state ξ = − K = K − ). Since s >
0, species S has a selective disadvantage with respect to strain F and φ α is therefore a decreasing function of δ when all the other parameters are fixed, see Fig. S2(c,d). S6. MEAN FIXATION TIME AND AVERAGE NUMBER OF SWITCHES
In addition to the fixation probability, we have also computed the MFT, T ( α ) ( x ) – the unconditional mean timeuntil the fixation of either species S or F , starting from a initial fraction x of individuals of type S . As in thecase δ = 0, T ( α ) ( x ) is obtained by averaging the unconditional MFT, T ( x ) | N , obtained in a population of constantsize N over P ( α ) ν/s ( N ) with a rescaled switching rate ν → ν/s [44, 45]. In the limits of slow and fast switching, wehave T ( α ) ( x ) = [(1 + δ ) T ( x ) | K + + (1 − δ ) T ( x ) | K − ] / ν/s (cid:28) T ( α ) ( x ) = T ( x ) | K when ν/s (cid:29)
1, seeFig. S3(a). When 1 /K (cid:28) s (cid:28) T ( x ) | N ∼ O (1 /s ) [48, 66], and the MFT under random switching also scalesas 1 /s , i.e., T ( α ) ( x ) = O (1 /s ); this result is evident since x deterministically relaxes on a time scale O (1 /s ). Asthe average population size (cid:104) N (cid:105) decreases with ν , see Fig. S3(b), environmental variability reduces the subleadingprefactor of T ( α ) ( x ) [45]. As a consequence, on average the population experiences O ( ν/s ) switches prior to fixationwhen 1 /K (cid:28) s (cid:28)
1. Fig. S3(c) confirms that in this regime the average number of switches prior to fixation scalesas 1 /s and increases linearly with ν to leading order. Since the PSD greatly varies when ν and δ change, see Figs. 2and S1, the fact that the average number of switches increases linearly with ν shows that it is essentially independentof the population size and supports the rescaling ν → ν/s in the approximations of Eqs. (3) and (S38).1 FIG. S4: (a) Effective parameter q versus s for δ = − . , . s = 0 . , .
05 (squares, circles). Dependenceof q on b is approximately linear while q depends weakly on δ and s (solid lines are eyeguides). (b,c,d,e) φ α versus ν for( K , γ, s, δ, x ) = (250 , . , . , . , .
6) in (b,d) and (250 , . , . , . , .
6) in (c,e). Here (b,c) and (d,e) show results forrandom ( α = r ) and periodic ( α = p ) switching, respectively, with the same parameters. In (b,c,d,e) b = (0 , . , . , . , φ PG r ( ν ) from Eq. (S42) in (b,c) and φ PG p ( ν ) from Eq. (S43) in (d,e). In (b,d), φ PG α ( ν ) is an increasing function of ν at low valuesof b , and varies nonmonotonically with ν for intermediate b ’s. In (c,e), φ PG α ( ν ) is a nonmonotonic function of ν at low b ’s andbecomes a decreasing function of ν as b increases. S7. ECO-EVOLUTIONARY DYNAMICS & FIXATION PROBABILITY IN A PUBLIC GOODSCENARIO
The model studied in the main text describes the competition for resources of the slow and fast growing strains S and F without assuming any explicit interactions between them. Yet, as discussed in Sec. S1.1 of this SM, themodel can be generalized to describe the situation where strain S is a public good (PG) producer. Here, we considerthe situation where S produces a PG that benefits the entire population that is subject to a time-varying carryingcapacity.A simple way to describe a PG scenario in the general framework outlined in Section S1.1 is to multiply the percapita birth rate by the global term g ( x ) = 1 + bx , where b ≥ N S/F T − S/F −−−→ N S/F −
1, and N S/F T + S/F −−−→ N S/F + 1,with the modified transition rates T + S = g ( x ) − s ¯ f N S , T + F = g ( x )¯ f N F , and T − S = NK ( t ) N S , T − F = NK ( t ) N F , where b = O (1) while K ( t ) is given by Eq. (1) of the main text. To discuss how the properties of this model can be studiedby extending the analysis carried out in the main text, it is convenient to first consider specifically the case of randomswitching ( α = r ). When demographic noise is neglected and the only source of randomness stems from the randomlyswitching K ( t ), the population’s mean-field dynamics obeys [44, 45] (see also Sec. S1.2) dxdt = − sg ( x ) x (1 − x )1 − sx and dNdt = N (cid:20) g ( x ) − N K (cid:18) − γξ r ( t )1 − γδ (cid:19)(cid:21) , (S40)with K ≡ K (1 − γ ) / (1 − γδ ). N and x are thus explicitly coupled, which breaks the time separation and yields anexplicit form of eco-evolutionary dynamics. Analytical progress can be made by using the effective theory devised inRefs. [44, 45]. Since the model’s dynamics under a constant carrying capacity is well described in terms of a populationof an effective size, as in the case δ = 0 [44, 45], we introduce a suitable parameter q (with 0 ≤ q ≤ b ) and replace g ( x ) by 1 + q in (S40). This decouples N and x , and one can thus perform a similar PDMP-based approximation as2before, yielding P PDMP ν,q ( N ) ∝ N (cid:34)(cid:18) (1 + q ) K + N − (cid:19) ν +(1+ q ) − (cid:18) − (1 + q ) K − N (cid:19) ν − (1+ q ) − (cid:35) , (S41)where we have omitted the normalization constant. As in Refs. [44, 45], the parameter q is obtained by matching thesimulation results for the S fixation probability in the fast switching limit (i.e., when ν/s (cid:29)
1) with φ ( x ) | (1+ q ) K .Results reported in Fig. S4(a), obtained using the diffusion approximation [see below and Eq. (S38)], show that q = q ( b )increases almost linearly with b , and depends only weakly on s and δ , with q (0) = 0 when b = 0. From P PDMP ν,q itis clear that the effect of increasing b , and therefore the effective parameter q ( b ), results in effectively increasing thecarrying capacity K ± → (1+ q ( b )) K ± and reducing the switching rates ν ± → ν ± / (1+ q ) = ν (1 ∓ δ ) / (1+ q ). Proceedingas in the case b = q = 0 and δ = 0 [44, 45], the fixation probability is obtained by averaging φ ( s, x ) | N over the PSDin Eq. (S41) with ν → ν/s [88]. Furthermore, by changing the variable of integration to N (cid:48) = N/ (1 + q ) we find thatthis is equivalent to rescaling the selection strength to s eff = (1 + q ) s in the model without a PG φ PG r ( ν, q ) = (cid:90) (1+ q ) K + (1+ q ) K − φ ( s, x ) | N P PDMP ν/s,q ( N ) dN = (cid:90) K + K − φ ( s eff , x ) | N (cid:48) P PDMP ν/s eff , ( N (cid:48) ) dN (cid:48) , (S42)where N (cid:48) = N/ (1 + q ) and we have used φ ( s, x ) | N (cid:39) ( e − Ns (1 − x ) − e − Ns ) / (1 − e − Ns ).According to Eq. (S42), the effect of increasing b results in raising the value of the corresponding value of q , seeFig. S4(a), which in turn results in a carrying capacity switching between (1 + q ( b )) K ± . Thus, as in the case δ = 0,one can transform the expression of the S fixation probability, φ PG r ( ν, q ), to the (approximate) fixation probabilityin the absence of PG but under an effective (increased) selection intensity s eff = [1 + q ( b )] s . This results in φ PG r decaying approximately exponentially with b [45].Equation (S42) is an approximation of the actual fixation probability φ r that is valid over a broad range of frequencies ν and gives an accurate description of φ r when δ = 0 [44, 45] and | δ | (cid:28) | δ | and b increase. Here, we are chiefly interested in the qualitative dependence of the fixation probabilityon ν when b changes and γ = O (1), δ = O (1) (see Figs. 3(e) and S2(e)). With Fig. S2(a) in mind, we can understandhow raising b changes the diagram of Fig. S2(a): As b is increased, the triangular-like region is squashed since γ c increases under the effect of s → s eff = (1 + q ( b )) s . This allows us to qualitatively explain how the fixation probability φ PG r varies with ν under intermediate switching at γ, δ, s fixed. In the case of Fig. S4(b), δ < δ c at low b and therefore φ PG r ( ν ) increases monotonically; then as γ c increases together with b , ( γ, δ ) enter the triangular-shaped region (i.e., γ > γ c , δ > δ c ) of Fig. S2(a) where φ PG r ( ν ) varies non-monotonically with ν . In the example of Fig. S4 (c), γ > γ c and δ > δ c at low b implying that φ PG r ( ν ) is a nonmonotonic function of ν ; then γ c increases along with b and attainsa value such that γ < γ c with δ > δ c , and in this case φ PG r ( ν ) decreases monotonically with ν . Hence, while Eq. (S42)cannot accurately predict the full ν dependence of φ r , it qualitatively captures the emergence of a peak in S4(b) atsome nontrivial intermediate switching rate, and the disappearance of the peak in S4(c), when b is increased. Theseare examples of the rich and complex behavior that eco-evolutionary loops can generate.The results of this section have so far focused on the case of random switching, but we have again obtained asimilar qualitative behavior with periodic switching, as shown in Fig. S4(d,e). This can be explained in terms ofa PPP-based approximation in the realm of an effective theory as in the random switching case. In fact, we haveverified that the effective parameter q ( b ) allows us to obtain a suitable approximation of the fixation probability underfast periodic switching, i.e. φ p (cid:39) φ | (1+ q ) K . This suggests to use the PPP-based approximation P PPP ν,q as an effectiveapproximate PSD in the PG scenario with periodic switching, where P PPP ν/s,q is obtained from Eqs. (S19)-(S21) byrescaling ν ± → ν ± / [(1 + q ) s ] and K ± → (1 + q ) K ± . This rescaling of ν ± and K ± results in a support of P PPP ν/s,q thatis now denoted by [ N min ( q ) , N max ( q )]. In the same vein as in the random switching case, we thus write φ PG p ( ν, q ) = (cid:90) N max ( q ) N min ( q ) φ ( s, x ) | N P PPP ν/s,q ( N ) dN, (S43)which is expected to be a suitable approximation of φ p when | δ | (cid:28)
1, and to qualitatively capture the ν dependence of φ p when δ = O (1). The results of Fig. S4(d,e) indeed show that φ PG p provides the same qualitative description of thefixation probability as Eq. (S42) in the random switching case. In particular, Eq. (S43) qualitatively reproduces theemergence of a peak at a nontrivial frequency in S4(d), and the disappearance of the peak in S4(e), as bb