Noise-induced stabilization and fixation in fluctuating environment

The dynamics of a two-species community of N competing individuals are considered, with an emphasis on the role of environmental variations that affect coherently the fitness of entire populations. The chance of fixation of a mutant (or invading) population is calculated as a function of its mean relative fitness, the amplitude of fitness variations and their typical duration. We emphasize the distinction between the case of pairwise competition and the case of global competition; in the latter a noise-induced stabilization mechanism yields a higher chance of fixation for a single mutant. This distinction becomes dramatic in the weak selection regime, where the chance of fixation for a single deleterious mutant is an N-independent constant for global competition and decays like (ln N)−1 in the pairwise competition case. A Wentzel-Kramers-Brillouin (WKB) technique yields a general formula for the chance of fixation of a deleterious mutant in the strong selection regime. The possibility of long-term persistence of large [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr{O}}$$\end{document}O(N)] suboptimal (and extinction-prone) populations is discussed, as well as its relevance to stochastic tunneling between fitness peaks.

In the weak selection regime, | |  N s 1 and | |  s 1, the effect of selection is negligible and the neutral result reemerges. When the mutation is beneficial (s > 0, but still  s 1) and  Ns 1 (strong selection regime), Π n=1 is N-independent and converges, at large N, to s. A simple and intuitive argument for this result relies on the distinction between the region 1 ≤ n ≤ 1/s, where demographic fluctuations are dominant and the dynamic is more or less neutral, and the region 1/s < n, where selection dominates and fixation is almost assured. The chance of fixation is thus determined by the chance to cross a neutral region of length 1/s, which is exactly s 3 .
When the mutation is deleterious (s < 0), Π n=1 decays exponentially with N|s| in the strong selection regime, since now fixation is reached through a series of unlikely events against a constant bias. Accordingly, for any practical purpose one may neglect the chance of a deleterious mutation to reach fixation when N is large. This observation poses a serious question to the standard theory of species evolution. If genotypes of existing species are associated with local maxima in the fitness landscape, evolutionary pathways must cross fitness valleys. Because the chance of such "tunneling" events is vanishingly small, the timescales associated with it turn out to be unrealistically high 4 .
This set of results was obtained for a system with pure demographic noise, where the stochastic component in the reproductive success of each individual is independent of the success of its conspecific. As a result, the per-generation noise-induced abundance variations scale like the square root of the population size [i.e., are n ( )  ]. Environmental changes that affect coherently the fitness of entire populations lead to much stronger, Scientific REPORtS | (2018) 8:9726 | DOI: 10.1038/s41598-018-27982-1 n ( )  , abundance variations 5 , therefore one would expect a substantial impact of these fluctuations on the chance of fixation. Recently, many empirical studies showed that the effect of coherent fitness fluctuations is indeed much more pronounced than that of the demographic noise [6][7][8][9] , and that selection changes its amplitude and sign through time 10,11 . Consequently, the study of temporal environmental stochasticity received considerable attention [12][13][14][15][16][17][18][19][20][21] . In parallel, a few recent experimental studies have considered the response of microorganism communities (and their evolutionary pathways) to fitness fluctuations [22][23][24] .
In some scenarios fluctuating selection may activate a noise-induced stabilizing mechanism 25 . The main aim of this work is to calculate fixation probabilities in that case and to contrast them with the known results that were obtained in the absence of such a stabilizing mechanism. As we shall see, the fixation probabilities change dramatically under the influence of this mechanism. Moreover, in that case the system appears to allow for long-term persistence of suboptimal mutant populations, a phenomenon that may facilitate stochastic tunneling through fitness valleys (see discussion section).
In the next section we discuss the two models considered throughout this paper and clarify the conditions under which noise-induced stabilization occurs. In the third section we define mathematically these two models and explain the methodology used to obtain the chance of fixation when the environment fluctuates. The results are presented in the forth and the fifth sections and the possibility of stochastic tunneling is considered in the discussion. A glossary is provided in Table 1 and the main results of this paper are summarized in Table 2.

