Noise-induced stabilization and fixation in fluctuating environment
The dynamics of a two-species community of competing individuals is 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 -independent constant for global competition and decays like in the pairwise competition case. A 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 suboptimal (and extinction-prone) populations is discussed, as well as its relevance to stochastic tunneling between fitness peaks.
A fundamental problem in the fields of population genetics, evolution, and community ecology, is the fate of a single mutant (or invader) introduced in a finite population of wild types. For a fixed-size community of individuals, with Markovian, zero-sum dynamics driven by stochastic birth-death events, the mutant population eventually reaches either fixation or extinction. The classical analysis, provided by Kimura and his successors Crow et al. (1970); Ewens (2012), is focused on the neutral case, (where the dynamics is only due to demographic stochasticity, i.e., the noise inherent in the birth-death dynamics), and on time-independent selective forces (deleterious/beneficial mutation).
When the system is neutral (no fitness differences, all individuals are demographically equivalent) the chance of a single mutant to reach fixation is, by symmetry, . In general when the mutant population has abundance , .
Under fixed selection the fixation probability is,
In the weak selection regime, and , the effect of selection is negligible and the neutral result reemerges. When the mutation is beneficial (, but still ) and (strong selection regime), is -independent and converges, at large , to . A simple and intuitive argument for this result relays on the distinction between the region , where demographic fluctuations are dominant and the dynamics is more or less neutral, and the region , 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 , which is exactly Desai et al. (2007).
When the mutation is deleterious (), decays exponentially with in the strong selection regime, since now fixation involves a stochastic tunneling against a constant bias. Accordingly, for any practical purpose one may neglect the chance of a deleterious mutation to reach fixation when 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 Gavrilets (2010).
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 are . Environmental changes that affect coherently the fitness of entire populations lead to much stronger, , abundance variations Lande et al. (2003), 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 Leigh (2007); Bell (2010); Kalyuzhny et al. (2014a, b); Chisholm et al. (2014). Consequently, the study of temporal environmental stochasticity received a considerable attention Kessler and Shnerb (2014); Kessler et al. (2015); Sæther and Engen (2015); Cvijović et al. (2015); Kalyuzhny et al. (2015); Danino et al. (2016, 2018); Fung et al. (2016); Hidalgo et al. (2017); Danino et al. (2017). In parallel, a few recent experimental studies have considered the response of microorganism communities and their evolutionary pathways to fitness fluctuations Dean et al. (2017); Steinberg and Ostermeier (2016); Taute et al. (2014).
In some scenarios selection may activate a noise-induced stabilizing mechanism that could change dramatically the chance of fixation. The main aim of this work is do provide the results for these cases and to contrast them with the known results that were obtained in the absence of such a stabilizing mechanism. As we shall see, our analysis suggest a new mechanism that allows for long-term persistence of suboptimal mutant populations, a phenomenon that may facilitate stochastic tunneling through fitness valleys (see discussion section).
To begin, let us define two zero-sum competition models, one that does not allow for noise-induced stabilization (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 is the population fraction of the mutant, the chance of an interspecific ”duel” is and the chance of the mutant individual to win such a duel is defined to be . Accordingly, the deterministic growth/decay of (when time is measured in generations, elementary duels in each generation) satisfies the logistic equation,
If fluctuates in time such that its mean value is zero (a time-averaged neutral model Kalyuzhny et al. (2015)) the system performs an unbiased random walk along the axis. When the mean of , , 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 species with relative log-fitness to gain one individual in an elementary death-birth event is , while its probability to lose one is . In contrast with model A, here the fitness dependence is nonlinear. As a result, the deterministic dynamics satisfies, to second order in ,
While the linear () term in Eq. (3) gives, as in Eq. (2), a flow towards zero or one, the term has an attractive fixed point at . When is fixed in time this second term is negligible, but when fluctuates the term may dominate. Therefore, in model B environmental variations may induce stability through nonlinear fitness dependence Chesson and Warner (1981).
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 Danino et al. (2018).
Ii Model definitions and the Backward Kolmogorov Equation
For both model A and model B we assume that , where is the mean log-fitness and may take two values, (telegraphic, or dichotomous, noise).
The chance of the environment to switch (from to ) is per elementary birth-death event, so the sojourn time of the environment is taken from a geometric distribution with mean elementary birth-death steps, or generations. The chance of fixation when the mutant is in the plus (minus) state and its abundance is , (), satisfies the discrete backward Kolmogorov equation (BKE),
Where s are the transition probabilities in the state and the s are the probabilities in the states. The transition probabilities of model A are (we replaced by ),
and the corresponding probabilities for model B are,
The exact difference equations (II), with the appropriate set of s and s and with the boundary conditions , , may be solved numerically as a Markov chain Ashcroft et al. (2014) or as a linear system Danino et al. (2018), and these solutions are compared below with the analytic formulas.
One can translate this pair of discrete BKEs for into an equivalent set for and . Taking the continuum limit where is replaced by and functions of are expanded to second order in , a pair of coupled, second order differential equations for and emerges. In Danino et al. (2018) we have analyzed these equations in the limit of large and small . Using a dominant balance argument we showed that the dynamics is governed by a single second-order equation (in Danino et al. (2018) only the time to absorption was discussed, but the equation for is the homogenous version of the corresponding BKE, see Redner (2001)). This equation is,
with the boundary conditions and . Here is the strength of the environmental noise and is this strength divided by the strength of demographic stochasticity . 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 increases such that .
Eq. (7) may be solved using integrating factors, but this leads to complicated and hard to interpret nested integrals expressions. Instead one may analyze this equation in the inner (), middle () and outer () regimes and then match asymptotically the solutions in the large (more precisely, large ) limit. In the next section we present briefly the results for model A, following Danino and Shnerb (2017). Our purpose is to contrast these result with the outcomes of model B and to emphasize the effects of the noise-induced stabilizing mechanism.
Iii Model A: local competition and linear selection
The solutions of model A in the inner, middle and outer regimes are given by,
For one may use the uniform approximation solution for an arbitrary ,
For , if is negligible, i.e., if , the uniform approximation takes the form,
The chance of a single mutant to reach fixation is obtained by plugging into the inner solution,
The most important moral from the comparison between Eqs. (9, III) and Eq. (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 independent. This happens when where 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 , here the criteria for strong selection is , i.e.,
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 leaving systems. Accordingly, under environmental stochasticity systems may be in the weak selection regime even if is very large.
Two other characteristic scales in this system are and . For , , so . The chance of fixation of the mutant population becomes large when it reaches such that , thus,
This scale has been identified in Cvijović et al. (2015). Clearly, all of the three scales have similar features and they grow exponentially when is much smaller than . For a system with pure demographic noise and selection (Eq. 1), the chance of fixation becomes independent above and the condition for yields as well. While of (15) converges to this limit when , does not, and this reflect the inadequacy of our expression for in the limit Danino and Shnerb (2017).
While is small, may be large, so in the weak selection regime, and in particular in the time-averaged neutral scenario where , the chance of fixation decays logarithmically with system’s size as demonstrated in Figure 2.
When selection is weak an increase in increases the chance of fixation. In the purely demographic neutral case is determined by abundance so . 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 Cvijović et al. (2015). The situation is completely different in the strong selection regime Danino and Shnerb (2017), where is a monotonously decreasing function of . Here the reason is the divergence of when increases, meaning that the chance of the beneficial mutant population to enter the region of deterministic selective growth is much smaller.
Iv Model B: noise induced stabilization
In model B the dynamics is affected by the noise-induced stabilizing mechanism that facilitates the invasion 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 , as described above, and the environmental fluctuations are rapid, Eq. (3) for the deterministic dynamics of the population takes the form
This equation supports an attractive fixed point at
and is between zero and one if
Therefore, , the ratio between mean selection and environmental fluctuations, determines the qualitative behavior of the system. When the noise induces a stable coexistence point and the dynamics of model B differs substantially from the dynamics of model A. When 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 Danino et al. (2018) the time to absorption (either fixation or extinction) for model B was found to diverge like when is large and , while for the scaling is sublinear. Therefore, in this section we consider only the case.
Implementing the technique of asymptotic matching to model B equation when , the solutions for Eq. (7) with are,
where as before and the constants are given by,
The main difference between Eqs. (20) and (IV) and their model A counterparts, Eqs. (9) and (III), is the different scaling of . In model B, goes to zero when . Above this -independent point the chance of fixation in the middle region is fixed, , as demonstrated in Figure 3. A wide plateau appears in the middle region due to the force towards which is so strong that becomes almost -independent.
The uniform solution when has a relatively simple form (see Figure 3),
and the chance of fixation of a single mutant () is,
Amazingly, turns out to be an increasing function of , a behaviour that manifests itself in Figure 4. This phenomenon reflects the stabilizing effect of the nonlinear mechanism: the chance to reach the plateau does not depend on because the plateau occurs at values of that scale like . For example, in the large limit sticks to for any . For the chance of fixation decays like since this regime (in which our expressions fail) is dominated by demographic stochasticity. When increases the stabilizing mechanism wins against the demographic noise and leads to an increase of the chance of fixation.
in Eq. (24) is a multiplication of two factors: its numerator is the chance of establishment , which is the probability that the mutant population will reach the basin of attraction of the coexistence fixed point (the plateau). , that determines the denominator, is the chance that the mutant population will reach fixation given establishment. Accordingly, for a single beneficial mutant,
If the mutant is advantageous () and the system is in its strong selection regime ( or ), and . In terms of ,
This is a monotonously decreasing function of . When increases, the stochasticity becomes stronger and the stabilizing mechanism weakens Danino et al. (2016), 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 is the dominant effect and decreases.
Model A and model B differ even more dramatically in the weak selection regime , where . For model B, the chance of fixation becomes,
Unlike model A, where the chance of fixation decays logarithmically with in the weak selection regime, here the denominator of (27) is an constant. When , for example, the chance of fixation is one half of the chance of establishment: the effective strength of the selection bias is zero, so once the mutant population reaches the plateau its odds to win or to lose are equal. For nonzero there is a linear increase or decrease of as a function of . This relatively weak effect is demonstrated in panel (B) of Fig. 4.
Before concluding this section we would like to add a technical comment. is obtained from Eq. (II) by solving a linear problem (dividing a matrix by a vector). Using the sparsity of the matrix we were able to analyze systems with up to 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 solve that we have used quadrapole precision and this makes the numerics much slower and limits available system sizes to values up to .
V 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 () and deleterious () mutants, but the strong selection regime () 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 and . In the latter case, the growth rate of a deleterious mutant is still positive () 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 . Therefore, the chance of fixation, namely the chance to pick such a lucky sequence, decays like a power-law in . On the other hand, when , the route to fixation involves an improbable series of successes in consecutive elementary duels (reflecting demographic stochasticity) and in such a case decays exponentially in , like in the purely demographic case (1).
Looking at the equations above, one notices that the decay of when (and ) is due to the divergence of the term in the denominator of when the selection is strong. Accordingly, our theory predicts in that regime (as suggested in Assaf et al. (2013) for model A with ) a power-law decay, . The exponent of grows with and shrinks when the environmental stochasticity become stronger, as expected. However it does not show any qualitative shift when .
This difficulty turned out to be related to the failure of the continuum approximation that has been used when we have translated the difference equations (II) to the differential equation (7). As explained in Kessler and Shnerb (2007), this procedure fails when the differences between neighboring points (say, ) are too large and cannot be approximated using first and second derivatives. To overcome this obstacle a WKB approach was suggested Kessler and Shnerb (2007); 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 and , correspondingly.
As explained in the introduction [following Eq. (2)], in the absence of demographic noise and under model A dynamics , where . Accordingly, during generations the dynamics of satisfies,
so one may consider the stochastic process as a biased random walk along the axis. The random walker picks a left or a right move with equal chance , but left moves towards extinction () are longer than right moves towards fixation (). The backward Kolmogorov equation is,
with the boundary conditions and .
To implement the WKB technique, one writes and , where . In this WKB formalism we implement the continuum approximation to , i.e., for the logarithm of . Eq. (29) then takes the form,
and yields the transcendental equation,
Since is independent, , and given the boundary conditions one obtains,
If a ”single mutant” is associated with (since we impose the boundary condition at , we have to define it that way, but the results must be independent of this choice) the chance of fixation decays like a power-law in ,
The value of is given by (31). For small ,
in agreement with the definition of above. On the other hand, if is large,
and this expression diverges when , as required, to mark the transition to the exponential phase. Between these two limits, the expression
provides a decent approximation. The accuracy of this WKB argument is demonstrated in Figure 5.
This WKB argument gives the decay of in the power-law regime. When it fails, of course, since stochastic tunneling in this regime occurs only due to demographic stochasticity that was neglected in (29). For further discussion of the exponential phase (in the context of extinction times) see Ovaskainen and Meerson (2010); Kamenev et al. (2008).
The two models considered in this paper have one important feature in common: the abundance scale , below which the mutant population dynamics is dominated by fluctuations and above it by selection. This scale may become extremely large when the differences in the mean fitness are much smaller that the amplitude of the temporal fitness fluctuations, and one may easily imagine a situation where it becomes comparable or even larger than the effective size of an empirical community, meaning that the ecological or the evolutionary process takes place in the weak selection regime.
In the opposite, strong selection phase, the qualitative features of the chance of fixation are not much different between model A and model B. In both cases the chance of fixation for a beneficial mutant is -independent while the chance of fixation of a deleterious mutant falls like when and exponentially with if . 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., . Otherwise, the behavior of model B dynamics is qualitatively the same as model A).
On the other hand, in the weak selection regime 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 and under model A dynamics .
In sharp contrast with this result, in model B the system supports an attractive fixed point at . The plateau that characterizes in that case (Figure 3) reflects the effect of this attractive fixed point, marking the range of values which lies in its basin of attraction. The attractiveness of this coexistence point grows with and leads to an apparently paradoxical behavior: an increase of the chance of fixation with . Once this fixed point becomes dominant, the fate of the mutant population depends on its chance of establishment , i.e., of reaching the plateau, and on , the probability to jump from the plateau region to fixation. Since the plateau is wide, these two probabilities are independent, and so is the chance of fixation itself. Even if is huge, as long as it is much smaller than , model B yields an -independent value for 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 dynamics still supports an attractive fixed point at as long as . The chance of establishment is still -independent; it is the chance of fixation conditioned on establishment, , which goes to zero in that case. As shown in Danino et al. (2018), 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, , extinction-prone population with negative fitness.
This phenomenon may provide a plausible 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 bridge 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 Gavrilets (2010). To overcome this problem modern theoretical studies consider evolution on a neutral or nearly neutral (holely, 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 , 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 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.
This research was supported by the ISF-NRF Singapore joint research program (grant number 2669/17).
- Crow et al. (1970) J. F. Crow, M. Kimura, et al., An introduction to population genetics theory. (1970).
- Ewens (2012) W. J. Ewens, Mathematical population genetics 1: theoretical introduction, Vol. 27 (Springer Science & Business Media, 2012).
- Desai et al. (2007) M. M. Desai, D. S. Fisher, and A. W. Murray, Current biology 17, 385 (2007).
- Gavrilets (2010) S. Gavrilets, Evolution: the extended synthesis (eds M Pigliucci, GB Müller) , 45 (2010).
- Lande et al. (2003) R. Lande, S. Engen, and B.-E. Saether, Stochastic population dynamics in ecology and conservation (Oxford University Press, 2003).
- Leigh (2007) E. G. Leigh, Journal of Evolutionary Biology 20, 2075 (2007).
- Bell (2010) G. Bell, Philosophical Transactions of the Royal Society of London B: Biological Sciences 365, 87 (2010).
- Kalyuzhny et al. (2014a) M. Kalyuzhny, E. Seri, R. Chocron, C. H. Flather, R. Kadmon, and N. M. Shnerb, The American Naturalist 184, 439 (2014a).
- Kalyuzhny et al. (2014b) M. Kalyuzhny, Y. Schreiber, R. Chocron, C. H. Flather, R. Kadmon, D. A. Kessler, and N. M. Shnerb, Ecology 95, 1701 (2014b).
- Chisholm et al. (2014) R. A. Chisholm, R. Condit, K. A. Rahman, P. J. Baker, S. Bunyavejchewin, Y.-Y. Chen, G. Chuyong, H. Dattaraja, S. Davies, C. E. Ewango, et al., Ecology letters 17, 855 (2014).
- Kessler and Shnerb (2014) D. A. Kessler and N. M. Shnerb, Journal of Theoretical Biology 345, 1 (2014).
- Kessler et al. (2015) D. Kessler, S. Suweis, M. Formentin, and N. M. Shnerb, Physical Review E 92, 022722 (2015).
- Sæther and Engen (2015) B.-E. Sæther and S. Engen, Trends in ecology & evolution 30, 273 (2015).
- Cvijović et al. (2015) I. Cvijović, B. H. Good, E. R. Jerison, and M. M. Desai, Proceedings of the National Academy of Sciences 112, E5021 (2015).
- Kalyuzhny et al. (2015) M. Kalyuzhny, R. Kadmon, and N. M. Shnerb, Ecology letters 18, 572 (2015).
- Danino et al. (2016) M. Danino, N. M. Shnerb, S. Azaele, W. E. Kunin, and D. A. Kessler, Journal of theoretical biology 409, 155 (2016).
- Danino et al. (2018) M. Danino, D. A. Kessler, and N. M. Shnerb, Theoretical Population Biology 119, 57 (2018).
- Fung et al. (2016) T. Fung, J. P. O’Dwyer, K. A. Rahman, C. D. Fletcher, and R. A. Chisholm, Ecology 97, 1207 (2016).
- Hidalgo et al. (2017) J. Hidalgo, S. Suweis, and A. Maritan, Journal of theoretical biology 413, 1 (2017).
- Danino et al. (2017) M. Danino, D. A. Kessler, and N. M. Shnerb, arXiv preprint arXiv:1707.05089 (2017).
- Dean et al. (2017) A. M. Dean, C. Lehman, and X. Yi, Genetics 205, 1271 (2017).
- Steinberg and Ostermeier (2016) B. Steinberg and M. Ostermeier, Science advances 2, e1500921 (2016).
- Taute et al. (2014) K. M. Taute, S. Gude, P. Nghe, and S. J. Tans, Trends in Genetics 30, 192 (2014).
- Chesson and Warner (1981) P. L. Chesson and R. R. Warner, American Naturalist , 923 (1981).
- Ashcroft et al. (2014) P. Ashcroft, P. M. Altrock, and T. Galla, Journal of The Royal Society Interface 11, 20140663 (2014).
- Redner (2001) S. Redner, A guide to first-passage processes (Cambridge University Press, 2001).
- Danino and Shnerb (2017) M. Danino and N. M. Shnerb, arXiv preprint arXiv:1710.08807, , Journal of Theoretical Biology (in press) (2017).
- Assaf et al. (2013) M. Assaf, M. Mobilia, and E. Roberts, Physical review letters 111, 238101 (2013).
- Kessler and Shnerb (2007) D. A. Kessler and N. M. Shnerb, Journal of Statistical Physics 127, 861 (2007).
- Ovaskainen and Meerson (2010) O. Ovaskainen and B. Meerson, Trends in ecology & evolution 25, 643 (2010).
- Kamenev et al. (2008) A. Kamenev, B. Meerson, and B. Shklovskii, Physical review letters 101, 268103 (2008).