Bet-hedging strategies in expanding populations

In ecology, species can mitigate their extinction risks in uncertain environments by diversifying individual phenotypes. This observation is quantified by the theory of bet-hedging, which provides a reason for the degree of phenotypic diversity observed even in clonal populations. Bet-hedging in well-mixed populations is rather well understood. However, many species underwent range expansions during their evolutionary history, and the importance of phenotypic diversity in such scenarios still needs to be understood. In this paper, we develop a theory of bet-hedging for populations colonizing new, unknown environments that fluctuate either in space or time. In this case, we find that bet-hedging is a more favorable strategy than in well-mixed populations. For slow rates of variation, temporal and spatial fluctuations lead to different outcomes. In spatially fluctuating environments, bet-hedging is favored compared to temporally fluctuating environments. In the limit of frequent environmental variation, no opportunity for bet-hedging exists, regardless of the nature of the environmental fluctuations. For the same model, bet-hedging is never an advantageous strategy in the well-mixed case, supporting the view that range expansions strongly promote diversification. These conclusions are robust against stochasticity induced by finite population sizes. Our findings shed light on the importance of phenotypic heterogeneity in range expansions, paving the way to novel approaches to understand how biodiversity emerges and is maintained.


Introduction
The dynamics and evolutionary history of many biological species, from bacteria to humans, are characterized by invasions and expansions into new territory. The effectiveness of such expansions is crucial in determining the ecological range and therefore the success of a species. A large body of observational [1,2] and experimental [3][4][5][6] literature indicates that evolution and selection of species undergoing range expansions can be dramatically different from that of other species resident in a fixed habitat. Theoretical studies of range expansions based on the Fisher-Kolmogorov equation [7,8] or variants [9][10][11] also support this idea. Adaptive dispersal strategies [2] and small population sizes at the edges of expanding fronts [12,13] are among the main reasons for such difference.
Range expansions often occur in non-homogeneous and fluctuating environments. Under such conditions, it is possible to mathematically predict the expansion velocity of a community of phenotypically identical individuals [14][15][16][17][18][19]. However, diversity among individuals is expected to play an important positive role when populations expand in fluctuating environments. For instance, diverse behavioral strategies help animal populations to overcome different invasion stages and conditions [20][21][22][23]. Analyses of phenotypic diversity in motile cells suggest that it also may lead to a selective advantage at a population level [24][25][26]. Although several studies have tackled the problem of how individual variability affects population expansion [6,9,10,[27][28][29][30][31], systematic and predictive theory is still lacking [23].
In this paper, we study how bet-hedging strategies can aid populations in invading new territories characterized by fluctuating environments. In particular, we analyze the effect of spatial expansion, different types of environmental heterogeneity, and demographic stochasticity on development of bet-hedging strategies for a population front evolving according to a Fisher wave.
By employing mathematical as well as extensive computational analyses, we find that the advantage of bet-hedging in range expansions depends on the rate of environmental variation.
In particular, bet-hedging is more convenient for infrequently varying environments, whereas its advantages vanish for frequent environmental variation. For the same model, bet-hedging is never an advantageous strategy in the well-mixed case, supporting the view that range expansions strongly promote diversification. We further find that spatial environmental variations provide more opportunities for bet-hedging than temporal fluctuations. Finally, we show that our conclusions still hold when considering stochastic effects on the front propagation induced by a finite population size. The paper is organized as follows. We introduce a general population model and an example with two available phenotypes and two environmental states. We present an extensive study of this example. We then generalize the main conclusions obtained for the example for a case with an arbitrary number of environmental states, and then with also an arbitrary number of strategies. We conclude with a discussion and future perspectives. i = 1. . .N of the population assuming each phenotype i with ∑ i α i = 1 and 0 � α i � 1 8i ( Fig  1A). As customary in game theory, we say that a strategy is a "pure strategy" if α i = δ ik for some phenotype k, and a "mixed strategy" otherwise. We assume that the α i 's remain constant in time within the population.
The environment can be found in one of M different states, which can randomly alternate either in time or in space. We call p i the probability of encountering environment i. We further define the growth rate s ij � 0 of phenotype j in environment i (Fig 1A). When the population size is sufficiently large, so that demographic stochasticity can be neglected, the populationaveraged growth rate given the state i = i(x, t) of the environment at position x and time t is Since Eq (1) is linear in the α j 's, the population-averaged growth rate in a given environment is always maximized by the pure strategy with the highest growth rate. However, in the presence of uncertainty about the environment, the population might choose other strategies. One possibility is to select a different pure strategy, less risky than the optimal one. This case is often termed "conservative bet-hedging" in the ecological literature [41]. Another option is to adopt a mixed strategy, with different phenotypes more adapted to different environments. This case is termed "diversifying bet-hedging" in the literature [41,56]. Since our interest is in diversification, the term "bet-hedging" will refer herein to diversifying bet-hedging.
Before presenting our results in full generality, weconsider a simple, yet ecologically relevant instance of the model with only two phenotypes: "safe" and "risky" and two environmental states: "adverse" (a) and "favorable" (b). The safe phenotype is characterized by a growth rate s s both in the adverse and favorable environments. The growth rate of the risky phenotype is s a in environment (a) and s b in environment (b) (Fig 1B) [57]. The two environments occur with the same probability, p a = p b = 1/2. A fraction of individuals α adopts the risky phenotype and the complementary fraction (1 − α) adopts the safe phenotype ( Fig 1B). For this model, the population-averaged growth rate reads sðx; tÞ ¼ s a ¼ ð1 À aÞs s þ as a ; in env: a Note that, with a slight abuse of notation, we use equivalently σ i or σ(x, t) to denote the population-averaged growth rate in the environment i(x, t). For pure strategies, α = 0 or α = 1, the population-averaged growth rate σ reduces to the growth rate of the safe or risky phenotype, respectively.

Two-phenotype, two-environment model
We seek to understand those conditions under which bet-hedging is advantageous for the population. To this end, we shall compare three situations: i) well-mixed populations, ii) range expansions in environments that fluctuate temporally, but that are homogeneous in space ( Fig  1C), and iii) range expansions in spatially fluctuating environments that are homogeneous in time ( Fig 1D).
Well-mixed case. We start by analyzing the well-mixed case, where the spatial coordinates of individuals can be ignored. The total population density f(t) evolves according to the The fitness of an individual with phenotype j in an environment i is given by s ij . B) Two-phenotypes model: Individuals can adopt either a "risky" or a "safe" phenotype with probabilities α, and 1 − α respectively. The safe phenotype is characterized by an environment-independent growth rate s s . The growth rate of the risky phenotype is s a or s b , depending on whether the current environment is "adverse" (a) or "favorable" (b). C) and D) Sketch of range expansion in a population having 0 � α � 1 for temporally varying C) and spatially varying D) environments, respectively.
In writing Eq (3), we used the assumption that the fraction α of the population adopting the risky phenotype remains constant in time (see [58,59] for cases in which this assumption is relaxed). Eq (3) can be readily integrated, obtaining where hσ i i = ∑ i p i σ i denotes an average over the environmental states. For Eq (4) to hold, we do not need to make strong assumptions about the statistics of the environmental states, other than it should be stationary, ergodic, and with a finite correlation time.
The optimal strategy α � is obtained by maximizing the right-hand side of Eq (4) respect to the strategy α. Since hσ i i is a linear function of α, its maximum is always reached at the extremes of the interval (α 2 [0, 1]). In particular, defining the normalized growth rates s a � s a =s s ands b � s b =s s , we find that the optimal strategy is α � = 1 whens b > 2 Às a and α � = 0 otherwise. This means that no bet-hedging strategy is possible in this model in the wellmixed case [57].
This simple result illustrates an aspect of bet-hedging that is sometimes under-appreciated. In well-mixed systems, bet-hedging optimal strategies appear when the model includes at least one of the following ingredients: a) discrete generations, as in the seminal work of Kelly [42], b) finite switching rates among strategies [33,59], or c) a delta-correlated environment [53]. Any of these ingredients can lead to nonlinearities in the average exponential growth rate, therefore opening the way for a non-trivial optimal strategy.
Note that, in this model, the frequency of environmental change does not play a role, as far as it is finite [53]. The physical reason can be understood from the right-hand side of Eq (4): the optimal strategy depends on the frequency of different environmental states but not on the switching rates. This feature is also shared by other well-mixed models that do allow for optimal bet-hedging strategies, such as the classic model by Kelly [42]. We shall see in the following that, on the contrary, the rate of environmental change plays an important role for expanding populations.
Range expansion in fluctuating environments. We now consider a population expanding into an unoccupied, one-dimensional space under the influence of a stochastically changing environment. Its population dynamics are described by the Fisher equation [7,60]: where f(x, t) is the population density at spatial coordinate x and time t, and D is the diffusion constant, which characterizes the motility of individuals. For a constant growth rate σ, the stationary solution of Eq (5) is characterized by a front advancing in space with velocity v F ¼ 2 ffi ffi ffi ffi ffi ffiffi Ds p . Instead, we consider a fluctuating case in which the growth rate σ(x, t) depends on the population strategy α and on environmental conditions according to Eq (2). In such case, we define an asymptotic mean velocity of the front as v M ¼ lim In what follows, we take v M as a proxy of the long-term population fitness and maximize it with respect to α to determine the optimal strategy.