Noise-Induced Stabilization
To begin, let us define two zero-sum competition models, one that does not allow for noise-induced stabilization when the environment fluctuates (model A) and one that admits this phenomenon (model B).
In model A competition is pairwise and selection acts linearly. As an example one may envisage a population of competing animals, where a random encounter between two of them may end up in a struggle over, say, a piece of food, a mate or a territory. To model such a system we assume that these "duels" occur at a constant rate between two randomly picked individuals and in each duel the loser dies and the winner produces a single offspring. If x = n/N is the population fraction of the mutant, the chance of an interspecific "duel" is 2x(1 − x) and the chance of the mutant individual winning such a duel is defined to be 1/2 + s/4. Accordingly, the deterministic growth/ decay of x (when time is measured in generations, N elementary duels in each generation) satisfies the logistic equation,  (11) and (12) Eq. (23) If s fluctuates in time such that its mean value is zero (a time-averaged neutral model 16 ) the system performs an unbiased random walk along the z = ln[x/(1 − x)] axis. When the mean of s, s 0 , is nonzero, the random walk is biased towards either fixation or extinction.
In model B, on the other hand, the competition is global. In a forest, for example, following the death of an adult tree local seeds or seedlings are competing for the opened gap. If the seed dispersal length is larger than the linear size of the forest, the seed bank at any given location reflects the composition of the whole community. Death events are assumed to be random and fitness-independent. Accordingly, the chance of a mutant species with relative log-fitness s to gain one individual in an elementary death-birth event is a multiplication of the chance that a wild-type was chosen to die, 1 − x, by the chance of the mutant to win the gap, xe s /(1 − x + xe s ) (a Moran process). The probability of the mutant species to lose one individual is, accordingly, (1 − x)x/(1 − x + xe s ). In contrast with model A, here the fitness dependence is nonlinear. As a result, the deterministic dynamic satisfies, (2), a flow towards zero or one, the s 2 term has an attractive fixed point at x = 1/2. When s is fixed in time this second term is negligible, but when s fluctuates the s 2 term may dominate. Therefore, in model B environmental variations may induce stability through nonlinear fitness dependence 25 .
While for pairwise competition (model A) the effect of environmental fluctuations weakened when their correlation time decreases, the stabilizing effect of the global competition model reaches its maximum when the correlation time is minimal 18 .
The ability of environmental fluctuations to stabilize a coexistence state was pointed out by Chesson and coworkers 25,26 and is known in the ecological literature as the storage effect. The storage effect stabilizes a coexistence state when the fitness affects recruitment but death occurs at random. This is the situation in our model B, which is an individual based version of the lottery game considered by Chesson and Warner 25 , see a detailed discussion in 17 . On the other hand, in model A fitness affects both birth and death in an anticorrelated manner. As a result there is no storage effect stabilization in that case. Models A and B are thus the two extreme scenarios; in general, as long as the effect of fitness on recruitment is larger than its effect on death, one should expect the stabilizing mechanism to affect the system.

Model Definitions and the Backward Kolmogorov Equation
For both model A and model B we assume that where s 0 is the mean log-fitness and ζ(t) may take two values, ±γ (telegraphic, or dichotomous, noise). Both white Gaussian noise and white Poisson noise can be recovered from the dichotomous noise by taking suitable limits 27 , so the results obtained here are quite generic. The chance of the environment to switch (from ±γ to γ) is 1/(δN) per elementary birth-death event, so the sojourn time of the environment is taken from a geometric distribution with mean δN elementary birth-death steps, or δ generations. The chance of fixation when the mutant is in the plus (minus) state and its abundance is n, Π +,n (Π −,n ), satisfies the discrete backward Kolmogorov equation (BKE), where W ± s are the transition probabilities in the +γ state and the Q ± s are the corresponding probabilities in the (−γ) states. The transition probabilities of model A are (we replaced n/N by x),  The exact difference equation (4), with the appropriate set of W's and Q's and with the boundary conditions Π +,0 = Π −,0 = 0, Π +,N = Π −,N = 1, may be solved numerically as a Markov chain 28 or as a linear system 18 , and these solutions are compared below with the analytic formulas.
One can translate this pair of discrete BKEs for Π ±,n into an equivalent set for Π n ≡ (Π +,n + Π −,n )/2 and Δ n ≡ (Π +,n − Π −,n )/2. Taking the continuum limit where n is replaced by Nx and functions of x ± 1/N are expanded to second order in 1/N, a pair of coupled, second order differential equations for Π(x) and Δ(x) emerges. In 18 we have analyzed these equations in the limit of large N and small s 0 . Using a dominant balance argument we showed that the dynamic is governed by a single second-order equation (in 18 only the time to absorption was discussed, but the equation for Π is the homogenous version of the corresponding BKE, see 29 ). This equation is, 0 with the boundary conditions Π(0) = 0 and Π(1) = 1.
is the strength of the environmental noise and ≡ G Ng is this strength divided by the strength of demographic stochasticity, 1/N. The differences between model A and B are encapsulated in the parameter η: As δ grows, model B becomes closer to model A. However, the derivation of Eq. (7) assumes that fixation cannot occur during a single sweep of the environment, so an increase in δ is legal only if N increases such that δ γ Eq. (7) may be solved using integrating factors, but this leads to complicated and hard to interpret nested integral expressions. Instead one may analyze this equation in the inner 1) regimes and then match asymptotically the solutions in the large N (more precisely, large G) limit. In the next section we present briefly the results for model A, following 30 . Our purpose is to contrast these result with the outcomes of model B and to emphasize the effects of the noise-induced stabilizing mechanism.

