The Impact of Environmental Fluctuations on Evolutionary Fitness Functions

The concept of fitness as a measure for a species’ success in natural selection is central to the theory of evolution. We here investigate how reproduction rates which are not constant but vary in response to environmental fluctuations, influence a species’ prosperity and thereby its fitness. Interestingly, we find that not only larger growth rates but also reduced sensitivities to environmental changes substantially increase the fitness. Thereby, depending on the noise level of the environment, it might be an evolutionary successful strategy to minimize this sensitivity rather than to optimize the reproduction speed. Also for neutral evolution, where species with exactly the same properties compete, variability in the growth rates plays a crucial role. The time for one species to fixate is strongly reduced in the presence of environmental noise. Hence, environmental fluctuations constitute a possible explanation for effective population sizes inferred from genetic data that often are much smaller than the census population size.

In this Supporting Material, we provide more details on the calculations leading to the one-dimensional Fokker-Planck equation (FPE), Eq. (3) main text. Moreover, we present calculations for the average fixation time and the extinction probability. We discuss the mapping of the individual based model (IBM) to the Langevin model. For the no-memory limit we present an analytic derivation of the mapping. For other parameter values we give heuristic arguments that are supported by additional data. Finally, we present results for non-neutral evolution investigating the regime in which the white noise approximation is an adequate description for the evolutionary process.

Derivation of the one-dimensional Fokker-Planck equation
In this Section, we provide details how the Langevin equations,Ṅ 1 = N 1 (ν 1 − γ N K ) + N 1 ν 1 + γ N K µ 1 +σ 1 N 1 ξ 1 N 2 = N 2 (ν 2 − γ N K ) + N 2 (ν 2 + γ N K )µ 2 +σ 2 N 2 ξ 2 (S1) can be transformed into the one-dimensional FPE presented in the main text. Both noise terms in Eq. (S1) are interpreted in the Ito sense ? because they only should act as noise, i.e. enter the variances but do not affect the means, as it can be seen in the Fokker-Planck equation (S4) [or (3) main text]. From general results on stochastic processes (see 1 ), it follows that the previous Langevin equation is associated to the following two-dimensional FPE : where ∂ i ≡ ∂ N i . The drift part is directly stemming from the non-fluctuating parts of the Langevin equations N S (ν S − γN/K). Diffusion depends on the correlation level of the noises experienced by the two species. In particular, we have introduced the correlation coefficient ε ≡ ξ 1 ξ 2 / ξ 2 1 ξ 2 2 .
The case when the two noises are the same is given by ε = 1, when they are independent is ε = 0 and when they are anticorrelated is ε = −1.
To study the evolutionary dynamics associated to Eq. (S2), the relative abundances are the natural choice of variables. Therefore, we transform the absolute abundances N 1 and N 2 to x = N 1 N 1 +N 2 and N = N 1 + N 2 . To perform the change of variables, not only N 1 = xN and N 2 = (1 − x)N have to be replaced, also the differential operators and the probability distribution have to be transformed. Ensuring that the latter is still normalized after change of variables, the Jacobian has to be introduced, P(N 1 , N 2 ) → 1 N P(x, N). The derivatives are given by, After the change of variables, the FPE for x and N can now be further simplified exploiting the fact that the time scale of selection, s = ν 1 − ν 2 , is much slower than the one of the population growth ν 1 x + ν 2 (1 − x). 2 Therefore, we marginalize the FPE with respect to the total population size N. Thereby, the integrals ∞ 0 dN of N-derivative terms such as where Q ≡ x(1 − x)P(x,t). The drift term in the second line stemming from demographic fluctuations can be neglected as N 1 holds. To finally arrive at the one-dimensional FPE employed in the main text, we compute the steady state population size N * . As the deterministic differential equation for N is given bẏ the fixed point for the populations size is Employing that relation and the aforementioned condition s ν 1 x + ν 2 (1 − x), the last term in Eq. (S3) can be simplified as ν 1 −sx 2N ≈ γ 2K , which finally leads to the onedimensional FPE in the main text:

Fixation probability
In the following, we derive a general expression for the fixation probability. The calculations are analogous to the procedure for the neutral case described in the body of the paper.
To determine the fixation probability the following backward equation has to be solved, Boundary conditions are P fix (0) = 0 and P fix (1) = 1. The solution to Eq. (S5) for the fixation probability is (S7) The solution (S6) is obtained by integrating (S5) once, to find the gradient The expression (S8) is verified to be proportional to the deriva- and boundary conditions are then imposed to fix the two constants of integration. The resulting expression is finally transformed into Eq. (S6) by using the elementary identity: 2 Tanh −1 (x) = log [(1 + x)/(1 − x)]. It is verified that in the limit ζ → 0, one recovers the expression given in the main text.
All in all, the behavior we discussed in the main text is validated by analyzing the fixation probability: Both a higher growth rate and a smaller variability are beneficial for an individual.