Range expansion in temporally varying environments.
We first consider the case in which environmental conditions change randomly with time, but are homogeneous across space, σ(x, t) = σ(t) (see Fig 1C). Switching rates between adverse and favorable environments are k a ! b = k b ! a = k. We first estimate the asymptotic mean velocity defined in Eq (6) in the limiting cases of k ! 0 and k ! 1.
When the environment changes very infrequently, k ! 0, the population front has the time to relax to the asymptotic shape characterized by its corresponding Fisher velocity, v a ¼ 2 ffi ffi ffi ffi ffi ffi ffi ffi Ds a p or v b ¼ 2 ffi ffi ffi ffi ffi ffi ffi ffi Ds b p depending on the environment [7,61]. Thus, the asymptotic mean velocity can be estimated as v M = (v a + v b )/2. Maximizing v M with respect to α, we find that in this case, a bet-hedging optimal strategy exists under the conditions (Fig 2A): In the opposite limiting case of a rapidly fluctuating environment, k ! 1, the population effectively experiences the average of the two growth rates, so that the velocity can be estimated as v M � 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Dhsi p , where h. . .i denotes an average over the environmental states. In this case, the optimal strategy α � is achieved by maximizing the average growth rate hσi with respect to α. Since hσi is linear in α, the maximum always lies at the extremes of the interval [0, 1]. In particular, we find α � = 1 whens b > 2 Às a and α � = 0 otherwise, as in the well-mixed case. This implies that no bet-hedging regime exists in this limit, similarly to the well-mixed case ( Fig  2B).
To explore the intermediate regimes of finite k, it is necessary to resort to numerical simulations of Eq (5). For a set of parameters such that the optimal strategy is α � = 1 for k ! 0, the optimal strategy remains α � = 1 for all values of k, see Fig 3A. Instead, in a case where the optimal solution is in the bet-hedging region for k ! 0, the optimal strategy α � increases with the switching rate, so that for large k the optimal strategy is outside the bet-hedging region, α � = 1, see Fig 3B. These results support our analytical estimates of limiting values and suggest that the asymptotic mean velocity is a monotonically increasing function of the switching rate k in Bet-hedging strategies in expanding populations this case. Note that, in the example of Fig 3B, the velocity corresponding to the optimal bethedging strategy is only a few percent larger than the velocity for α = 0. For other parameters values, we found velocities up to 15% larger than for pure strategies.
Range expansion in spatially varying environments. We now consider the case in which environmental conditions are constant in time, but depend on the spatial coordinate x. The dynamics are described by the Fisher Eq (5) with two types of environment randomly alternating in space, σ(x, t) = σ(x). We call k S the spatial rate of environmental switch, so that the probability of encountering an environmental shift within an infinitesimal spatial interval dx is equal to k S dx. The switching rates from environment a to b and vice-versa are both equal to k S . As above, we first analyze the two limits k S ! 0 and k S ! 1.
In the limit k S ! 0, the population front traverses large regions of space characterized by a constant environment, either a or b, thus being able to reach the corresponding Fisher velocity, v a or v b , respectively. The mean traversed lengths Δx a and Δx b are equal for the two environments. On the other hand, the mean times spent in each of them, t a and t b , are different, and satisfy the relation Therefore, in this case, the asymptotic mean velocity is given by the harmonic mean of the velocities in the two environments Here, for k S ! 0 the bet-hedging region is broader with respect to the temporally fluctuating environment for k ! 0, see Fig 4A. At the opposite limit of large k S , the environment is characterized by frequent spatial variations. In this case, the population front occupies multiple a and b sectors with an effective Bet-hedging strategies in expanding populations growth rate hσi. As in the time-varying case, the asymptotic mean velocity in this limit is v M ðk S ! 1Þ ¼ 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Dhsi p , see also [15,16]. The corresponding optimal strategy is the same as in Fig 2C, so that there is no bet-hedging regime.
We numerically solved Eq (5) for intermediate values of k S and obtained the mean asymptotic velocities as a function of α, see Fig 4B. Results support theoretical predictions in the limiting cases k S ! 0 and k S ! 1. In this case, we observe a non-monotonic behavior of v M as a function of k S , so that the maximum mean velocity is attained at an intermediate switching rate. An analytical explanation of this non-trivial effect goes beyond the scope of this work.
Effect of finite population size. The deterministic Fisher Eq (5) is rigorously valid only in the limit of infinite local population sizes. We now explore the robustness of our results when considering stochasticity induced by the finite size of populations, i.e. "demographic noise". We focus on the case of a front propagating in a temporally varying environment. To study finite population size, we solve numerically a stochastic counterpart of the Fisher equation ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi see e.g. [62]. In Eq (10), ξ(x, t) is Gaussian white noise with hξ(x, . The parameter N represents the number of individuals per unit length corresponding to f(x, t) = 1. For large population sizes, N � 1, Eq (10) reduces to Eq (5). Numerical integration of Eq (10) requires some care due to the fact that both noise and the deterministic terms go to zero as the absorbing states f(x, t) = 0 and f(x, t) = 1 are approached [63][64][65]. A detailed description of our integration scheme is presented in the Supporting S1 Appendix. , respectively. The red curve is the analytical solution for a spatially fluctuating environment with k S ! 0, see Eq (9). Note that in this case, the asymptotic mean velocity does not increase monotonically with k S but is maximal at k S � 0.
where C is a constant, N is the maximum population size per unit length, and v F ¼ 2 ffi ffi ffi ffi ffi ffiffi Ds p is the Fisher velocity in the absence of demographic noise. Eq (11) is valid in the weak noise limit; for the corresponding strong noise expression, see [68]. Asymptotic mean velocities for stochastic waves in temporally varying environments are shown in Fig 5. Also in this case, small populations, subject to relatively strong demographic noise, propagate more slowly than large populations. In particular, curves at different values of N can be approximately rescaled using Eq (11), assuming that C does not depend on α (insets of Fig 5). These results imply that the optimal strategy α � is robust with respect to demographic noise, at least for moderately to relatively large values of N. The same scaling holds for spatially varying environments, but with mild deviations that seem to expand the bet-hedging region even further, compared with the infinite population size limit (see Supporting S2 Appendix). Finally, we remark that the effect of finite population size on well-mixed bet-hedging populations has been studied in the literature [33, 57-59, 69, 70].

Two-phenotype, multiple-environment model
In this section, we generalize our results to a model with two strategies, but an arbitrary number i = 1. . .N of environmental states. Let us start with the temporally varying case. Following the usual logic, the mean velocity for k ! 0 reads ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where s i1 and s i2 are the growth rates of the two strategies in environment i. The first derivative of the mean velocity respect to α reads @v m @a ¼ ffi ffi ffi ffi D p X i p i s i1 À s i2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi as i1 þ ð1 À aÞs i2 p ð13Þ Since v(α) is a concave function, the condition for having a bet-hedging strategy, i.e. a maximum in the interior of the interval (0, 1) is These conditions reduce to the Eq (7) in the limiting case of the two-environment model. With a similar strategy we can compute the limits of the bet-hedging region also for the spatially varying case. In this case we have To determine the bet-hedging region we follow the same logic as in the temporally varying case, yielding so that the condition in this case reads Even in this case, the bet-hedging region is broader in the spatially-fluctuating than in the temporally-fluctuating case. This fact is proven in full generality in the next subsection.

General bet-hedging model
In this Section, we demonstrate that our main conclusions hold in full generality for arbitrary numbers of phenotypes N and environmental states M (see Section Model). In particular, for a temporally fluctuating environment in the limit of very slow switching rates, the bet-hedging regime occupies a reduced region of parameter space compared to temporally constant environments fluctuating slowly in space. Also in this case, we find that for frequent environmental change, the propagation velocity tends to v M � 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Dhsi p , regardless of whether the environmental fluctuations depend on time or space. Therefore, the optimal strategy maximizes the linear function of the α i s hσi and is therefore a pure strategy as discussed after Eq (1).
We consider a range expansion where the environment fluctuates in time and the stochastic switching rates among the M environmental states are small. Following the same line of thought of the two-strategy, two-environment model, the optimal strategy maximizes where σ i = ∑ j s ij α j . For spatially varying environments, the optimal strategy maximizes the harmonic mean Both for Eqs (19) and (20), maximization has to be performed with the constraint ∑ j α j = 1 and 0 � α j � 1 8j. We recall that the bet-hedging regime is the region of parameter space where the optimal solution is a mixture of all phenotypes, α i > 0 8i. Here we show that if, for a given choice of the s ij 's and p i 's, a population advancing in a temporally varying environment is in a bet-hedging regime, then the same holds for spatially varying environments. For the demonstration, we borrow a mathematical tool from evolutionary game theory [71]. We introduce the gradients where hxi = ∑ i p i x i is the average over environments. We now associate replicator equations to Eqs (19) and (20): The system is in a bet-hedging regime when the replicator equations admit a stable fixed point in the interior of the unit simplex, 0 < α i < 1. Instead of computing the fixed point explicitly, we check whether each phenotype l has a positive growth rate for α l � 1. Brouwer's fixed point theorem ensures that, under this condition, there must be a fixed point in the interior (see [71], chapter 13). For our aims, it is therefore sufficient to prove that, for small α l , if ðF T l À � F T Þ is positive, then ðF S l À � F S Þ must be positive as well. Note that for α l � 1, the average σ = ∑ j s ij α j does not depend on α l , and therefore, σ and s l are uncorrelated random variables respect to the average over the environment. Since ffi ffi ffi s p > 0, this means that the sign of ðF T l À � F T Þ is the same than the quantity Following the same logic, the sign of ðF S l À � F S Þ is the same than This means that, in the general case, the bet-hedging region is defined by the conditions temporally varying case : spatially varying case : We now turn to the demonstration that the bet-hedging region in the spatially varying case is always broader than in the temporally varying case. Since hs l i > 0, we need to demonstrate that the following inequality always holds This can be proven from the chain of inequalities In Eq (28), the second and third inequalities are consequences of Jensen's inequality, since both x 2 and 1/x are convex functions. For the first inequality in Eq (28), since s > 0, we can use the result hx i i � hx j i i/j proved for i > j in [72]. Combining this result for (i = 3, j = 2) and (i = 2, j = 1), we obtain hx 3 i � hx 2 ihxi. Taking hxi ¼ h1= ffi ffi ffi s p i we finally prove Eq (28). Therefore, in the limit of small switching rates of the environment, the bet-hedging region is wider in the spatially varying case than in the temporally varying case.
In the opposite limit of high rates of environmental switch, the function to be optimized is linear, and the optimal strategy is a pure strategy, i.e. the bet-hedging region shrinks to a set of measure zero. In this case, the particular phenotype l adopted by the whole population is that maximizing ∑ i p i s il . This conclusion holds both for temporally and spatially varying environments.

Discussion
Understanding the precise mechanisms of population expansions is of utmost importance, not only for understanding species diversity, but also to cope with invasive species in new habitats [20][21][22][23], bacterial infections [24][25][26]73], and cell migration, such as those occurring during tissue renewal or cancer metastasis [5]. Phenotypic diversity is a convenient strategy for the success of population expansions in a broad range of contexts [20][21][22][23][24][25][26]. Although precise experimental measures are not easy to obtain, a recent study shows that populations with increased variability in individual risk-taking can colonize wider ranges of territories [74].
In this work, we proposed a general mathematical and computational framework to analyze such scenarios. In particular, we introduced a population model with diverse phenotypes that perform differently depending on the type of environment. We focused on the "optimal" degree of diversity leading to the fastest average population expansion in an environment fluctuating either in space or in time. We found that, contrarily to the well-mixed case, bet-hedging can be convenient in expanding populations. This result complements the study in [53] for a fixed habitat and supports the view that diversification is of broad importance for spatially-structured populations. For environments varying slowly in time, the expansion is relatively slow, and diverse communities can be optimal depending on the parameters. On the contrary, for fast environmental changes, the optimal population always adopts a unique strategy.
A remarkable outcome of our analysis is that spatial fluctuations create more opportunities for bet-hedging than temporal fluctuations, in that the region of parameter space where the optimal population is diverse, is always larger in the former case. One intuitive explanation is that in the case of spatial fluctuations, the population spends less time traversing favorable patches than adverse ones. This means that the beneficial effect of favorable patches is reduced with respect to the case of temporal fluctuations. Therefore, a pure risky strategy is less efficient in the case of spatial variability and can be more easily outcompeted by a diversified bet-hedging strategy.
The framework presented here can be extended to accommodate other scenarios. We have assumed that the fraction of individuals adopting each phenotype is fixed by the phenotypic switching rates. To understand the evolution of bet-hedging, it could be interesting to study scenarios in which the phenotypic switching rates are slower, so that phenotypes can be selected, and/or are themselves subject to evolution and selection [57,70]. Another potentially relevant extension would be to consider two-dimensional habitats. Although the classic theory for Fisher waves [7,8] is unaffected in higher dimensions, in the presence of spatial heterogeneity the front shape can become anisotropic, potentially affecting the results. Similarly, it would be interesting to analyze the combined effect of spatial and temporal variability. We also limited ourselves to the case where the different environments affect individual growth rates, whereas in general, one could also expect them to have an effect on motility [14,15,[75][76][77], opening the way for different forms of bet-hedging. Finally, the present study was limited to pulled waves. It would be interesting to study the effect of bet-hedging on pushed waves, for example to describe population expansion in the presence of an Allee effect [78,79].
It would be also interesting to experimentally test our results. Experiments of expanding bacterial colonies in non-homogeneous environments have already been performed and shed light, for example, on the evolution of antibiotic resistance in spatially-structured populations [80]. To perform experiments within the limits of our theory, a challenge can be to maintain the environmental variability sufficiently low to avoid exposing the population to an excessive evolutionary pressure. Similar problems appear, for example, in studies of range expansion of mutualistic bacteria [81]. An extension of the theory including both phenotypic and genetic diversity could account for these scenarios.
In summary, we have introduced a model to understand conditions favoring diversification of an expanding population. Our work provides a bridge between the theory of bet-hedging and that of ecological range expansion described by reaction-diffusion equations. The results of the model highlight the relation between population diversity and fluctuations of the environment encountered during range expansion. The flexibility and generality of our framework make it a useful starting point for applications to a wide range of ecological scenarios.