Model A: Local Competition and Linear Selection
The solutions of model A in the inner, middle and outer regimes are given by, where α ≡ s g / 0 and the constants, that were determined using asymptotic matching, are, For |α| < 1 one may use the uniform approximation solution for an arbitrary x, 1, the uniform approximation takes the form, The agreement between Π unif and the outcomes of the numerical solutions of Eq. (4) is demonstrated in Fig. 1. The theory and the numerics become closer and closer as N increases.
The chance of a single mutant to reach fixation is obtained by plugging x = 1/N into the inner solution, The most important conclusion from the comparison between Eqs (9 and 10) and (1) is the modification of the criteria for strong selection. We define the strong selection sector as the parameter regime where the chance of fixation of a single beneficial mutant becomes N independent. This happens when  N n c where n c marked the point where the deterministic effect of selection dominates against the stochastic effects of fluctuations. While for system with selection and pure demographic noise n c = 1/s, here the criteria for strong selection is C 1 = C 3 = 1, i.e., c 0 This scale diverges exponentially when the mean selection is much smaller than the effective strength of fitness fluctuations, which might be the generic situation in living systems. Accordingly, under environmental stochasticity systems may be in the weak selection regime even if N is very large.
Two other characteristic scales in this system are N 1 and n 2 . For  N N 1 , C 2 = 0, so Ñ gs g exp( / )/ 1 0 . The chance of fixation of the mutant population becomes large when it reaches n 2 such that Π in (n 2 /N) > 1 − e −1 , thus,  This scale has been identified in 15 . Clearly, all three scales have similar features as they grow exponentially when s 0 is much smaller than g. For a system with pure demographic noise and selection (Eq. 1), the chance of fixation becomes N independent above n c = 1/s 0 and the condition for n 2 yields 1/s 0 as well. While n 2 of (15) converges to this limit when g → 0, n c does not, and this reflect the inadequacy of our expression for C 1 in the γ  s 0 limit 30 .
In the weak selection regime, N < n c , the chance of fixation is N-dependent. When α  G ln 1 one can expand the inner solution in small α to obtain 15,30 , While α  G ln 1 is small, ln G may be large, so in the weak selection regime, and in particular in the time-averaged neutral scenario where s 0 = 0, the chance of fixation decays logarithmically with system's size as demonstrated in Fig. 2.
When selection is weak an increase in g increases the chance of fixation. In the purely demographic neutral case Π is determined by abundance so Π n=1 = 1/N. In the limit of infinitely strong environmental variations a mutant will reach fixation for certainty if it was born in the right time, so the chance of fixation will grow to one half. In general the transition from abundance-dependence to environment dependence facilitates the chance of low-abundance populations to win 15 .
The situation is completely different in the strong selection regime 30 , where Π n=1 is a monotonously decreasing function of g. Here the reason is the divergence of n c when g increases, meaning that the chance of the beneficial mutant population to enter the region of deterministic selective growth is much smaller.