Neutral case
The expression for the time of fixation in the neutral case that we presented in the body of the paper is derived as follows. The average time for fixation obeys the following backward equation, Integrating Eq. (S9) and by variation of constants, we obtain: where A is a constant to be fixed by the boundary conditions. The integrals x 0 of Eq. (S10) needed for T (x) are performed by decomposing the rational function at the prefactor and using the formula : that follows from the very definition of the dilogarithm . The formula (S11) is used four times either directly (with a simple change of variables) or first integrating by parts to satisfy the condition a > 0 in (S11). The resulting expression is then transformed to the form given in the main text (which is the one given by Mathematica) by using the reflection property,

General case
In the general case when selection is present, the expression for the average fixation time cannot be found explicitly but is reducible to quadratures as follows. The fixation time obeys the backward equation (S5) with the left-hand side replaced by −1. Using the definitions (S7), we obtain Boundary conditions are T (0) = T (1) = 0. The homogeneous solution was already found following (S8) and reads where C 1 and C 2 are constants and we defined to simplify notation. The non-homogeneous solution for the gradient of T is obtained by varying the constant in (S8), and integrating the resulting first-order differential equation to obtain where F 1 is the hypergeometric function of two variables. 4 The solution for T involves the integral x 0 ∂ y T part (y) dy of the expression above (for which a closed form does not seem to be available), and the two constants in (S13) are fixed by It is verified from the expression above or directly from the original equation (S12) that in the two limits x → 0 and x → 1 the solution behaves like in the neutral case, i.e. −K/γx log x and −K/γ(1 − x) log(1 − x). Selection and the rest of the parameters affect of course the solution in the rest of the interval of definition x ∈ [0, 1].

Coexistence time
Depending on the position of the stable fixed point, coexistence between two species (one with a larger growth rate, one with a smaller variability, ν 1 > ν 2 and σ 1 = ∆, σ 2 = 0) is possible. In this section we present some additional data demonstrating this. In Fig. S1, the extinction time which corresponds to the time of coexistence is shown depending on ∆ is shown for different values of s. Dots correspond to solutions of Eqs. (S1) and black lines are numerical solutions of Eq. (S12). The extinction time has a maximum which exactly coincides with the parameter values of a fixed point x * = 0.5. The dependence of this maximal extinction time on the selection strength s is shown in Fig. S2.

Mapping individual-based models onto the Langevin dynamics
The aim of this Section is to show that individual-based models are described by the Langevin equations, Eqs. (S1), discussed in the main text and to analyze the mapping between the parameters of the two models. The environmental conditions change stochastically at the rate 1/τ and are distributed according to a distribution, p(E), with mean E and variance Var [E]. The dependency of the instantaneous reproduction rate λ S (E) on E is given by the sigmoidal function : which reduces to φ S ± ω S in the limit of large variances Var [E]. Birth rates are defined as, In the no-memory limit m = 0, the growth rate is therefore given by the instantaneous growth rate λ S (E), while for m → 1 the current growth rate is the arithmetic mean of all previously experienced environments. Death rates are given by Γ i death,S = γN/K.

No memory, m = 0
We discuss first the model without memory, where the memory parameter, m, is zero : Individuals reproduce with the instantaneous reproduction rates [Eq. (S15)], which reduce to φ S ± ω S in the limit of large environmental variance. We consider an interval of length δt τ such that the probability for 3/7 an individual to reproduce or die is small, yet the total number of events occurring over the whole population (∼ K 1) is large. Neglecting the standard demographic noise term, 5 the variation of the S-type population is given by whereσ (s) is the environmental Boolean random variable that takes values ±1 and switches with characteristic time τ.
The last term of Eq. (S17) is estimated as follows where G e is a Gaussian random variable having zero mean and variance Here, we used that σ (t)σ (t ) = e −|t−t |/τ and δt τ. The second term in Eq. (S18) is evaluated at the order δt using the same integral, Eq. (S19), and gives N S (t)ω S τδt. Combining back all the terms, we conclude that the equation (S17) is equivalent to the Langevin equation (S1) with the mapping of the parameters Note that the standard demographic noise term in Eq. (S1) should a priori include the fluctuating environmental term N S (t)ω S G e in the sum of the rates. In fact, it can be safely ignored as φ S N S δt ω S N S √ 2τδt due to φ S ≥ ω S and δt τ. Finally, the factor 2 appearing in σ 2 S in (S20) depends on the Poisson statistics of the environmental fluctuations. If the duration is fixed and equal to τ, Eq. (S19) becomes τδt. In that case, the corresponding mappings are ν S = φ S + ω 2 S τ/2 and σ 2 S = ω 2 S τ. This is confirmed numerically in Fig. S3 where we show data for exponentially distributed (black) and fixed duration (red) environments. Solid lines are analytic solutions of the fixation time [Eq. (6) main text] employing the respective mappings.

