Exclusion of the fittest predicts microbial community diversity in fluctuating environments
MMicrobial species interactions determine community diversity influctuating environments
Shota Shibasaki , Mauro Mobilia , and Sara Mitri Department of Fundamental Microbiology, University of Lausanne, Switzerland Department of Applied Mathematics, School of Mathematics, University of Leeds, United Kingdom
E-mail addresses of all authors.SS: [email protected]: [email protected]: [email protected] running title: Microbial diversity under fluctuationKey words: competitive exclusion, environmental switching, demographic noise, chemostat, beta diversity,intermediate disturbance hypothesisType of article: articleNumbers of words: 150 words in abstract, and 5052 words in main text.Number of references: 53.Numbers of figures: 4, tables: 1, text box 0.Corresponding author: Sara Mitri: Quartier UNIL-Sorge, Batiment Biophore, Office 2414, CH-1015 Lau-sanne. Phone : +41 21 692 56 12. E-mail: [email protected] of authorship: SS, MM and SM designed the study, SS performed simulations, analyzed data,and wrote the first draft, and all authors contributed to revisions.1 a r X i v : . [ q - b i o . P E ] S e p bstract Microorganisms often live in environments that fluctuate between mild and harsh conditions. Althoughsuch fluctuations are bound to cause local extinctions and affect species diversity, it is unknown how diversitychanges at different fluctuation rates and how this relates to changes in species interactions. Here, we usea mathematical model describing the dynamics of resources, toxins, and microbial species in a chemostatwhere resource supplies switch. Over most of the explored parameter space, species competed, but thestrength of competition peaked at either low, high or intermediate switching rates depending on the species’sensitivity to toxins. Importantly, however, the strength of competition in species pairs was a good predictorfor how community diversity changed over the switching rate. In sum, predicting the effect of environmentalswitching on competition and community diversity is difficult, as species’ properties matter. This mayexplain contradicting results of earlier studies on the intermediate disturbance hypothesis. Introduction
Natural environments are not static: temperature, pH, or availability of resources change over time. Manystudies in microbiology, ecology and evolution have focused on responses to fluctuations in resource abundancein the regime of feast and famine periods (Hengge-Aronis, 1993; Vasi et al., 1994; Srinivasan and Kjelleberg,1998; Xavier et al., 2005; Merritt and Kuehn, 2018; Himeoka and Mitarai, 2019). These models capture thedynamics within many natural ecosystems. For example, the gut microbiome of a host is exposed to fluctuatingresources that depend on its host’s feeding rhythm, which may affect microbiota diversity (Cignarella et al.,2018; Li et al., 2017; Thaiss et al., 2014). In addition to their magnitude, environmental fluctuations can alsodiffer in their time scales: for the gut microbiota, a host’s feeding rhythm may vary from hourly to daily, oreven monthly if feeding depends on seasonal changes (Davenport et al., 2014; Smits et al., 2017). Finally, thehost does not always feed its gut flora, but may also expose it to lethal toxins through the consumption ofantimicrobial treatments.How fluctuations affect species diversity has been a highly contested topic in ecology. The intermediatedisturbance hypothesis argues that intermediate intensity and frequency of disturbance maximize species di-versity (Connell, 1978; Grime, 1973) because species with high competitive ability dominate at a low level ofdisturbance while species that adapted to the disturbance dominate at high disturbance (Grime, 1977). Thishypothesis is controversial (Fox, 2013) and other relationships between disturbance and species diversity havebeen reported both empirically and theoretically (Mackey and Currie, 2001; Miller et al., 2011).Another body of theory that considers how environmental fluctuations can affect species diversity is moderncoexistence theory (Chesson, 1994). In modern coexistence theory, spatial and temporal differences in theenvironment mediate species coexistence when species differ in their preferences for environmental conditionsand/or in their response to fluctuations in factors like resources and toxins (Amarasekare, 2019; Chesson,2000a,b; Letten et al., 2018b; Barab´as et al., 2018; Ellner et al., 2019). Resources and toxins have oppositeeffects on growth rates: consumption of resources by other species decreases the growth rate of a focal species(competition/exploitation) while absorption of toxin by others can increase its growth rate (facilitation). Thisimplies that environmental conditions, especially amounts of resources and toxins, can affect the sign and/ormagnitude of species interactions (Hoek et al., 2016; Piccardi et al., 2019; Zu˜niga et al., 2019). In turn,how pairs of species in a community interact (i.e., competition or facilitation) affects species diversity andcommunity stability (Mougi and Kondoh, 2012; Coyte et al., 2015). Environmental fluctuations in nutrient andtoxin concentrations are therefore expected to affect species interactions, which should then affect communitydiversity.Recently, Rodr´ıguez-Verdugo et al. (2019) used theory and experiments to show that the rate of environ-mental fluctuations affects the coexistence of two species when the environment switches between two carbonsources that promote either exploitation or competition between the species. However, it is unclear whetherand how environmental fluctuations in amounts of – rather than types of – resources affect species coexistenceand diversity because we do not know how species interactions change a priori . More generally, we ask: can3ne predict patterns of diversity under fluctuating environments in larger communities composed of more thantwo species? How do these changes in diversity relate to changes in inter-species interactions? And how doenvironmental conditions (e.g., resource or toxin concentrations) affect these changes?To shed light on these questions, we develop a mathematical model to investigate how the rate of environmen-tal fluctuations affects species interactions and diversity. We start with a simple scenario with just two speciesin an environment that switches between scarce and abundant resource supplies. In addition to resources, ourmodel includes the dynamics of toxins, as we expect them to affect inter-species interactions (Hoek et al., 2016;Piccardi et al., 2019; Zu˜niga et al., 2019). Accordingly, the two species can both compete for resources andfacilitate each other by absorbing toxic compounds. We focus on the faster grower and ask: (i) How often doesthis species go extinct due to the presence of a slower-growing species, and (ii) how does such extinction dependon the environmental switching rate and other properties of the environment? For most parameter values, thefaster grower’s extinction probability increases because of competitive exclusion by the slower grower. To oursurprise, we found that this competitive effect peaked at different switching rates depending on environmentaltoxicity. Second, to verify whether these findings generalize to, and predict the dynamics of larger communitydynamics, we simulated communities of up to 10 species and found a similar pattern: beta diversity was highestat different switching rates in different environments. Although various relationships between disturbance anddiversity have been previously reported (Mackey and Currie, 2001; Miller et al., 2011), our model associates thestrength of competition between two species with the loss of beta diversity in a larger community in a fluctuatingenvironment. In other words, to estimate how beta diversity changes over the environmental switching rate,it may be sufficient to analyze the sign and strength of interactions between two representative species underthose same conditions.
In this paper, we analyze how environmental switching rate ν affects species interactions and diversity in achemostat model. For simplicity, we assume symmetric random switching (but see Taitelbaum et al. (2020)for asymmetric, random and periodic switching). This model combines two sources of noise: birth and deathevents occur randomly yielding demographic noise, while the environment randomly switches from being harshto mild and vice versa and drives the population size (Wienand et al., 2017, 2018; West and Mobilia, 2020;Taitelbaum et al., 2020). When the environment becomes harsh (e.g., abundant toxins or scarce resources), thepopulation shrinks and demographic noise can lead to extinction. As in previous work, a distinctive feature ofthis model is therefore the coupling of demographic and environmental noise: environmental switching changesthe population size, which in turn modulates demographic fluctuations, which results in feedback loops. Here,we implement environmental switching by modeling time-varying resource and/or toxin supplies, as we areinterested in how these environmental conditions mediate species interactions and diversity.Our model includes the dynamics of the amount r i of resource i ( i = 1 , . . . , N/
2, where N is an even number)that allows all species to grow, the amount t j of toxic compound j ( j = 1 , . . . , N/
2) that kills all species, and the4bundance s k of species k ( k = 1 , . . . , N ). When N = 2 species, we recapitulate the environmentally mediatedspecies interactions shown in the model of Piccardi et al. (2019), where species affect each other positively intoxic environments, with increasing competition in more benign environments. The dynamics are modelled interms of a multivariate continuous-time birth-and-death process (see Novozhilov et al. (2006); Allen (2010) forexample) defined by the following “birth” and “death” reactions: r i τ + ri −−→ r i + 1 , (1a) r i τ − ri −−→ r i − , (1b) t j τ + tj −−→ t j + 1 , (1c) t j τ − tj −−→ t j − , (1d) s k τ + sk −−→ s k + 1 , (1e) s k τ − sk −−→ s k − , (1f)occurring with transition rates: τ + r i = αR i ( ξ ) , (2a) τ − r i = N (cid:88) k =1 µ ik Y rik r i r i + K rik s k + αr i , (2b) τ + t j = αT j ( ξ ) , (2c) τ − t j = N (cid:88) k =1 δ jk Y tjk t j t j + K tjk s k + αt j , (2d) τ + s k = N/ (cid:88) i =1 µ ik r i r i + K rik s k , (2e) τ − s k = N/ (cid:88) j =1 δ jk t j t j + K tjk s k + αs j , (2f)where α is the dilution rate of the chemostat, ξ = ± R i ( ξ ) ( T j ( ξ )) is resource i ’s (toxin j ’s) supply underthe environmental condition ξ , Y rik ( Y tjk ) is species k ’s biomass yield for resource i (toxin j ), µ ik is the maximumgrowth rate of species k by resource i , δ jk is the maximum death rate of species k by toxin j (species k ’ssensitivity to toxin j ), and K rik ( K tjk ) is the amount of resource i (toxin j ) that gives the half-maximum growth(death) rate of species k . In this model, we assume that growth and death rates depend on the amounts ofresource and toxin in Monod function forms, respectively, and that the effects are additive.Environmental switching changes resource and/or toxin supplies at rate ν according to ξ ν −→ − ξ, (3)5here ξ ∈ {− , +1 } is a coloured dichotomous (telegraph) noise assumed to be symmetric and stationary (Bena,2006; Horsthemke and Lefever, 2006). Here, the mean of ξ ( t ) is therefore always zero and its finite correlationtime is 1 /ν . One can imagine three environmental switching scenarios, where either or both resource and toxinsupplies change (Table 1). In the main text, we analyze only scenario 1 where only resource supplies switch,but see Appendix 4 for the remaining scenarios. When environmental switching affects the resource supply R i ,it switches between abundant and scarce: R i ( ξ ) = R + i ξ = 1 (abundant) R − i ξ = − R + i > R − i > , ∀ i. (4)It is clear that (2a)-(4) define a multivariate birth-and-death process coupled with environmental switching. Asenvironmental switching does not affect toxin supply in this scenario, the mean amount of toxins is supplied (cid:104) T j (cid:105) , where (cid:104)·(cid:105) represents the mean of abundant and scarce supplies of a focal resource or toxin (e.g., (cid:104) T j (cid:105) ≡ (cid:0) T + j + T − j (cid:1) / T + j > T − j > , ∀ j ).Note that when demographic noise is ignored, the sole source of randomness stems from the dichotomousMarkov noise, which clearly results in the waiting time between two environmental switches to be exponen-tially distributed with a mean 1 /ν (Bena, 2006). In the joint presence of demographic noise and environmentalswitching, see, e.g., (Hufton et al., 2016; West et al., 2018), the waiting time approximately follows the expo-nential distribution with mean 1 /ν , with the larger the switching rate the shorter the sojourn in either of theenvironmental states.The master equation for this model is defined by combining the dynamics of amount of resources and toxins,abundances of species, and environmental switching (Eqs (1a) - (3)):˙ P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1) = N/ (cid:88) i =1 (cid:0) E − r i − (cid:1) (cid:8) τ + r i P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1)(cid:9) + N/ (cid:88) i =1 (cid:0) E + r i − (cid:1) (cid:8) τ − r i P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1)(cid:9) + N/ (cid:88) j =1 (cid:16) E − t j − (cid:17) (cid:110) τ + t j P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1)(cid:111) + N/ (cid:88) j =1 (cid:16) E + t j − (cid:17) (cid:110) τ − t j P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1)(cid:111) + N (cid:88) k =1 (cid:0) E − s k − (cid:1) (cid:8) τ + s k P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1)(cid:9) + N (cid:88) k =1 (cid:0) E + s k − (cid:1) (cid:8) τ − s k P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1)(cid:9) + ν (cid:8) P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, − ξ (cid:1) − P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1)(cid:9) (5)where P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1) represents the probability density of (cid:126)r = ( r i ), (cid:126)t = ( t j ), (cid:126)s = ( s k ), and ξ , the dot denotes the6ime derivative, and E ± r i is a shift operator such that E ± r i P (cid:0) (cid:126)r, (cid:126)t, (cid:126)s, ξ (cid:1) = P (cid:0) r , . . . , r i ± , . . . , r N/ , (cid:126)t, (cid:126)s, ξ (cid:1) , (6)and E ± t j and E ± s k are the equivalent shift operators for t j and s k , respectively. The last line on the right-hand-side of (5) accounts for the random switching reaction (3). The model given by Eq. (5) was implementedusing the Gillespie algorithm (Gillespie, 1977) until, after a sufficiently large time T end , the population sizedistribution converges to a quasi-stationary distribution (i.e., where distributions of species abundances appearto be stationary for a long time; the real equilibrium state corresponding to the extinction of all species beingpractically unobservable). In two species scenarios, either or both two species go extinct before time T end inmost cases (see Fig. 4 column C.). For details on parameter values and the simulations, see Appendix 2. Aschematic illustration of this model with N = 2 is summarized in Fig. 1. We begin by analyzing interactions between two species ( N = 2) where the sign and magnitude of speciesinteractions can change. We used parameter values such that species 1 always outcompetes species 2 in adeterministic and fixed environment setting (i.e., species 1 grows faster than species 2, see Appendix 1 for analysisand Table A.1 for exact parameter values) to clarify the effects of demographic noise. Under demographic noisecoupled with environmental switching, either of the two species or both species tend to go extinct. As a proxyfor interactions, we focus on the net effect of species 2 on species 1, which is defined by the extinction probabilityof species 1 in mono-culture minus that in co-culture with species 2:∆ P ( s ( T end ) = 0) ≡ P ( s ( T end ) = 0; s (0) = 0) − P ( s ( T end ) = 0; s (0) > . (7)If ∆ P ( s ( T end ) = 0) is < − .
01, we consider that species 2 increases the extinction probability of species 1(negative effect), whereas ∆ P ( s ( T end ) = 0) > .
01 implies that species 2 has a positive effect on species 1. Weignore small values of | ∆ P ( s ( T end ) = 0) | , as they are biologically meaningless and may be due to the finitenumber of simulations.In the analysis of species diversity (see below) with N = 2, we also analyzed the species interactions in 100pairs of species that differ in µ ik , K rik , δ jk , and K tjk in the presence of demographic noise and environmentalswitching. Since pairs of species may coexist depending on their parameter values in the diversity analysis,species 1 (or 2) is the one that has larger (or smaller) population size than the other in the absence of any noise.If both species go extinct in the absence of noise, species 1 is chosen randomly. To explore how species diversity changes over switching rates, we ran simulations at different community sizesranging from N = 2 to N = 10 and different mean toxin sensitivities, from ¯ δ = 0 . δ = 1. For each7ondition (one N and one ¯ δ ), we sampled 100 sets of parameters µ ik , K rik , δ jk , and K tjk , which represented 100meta-communities composed of N species whose differences were reflected in the parameter values. Each meta-community had 100 independent communities (no inter-community migration and independent environmentalswitching) with equal initial species abundances. In short, we replicated 100 simulations for each of 100 sets ofparameter values, and we had in total 10,000 simulations for each combination of community size, mean toxinsensitivity, and switching rate.Because demographic noise and environmental switching affect species composition, we investigated theheterogeneity of communities at the end of each simulation run. We measured beta diversity (Jost, 2007; Chaoet al., 2012) of each meta-community and species richness (number of surviving species) for each community ata quasi-stationary state ( T end ): D β ( T end ) ≡ D γ ( T end ) D α ( T end ) , (8)with alpha and gamma diversities defined as below: D α ( T end ) ≡ exp (cid:32) − (cid:88) l =1 N (cid:88) k =1 w l p lk ( T end ) ln p ln ( T end ) (cid:33) , (9) D γ ( T end ) ≡ exp (cid:32) − N (cid:88) k =1 ¯ p k ln ¯ p k ( T end ) (cid:33) . (10) w l is a weight for community l calculated by size of community l (sum of species abundances in community l relative to the sum of community sizes over l ), p lk is the relative abundance of species k in community l (i.e.,in community l , p lk ( T end ) = s k ( T end ) / (cid:80) k s k ( T end )), and ¯ p k = (cid:80) l w l p lk is the mean relative abundance ofspecies k among communities l = 1 , . . . , l , it does not affect alpha,beta and gamma diversities as w l = 0. If all species go extinct in all communities, beta diversity becomes D β ( T end ) = 1. For more detail, see Appendix 2. Statistical analysis was performed with Python 3.7.6 and Scipy 1.4.1. For statistical tests of Spearman’s rank-order correlation, scipy.stats.spearmanr was used.
Our aim is to investigate how species interactions and diversity change over the environmental switching rate.Rather than measuring interactions through the effect of each species on the other’s biomass, we focus on afast-growing species (that we call species 1) and analyze how its extinction probability is affected by the presenceof a second slower-growing species (species 2). Our reasoning here is that in the absence of noise, species 18hould always out-compete species 2, so measuring any deviation from this outcome allows us to quantify theeffect of coupled demographic noise and environmental fluctuations. We explored how this proxy for speciesinteractions was affected by different environmental switching rates as well as different toxin sensitivities, whichwe varied simultaneously for both species, such that the species were always equally sensitive to the toxin.When both species were highly sensitive to the toxin (large δ ), species 2 has a positive effect on species 1. Thisoccurs because species 2 also degrades the toxin, which outweighs the competition for nutrients (Piccardi et al.,2019). However, for most parameter values in Fig. 2A, species 2 has a negative effect on species 1 by increasingits extinction probability. We therefore focus on competitive interactions (i.e., ∆ P ( s ( T end ) = 0) <
0) for theremainder of the main part of the study and consider positive interactions in Appendix 5.As we varied the switching rate, we observed that the strength of the competitive effect was highly de-pendent on the toxin sensitivity δ of the two species: monotonically increasing, monotonically decreasing, ornon-monotonically changing with a minimum or maximum value at an intermediate switching rate (Fig. 2B).How our system behaves under the two extreme switching rates becomes clear if one considers what happens inthe absence of environmental switching (Fig. A.3). At very low rate ν →
0, there are almost no switches: theenvironment remains either poor or rich in nutrient (probability 1/2) and microbes are as likely to experienceeither scarce or abundant resources from t = 0 until t = T end . The outcome at this switching rate thereforecorresponds to the mean outcome of those two environments. On the other hand, a very fast environmentalswitching rate ( ν → ∞ ) corresponds to a scenario where the resource supply is at mean concentration (envi-ronmental noise self averages, see, e.g.,Wienand et al. (2017, 2018); West and Mobilia (2020); Taitelbaum et al.(2020)). How the strength of competition varies at an intermediate switching rate is, however, less intuitive.We explore this next. To better understand the unexpected changes in the effect of species 2 on species 1 (Fig. 2B), we analyseour proxy for species interactions in more detail. ∆ P ( s ( T end ) = 0) represents the difference between theprobability of species 1 going extinct in mono-culture ( P ( s ( T end ) = 0; s (0) = 0)) versus in co-culture withspecies 2 ( P ( s ( T end ) = 0; s (0) > P ( s ( T end ) = 0) = − P ( s ( T end ) = 0 , s ( T end ) > s (0) > (cid:124) (cid:123)(cid:122) (cid:125) competitive exclusion+ P ( s ( T end ) = 0; s (0) = 0) (cid:124) (cid:123)(cid:122) (cid:125) sp 1 goes extinct in mono-culture − P ( s ( T end ) = 0 , s ( T end ) = 0; s (0) > (cid:124) (cid:123)(cid:122) (cid:125) both species go extinct . (11)The first line of Eq (11) represents the probability that species 2 excludes species 1 and survives (competitiveexclusion), while the second line is species 1’s extinction probability in mono-culture minus the probability of9oth species going extinct.Fig. A.2 shows that the overall interaction strength is not affected much by the second line of Eq (11), andthat the competitive exclusion probability explains the variation in species interaction strength in most cases(compare Fig. 2A and C). Intuitively, this is because under the environmental conditions where both species arelikely to go extinct (i.e., small resource supplies and high toxin sensitivity), species 1 is also likely to go extinctin mono-culture, resulting in a small magnitude of the second line of Eq (11). For this reason, the competitiveexclusion probability changes similarly to the effect of species 2 on species 1 over environmental switching ratesand toxin sensitivities (Fig. 2C). In other words, we can assume that the second line of Eq (11) is negligibleand focus instead on the competitive exclusion probability. Following the logic outlined above, we now use the competitive exclusion probability to investigate why thestrength of negative interactions changes over the switching rate in various ways, depending on toxin sensitivity.In the absence of environmental switching, the competitive exclusion probability is uni-modal over toxinsensitivity δ , regardless of the resource supply (Fig. 3A). As a reminder, competitive exclusion implies theextinction of species 1 and the survival of species 2. Since toxin sensitivity is identical for both species, whentoxin sensitivity is too high, both species are likely to go extinct (Fig. A.4), which does not count as competitiveexclusion. Instead, competitive exclusion is most likely at lower toxin sensitivities whose value depends on theamount of resource supply: when more resources are supplied, the peak moves to a higher toxin sensitivitybecause species are more likely to survive (Fig. 3A). Hereafter, we refer to the toxin sensitivity that maximizesthe competitive exclusion probability in the absence of environmental switching as the “critical toxin sensitivity”.This analysis clarifies what happens at the two extreme switching rates. At a very slow switching rate ( ν → δ = 0 . . ν → ∞ ), where resources remain at mean abundance, competitiveexclusion has one peak (at δ = 0 . δ = 0 . , .
4, respectively). By increasingenvironmental switching from a very slow rate, the peak at δ = 0 . δ = 0 . δ = 0 . δ = 0 . , . δ = 0 .
6, the environmentalswitching rate also non-monotonically affects the competitive exclusion probability, but for a different reason.Toxin sensitivities between the critical values under mean and abundant resource supplies can have a “valley”over the switching rate: at such toxin sensitivities, competitive exclusion probability changes in a U shape overthe switching rate (species interactions change in a humped shape, see Fig. 2B). At toxin sensitivities largerthan the critical value under the abundant resource supply, competitive exclusion probability is very small anddoes not change non-monotonically over the switching rate because both species frequently go extinct (Fig A.4).
We have shown that using a given set of parameters, the rugged landscape shown in Fig. 2C causes competitionto either increase, decrease or vary non-monotonically across switching rates, depending on toxin sensitivity.We next explore the generality of this finding and ask whether it will hold under different scenarios. In theappendix, we explore scenarios where (i) switching occurs in toxin rather than resource supplies, where (ii) bothresource and toxin supplies switch (Table 1, see Appendix 4), or where (iii) we change the amounts of scarceand abundant resource supplies (Appendix 5).In all these scenarios, the landscapes of competitive exclusion probability also contain two “mountain ranges”(Figs. A.5, A.7), leading, as before, to a variety of patterns as switching rates and toxin sensitivities change.In each scenario, the distances between the three critical toxin sensitivities (mild, mean and harsh, e.g. blackarrow in Fig. 3A) change (see Table A.2). Do these distances predict the shape of the landscape and theprobability of observing non-monotonic behavior? In Table A.2 and Fig. 3B we show that the distance betweencritical sensitivities under harsh and mean environments (i.e., very fast environmental switching) correlatespositively with the likelihood of observing non-monotonic effects of the switching rate on competition (Fig.3B, black circles; Spearman’s ρ = 0 .
77, P-value: 0 . ρ = − .
22, P-value: 0 .
64, and ρ = 0 .
42, P-value: 0 . .5 Beta diversity changes similarly to competitive exclusion In the previous sections, we have focused on interactions between two species and the conditions under whichthey may drive each other extinct. Ultimately, however, our interest is to predict how species diversity is affectedby the environment in communities comprised of tens, hundreds or even thousands of species. In this section,we ask whether the same rules that we have uncovered between two species apply for larger communities of upto 10 species. Will the diversity of a community similarly depend on environmental switching rates and toxinsensitivities?We address this question by extending our model to simulate communities of between 2 and 10 species.Simulating larger communities with this model would be too computationally expensive. As our model includesstochastic environmental switching and demographic noise, we ran 10’000 simulations for each community sizeand for each of 5 different mean toxin sensitivities (here species in a community could differ in their toxinsensitivity) ranging from ¯ δ = 0 . δ = 1 (see Methods). At the end of each simulation run, we measured thebeta diversity and species richness (number of surviving species). These 10’000 (100 × δ = 0 . . δ = 0 .
2) or minimum (¯ δ = 0 .
6) values at intermediate switching rates. At one mean toxin sensitivity ¯ δ = 0 . δ = 0 .
4, both species cansometimes co-exist, which keeps beta-diversity high. Although two species can also coexist when ¯ δ = 0 . . δ = 0 . . in Dand E columns in Fig. 4). In such cases, different species fixate in each simulation by competitive exclusion. Insum, estimating the competitive exclusion probability in two species in a given environment is a good predictorfor the beta diversity of larger communities under those same environmental conditions. Understanding how species diversity in microbial communities arises and is maintained is a central questionin microbial ecology and evolution. Diversity is affected by biotic interactions between species, but also byabiotic properties of the environment, such as its harshness (availability of nutrients and/or toxicity), and howconditions fluctuate over time. Here we have used a mathematical model to explore the relationship between12hese three factors and how they affect the diversity of small communities.Our study is centred on two main findings. First, we show that the rate at which resource abundancefluctuates in the environment changes the strength of negative interactions (competition) – quantified as theability of a slower-growing species to drive a fast-growing one extinct. However, these changes depend stronglyon how toxic the environment is for both species. As the switching rate increases, competition can monotonicallyincrease, monotonically decrease, or change non-monotonically (Fig. 2B), depending on toxin sensitivity. Thesenegative interactions mostly manifest themselves as competitive exclusion, rather than the extinction of bothspecies (Figs. 2C and A.2). By calculating the critical toxin sensitivities that maximize competitive exclusionin the absence of environmental switching, one can predict how the switching rate will affect the strength ofcompetition at a given toxin sensitivity. Our second main finding is that studying changes in the probability oftwo species competitively excluding one another is a good indicator for how the beta diversity of a communitycomposed of up to 10 species will change with environmental switching (Fig. 4).This brings us to a hypothesis that has been debated at length in ecology: the intermediate disturbancehypothesis (Connell, 1978; Grime, 1973), which states that intermediate intensity and frequency of disturbancemaximize species diversity. Fox (2013) argues that the intermediate disturbance hypothesis should be abandonedbecause many examples disagree with it (Mackey and Currie, 2001; Miller et al., 2011). In our model, fluctuationsin resource and toxin concentrations can be regarded as disturbances. In agreement with Mackey and Currie(2001) and Miller et al. (2011) then, an intermediate intensity (i.e., toxin sensitivity) or disturbance frequency(environmental switching rate) does not always maximize beta diversity: our analysis shows that intermediatefrequencies of disturbance maximize beta diversity only when mean toxin sensitivity is within a certain range.Mean toxin sensitivities at the two thresholds of this range show that beta diversity monotonically decreases orincreases over the switching rate. These thresholds depend on scenarios of environmental switching and amountsof resource supplies because these parameters change probability of competitive exclusion (see Appendix 4 andAppendix 5). High diversity at intermediate disturbances is then a consequence of a change in environmentalconditions and not expected to apply generally.The relationship between environmental fluctuations and species diversity is also an important question inmodern coexistence theory, which predicts that fluctuations will affect species coexistence by changing speciesgrowth rates when rare (Chesson, 2000a,b; Barab´as et al., 2018; Ellner et al., 2019; Letten et al., 2018a).Compared to the approach taken in this study – where we ask how many and which species persist at the end ofa long but fixed time frame (i.e., for a quasi-stationary distribution, Fig. 4) – modern coexistence theory allowsone to analyze whether or for how long a set of species will all coexist (Schreiber et al., 2020). An interestingfuture direction would be to apply modern coexistence theory to investigate how environmental fluctuationrates and toxin sensitivities affect the duration of all-species coexistence. This approach would help to proposebiological mechanisms behind species coexistence in our setup.Of course, our model makes some simplifying assumptions and has some limitations. First, we used arbitrarytime units, which in practice can be considered to be hours, corresponding to typical bacterial growth rates inrelevant experiments (Novick and Szilard, 1950; Lin et al., 2002; Zhao and Lin, 2003). This implies that species13nteractions and beta diversity will vary when environmental switching ranges from hourly ( ν = 10 ) to aboutonce every four days ( ν = 10 − ) on average, which is shorter than in some experimental studies (Benneir andLenski, 1999; Rodr´ıguez-Verdugo et al., 2019; Chen and Zhang, 2020) but not impractical. That said, under thisassumption, changing an hourly to a monthly scale, for example, would have different effects on ecological dy-namics. Second, the switching rate is assumed to be symmetric between abundant and scarce resources (and/ortoxin supplies, see Appendix 4). One may include more than two environmental conditions or assume asymmet-ric environmental switching, as in work by Taitelbaum et al. (2020), who show that asymmetric switching cannon-trivially change the effects of environmental switching. Third, our model focuses on competitive exclusionbut other types of interactions can also affect diversity (Rodr´ıguez-Verdugo et al., 2019). Positive interactionsbetween pairs of species (e.g., cross-feeding), for example, might increase alpha and gamma diversities, becausesuch interactions enable species to coexist (Sun et al., 2019). This could result in an increase in beta diversitybecause the extinction of one species by demographic noise, for example, increases its partner species’ extinctionprobability. Finally, our community analysis considers up to ten microbial species, which is orders of magnitudebelow the size of natural microbial communities, according to genomic sampling (Gans, 2005; Roesch et al.,2007). However, it may also be reasonable to assume that species live in structured environments where theycannot possibly interact with more than a handful of other genotypes (Tecon et al., 2019). This suggests thata 10-species community may already be biologically meaningful.In conclusion, the time scale of environmental switching affects the strength of species interactions, whichresults in changing beta diversity via competitive exclusion. In addition, how species interactions and betadiversity change with respect to the environmental switching rate varies and depends on how the species aresensitive species to the environmental toxicity. This may be one mechanism which explains why the intermediatedisturbance hypothesis does not always hold. The variety with which species interactions and beta diversitychange over the switching rate means that it will be very difficult to predict how a given ecosystem, such asthe gut microbiome, will behave under environmental fluctuations. Nevertheless, the similarity between howcompetitive exclusion plays out between two species and beta diversity at the community level means thatanalysing how the former changes over switching rates (as in Rodr´ıguez-Verdugo et al. (2019), for example)may be sufficient to predict the latter. More generally, our study, along with others (Mackey and Currie, 2001;Miller et al., 2011) show that predicting how the environment shapes diversity in ecology remains a challengingproblem. Nevertheless, models like the one presented here can help to disentangle the environmental factorscontributing to species diversity. Acknowledgement
S.S. is funded by University of Lausanne and Nakajima foundation. S.M. is funded by European ResearchCouncil Starting Grant 715097 and the University of Lausanne. The authors declare no conflict of interest.14 ata accessibility
The programming codes and csv files for this manuscript are available in Github.15 eferences
Allen, L. J. S.
An Introduction to Stochastic Processes with Applications to Biology . Chapman andHall/CRC, New York, 2nd edition, 2010. ISBN 9780429184604. doi: 10.1201/b12537. URL .Amarasekare, P. The evolution of coexistence theory.
Theoretical Population Biology , 133:49–51, 2019. ISSN10960325. doi: 10.1016/j.tpb.2019.09.005. URL https://doi.org/10.1016/j.tpb.2019.09.005 .Barab´as, G., D’Andrea, R., and Stump, S. M. Chesson’s coexistence theory.
Ecological Monographs , 0(0):1–27,2018. ISSN 00129615. doi: 10.1002/ecm.1302. URL http://doi.wiley.com/10.1002/ecm.1302 .Bena, I. Dichotomous Markov noise: Exact results for out-of-equilibrium systems. A review.
InternationalJournal of Modern Physics B , 20(20):2825–2888, 6 2006. ISSN 02179792. doi: 10.1142/S0217979206034881.URL http://arxiv.org/abs/cond-mat/0606116 .Benneir, A. F. and Lenski, R. E. Experimental evolution and its role in evolutionary physiology.
AmericanZoologist , 39(2):346–362, 1999. ISSN 00031569. doi: 10.1093/icb/39.2.346.Chao, A., Chiu, C. H., Hsieh, T. C., and Inouye, B. D. Proposing a resolution to debates on diversity partitioning.
Ecology , 93(9):2037–2051, 9 2012. ISSN 00129658. doi: 10.1890/11-1817.1. URL http://doi.wiley.com/10.1890/11-1817.1 .Chen, P. and Zhang, J. Antagonistic pleiotropy conceals molecular adaptations in changing environments.
Nature Ecology & Evolution , 2 2020. ISSN 2397-334X. doi: 10.1038/s41559-020-1107-8. URL .Chesson, P. Multispecies Competition in Variable Environments.
Theoretical Population Biology , 45:227–276,1994.Chesson, P. Mechanisms of Maintenance of Species Diversity.
Annual Review of Ecology and System-atics , 31(1):343–366, 11 2000a. ISSN 0066-4162. doi: 10.1146/annurev.ecolsys.31.1.343. URL .Chesson, P. General theory of competitive coexistence in spatially-varying environments.
Theoretical PopulationBiology , 58(3):211–237, 2000b. ISSN 00405809. doi: 10.1006/tpbi.2000.1486.Cignarella, F., Cantoni, C., Ghezzi, L., Salter, A., Dorsett, Y., Chen, L., Phillips, D., Weinstock, G. M., Fontana,L., Cross, A. H., Zhou, Y., and Piccio, L. Intermittent Fasting Confers Protection in CNS Autoimmunity byAltering the Gut Microbiota.
Cell Metabolism , 27(6):1222–1235, 2018. ISSN 19327420. doi: 10.1016/j.cmet.2018.05.006. URL https://doi.org/10.1016/j.cmet.2018.05.006 .Connell, J. H. Diversity in Tropical Rain Forests and Coral Reefs.
Science , 199(4335):1302–1309, 1978.16oyte, K. Z., Schluter, J., and Foster, K. R. The ecology of the microbiome: Networks, competition, andstability.
Science , 350(6261):663–666, 2015. ISSN 0036-8075. doi: 10.1126/science.aad2602. URL .Davenport, E. R., Mizrahi-Man, O., Michelini, K., Barreiro, L. B., Ober, C., and Gilad, Y. Seasonal variationin human gut microbiome composition.
PLoS ONE , 9(3), 2014. ISSN 19326203. doi: 10.1371/journal.pone.0090731.Ellner, S. P., Snyder, R. E., Adler, P. B., and Hooker, G. An expanded modern coexistence theory for empiricalapplications.
Ecology Letters , 22(1):3–18, 1 2019. ISSN 1461-023X. doi: 10.1111/ele.13159. URL https://onlinelibrary.wiley.com/doi/abs/10.1111/ele.13159 .Fox, J. W. The intermediate disturbance hypothesis should be abandoned.
Trends in Ecology & Evolution , 28(2):86–92, 2 2013. ISSN 01695347. doi: 10.1016/j.tree.2012.08.014. URL https://linkinghub.elsevier.com/retrieve/pii/S0169534712002091 .Gans, J. Computational Improvements Reveal Great Bacterial Diversity and High Metal Toxicity in Soil.
Science , 309(5739):1387–1390, 8 2005. ISSN 0036-8075. doi: 10.1126/science.1112665. URL .Gillespie, D. T. Exact stochastic simulation of coupled chemical reactions.
Journal of Physical Chemistry , 81(25):2340–2361, 1977. ISSN 00223654. doi: 10.1021/j100540a008.Grime, J. P. Competitive Exclusion in Herbaceous Vegetation.
Nature , 242(5396):344–347, 3 1973. ISSN0028-0836. doi: 10.1038/242344a0. URL .Grime, J. P. Evidence for the Existence of Three Primary Strategies in Plants and Its Relevance to Ecologicaland Evolutionary Theory.
The American Naturalist , 111(982):1169–1194, 1977. ISSN 0003-0147. doi: 10.1086/283244.Hengge-Aronis, R. Survival of hunger and stress: The role of rpoS in early stationary phase gene regulationin E. coli.
Cell , 72(2):165–168, 1 1993. ISSN 00928674. doi: 10.1016/0092-8674(93)90655-A. URL https://linkinghub.elsevier.com/retrieve/pii/009286749390655A .Himeoka, Y. and Mitarai, N. Dynamics of bacterial populations under the feast-famine cycles. arXiv , pages1–39, 10 2019. URL http://arxiv.org/abs/1910.05673 .Hoek, T. A., Axelrod, K., Biancalani, T., Yurtsev, E. A., Liu, J., and Gore, J. Resource Availability Modulatesthe Cooperative and Competitive Nature of a Microbial Cross-Feeding Mutualism.
PLOS Biology , 14(8):e1002540, 8 2016. ISSN 1545-7885. doi: 10.1371/journal.pbio.1002540. URL http://dx.plos.org/10.1371/journal.pbio.1002540 . 17orsthemke, W. and Lefever, R.
Noise-Induced Transitions , volume 15 of
Springer Series in Synergetics .Springer Berlin Heidelberg, 2nd edition, 4 2006. ISBN 978-3-540-11359-1. doi: 10.1007/3-540-36852-3. URL http://link.springer.com/10.1007/3-540-36852-3 .Hufton, P. G., Lin, Y. T., Galla, T., and McKane, A. J. Intrinsic noise in systems with switching environments.
Physical Review E , 93(5):1–13, 2016. ISSN 24700053. doi: 10.1103/PhysRevE.93.052119.Jost, L. PARTITIONING DIVERSITY INTO INDEPENDENT ALPHA AND BETA COMPONENTS.
Ecol-ogy , 88(10):2427–2439, 10 2007. ISSN 0012-9658. doi: 10.1890/06-1736.1. URL http://doi.wiley.com/10.1890/06-1736.1 .Letten, A. D., Dhami, M. K., Ke, P. J., and Fukami, T. Species coexistence through simultaneous fluctuation-dependent mechanisms.
Proceedings of the National Academy of Sciences of the United States of America ,115(26):6745–6750, 2018a. ISSN 10916490. doi: 10.1073/pnas.1801846115.Letten, A. D., Dhami, M. K., Ke, P. J., and Fukami, T. Species coexistence through simultaneous fluctuation-dependent mechanisms.
Proceedings of the National Academy of Sciences of the United States of America ,115(26):6745–6750, 2018b. ISSN 10916490. doi: 10.1073/pnas.1801846115.Li, G., Xie, C., Lu, S., Nichols, R. G., Tian, Y., Li, L., Patel, D., Ma, Y., Brocker, C. N., Yan, T., Krausz, K. W.,Xiang, R., Gavrilova, O., Patterson, A. D., and Gonzalez, F. J. Intermittent Fasting Promotes White AdiposeBrowning and Decreases Obesity by Shaping the Gut Microbiota.
Cell Metabolism , 26(4):672–685, 2017. ISSN19327420. doi: 10.1016/j.cmet.2017.08.019. URL http://dx.doi.org/10.1016/j.cmet.2017.08.019 .Lin, Y. H., Bayrock, D. P., and Ingledew, W. M. Evaluation of Saccharomyces cerevisiae grown in a multistagechemostat environment under increasing levels of glucose.
Biotechnology Letters , 24(6):449–453, 2002. doi:10.1023/A:1014501125355.Mackey, R. L. and Currie, D. J. The diversity-disturbance relationship: Is it generally strong and peaked?
Ecology , 82(12):3479–3492, 2001. ISSN 00129658. doi: 10.1890/0012-9658(2001)082[3479:TDDRII]2.0.CO;2.Merritt, J. and Kuehn, S. Frequency- and Amplitude-Dependent Microbial Population Dynamics during Cyclesof Feast and Famine.
Physical Review Letters , 121(9):098101, 8 2018. ISSN 0031-9007. doi: 10.1103/PhysRevLett.121.098101. URL https://link.aps.org/doi/10.1103/PhysRevLett.121.098101 .Miller, A. D., Roxburgh, S. H., and Shea, K. How frequency and intensity shape diversity-disturbance relation-ships.
Proceedings of the National Academy of Sciences of the United States of America , 108(14):5643–5648,2011. ISSN 10916490. doi: 10.1073/pnas.1018594108.Mougi, A. and Kondoh, M. Diversity of Interaction Types and Ecological Community Stability.
Science , 337(6092):349–351, 2012. ISSN 0036-8075. doi: 10.1126/science.1220529. URL . 18ovick, A. and Szilard, L. Experiments with the Chemostat on Spontaneous Mutations of Bacteria.
Proceedingsof the National Academy of Sciences , 36(12):708–719, 12 1950. ISSN 0027-8424. doi: 10.1073/pnas.36.12.708.URL .Novozhilov, A. S., Karev, G. P., and Koonin, E. V. Biological applications of the theory of birth-and-deathprocesses.
Briefings in Bioinformatics , 7(1):70–85, 3 2006. ISSN 1467-5463. doi: 10.1093/bib/bbk006. URL https://academic.oup.com/bib/article/7/1/70/263777 .Piccardi, P., Vessman, B., and Mitri, S. Toxicity drives facilitation between 4 bacterial species.
Proceedings of theNational Academy of Sciences , 116(32):15979–15984, 2019. ISSN 0027-8424. doi: 10.1073/pnas.1906172116.Rodr´ıguez-Verdugo, A., Vulin, C., and Ackermann, M. The rate of environmental fluctuations shapes ecologicaldynamics in a two-species microbial system.
Ecology Letters , 22(5):838–846, 5 2019. ISSN 1461-023X. doi:10.1111/ele.13241. URL http://doi.wiley.com/10.1111/ele.13241 .Roesch, L. F., Fulthorpe, R. R., Riva, A., Casella, G., Hadwin, A. K., Kent, A. D., Daroub, S. H., Camargo,F. A., Farmerie, W. G., and Triplett, E. W. Pyrosequencing enumerates and contrasts soil microbial diversity.
ISME Journal , 1(4):283–290, 2007. ISSN 17517362. doi: 10.1038/ismej.2007.53.Schreiber, S., Levine, J., Godoy, O., Kraft, N., and Hart, S. Does deterministic coexistence theory matter in afinite world? bioRxiv , 2020. doi: 10.1101/290882.Smits, S. A., Leach, J., Sonnenburg, E. D., Gonzalez, C. G., Lichtman, J. S., Reid, G., Knight, R., Manjurano,A., Changalucha, J., Elias, J. E., Dominguez-Bello, M. G., and Sonnenburg, J. L. Seasonal cycling in the gutmicrobiome of the Hadza hunter-gatherers of Tanzania.
Science , 357(6353):802–805, 2017. ISSN 10959203.doi: 10.1126/science.aan4834.Srinivasan, S. and Kjelleberg, S. Cycles of famine and feast: The starvation and outgrowth strategies of amarine Vibrio.
Journal of Biosciences , 23(4):501–511, 1998. ISSN 02505991. doi: 10.1007/BF02936144.Sun, Z., Koffel, T., Stump, S. M., Grimaud, G. M., and Klausmeier, C. A. Microbial cross-feeding promotesmultiple stable states and species coexistence, but also susceptibility to cheaters.
Journal of TheoreticalBiology , 465:63–77, 2019. ISSN 00225193. doi: 10.1016/j.jtbi.2019.01.009. URL https://linkinghub.elsevier.com/retrieve/pii/S0022519319300104 .Taitelbaum, A., West, R., Assaf, M., and Mobilia, M. Population Dynamics in a Changing Environment:Random versus Periodic Switching.
Physical Review Letters , 125(4):048105, 7 2020. ISSN 0031-9007. doi:10.1103/PhysRevLett.125.048105. URL https://link.aps.org/doi/10.1103/PhysRevLett.125.048105 .Tecon, R., Mitri, S., Ciccarese, D., Or, D., van der Meer, J. R., and Johnson, D. R. Bridging the Holistic-Reductionist Divide in Microbial Ecology. mSystems , 4(1):17–21, 2019. ISSN 2379-5077. doi: 10.1128/msystems.00265-18. 19haiss, C. A., Zeevi, D., Levy, M., Zilberman-Schapira, G., Suez, J., Tengeler, A. C., Abramson, L., Katz,M. N., Korem, T., Zmora, N., Kuperman, Y., Biton, I., Gilad, S., Harmelin, A., Shapiro, H., Halpern,Z., Segal, E., and Elinav, E. Transkingdom control of microbiota diurnal oscillations promotes metabolichomeostasis.
Cell , 159(3):514–529, 2014. ISSN 10974172. doi: 10.1016/j.cell.2014.09.048. URL http://dx.doi.org/10.1016/j.cell.2014.09.048 .Vasi, F., Travisano, M., and Lenski, R. E. Long-term experimental evolution in Escherichia coli. II.Changes inlife-history traits during adaptation to a seasonal environment.
American Naturalist , 144(3):432–456, 1994.ISSN 00030147. doi: 10.1086/285685. URL .West, R. and Mobilia, M. Fixation properties of rock-paper-scissors games in fluctuating populations.
Journalof Theoretical Biology , 491:110135, 4 2020. ISSN 00225193. doi: 10.1016/j.jtbi.2019.110135. URL http://dx.doi.org/10.1016/j.jtbi.2019.110135 .West, R., Mobilia, M., and Rucklidge, A. M. Survival behavior in the cyclic Lotka-Volterra model with arandomly switching reaction rate.
Physical Review E , 97(2):022406, 2 2018. ISSN 2470-0045. doi: 10.1103/PhysRevE.97.022406. URL https://link.aps.org/doi/10.1103/PhysRevE.97.022406 .Wienand, K., Frey, E., and Mobilia, M. Evolution of a Fluctuating Population in a Randomly SwitchingEnvironment.
Physical Review Letters , 119(15):158301, 10 2017. ISSN 0031-9007. doi: 10.1103/PhysRevLett.119.158301. URL https://link.aps.org/doi/10.1103/PhysRevLett.119.158301 .Wienand, K., Frey, E., and Mobilia, M. Eco-evolutionary dynamics of a population with randomly switchingcarrying capacity.
Journal of the Royal Society Interface , 15(145):20180343, 8 2018. ISSN 17425662. doi:10.1098/rsif.2018.0343. URL https://royalsocietypublishing.org/doi/10.1098/rsif.2018.0343 .Xavier, J. B., Picioreanu, C., and Van Loosdrecht, M. C. M. A framework for multidimensional modellingof activity and structure of multispecies biofilms.
Environmental Microbiology , 7(8):1085–1103, 2005. ISSN14622912. doi: 10.1111/j.1462-2920.2005.00787.x.Zhao, Y. and Lin, Y. H. Growth of Saccharomyces cerevisiae in a chemostat under high glucose conditions.
Biotechnology Letters , 25(14):1151–1154, 2003. ISSN 01415492. doi: 10.1023/A:1024577414157.Zu˜niga, C., Li, C.-T., Yu, G., Al-Bassam, M. M., Li, T., Jiang, L., Zaramela, L. S., Guarnieri, M., Betenbaugh,M. J., and Zengler, K. Environmental stimuli drive a transition from cooperation to competition in syntheticphototrophic communities.
Nature Microbiology , 4(12):2184–2191, 12 2019. ISSN 2058-5276. doi: 10.1038/s41564-019-0567-6. URL .20
TTT TT abundant R scarce R ¯ ? TR R S1 S2
R R R R R
A B ν R Figure 1: Schematic illustration of the model
A chemostat model with environmental switching. A: Interaction network when N = 2. A → B represents that Aincreases B while A ⊥ B represents that A decreases B. Two species compete for the same resource (R in a circle) butare killed by the same toxic compound (T in a triangle). As a proxy for species interactions, we follow the net effect ofspecies 2 on species 1 (large arrow from species 2 to 1) – due to resource competition plus facilitation by detoxification– which can be environmentally mediated. B: An example of a chemostat model with environmental switching and N = 2. Environmental switching is realized by changing the media flowing into a chemostat. In this example, thecurrent environmental condition is ξ = 1 (abundant resource supply R +1 ). BC D
Figure 2: Species interaction strength changes differently over the switching rate
A: The effect of species 2 on species 1’s extinction probability ∆ P ( s ( T end ) = 0) (Eq 7) changes over the switching rate ν and the toxin sensitivity δ = δ = δ . In most of the parameter space, species 2 has a negative effect on species 1.B: Some illustrative examples from panel A plotted differently to show how ∆ P changes over the switching rate atgiven toxin sensitivities. The difference in extinction probability can monotonically increase ( δ = 0 .
1, green),monotonically decrease ( δ = 0 .
4, red), or non-monotonically change with a minimum ( δ = 0 .
2, purple) or a maximum( δ = 0 .
6, blue) value at an intermediate switching rate. C: Probability that species 2 persists but species 1 goes extinct(i.e., competitive exclusion) over switching rate and toxin sensitivity. D: competitive exclusion probabilities over thetoxin sensitivity, when the environmental switching rate is slow ( ν = 10 − ), intermediate ( ν = 10 − ), or fast ( ν = 10 ). A Figure 3: Competitive exclusion explains changes in species interactions
Analysis on competitive exclusion predicts at which toxin sensitivity the species interaction can changenon-monotonically. A: In the absence of environmental switching, the competitive exclusion probability (i.e.,probability that species 2 excludes species 1) is uni-modal over the switching rate. The toxin sensitivities giving thepeak values (critical toxin sensitivities) depend on the resource supply: scarce R = R − (green), mean R = (cid:104) R (cid:105) (red),or abundant R = R +1 (blue). B: The number of times we observe non-monotonic species interactions across theexplored parameter range changes with the distance between the two critical toxin sensitivities, depending on wheredistances are measured (between harsh and mean environments: dots, between mean and mild: environment diamonds,and between harsh and mild environments: crosses). These three distances were measured in each of the followingseven scenarios: three different scenarios of environmental switching (Table 1 and Appendix 4) and four environmentalswitching scenario 1s with changing amounts of resource supplies (Appendix 5). The correlation is only significantlypositive for the distance between scarce resource or abundant toxin supplies (i.e., harsh environments) and meanresource/toxin supplies (Spearman’s ρ = 0 .
77, P-value: 0 . o x i n s e n s i t i v i t y . . . . . b e t a d i v . r i c h n e ss c o m p . e x c l . s p e c i e s r i c h n e ss b e t a d i v . r i c h n e ss t e n - s p e c i e s t w o - s p e c i e s A E D C B F i g u r e : C o m p a r i s o n o f c o m p e t i t i v ee x c l u s i o n a ndb e t a d i v e r s i t y Competitive exclusion probabilities (column A), beta diversities (column B), and species richness (column C) over theswitching rate in two-species communities. Here, competitive exclusion refers the event where the slower-growingspecies excludes the faster-growing species. Columns D and E show the beta diversity and species richness inten-species communities, respectively. In the plots of competitive exclusion probability and beta diversity, the blacklines show the means and blue areas represent the probability distributions calculated from 10’000 simulations (100beta diversity and each of them from 100 replicate runs). Each color in the species richness plots represents theproportion of 10’000 runs where at the end of the run there were that number of species surviving. Each row representsthe mean toxin sensitivity ¯ δ of communities in those runs. R i ( ξ = 1) R i ( ξ = − T j ( ξ = 1) T j ( ξ = − R + i R − i (cid:104) T j (cid:105) (cid:104) T j (cid:105) (cid:104) R i (cid:105) (cid:104) R i (cid:105) T − j T + j R + i R − i T − j T + j upplementary information Appendix 1 Analysis excluding noise
In this section, we analyse the equilibrium state of a two-species, one-resource, and one-toxin model, neglectingenvironmental switching and demographic noise. Although an equilibrium state cannot be analytically obtained,we shall see that there exists at most one stable and feasible (species abundance s k >
0) equilibrium state. Byneglecting demographic noise and environmental switching from Eq (5) in the main text, the dynamics aregoverned by the following ordinary differential equations:˙ r = α ( R − r ) − (cid:88) k =1 , µ k Y r k r r + K r k s k , (A.1a)˙ t = α ( T − t ) − (cid:88) k =1 , δ k Y t k t t + K t k s k , (A.1b)˙ s k = (cid:18) µ k r r + K r k − δ k t t + K t k − α (cid:19) s k , (A.1c)When only species k persists in the system, a feasible equilibrium state of Eqs (A.1a) -(A.1c), ( r , t , s k ) =( r ∗ k , t ∗ k , s ∗ k ), should satisfy below: α ( R − r ∗ k ) = µ k Y r k r ∗ k r ∗ k + K r k s ∗ k , (A.2a) α ( T − t ∗ k ) = δ k Y t k t ∗ k t ∗ k + K t k t ∗ k , (A.2b) µ j r ∗ k r ∗ k + K r k = δ k t ∗ k t ∗ k + K t k + α. (A.2c)By rearranging the first and second equations, we see that they represent quadratic functions of r ∗ k and t ∗ k : − r ∗ k + { R − K r k − µ k s ∗ k / ( αY r k ) } r ∗ k + K r k R = 0 , (A.3) − t ∗ k + (cid:8) T − K t k − δ k s ∗ k / ( αY t k ) (cid:9) t ∗ k + K t k T = 0 . (A.4)Notice that K r k R and K t k T are positive. This implies that the above equations always have a unique positiveroot. In other words, we can obtain unique solutions for r ∗ k and t ∗ k once we obtain s ∗ k .By substituting the positive roots of Eqs (A.3) and (A.4) into Eq (A.2c), we obtain the following equationwhose positive root is s ∗ k : f ( s ) = 12 α (cid:110) ( − µ k + 2 α + δ k ) s + (cid:112) Q ( s ) − (cid:112) Q ( s ) + α (cid:0) − K r k Y r k + K t k Y t k − Y r k R + Y t k T (cid:1)(cid:111) (A.5) Here, feasibility means r ∗ k , t ∗ k , s ∗ k > Q ( s ) and Q ( s ) are quadratic functions of s : Q ( s ) = { µ k s + αY r k ( K r k − R ) } + 4 αY r k K r k R > Q ( s ) = (cid:8) δ k s + αY t k (cid:0) K t k − T (cid:1)(cid:9) + 4 αY t k K t k T > . (A.6b)Notice that f ( s ) always has a root s = 0 because f (0) = 12 α (cid:110)(cid:112) Q (0) − (cid:112) Q (0) + α (cid:0) − K r k Y r k + K t k Y t k − Y r k R + Y t k T (cid:1)(cid:111) = 12 α (cid:8) αY r k ( K r k + R ) − αY t k (cid:0) K t k + T (cid:1) + α (cid:0) − K r k Y r k + K t k Y t k − Y r k T + Y t k T (cid:1)(cid:9) = 0 . (A.7)Although Newton’s method numerically provides a root of f ( s ), this root depends on the initial value used inNewton’s method. In addition, as f ( s ) has root s = 0, Newton’s method may provide this root with variousinitial values, which does not always mean that f ( s ) does not have positive roots (i.e., s ∗ k ). In other words, itis recommended to investigate how many positive roots f ( s ) has before using Newton’s method.To investigate the number of f ( s )’s positive roots, it is useful to obtain df /ds : dfds = 12 α (cid:40) ( − µ k + 2 α + δ k ) + dQ /ds (cid:112) Q ( s ) − dQ /ds (cid:112) Q ( s ) (cid:41) . (A.8)Although it is analytically difficult to obtain the solution(s) of df /ds , we can obtain the maximum number ofpositive roots of f ( s ) by analyzing the number of df /ds ’s sign changes. Notice that dQ /ds and dQ /ds arelinear functions of s : dQ ds = 2 µ k { µ k s + αY r k ( K r k − R ) } , (A.9a) dQ ds = 2 δ k (cid:8) δ k s + αY t k (cid:0) K t k − T (cid:1)(cid:9) . (A.9b)In addition, Q ( s ) and Q ( s ) are always positive. Then, the second and third terms of Eq (A.8) change theirsign at most once by increasing s . The maximum number of df /ds ’s sign change is, therefore, two. This impliesthat the maximum number of positive roots of f ( s ) is also two. To obtain the exact number of positive rootsof f ( s ), it is necessary to numerically calculate the root(s) of df /ds . Substituting the root(s) into Eq (A.5) andcalculating the sign of f ( s ), it is possible to obtain the exact value of s ∗ k .Once we have obtained a feasible equilibrium state ( r ∗ k , s ∗ k , t ∗ k ), it is necessary to analyze the stability ofthis equilibrium state. Although the stability analysis requires the evaluation of the Jacobian matrix at eachequilibrium state, we can see that there exists at most one feasible and stable equilibrium state without such2tability analysis. Notice that: ˙ s j ≶ ⇔ µ k r ∗ k r ∗ k + K r k ≶ δ k t ∗ k t ∗ k + K t k + α ⇔ f ( s ) ≷ . (A.10)These inequalities implies that a stable equilibrium state satisfies the following inequality: dfds (cid:12)(cid:12)(cid:12)(cid:12) s = s ∗ k > . (A.11)Although there can exist at most two feasible equilibria, only one of them satisfies the above inequality (Fig.A.1). The number of feasible and stable equilibria is, therefore, one at most. Appendix 2 Details of simulations
Appendix 2.1 Species interaction analysis
Here, we summarize the details of the simulations and parameter values used for species interaction analysis inthe main text. In the analysis of species interactions, we used the minimal model ( N = 2): two species competefor one resource and absorb and are killed by one toxin. As the model represents a stochastic process in acontinuous-time scale, the simulation is implemented by the Gillespie algorithm (Gillespie, 1977). In each run,the initial condition is either ( r (0) , t (0) , s (0) , s (0)) = (150 , , ,
0) or (150 , , , ξ = 1 with a probability of 0.5; otherwise ξ = −
1. Each simulation continues until at time T end = 200 when the distributions of species abundancesconverge to a quasi-stationary distribution and do not change for a long time. The extinction probability isestimated by running 10 simulations for each condition.Table A.1 summarizes the parameter values used in the species interaction analysis. In the main text, species2 is assumed to grow slower than species 1, meaning that species 2 goes extinct in the absence of demographicnoise. For simplicity, we assumed that species 2’s parameter values are the same as those of species 1 exceptfor the maximum growth rate.In the main text, to investigate how species 2’s effect on species 1 changes over the switching rate, we beganby analyzing the difference in the extinction probability of species 1 in the presence/absence of species 2 as aproxy for species interactions, but then moved to using the competitive exclusion probability instead. AlthoughFig. 2A and C shows that these two measures give similar results, Fig. A.2 confirms this conclusion as the ratioof these two measures is around 1. 3 ppendix 2.2 Species diversity analysis In the community diversity analysis, we changed the number of species from N = 2 to N = 10. Some parametervalues were not fixed in this analysis, and we sampled them from the following probability density functions: µ ik ∼ N (cid:0) , . (cid:1) , (A.12a) δ jk ∼ Beta (cid:0) δ, (cid:0) − ¯ δ (cid:1)(cid:1) , (A.12b) K rik , K tjk ∼ N (cid:0) , (cid:1) . (A.12c)Here, each function is uni-modal with a mean of 1 . µ ik , ¯ δ for δ jk , and 100 for K rik , K tjk . For µ ik , K rik , and K tjk , the mean values are the same as in Table A.1 and they are rarely negative due to the small variances. Wesampled δ jk from a beta distribution so that 0 ≤ δ jk ≤ δ jk should be non-negative by definition and not belarger than 1 because a large δ jk is likely to drive species k extinct. We systematically vary the value of themean toxin sensitivity ¯ δ to be ¯ δ = 0 . , . , . , . .
99 in each set of simulations. We used ¯ δ = 0 .
99 insteadof ¯ δ = 1 . δ jk with ¯ δ = 1 .
0. For each numberof species N (2, 4, 6, 8 or 10) and ¯ δ , 100 sets of parameter values are sampled. With each parameter set, weperformed the simulation 100 times until T end = 200 at each value of the environmental switching rate ν . Then,we calculated the species richness and beta diversity as described in the main text.In each run, initial resource amounts, species abundances, and toxin amounts are given by ( r i (0) , t j (0) , s k (0)) =(150 , ,
10) for any i, j, k . As an environmental switching scenario, we chose scenario 1 (Table 1) with R + i = 200, R − i = 50, and (cid:104) T k (cid:105) = 125 for any i and k . The initial environmental condition is ξ = 1 withprobability of 0.5; otherwise ξ = − Appendix 3 Analysis in the absence of environmental switching
The analysis with an environmental switching rate of 0 facilitates the understanding of the results in the twoextreme switching rates ( ν → ν → ∞ ) in the main text. In this section, we analyze species 2’s effect onspecies 1 when the resource supply is fixed to be scarce, mean, or abundant.In the absence of environmental switching, resource supply and toxin sensitivity change species 2’s effect onspecies 1 (Fig. A.3): each resource supply has a unique maximum negative interaction, which is reached at acertain toxin sensitivity. These toxin sensitivities correspond to the peak sensitivities of competitive exclusionprobabilities (Fig. 3A). Species 2’s effect on species 1 converges to zero when the toxin sensitivities exceed thepeak sensitivities as competitive exclusion probabilities decrease (see section 3.3). This is because both species1 and species 2 go extinct at too high toxin sensitivity (Fig. A.4).Note that species 2 has a positive effect on species 1 when the resource supply is abundant and toxinsensitivity is 1 .
0. In this case, the benefits of absorbing toxic compounds by species 2 overcome the costs ofcompeting for the resource. This is why in Fig. 2A, species 2 has a positive effect on species 1 when theswitching rate is very small ( ν ≤ − ) and the toxin sensitivity is 1 . ppendix 4 Alternative environmental switching scenarios In the main text, environmental switching affects only the resource supply (scenario 1) while the amount oftoxin supply is fixed. Here, the results of other environmental switching scenarios are shown: In scenario2, environmental switching affects only the toxin supply, while both resource and toxin supplies change andcorrelate negatively in scenario 3 (see Table 1).In both scenarios 2 and 3, the difference in species 1’s extinction probability ∆ P ( s ( T end ) = 0) is verysimilar to the negative competitive exclusion probability when the sign of ∆ P ( s ( T end ) = 0) is negative (Fig.A.5). This once again confirms that we can use one measure or the other.In addition, in both scenarios, the probabilities of competitive exclusion are bi-modal across toxin sensitivitiesat very slow environmental switching ν = 10 − , but uni-modal at very fast environmental switching ( ν = 10 )(Fig. A.5C, D). As explained in the main text, when ν →
0, there are no switches and the environmental stateis randomly allocated to harsh or mild conditions at t = 0 with the same probability (the mean of ξ is zero),yielding a bi-modal distribution at low switching rate. In the limit ν → ∞ , there are so many switches thatenvironmental noise averages out, i.e. ξ is replaced by its mean (that is zero). This results in a uni-modaldistribution when ν (cid:29) δ = 0 .
3) is close to that under mean toxin supply ( δ = 0 . δ = 0 .
2, an intermediate switching rate ( ν = 10 − ) shows the minimum difference in extinction probability, andthe maximum probability of competitive exclusion. Although in scenario 3 the same intermediate switching rateminimizes the probability of competitive exclusion at δ = 0 .
6, the same non-monotonic effect is not observed inthe difference in extinction probabilities (Fig. A.5B). Note that the critical toxin sensitivities under the mildenvironments in scenarios 2 (scarce toxin supply) and 3 (abundant resource supply and scarce toxin supply) areslightly larger than 1.0 (Fig. A.6) and therefore not visible in Fig. A.5D. Table A.2 shows the critical toxinsensitivities in each scenario.
Appendix 5 Effects of resource supply
In this section, we again focus on environmental switching scenario 1 (changing only resource supply). We willsee that the amount of resource supplies R ( ξ ) changes the critical toxin sensitivities under scarce, mean, andabundant, resource supplies, affecting the likelihood of the non-monotonic effect of the environmental switchingrate.By increasing the abundant or decreasing the scarce resource supply (Figs. A.7A and D, respectively), thedistance between the critical toxin sensitivities under scarce and mean (or mean and abundant) resource suppliesbecomes larger (Table A.2). Conversely, decreasing abundant or increasing scarce resource supply (Figs. A.7B5nd C, respectively), decreases the distance between the critical toxin sensitivities under scarce and mean (ormean and abundant) resource supplies (Table A.2). Once again, changes in competitive exclusion probability(A.7) match changes in species 2’s effect on species 1 (Fig. A.8).Analyzing competitive exclusion probability instead of the difference in extinction probability is valid onlywhen the difference in extinction probability is negative: if it is positive, the second line in (11) cannot beignored. The sign of the difference can, however, become positive at δ = 1 .
0. Indeed, when R +1 = 400 and δ = 1 .
0, species 2 has a positive effect on species 1 whose strength varies non-monotonically with the rate ofenvironmental switching (A.9A). In this case, we analyze both (i) the competitive exclusion probability (thefirst line in Eq (11), Fig. A.9B) and (ii) the difference in species 1’s extinction probability in mono-culture andboth species extinction in co-culture (the second line in Eq (11), Fig. A.9C). The effects of the environmentalswitching rate on (i) and (ii) are similar, leading to similar non-monotonic effects of species 2 on species 1. Insum, non-monotonic effects of environmental switching rates on species interactions can be observed whetherthese interactions are positive or negative. Although the main text explains why non-monotonic effects ofenvironmental switching rates on species interactions happen when interactions are negative, it remains unclearwhy such non-monotonic changes happen when species interactions are positive.
Appendix 6 Diversity in communities of increasing species number
In the main text, we show the distributions of beta diversity and species richness in communities when theinitial number of species N is two or ten. This section shows the results of intermediates values of N = 4 , , δ = 0 . N in a community affects the maximum values of species richness (Fig. A.11).However, how the switching rate and the mean toxin sensitivity affect the distribution of species richnesswas consistent for different values of N . In particular, increasing the mean toxin sensitivity decreases speciesrichness. The effects of the environmental switching rate also depend on the mean toxin sensitivity: at meantoxin sensitivity 1 .
0, species in all community sizes are more likely to go extinct as the switching rate increases,while the likelihood of all species going extinct consistently shows a humped shape at toxin sensitivity 0 . .
6. 6igure A.1: Feasible and stable equilibrium state
An example of f ( s ) and equilibrium states. In this example, f ( s ) has three roots: s = 0, and two feasible equilibria.The blue arrows indicate that s decreases or increases when f ( s ) is positive or negative, respectively. As df/ds isnegative at the smaller feasible equilibrium state, this equilibrium state is unstable. On the other hand, the largerfeasible equilibrium has a positive df/ds . Note that the equilibrium state corresponding to s = 0 can be also stable.Parameter values are δ k = 1 . R = 200, T = 125 and as in Table A.1 otherwise. -competitive exclusion/ P ( s ( T end ) = 0) f r e q u e n c y Figure A.2: Competitive exclusion probability relative to the difference in extinction probability of species 1 inthe presence/absence of species 2
A histogram of the ratio of the competitive exclusion probability and the difference in species 1’s extinction probabilitiesalone versus in the presence of species 1, showing the results of 81 sets of the environmental switching rate ν and thetoxin sensitivity δ (9 values of ν = 10 − , . . . , and 9 values of δ = 0 . , . . . , . simulations were run to calculate the competitive exclusion probability and the difference in species 1’s extinctionprobabilities in the presence/absence of species 2. In many of these 81 parameter sets, this ratio is close to 1, indicatingthat both measures yield similar results. As in the manuscript we focus on conditions leading to competition betweenthe two species, we ignore toxin sensitivity δ = 1 . o x i n s e n s i t i v i t y scarce mean abundant Figure A.3: Species 2’s effect on species 1 in the absence of environmental switching
Species 2’s effect on species 1 when the resource supply is fixed to be scarce ( R − ), mean ( (cid:104) R (cid:105) ), or abundant ( R +1 ).The toxin sensitivities that minimize species 2’s effect on species 1 correspond to the peak sensitivities in Fig. 3A. toxin sensitivity p r o b . o f b o t h e x t i n c t i o n scarcemeanabundant Figure A.4: Probability of both species going extinct in the absence of environmental switching
Probabilities that both species 1 and species 2 go extinct when resource supply is fixed to be scarce ( R − ), mean ( (cid:104) R (cid:105) ),or abundant ( R +1 ). o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y A BC D t o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y Figure A.5: Effects of the environmental switching rate in alternative scenarios
Examples of the effect of switching rate in alternative scenarios. In the left column (A and C), toxin supply isswitching (scenario 2), while both resource and toxin supplies switch and are negatively correlated (scenario 3) in theright column (B, D). A and B: difference between extinction probabilities in absence and presence of species 2. C andD: competitive exclusion probability. Parameter values: R +1 = 200, R − = 50, T +1 = 200, and T − = 50. toxin sensitivity p r o b . o f c o m p e t i t i v e e x c l u s i o n scenario 2 scenario 3 Figure A.6: Critical toxin sensitivities under mild environments
The critical toxin sensitivities (i.e., toxin sensitivity that maximizes the competitive exclusion probability in theabsence of environmental switching) under the mild environments (scenario 2 :scarce toxin supply T − = 50, andscenario 3: abundant resource supply R +1 = 200 and scarce toxin supply) are > C BD t o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y Figure A.7: Effects of resource supplies on competitive exclusion
Top: abundant resource supply becomes twice (A) or half (B) of R +1 = 200, i.e. R +1 = 400 in panel (A) and R +1 = 100in panel (B). Bottom: scarce resource supply becomes twice (C) or half (D) of R − = 50, i.e. R − = 100 in panel (C) and R − = 25 in panel (D). C BD t o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y t o x i n s e n s i t i v i t y Figure A.8: Effects of resource supplies on species interactions
Similar to Fig. A.7 but showing species 2’s effect on species 1’s extinction probability. Top: abundant resource supplybecomes twice (A) or half (B) of R +1 = 200, i.e. R +1 = 400 in panel (A) and R +1 = 100 in panel (B). Bottom: the scarceresource supply becomes twice (C) or half (D) of R − = 50, i.e. R − = 100 in panel (C) and R − = 25 in panel (D). Weplotted toxin sensitivity from 0.1 to 2.0 in panel A to see non-monotonic positive species interactions (see also Fig.A.9 ) A B C p r o b . o f c o m p e t i t i v e e x c l u s i o n m o n o – b o t h e x t i n c t i o n p r o b . Figure A.9: Positive interaction strength varies non-monotonically with switching rate
A: Difference in species 1’s extinction probability with positive sign at δ = 1 . R +1 = 400, showing a non-monotoniceffect of the environmental switching rate. B: Probability of competitive exclusion. C: Difference between theprobability of species 1 going extinct alone and both species going extinct. Parameter values are + s = 400, δ = 1 . o x i n s e n s i t i v i t y . . . . . f o u r - s p e c i e ss i x - s p e c i e s t w o - s p e c i e s b e t a d i v e r s i t y e i g h t - s p e c i e s t e n - s p e c i e s F i g u r e A . : B e t a d i v e r s i t y w i t h i n c r e a s i n g i n i t i a l nu m b e r o f s p ec i e s i n a c o mm un i t y Beta diversities with increasing initial numbers of species N and mean toxin sensitivities ¯ δ . The black lines show themeans and blue areas represent the probability distributions calculated by 10’000 simulations (100 beta diversitymeasurements using different parameter sets, each from 100 replicate runs). o x i n s e n s i t i v i t y . . . . . f o u r - s p e c i e ss i x - s p e c i e s t w o - s p e c i e s s p e c i e s r i c h n e sss p e c i e s r i c h n e ss e i g h t - s p e c i e s t e n - s p e c i e s F i g u r e A . : Sp ec i e s r i c hn e ss w i t h i n c r e a s i n g i n i t i a l nu m b e r o f s p ec i e s i n a c o mm un i t y Species richness with increasing initial numbers of species N and mean toxin sensitivities ¯ δ . Each bar plot representsresults of 10’000 simulations. α . R ± R +1 = 200 , R − = 50 abundant or scarce resource supply concentration T ± T +1 = 200 , T − = 50 abundant or scarce toxin supply concentration Y r k Y r k = 1 for k = 1 , k ’s biomass yields of resource Y t k Y t k = 1 for k = 1 , k ’s biomass yields of toxin µ . µ .
91 maximum growth rate of species 2 on resource 1 δ k [0 . , . . . , .
0] and δ = δ = δ sensitivity of species k to toxin 1. K r k
100 amount of resource 1 that gives half-max growth rate of species kK t k
100 amount of toxin 1 that gives half-max death rate of species k Table A.2: Summary of critical toxin sensitivities and number of times non-monotonicity is observedSwitching scenario 1 2 3 1Amounts of supplies base line ∗ ↑ R +1 ↓ R +1 ↑ R − ↓ R − Critical toxin harsh 0.1 0.3 0.1 0.1 0.1 0.3 0.1sensitivities mean 0.4 0.4 0.4 0.9 0.2 0.5 0.3mild 0.8 1.1 1.3 1.2 0.3 0.8 0.8number of mean and harsh 2 0 2 7 0 1 3non-monotonic changes mean and mild 1 0 0 3 0 0 2observed between mild and harsh 3 0 2 10 0 1 4 † ∗ For exact parameter values, see Table A.1. † At the critical toxin sensitivity corresponding to the mean environment ( δ = 0 . −1 = 4.