Model B: Noise Induced Stabilization
In model B the dynamic is affected by the noise-induced stabilizing mechanism that facilitates the establishment of a mutant. Before we introduce the expressions for the chance of fixation, we would like to discuss the conditions under which this stabilizing mechanism takes place.
When s(t) = s 0 ± γ, as described above, and the environmental fluctuations are rapid, Eq. (3) for the deterministic dynamic of the population takes the form Assuming, γ | |  s 0 , this equation supports an attractive fixed point at Therefore,  s , the ratio between mean selection and environmental fluctuations, determines the qualitative behavior of the system. When | | <  s 1 the noise induces a stable coexistence point and the dynamic of model B differs substantially from the dynamic of model A. When | | >  s 1 the deterministic force does not change its sign in the region between fixation and extinction, so the behaviors of model A and model B are qualitatively similar. In agreement with this observation, in 18   where as before α ≡ s g / 0 and the constants are given by, The main difference between Eqs (20) and (21) and their model A counterparts, Eqs (9) and (10), is the different scaling of C 2 . In model B, C 2 goes to zero when δ  N g exp( )/ . Above this s 0 -independent point, the chance of fixation in the middle region is fixed, C 3 = C 1 , as demonstrated in Fig. 3. A wide plateau emerges in the middle region due to the force towards x * . This force is strong, hence Π becomes almost x-independent.
The uniform solution when C 2 → 0 has a relatively simple form (see Fig. 3), Amazingly, Π n=1 (for a beneficial mutant) turns out to be an increasing function of N, a behaviour that manifests itself in Fig. 4. This phenomenon reflects the stabilizing effect of the nonlinear mechanism: the chance to reach the plateau does not depend on N because the plateau occurs at values of x that scale like 1/N. For example, in the large N limit Π = x ( ) s 0 0 sticks to 1/2 for any x(1 − x) > 2/(Nγ 2 ). For =  G Ng 1 the chance of fixation decays like 1/N since this regime (in which our expressions fail) is dominated by demographic stochasticity. When N increases the stabilizing mechanism wins against the demographic noise and leads to an increase of the chance of fixation.
Π n=1 in Eq. (24) is a multiplication of two factors: its numerator is the chance of establishment Π est , which is the probability that the mutant population will reach the basin of attraction of the coexistence fixed point (the plateau). C 1 , that determines the denominator, is the chance that the mutant population will reach fixation given establishment.
If the mutant is advantageous (s 0 > 0) and the system is in its strong selection regime (G −2α → 0 or  N n c ), C 1 = 1 and, This is a monotonously decreasing function of δ. When δ increases, the stochasticity becomes stronger and the stabilizing mechanism weakens 17 , both effects tend to decrease the chance of fixation. The dependence on the amplitude of environmental fluctuations, γ, is more complicated. When γ is small, its increase facilitates the stabilizing mechanism that increases the chance of fixation, while for large γ the increase in n c is the dominant effect and Π n=1 decreases. Model A and model B differ even more dramatically in the weak selection regime  N n c , where G −α ≈ 1 − α ln G. For model B, the chance of fixation becomes, Unlike model A, where the chance of fixation decays logarithmically with N in the weak selection regime, for model B in the same case the chance of fixation is N-independent. As long as δ s Ng ln( )/ 1 the chance of fixation is simply one half of the chance of establishment: the effective strength of the selection bias,  s , in that case is zero, so once the mutant population reaches the plateau its odds to win or to lose are equal. For nonzero  s there is a linear increase or decrease of Π n=1 as a function of ln N. This relatively weak effect is demonstrated in panel (B) of Fig. 4.
Before concluding this section we would like to add a technical comment about our numerical calculations. Π(x) is obtained from Eq. (4) by solving a linear problem (dividing a matrix by a vector). Using the sparsity of the matrix, in model A we were able to analyze systems with up to N = 10 6 individuals. Because of the plateau that characterizes model B in the strong selection regime, this numerical solution becomes difficult; the plateau indicates that the matrix to be inverted is almost singular. To overcome this problem we have used quadrapole precision algorithm and this makes the numerics much slower and limits available system sizes to N values up to 20,000.