Finite memory, m > 0
We now turn to the scenario where memory extends over several environmental conditions that an individual previously experienced. Whilst for m = 0 an exact analytic mapping can be found in the limit of small τ, for finite memories the situation is more intricate. However, for the special case m = 1 the variability in the growth rate can be well approximated by the following argument. The reproduction rate Γ i repr,S of 0 400 800 0 0.2 0.4 0.6 0.8 1 Figure S3. Comparison of data obtained by simulations in the neutral case with fixed (red) and exponentially distributed (black) environmental changes. Both sets of data agree with our analytic calculations, where we used the mappings σ 2 = ω 2 τ for fixed times of environmental changes and σ 2 = ω 2 2τ for exponentially distributed switches. Thereby, the data confirms that the origin of the factor 2 in the mapping is solely the exponential distribution of the environmental changes. Parameters are φ 1 = φ 2 = 1 and ω 1 = ω 2 = 0.9, γ = 1, K = 5000, E = 0, Var[E] = 100 and τ = 0.01.
an individual, i, of type S at time t depends on the average of all instantaneous reproduction rates λ S (E) experienced by the individual : At the time of reproduction, we assume for simplicity that offspring looses its memory of past environments experienced by the progenitor. We consider again a time interval of length δt as in (S17). Fluctuations in the ratesΓ i S (t) decorrelate on timescales of the order of the lifetimes of individuals, t life , which are much longer than τ and δt. Therefore on the δt scale, noise is smooth, contrary to (S17). Conversely, timescales of several lifetimes are much smaller than those on which selection acts and much longer than the characteristic time of the noise. Therefore, to describe the dynamics of the fractions, the environmental noises are well approximated by a shortly correlated noise. An estimation of the amplitude of the noise is obtained by calculating the sum The durations ∆t of the environmental intervals are independent for different k and (and the contribution k = is negligible with respect to the rest of the sum) so that one can replace them by τ. In addition, the symmetry in the indices of the intervals allows us to further simplify the expression