The Chance of Fixation For A Deleterious Mutant Under Strong Selection -A WKB Approach
Until now we discussed the weak selection regime for both beneficial (s 0 > 0) and deleterious (s 0 < 0) mutants, but the strong selection regime (  N n c ) for deleterious mutant has not yet been considered. In this section we would like to provide a few basic insights for that case.
Quantitatively, one may guess that the chance of fixation in this regime behaves differently when γ < |s 0 | and γ > |s 0 |. In the latter case, the growth rate of a deleterious mutant is still positive (γ − |s 0 |) during half of the time, so the most probable (yet rare) route to fixation is based on picking a sequence of good years. During this series of lucky events the (on average) deleterious mutant plays the role of a beneficial one, and its time to fixation scales like ln N. Therefore, the chance of fixation, namely the chance to pick such a lucky sequence (which shrinks exponentially with the length of the sequence), decays like a power-law in N.
On the other hand, when γ < |s 0 |, the route to fixation involves an improbable series of successes in consecutive elementary duels (reflecting demographic stochasticity) and in such a case Π n=1 decays exponentially in N, like in the purely demographic case (1).
Looking at Eqs (13) and (24), one notices that the decay of Π n=1 when s 0 < 0 (and α < 0) is due to the divergence of the This difficulty turned out to be related to the failure of the continuum approximation that has been used when we have translated the difference equation (4) to the differential equation (7). As explained in 32 , this procedure fails when the differences between neighboring points (say, Π n+1 − Π n ) are too large and cannot be approximated using first and second derivatives only. To overcome this obstacle a WKB approach was suggested 32 ; here we would like to implement it for a model with environmental stochasticity. We shall neglect, for the moment, the effect of demographic noise and assume that extinction and fixation happen when the abundance reaches 1/N and 1 − 1/N, correspondingly.
As explained in the introduction [following Eq. (2)], in the absence of demographic noise and under model A dynamic =  z sz, where ≡ − z x x /(1 ). Accordingly, during δ generations the dynamic of z satisfies, In the opposite, strong selection phase, the qualitative features of the chance of fixation Π n=1 are not much different between model A and model B. In both cases the chance of fixation for a beneficial mutant is N-independent (Table 2, sixth line) while the chance of fixation of a deleterious mutant falls like N −2f ′ when γ > |s 0 | (seventh line) and exponentially with N if γ < |s 0 |. In model A the chance of fixation in the strong selection regime decreases as the environmental noise becomes stronger, while in model B it decreases with the correlation time δ but increases with noise amplitude γ. (Through this discussion, when the features of model B are contrasted with those of model A, model B is assumed to support a noise induced attractive fixed point, i.e., | | <  s 1. Otherwise, the behavior of model B dynamic is qualitatively the same as model A).
On the other hand, in the weak selection regime ( Table 2, fifth line) there are substential differences between the two scenarios. As required by its name, in this regime selection is a second order effect and the fate of the mutant population is determined by stochasticity. In a stochastic and balanced game, like the classical gambler's ruin problem, the chance to win is inversely proportional to the effective size of the community, so under purely demographic noise it is 1/N and under model A dynamic it is 1/ln N.
In sharp contrast with this result, in model B the system supports an attractive fixed point at x * . The plateau that characterizes Π(x) in that case (Fig. 3) reflects the effect of this attractive fixed point, marking the range of x values which lies in its basin of attraction. The attractiveness of this coexistence point grows with N and leads to an apparently paradoxical behavior: an increase of the chance of fixation with N. Once this fixed point becomes dominant, the fate of the mutant population depends on its chance of establishment Π est , i.e., of reaching the plateau, and on C 1 , the probability to jump from the plateau region to fixation. Since the plateau is wide, these two probabilities are N independent when the strength of selection s 0 is relatively weak, and so is the chance of fixation itself. Even if N is huge, as long as it is much smaller than n c , model B yields an N-independent value for Π n=1 for both beneficial and deleterious mutants.
Even in the strong selection regime, where the chance of fixation of a deleterious mutant are vanishingly small, model B dynamic still supports an attractive fixed point at x * as long as | | <  s 1. The chance of establishment is still N-independent; it is the chance of fixation conditioned on establishment, C 1 , which goes to zero in that case. As shown in 18 , the lifetime of the mutant population, once established, is = generations, a huge time for large communities. The stabilizing mechanism of model B thus allows for the long-term persistence of a macroscopic,  N ( ), extinction-prone population with negative fitness. This phenomenon may provide a mechanistic explanation to one of the the mysteries of evolutionary dynamics: the ability of evolutionary pathways to cross fitness valleys, i.e., to sustain a chain of suboptimal intermediate forms that bridges between two fitness peaks in a rugged fitness landscape. This stochastic tunneling has been recognized a while ago as a major theoretical problem, since the chances for a tunneling event are vanishingly small 4 . To overcome this problem modern theoretical studies consider evolution on a neutral or nearly neutral (holey, high-dimensional, connected) fitness landscape. In this neutral picture a separate mechanism has to be invoked to explain speciation, while on rugged landscape each species corresponds to a separate fitness pick and disruptive selection is guaranteed by the landscape itself. The long term existence of macroscopic suboptimal populations, considered here in context of model B, may allow for such a tunneling to occur with relatively high probability through a chain of mutation as long as the depth of the fitness valley is smaller than γ 2 /2, while keeping the intermediate forms extinction-prone. The relevance of the mechanisms suggested here to the development of natural communities depends on the amplitude of fitness variations with respect to their time-averaged value, on the typical correlation time of these fluctuations and on the range of competition -whether it is local/pairwise (model A) or global (model B). It is quite difficult to quantify s 0 and γ from field data, and in experimental systems the external conditions are usually kept fixed, as opposed to the intrinsic variability of natural environments. Still, we believe that the theory presented here, when applied to some experiments and to field data in population genetics and community ecology, may suggest many new insights into the processes that govern the composition and the evolution of natural communities.