4/7
To compute the average in Eq. (S22) three different orders of events have to be distinguished: a) If the birth of the j-th individual was prior to the one of the i-th individual, then it follows from Eq (S21) that the quantity to be averaged is where M j is the number of environmental switches since the birth of the i-th individual up to time t k . To derive Eq. (S22) we have used that the terms in the sum (S21) take independent values ±ω S with equal probability. Conversely, case (b) is when the birth of the j individual is posterior to the one of the i-th individual. Then the quantity to be averaged is ω 2 The integral over t is the continuous approximation of the sum over − k appearing in Eq. (S22) while the variables u and v refer to the variables M i k and δ b. The first and second term in the square parentheses of (S23) correspond to cases (a) and (b), respectively. By a series of change of variables and integrations by parts, it is shown that (S23) reduces to The validity of this approximation is confirmed for the neutral case in Fig. 2 in the main body of the paper where the fixation probability and time are compared to the analytic calculations employing Eq. (S24). Additional data for non-neutral evolution is presented in Fig. S4 where two species (one with finite variability, ω = 0.9, one with vanishing variability) are analyzed. Analytic solutions for the fixation time and probability are fitted to simulation data. The best fit deviates less than 1% from Eq. (S24). We now briefly discuss the dependence of σ S on the memory parameter m. In Fig. S5, the STD of the noise in the growth rate, σ ,depending on the memory parameter, m, is analyzed. Results were obtained by simulating the neutral evolution case, measuring the fixation time and calculating σ employing the analytic expression for the extinction time  To obtain the latter we fitted Eq. (S6) and the solution of Eq. (S12) to the IBM. We used A = s + σ 2 2 − σ 1 σ 2 ε andB = σ 2 1 + 2σ 1 σ 2 ε + σ 2 2 as fitting parameters and obtainedÃ = −3.04 × 10 −3 and B = 8.25 × 10 −3 . Other parameters are φ 1 = φ 2 = 1, ω 1 = 0.9, ω 2 = 0, τ = 1/100, γ = 1, K = 5000, E = 0 and Var[E].
Let us now analyze the mapping of the average reproduction rate ν S . Importantly, this mapping is very sensitive to model details which we will exemplify in the following. As results for neutral species do not dependent on the average reproduction rate, we have to turn to the evolution of non-equal individuals to understand the mapping of the growth rates. In Fig. S4, we show the fixation probability for two species with the same φ 1 = φ 2 = 1 but only the first species has a variable reproduction rate ω 1 = 0.9, ω 2 = 0 for m = 1. Red lines correspond to a fit with s = ν 1 − ν 2 ≈ −0.0030 and agree perfectly with simulation results. In other words, the first species does not only have a disadvantage due its sensitivity on environmental changes, σ 1 > σ 2 , but also has a smaller average growth rate. To study this effect in more detail, let us now analyze the fixation probability dependence on the memory parameter m, see Fig. S6 panel (a). Black dots correspond to the standard IBM (if not mentioned otherwise our discussion applies to this data), red dots to a slightly changed model which is going to be introduced in the following. For m = 0 both species are equally likely to fixate as the growth advantage of the more variable species ν 1 = φ + ω 2 1 τ > ν 2 = φ exactly compensates for its disadvantage due to the STD of the noise σ 1 > σ 2 . For increasing values of m first the more variable (m < 0.7) later the less variable species is favored (m > 0.7). Whether this behavior is caused by the STD of the noise or differences in the mean reproduction rates is not obvious as the influence of both fitness contributions is of comparable strength. Therefore, we estimate the selection coefficient, s = ν 1 − ν 2 , from the fixation probability data, see Fig. S6 panel (b). This is achieved by assuming that the variability of species 1 with ω 1 = 0 is zero (σ 1 = 0) and that the variability of species 2 with ω 2 = 8.237 is the same as in the neutral evolution scenario and thereby given by the data presented in Fig S5.  Note that this approximation might neglect some higher noise correlations arising due to the coupling of both species via the carrying capacity. For m = 0 the thereby obtained value of s agrees well with our analytic results, cf. Eq. (S20). With increasing m the growth rate of the more variable species is decreased till the selection coefficients becomes negative effectively favoring the more variable species. However the decrease of the selection coefficient with m is smaller than the reduction of σ shown in Fig. S5. Hence, for small m the more variable species is favored as its advantage due to a larger average reproduction rate is larger than its disadvantage due to its sensitivity on the environment. This advantage in the growth rate is more sensitive to details in the IBM in comparison to the variability discussed in the main text for the Langevin equation Eq. (1).We are going to illustrate this by analyzing a slightly modified version of the IBM. But before doing so, we present an intuitive argument explaining one factor influencing the average growth rate: When an individual is born it experiences the current environment shorter than the average length of an environment. However, the model weights all experienced environments equally, see Eq. (7) in the main text. As the first experienced environment is more likely to be a good environment [more reproduction events happen during more beneficial environments], higher growth rates have a larger weight in the average and the average growth rate of the variable species is effectively increased. To obtain a description including this factor, it would be best to perform a time average over all previously experienced environments. Unfortunately, such a procedure is computationally very expensive. We therefore, test our explanation for the bias by not including the very first environment, the one in which an individual is born, in the averaged reproduction rate. The average 25 different environments, the small modification of the model substantially changes the simulation results. While the modification almost has no impact on the STD of the noise, it alters the average reproduction rate. For instance the regime in which the more variable species is favored completely disappears, cf. red dots Fig. S6 where P fix ≤ 0.5 for all m. This example illustrates that on the first sight tiny details of an IBM might substantially influence the evolutionary outcome and that one should be cautious when drawing conclusions from them. Importantly, the mechanism discussed in the main text does not rely on specific assumptions of the microscopic models: A finite STD of the growth rate is always a disadvantage. It might be compensated for by a larger average reproduction rate but the same value of the growth rate without variability is always preferable.

DEPENDENCE ON THE SWITCHING RATE
In this Section, we present additional data for the non-neutral case. Fig. S7 shows the fixation probability depending on the environmental switching rate 1/τ. In particular, we investigate extinction for a species which is not sensitive to its environment (φ 1 = φ = 10, ω 1 = 0) competing with a sen-  Figure S7. Dependence of the fixation probability on the environmental switching rate. Both species have the same φ 1 = φ 2 = 10, but only the second species is sensitive to the environment (ω 1 = 0, ω 2 = 9). For no memory (m = 0), both species are equally likely to fixate, as the advantage in the average growth rate of species 2 exactly compensates for its disadvantage due to its sensitivity on the environment. For m > 0 those two effects do not cancel out anymore and the second species is favored. Other parameters are γ = 1, K = 100, α = 1, E = 0 and Var[E] = 100.
sitive species (φ 2 = φ = 10, ω 2 = 9) for different values of m. In the case of no memory m = 0 both species are equally likely to fixate as the advantage in the average reproduction rate ν 2 = φ + ω 2 τ exactly compensates for the disadvantage due to the STD of the noise in the growth rate [see Eq. (S20)].
For larger values of the memory parameter, a bias favoring the species with ω = 0 is present (the exact value of the fixation probability depends on mapping details as discussed above). Importantly, the bias is not only present for very quickly fluctuating environments, but already emerges if reproduction events happen on a time scale comparable to τ. This supports the conclusion, that we were already drawing in the body of this paper when discussing Fig. 4: the white noise approximation is an adequate description for such an evolutionary process in that parameter regime.