Stochastic delocalization of finite populations

Stochastic delocalization of finite populations

Lukas Geyrhofer, Oskar Hallatschek Biophysics and Evolutionary Dynamics Group, Max-Planck-Institute for Dynamics and Self-Organization, 37077 Göttingen, Germany,

The localization of populations of replicating bacteria, viruses or autocatalytic chemicals arises in various contexts, such as ecology, evolution, medicine or chemistry. Several deterministic mathematical models have been used to characterize the conditions under which localized states can form, and how they break down due to convective driving forces. It has been repeatedly found that populations remain localized unless the bias exceeds a critical threshold value, and that close to the transition the population is characterized by a diverging length scale. These results, however, have been obtained upon ignoring number fluctuations (“genetic drift”), which are inevitable given the discreteness of the replicating entities. Here, we study the localization/delocalization of a finite population in the presence of genetic drift. The population is modeled by a linear chain of subpopulations, or demes, which exchange migrants at a constant rate. Individuals in one particular deme, called “oasis”, receive a growth rate benefit, and the total population is regulated to have constant size . In this ecological setting, we find that any finite population delocalizes on sufficiently long time scales. Depending on parameters, however, populations may remain localized for a very long time. The typical waiting time to delocalization increases exponentially with both population size and distance to the critical wind speed of the deterministic approximation. We augment these simulation results by a mathematical analysis that treats the reproduction and migration of individuals as branching random walks subject to global constraints. For a particular constraint, different from a fixed population size constraint, this model yields a solvable first moment equation. We find that this solvable model approximates very well the fixed population size model for large populations, but starts to deviate as population sizes are small. Nevertheless, the qualitative behavior of the fixed population size model is properly reproduced. In particular, the analytical approach allows us to map out a phase diagram of an order parameter, characterizing the degree of localization, as a function of the two driving parameters, inverse population size and wind speed. Our results may be used to extend the analysis of delocalization transitions to different settings, such as the viral quasi-species scenario.

1 Introduction

Many population biological models simplify ecological complexity by assuming that habitats are homogeneous. Within these models, populations have a strong tendency to become evenly distributed in space as time proceeds. However, many biologically plausible habitats are strongly heterogeneous and provide gradients or localized spots of growth rather than homogeneous growth conditions. For instance, hydrothermal vents are highly localized environments that may have played a key role in the evolution of life itself [24]. Today chemo-synthetic bacteria using the heat, methane and sulfur of such vents form the basis of deep sea food chains, ranging from bacteria up to higher organisms as fish or octopuses [21, 22] (Figure (a)a). An oasis in the desert is another example, where member species shape the environment for one another, thereby feeding back on the stability of the oasis (Figure (b)b). When such heterogeneities are included in the population dynamics, the question arises whether populations will also become heterogeneously distributed, or localized, despite the mixing forces of diffusion or convection.

Figure 1: Examples for strongly heterogeneous environments. (a) Chemo-synthetic bacteria near hydrothermal vents. (Photo: NOAA) (b) Localized growth in a desert: an oasis. (Photo: Luca Galuzzi)

There have been several theoretical accounts to model the population dynamics at localized growth spots [2, 20, 3]. These studies focused in particular on the effect of a convection term, or wind, on the population localized at a growth spot, called the “oasis”. It has been found found that populations are localized as long as the convection forces remains below a certain critical value. As the convection speed is increased, the population becomes more widely spread around the oasis. The spread of the population is characterized by a correlation length scale that diverges at the critical convection velocity, signaling a continuous phase transition to extended states. This delocalization transition has also been investigated experimentally taking advantage of the fast growth rate of microbes [19, 16]. In these experiments, bacterial colonies were grown under hazardous UV light. Only a small localized (and moving) plate shielded the colonies, thus creating a preferred growth environment similar to the oasis considered in the theoretical models.

Closely related localization phenomena have been studied in evolutionary biology. In this context, populations evolve in a heterogeneous fitness landscape, rather than a geological landscape. Deleterious mutations act like a convection force that tends to push populations away from the optimal (i.e. most fit) genotype. For small mutation rates, natural selection is able to maintain a population cluster localized at the fitness optimum. For large mutation rates, however, populations drift away from the peak. The associated delocalization transition is called genetic meltdown in the population genetic literature [17, 26], respectively error threshold in the viral quasi-species literature [5, 6, 10, 27, 28].

The results on both the ecological delocalization transition and the error threshold, have been obtained largely by a linear stability analysis of the localized state. This approach ignores number fluctuations, which are however inevitable in populations of discrete individuals. The fluctuations are generated by variations in the offspring number, which itself depend on factors such as environment, lifespan or genetic background. The question thus arises how the standard view of the localization transition has to be modified to account for number fluctuations. We hypothesize that localized populations can be destabilized not only by convection but also by number fluctuations. Accordingly, we expect a modified phase diagram that depends on convection speed and population size, such as sketched in the tentative phase diagram (Figure 2).

Figure 2: Hypothesized phase diagram for finite populations localized at an “oasis”. Deterministic calculations of a population localized at a favorable growth spot (“oasis”) predict a sharp delocalization transition: When the convection speed exceeds a critical value , populations escape from the oasis. Our stochastic simulations of populations of fixed size suggest that this picture has to be modified to account for number fluctuations. i) Populations always delocalize on sufficiently long time scales, due to number fluctuations. The escape time, however, depends exponentially on population size. Thus, delocalization is effectively unobservable in sufficiently large populations. ii) In the presence of convection, the escape times are strongly reduced. In an exactly solvable model with fluctuations but finite populations, a sharp transition line can be derived (red dashed line). This line also describes qualitatively the crossover from small to very large escape times (compared to the coalescence time) in fixed population size models.

To test our hypothesis, we study quantitatively the localization transition in the presence of number fluctuations. We first simulate a one-dimensional population localized near a single lattice site with enhanced growth rate. This oasis is surrounded by a homogeneous habitat of lower growth rate. In order to fix the total population size, the excess population produced in each generation is removed by a (space-independent) death rate, which is adjusted in each time step. Our main goal is to reveal the impact of number fluctuations on the degree of localization, and the population distribution at steady state. We augment our simulation results by a mathematical analysis based on branching random walks subject to a tuned global constraint. This solvable model has been shown to be a good approximation to fixed population size models as far as the description of noisy traveling waves is concerned [9].

2 Simulations

2.1 Simulation model

Figure 3: Stochastic dynamics of the “oasis” model. A single time-step in our simulations consists of two sub-steps. (a) Individuals migrate between neighboring lattice sites, or demes, in general with a bias in one direction to account for convection. Number fluctuations are incorporated to account for the population turn over from one generation to the next. The subpopulation located at the oasis enjoys an increased reproduction rate . This first sub-step effectively describes a branching random walk with a branching rate that depends on the location. (b) In the second sub-step, the (global) population is regulated to comply with the constraint of fixed population size (see main text for details).

In our simulations, the population of reproducing individuals is distributed along a linear chain of sub-populations, called demes [14, 25]. These demes are considered to have a distance to their neighbors. The number of individuals in deme is denoted by . Each time step of duration (typically ) generations consists of the several sub-steps illustrated in Figure 3. First, we add to each occupancy variable the deterministic change given by


The first two terms on the right hand side account for biased migration: Individuals jump at rate from to and at rate from to . The last term represents deterministic growth at the oasis () at rate .

Next, we account for a stochastic change by adding given by


Here, is a random number drawn from a Poisson distribution with mean . The particular form of the noise in (2) ensures that i) occupancies do not drop below zero and ii) the generated variance in occupancy is given by . The variance produced thus equals the expected variance in a time for a (large) population of individuals that duplicate and die at the rate (in our computational unit of time). We note that the natural alternative of directly simulating individuals that duplicate and die is computational much more expensive if is much larger than .

After the previous sub-steps, the population size will typically be larger than the pre-set value . In order to restore the population size constraint, we add to each variable the change


with a (global) death rate


The rescaling of the occupancies by the factor of effectively represents a homogeneous culling of the population that is independent of the location. Notice that, although our simulations account for the fluctuations induced by the discreteness of the individuals, the occupancies are generally non-integer. The approximations involved are equivalent to the diffusion approach widely used in population genetics [13, 4]: As long as populations are large () and the temporal changes exhibit the same mean and variance as the discrete particle system, the deviations between both approaches should be small in the limit . This means in particular that the results become insensitive to the precise noise distribution used, as long as the first two moments are correct.

We impose periodic boundary conditions, as if demes were arranged on a circle like in the experiments of Ref. [19]. The number of demes is chosen large enough to ensure that population fits the simulation box. Time averages are taken after the simulations have become independent of their localized initial conditions.

2.2 Simulation results

Figure 4: Snapshots of simulations with fixed population size. Figure (a) depicts snapshots for our simulation model, in which individuals carry out an unbiased branching random walk subject to the global constraint of fixed population size (see Figure 3 for a description of the time step). All lattice sites (labeled by an index ) are equivalent, except for deme where the growth rate is increased by . The resulting density profiles are typically quite narrowly centered around the oasis (c.f. green and red lines). The average profile (black line) is much broader than typical profiles, due to large excursions of the center-of-mass. These excursions can be seen in the center-of-mass trajectory at the bottom part of the figure. Notice that the population occasionally escapes from the oasis, due to a rare large fluctuation. Such delocalizations are recorded when the occupancy at the the oasis drops below a critical value, . In the convection-less case, the population returns quickly after a random excursion (e.g. the first delocalization event in (a)). With convection, (b), the population usually traverses the (periodic) simulation box after a delocalization event until it reaches the oasis again. The escape time between two delocalizations is indicated by a black arrow. The green and red bullets indicate the spatio-temporal position of the center-of-mass corresponding to the red and green snapshots in the upper part of the figure. The convection speed in (b) is set to of the critical convection speed , at which delocalization is observed in the deterministic limit [2]. Notice that convection leads to more frequent escapes from the oasis. Simulation parameters were , , .

Snapshots of the population density are shown in Figure 4. For the chosen parameter values, fluctuations in the population density are clearly visible. The time traces of the center-of-mass trajectory (lower half of Figure 4) indicate that the population peak fluctuates wildly around the oasis, until it eventually escapes. The escape time (or delocalization time) is much greater in the case without convection (a) than with convection (b), but it eventually occurs. In fact, no sharp transition between localized and delocalized populations is detectable in our noisy system. If the population size is large enough, the delocalization is however unobservable within reasonable (simulation-) time. For definiteness, we say the population is on the oasis, when at least one individual is present there, , and off the oasis, when this condition is not met. We checked several alternative conditions, they yield similar or comparable results (see appendix A). Time averaged quantities (i.e. occupancies) are then obtained by averaging only over states with .

We measured in detail the time-averaged death rate , scaled by its value in the deterministic no-convection limit. If the overall death rate is large, a population can only be sustained in the vicinity of the oasis where net-growth is positive. Consequently, populations will be strongly localized at the oasis for . For smaller and smaller death rates, the population will be able to explore larger and larger parts of the habitat. When the death rate eventually changes sign (thus becoming a growth rate), growth becomes possible throughout the entire habitat. As a consequence, the population distribution becomes necessarily extended in space. Therefore, the death rate may be considered as an order parameter that measures the degree of localization of the population at the oasis. Consistent with this interpretation, the mean field theory predicts that continually decays from to as the convection speed is increased from to the critical value .

The death rates obtained from simulations are depicted in Figure 5, across a range of population sizes and convection velocities bounded by for which delocalization is observed in the deterministic limit. As expected from the mean field theory, we find that the death rate continually decays to zero with increasing convection . In addition, we find that death rates decrease when we decrease the population size. For any finite population, death rates go to at convection speeds significantly smaller than expected from the mean field predictions. Thus, increasing the demographic noise tends to destabilize the localized populations.

Figure 5: Mean death rates for varying population size and convection speed. In our simulations, the death rate is the factor with which the population has to be scaled down to keep its global size constant (c.f. section 2.1). Here, we display the scaled death rate normalized by its deterministic no-convection limit . According to the deterministic limit, the death rate should decay with increasing convection speed according to (solid black line) up to the critical convection speed  [2]. For finite populations the death rates are shifted to lower values roughly by the constant , as predicted by the analytical approach explained in section 3. As a consequence, smaller convection speeds are required for the death rates to cross , which signals the delocalization of the population. Our time average only includes states in which there is at least one individual at the oasis, . For convection speeds larger or equal than the deterministic critical velocity, this is just a tiny minority of all states. Hence, the points for large convection speeds (close to ) are slightly biased. Solid lines are obtained by our theoretical approach explained in section 3, in particular by the relation (48).
Figure 6: Distribution of escape times measured in simulations of fixed population sizes. (a) Distribution of escape times for different population sizes in the absence of convection. In all cases, the distribution is close to exponential. Deviations at early times (inset) are attributed to shape fluctuations, see main text. From the exponential tail of the distribution, a characteristic escape time can be extracted. This time strongly depends on the population size and convection velocity, as shown in figure (b). Solid lines in figure (b) represent the fit in (6). The three black arrows indicate the characteristic times of the distributions shown in (a). The simulation results were obtained using parameter values , .

Another more direct measure for the degree of localization is the escape time between a localization and the consecutive delocalization event, as indicated in the the center-of-mass trajectories in Figure 4. Histograms of escape times for various population sizes are shown in Figure (a)a. Notice that the escape time distribution is to a good approximation exponential,


indicating a memory-less delocalization process with a characteristic escape time . Deviations from the exponential distribution are discernible only on short times (inset in Figure (a)a). These deviations result from fluctuations in the density profile that temporarily satisfy the delocalization condition . After a relaxation process, the occupancy at the oasis is restored. Genuine delocalization is therefore described by only the exponential part of the distribution. The characteristic escape time can be extracted from the slope on a semi-logarithmic plot. This quantity is depicted in Figure (b)b as a function of population size for various convection velocities. Notice that decreases when either the population size is decreased (because the relative strength of fluctuations increases), or the convection speed is increased (because convection removes individuals from the oasis, thereby spreading the population). Numerically, our observations can be summarized by (Figure (b)b)


with and .

In summary, our numerical simulations indicate that, due to number fluctuations, populations inevitably escape the “pinning” at the oasis on sufficiently long time scales. The typical escape times, however, increase exponentially with population size, and thus delocalization is effectively unobservable for sufficiently large populations in the absence of convection. Our goal for the following is to develop an analytical framework to predict a critical population size in order to observe localization for a given convection velocity. This will be used to support the phase diagram in Figure 2.

3 Analytical approach

While our simulations are discrete in both space and time, they should be amenable to a continuous description if the time step and the reaction rate are small enough. Under these assumptions, the stochastic dynamics can be approximated by a reaction-diffusion equation of the form


Here, is a diffusion coefficient, which has to be set equal to deme per generation to match our simulations. The second term on the right hand side describes convection with velocity . The last “reaction” term summarizes all stochastic and deterministic changes due to the turn-over of the population. An example of the successful application of such reaction-diffusion equations in population biology is the Fisher-Kolmogorov wave equation, in which the reaction term describes logistic growth [7, 15]. This equation exhibits traveling wave solutions, which have been used to describe range expansions, epidemic outbreaks or the spreading of beneficial mutations  [12, 11, 9].

3.1 Deterministic limit

In the limit of very large populations, we can use a deterministic reaction term to describe the population dynamics in the oasis model,


The first term approximates the growth rate benefit provided at the oasis by a Dirac-delta-function . The pre-factor represents the product of the increased growth rate at the oasis and the linear size of the oasis, . The constant represents the death rate, which is tuned such that (7) admits a stable steady state.

The analysis of the resulting reaction-diffusion system has been carried out in Ref. [2]. It is found that the rescaled death rate in the absence of convection, , and continually decays to as the convection speed approaches the critical speed ,


It can be seen from Figure 5 that this mean field result clearly overestimates the death rates measured in the stochastic simulations.

The localized steady state profile is found to be given by


The density profile in equation (10) has characteristic width , given by


which diverges as a power law as the convection speed approaches the critical speed. This diverging length scale indicates that convection driven delocalization is a continuous phase transition in the deterministic limit.

Figure 7: Density profiles obtained for different models and approximations. (a) and (b) correspond to cases without and with convection, respectively. Symbols represent simulations with population size constraint. Green lines are predictions for the tuned constraint model, which we can derive analytically in the convection-less case (a). When including convection, only approximations are possible (b) (c.f. section 3.4). Blue lines represent the deterministic mean field approximations. For the parameters used (, , ), all approaches yield quite comparable density profiles. Significant differences arise however in the death rates even at these large populations, c.f. Figure 10.

The predicted density profiles are depicted in Figure 7, together with stochastic simulation data. The agreement is for the chosen parameter sets rather good, which is in contrast to the discrepancy between predicted and measured death rates.

3.2 General approach to finite population sizes

Modeling the spatial dynamics of finite populations typically leads to stochastic non-linear reaction diffusion equations: A non-linearity is necessary to fix the population size, and stochasticity arises due to the finiteness of the population. As a consequence, these models are generally hard to treat analytically. Systematic approximation schemes, such as perturbation expansions, often fail, as the example of noisy traveling waves demonstrates [23].

In a recent study, we have shown, however, that there is one particular form of the non-linearity that allows treating these and similar problems analytically [9]. We now review our methodology, and describe how the same method can be applied to diverse situations, including both ecological and evolutionary problems.

Our method applies to a whole class of models of branching random walks subject to a global constraint. A computational time-step of size of such models consists of two sub-steps, which are very similar to the sub-steps of our simulation model. The first sub-step evolves the population density according to a branching random walk, which is regulated in the second sub-step to comply with a global constraint. Specifically, the first sub-step takes the concentration field to an intermediate state


Here, the first term on the right-hand side describes the deterministic changes, which are due to net-growth and net-migration in our oasis model. In a time step , the associated change in density is given by , where the Liouville operator reads


Whereas the first two (particle conserving) terms describe biased diffusion with diffusivity and average bias , the last two reaction terms account for an overall negative growth rate except for the enhanced growth at the oasis. We note that, formally, we can reinterpret our model in an evolutionary context, if we consider to be a coordinate measuring either a quantitative trait with continuous characteristics, such as height, skin, or body mass, or a co-ordinate measuring the genetic distance. In both cases, the biased diffusion corresponds to the random change due to mutations that often tend to favor one direction. The delta-function growth spot would model an optimal phenotype which enjoys an increased fitness. The delocalization threshold would have to be considered as Eigen’s error threshold upon which selection is not able to maintain the optimal phenotype. More complicated quasi-species models can be easily devised, for instance, by considering a whole spectrum of mutations, recombination or more complex fitness landscape [8, 18, 4].

The stochastic second term on the left hand side in (12) accounts for the randomness in the reproduction process. The random function has vanishing mean and is uncorrelated across space and across the discrete time-steps,


The amplitude of the noise is chosen such that the produced variance represents the number fluctuations generated in the time span by a population of size , in which individuals branch and die at the same rate . The amplitude of the noise is thus consistent with the source of stochasticity in our simulations described in section 2.1.

In the second sub-step,


the population density is rescaled by a factor . Here, has to be chosen in such a way that the density at the end of the time step satisfies the constraint


The weighting function is at this point arbitrary. If one chooses one obtains a fixed population size model, as we have simulated above (up to the lattice discreteness). For any other choice of the population size is not fixed, but becomes a fluctuating quantity.

Notice that it is the constraint (16) that introduces a non-linearity in the problem, because depends on . The resulting noisy nonlinear reaction diffusion problems are intractable for general weighting functions . This can be seen from the first moment equation for the mean , for which one can derive the equation of motion


A detailed derivation of (17) can be found in [9]. Notice that the second term on the right hand side in general involves the second moment, . The equation of motion for the second moment, on the other hand, involves the third moment, and so on. The resulting infinite hierarchy of moment equations is intractable unless one is willing to introduce an ad-hoc truncation scheme.

However, it is possible to obtain a closed first moment equation for a special choice of the weighting function. Suppose, the solution of the nonlinear differential equation


exists, and is taken as our choice of the weighting function. Then, the non-linear second term of the moment equation (17) vanishes identically. As a consequence, the dynamics of the mean is controlled by a solvable linear differential equation


This equation is similar to the deterministic equation except for the negative last term that accounts for number fluctuations.

In summary, branching random walks occur whenever individuals proliferate, and are therefore a frequent element in both ecological and evolutionary models. A given ecological model often has an interpretation in an evolutionary context as well, if the dynamical variables are reinterpreted. It is usually difficult to analyze these models when populations are regulated in size, for instance by a logistic non-linearity. The algorithm described above provides a method to construct a solvable model of constrained branching random walks, in which the population size is fluctuating but finite. The algorithm to find and analyze this model consists of first solving equation (18) for the weighting function . Secondly, one has to determine the mean population density by solving the linear equation (19), which depends on the previously determined weighting function . It has been shown in [9] that this weighting function has a natural probabilistic interpretation: is the fixation probability that the descendants of an individual sampled from position will eventually take over the population. It has also been found in both Refs. [9, 8] that the tuned model is useful to determine those properties, which are only weakly dependent on the form of the constraint. In this way, it was possible to determine the speed of traveling waves in an primarily analytic way. In the following, we will use this approach of model tuning to derive the behavior of the death rate in our stochastic oasis model.

3.3 Stochastic oasis model

When the Liouville operator in (13) corresponding to our oasis problem is used, the steady state equations for and for take the form


To simplify the notation in the following, we introduce rescaled variables. We use to measure length in units of . We also introduce new variables for the population density and the weighting function . The equations of motion


together with the constraint


then form a closed set of equations that only depends on two control parameters, the rescaled velocity and the rescaled death rate 111The rescaled death rate was chosen such that, in the deterministic limit, we have in the convection less case. The rescaled velocity was chosen such that the deterministic delocalization threshold is at .

3.3.1 Stochastic oasis model without convection

It turns out that the convection terms in (22) and (23) considerably complicate the analysis. Therefore, we first focus on the case without convection, , and treat the windy situation in approximation later on, in section 3.4.

The equation for with restricted to can be mapped onto the mechanical problem of a point particle in a potential. When we interpret as time and as the location of the particle, its dynamics follows Newton’s equation


where the potential is given by


Using this analogy, we obtain the weighting function through the integral


where represents the sum of kinetic and potential energy of the moving particle. The value of is fixed by the requirement that the “trajectory” (here proportional to the density profile) approaches zero infinitely slowly in the limit . Thus, the total energy must equal the potential energy at . Hence, . The integral in (27) from to can then be performed analytically, which yields


The value has to be fixed by the boundary condition at . This boundary condition is obtained from the full equation for , equation (22), by integrating from to ,


which also holds for . In the absence of the convection term, we expect to be an even function. Therefore, we must demand the boundary condition


This boundary condition is satisfied by our solution (28) when we set


Due to the symmetry in , the full solution for is given by


Note that the form of the equations (2223) for and implies that is proportional to a stationary density profile,


where the proportionality constant is fixed by the constraint in (24).

Figure 8: Mean population size obtained for the analytical model with tuned constraint, in the absence of convection, . The population size diverges at , which corresponds to the deterministic limit. For , the population profile acquires a pronounced algebraic decay up to a length . This broad tail arises from large scale center-of-mass fluctuations. At the mean population size has a minimum.

By integrating over the density profile, , we obtain a relation between the mean population size and the rescaled death rate ,


The predicted dependence of the mean population size on the death rate is depicted in Figure 8. Notice that the (mean) population size reaches a minimum at . This minimum is a distinctive feature of the tuned model close to the delocalization transition. In this regime, the population frequently detaches from the oasis and takes large excursions away from the oasis until it is eventually pushed back to the oasis. During such excursions, the population explores the tail of the function . To obey the constraint (16) for a low value of , the population has to be rescaled to a large value, which leaves its footprint in an overall increased mean population size. Since this phenomenon cannot occur in a model with const., fixed population size models do not exhibit such a minimum in the population size.

3.3.2 Deterministic and stochastic limits

Our solutions for and are valid as long as


The bounds in (35) can be readily interpreted. For , we recover the exponential decay known from the deterministic limit (10),


The total population size scales as for . Since relative fluctuations in population size become small for large populations, we expect this leading order result to be independent of the form of the constraint. Indeed, the measured growth rates in the fixed model for agree very well with our prediction, as can be seen from Figure 10.

In the case of very small death rates, , we obtain an algebraic decay with distance,


This algebraic decay indicates that the center-of-mass of the population fluctuates strongly around the center of the oasis. Figure 9 depicts the density profile for small but finite , for which the decay is algebraic (as in (37)) up to a characteristic crossover scale (in original units)


For , the decay is exponential. We note that to clearly see the algebraic tail, the rescaled growth rate has to be finely tuned to a very small value. Yet, an algebraic decay of correlations and a diverging correlation length, is reminiscent of a second order phase transition. This raises the question of whether the underlying physics is generic for a wide range of systems.

Some more insight can be gained with the help of a simple scaling argument for the extinction dynamics close to the transition. For very small death rates, assume that the oasis sheds off subpopulations into the surrounding environment. Since this desert is only mildly deleterious, these trailing populations will survive up to a maximal time of roughly . Within this time, the subpopulations diffuse up to a length scale . Incidentally, this is precisely the correlation length, defined in (38), where the density profile crosses over from algebraic to a much faster exponential decay. The algebraic decay of the density profile within the correlation length can be explained by involving the statistics of branching random walks. Suppose, the populations are shed off from the oasis at a small, but constant rate . Each of these populations follows a critical branching process as long as the population is smaller then . (For larger populations, selection dominates and leads to a continual decrease in population size at rate .) First, the probability for a small population to reach a size larger than individuals through a critical branching process is proportional to . The time to reach the population size , conditional on survival, is roughly given by generations. During this time, the position of the center-of-mass will have diffused by a length scale . Combining these scaling arguments, we determine the probability that a sub-population reaches further than the distance as


Now, the mean population density at will be given by the probability density to reach , , times the typical density of a population that reaches . Such a population has size and a lateral spread of (acquired during a coalescence time of generations). Hence, the density of population that reaches scales as , and we obtain


The pre-factor clearly depends on the rate at which populations are shed off from the oasis.

The above scaling picture suggests, that the exponent is generic for a critical branching random walk with a boundary condition of a finite (source) population on one end. In this context, critical means that birth and death rate are equal. If the landscape is only nearly critical and characterized by a weak selective imbalance , then we expect to still see the power law up to the cutoff .

Figure 9: Density profile close to the stochastic threshold, . The plot depicts the density profile in (3332) predicted from the solvable model with tuned constraint. We used the parameters , such that the rescaled death rate is much smaller than one, . Notice the intermediate algebraic decay (with exponent ) up to the correlation length .
Figure 10: Relation between death rate and population size for different convection velocities . Symbols represent measured death rates in the fixed population size model for different population sizes and convection velocities (see legend). Dashed horizontal lines are the deterministic approximations (9), which strongly deviate from the simulations of finite populations. The solid red line is the predicted (mean) population size for a given death rate , c.f. (34), derived from the solvable model with tuned constrained in the convection-less case. Notice that these predictions also capture the quantitative trend of the simulation data for fixed population sizes (red symbols). The green and blue lines are approximations to the tuned model with convection, c.f. (48). Simulation parameters are , .

3.4 Stochastic oasis model with convection

A finite wind speed significantly complicates the mathematical analysis, because it results in a non-hermitian Liouville operator, defined in (13). As long as the population remains localized, it is however possible to return to a hermitian problem by the non-linear variable transformation,


In terms of the new functions and , the stationary equations of motion read


These equations cannot be solved in closed form, but a simple approximation is possible for small convection velocities where the exponential factors in (43,44) are close to . This approximation is justified if the exponent is smaller than on the considered length scale. If we consider the localization length of the deterministic problem (c.f. (11)) as the relevant scale, this implies that the convection speed should be much less than the deterministic threshold, .

After replacing the exponential factors in (43, 44) by , the equations of motion become identical to the convection free case with the replacement . Accordingly, the density profiles then obey


The resulting predictions for the population densities and growth rates have been used in Figure 7 to compare with our fixed simulations. The agreement is not perfect, due to the combined error generated by the model difference and our small velocity approximation, but reproduces well the qualitative behavior of the simulations.

3.4.1 Approximation for the death rate

Finally, we use the tuned model to derive a rather simple prediction for the death rate as a function of population size and convection speed. This approximation will help us better understand the delocalization mechanism in the presence of number fluctuations.

Integrating the governing equation (21) for the mean density profile over the space coordinate yields


Since the first term is proportional to the total population size , and the last term is fixed by the constraint (16), we obtain in general


This result is identical to the mean field prediction except for the negative last term. Hence, the death rate is lower than the mean-field prediction by the term . As a consequence, smaller convection speeds than in the mean-field approximation are required for the death rates to cross , which signals the delocalization of the population.

For using equation (47) predict the death rates for a given parameter set, we still need to find the occupancy at the oasis. This task requires the non-trivial solution of the set of equations in (43). By employing the small-velocity approximation (45) to express in (47), we obtain


This approximation (48) is used in Figures 5, 10 to explain fixed simulation data, and ultimately also to predict the fuzzy phase boundary in Figure 2, as implicit condition when the death rates reaches zero.

4 Discussion

We have studied localization of a finite population in an ecological context. Within our model, individuals enjoy an effective growth rate benefit in one site of a linear chain of subpopulations. We asked the question under which conditions this preferred growth spot, which has been called “oasis”, might be able to pin the population. We found that, under a model of globally regulated total population size, escape from the oasis inevitably occurs due to random number fluctuations. The escape time however sensitively depends on both the population size and the convection speed. The delocalization dynamics may be summarized in the phase diagram of Figure 2. When the convection speed is small enough and the population large enough, escape is extremely unlikely and the population can be considered as effectively localized. Due to the exponential dependence of escape times on and , the cross-over from long to short delocalization times is very sharp in the plane. The measured density profiles converge to the exponential decay for large population sizes, as predicted from deterministic mean field theories [2, 3]. The agreement with the deterministic theory generally becomes poorer as escape events become more frequent, either through an increase in convection velocity or through a reduction in population size. Deviations are pronounced in particular in the expected death rate that has to be used to control the population. In the absence of convection, this death rate is found to continually decrease as the population size is lowered. This is in contrast to the population size independent death rate predicted by the mean-field theory.

To analyze these findings, we have solved a variant of the model with a tuned population size constraint. In this model, a weighted sum over the occupancies of all demes is held strictly constant. The weighting function is chosen such that the first moment equation is linear. This model has the advantage of being analytically tractable. We found that the tuned model exhibited a sharp delocalization transition for all population sizes. The phase transition line in the plane could be determined in a simple analytical approximation, which is indicated by the dashed red line in the phase diagram Figure 2. Interestingly, we found that our tuned model exhibits a continuous delocalization transition driven solely by number fluctuations, even in the absence of convection. Close to the transition, the density profile exhibits a power law tail with exponent up to a correlation length that diverges at the transition.

It has to be remarked that the reason for the sharp phase transition in the tuned model is ultimately due to the inhomogeneity of the weighting function. Occupancy of the oasis is weighted more strongly than occupancy of lattice sites away from the oasis. This is in contrast to the fixed population size constraint, which has a spatially homogeneous weighting function: The occupancy of all demes are summed up, independent on their identity, and compared with the pre-fixed value. Therefore, the stochastic escape from the oasis is a characteristic effect of this degenerate weighting function. It remains to be seen, how one could explain analytically the exponential dependence of the typical escape time on population size and convection speed. Nevertheless, our simulations show that the narrow cross-over region from small to large escape times is located close to the sharp phase transition line of our tuned model in Figure 2. Furthermore, our tuned model predicts reasonably well the behavior of the death rate required to keep the population size finite. Compared with the mean-field approximation, the death rate is lowered by a value of roughly , in fair agreement with fixed simulations, c.f. Figure 5.

Our results also have a natural interpreted in an evolutionary context. Consider a population evolving in a fitness landscape characterized by an isolated fitness peak. Mutations cause random movements within this fitness landscape, typically biased in certain directions, and reproduction is often approximated by a branching process subject to the constraint of a fixed population size. The resulting quasi-species model of a finite population is a higher dimensional version of what we have considered in this paper, but otherwise analogous. We therefore expect our results could be useful for advancing the quasi-species theory, whose stochastic version has so far remained largely unexplored.


We would like to thank Jens Nullmeier and Paulo Pinto for many valuable discussions. We also acknowledge the anonymous referees for helpful comments on the manuscript. This work was supported by the Max Planck Society.

Appendix A Different choices for the definition “population on the oasis”

In our simulations, detailed in Section 2.1, the population always delocalizes from the growth spot if they are run long enough. Our analytically tractable model, on the other hand, describes populations that never fully detach from the oasis, provided the model parameter are within the regime in which the equations of motion, (2221), have non-trivial solutions. To compare both approaches, we needed to average observables only over those configurations in the simulation, which can reasonably be called “localized”. In Figure 11, we present the death rate for several alternative choices for the definition of localized states. We termed states localized if the center-of-mass (“com”) of the population was within a certain distance to the oasis and/or if the site of the oasis has an occupancy larger than a certain threshold value. The tested definitions gave very similar results. For definiteness, we chose for presenting averaged simulation data in the main text. Note that, with this definition, the death rate averaged over localized states cannot entirely decay to zero at the delocalization transition. This together with the model differences are reasons for the discrepancy between fixed simulations and our analysis of the tuned model, as seen in Figures 5, 10 and 11.

Figure 11: Different criteria for localization yield qualitatively similar results. In the main text, the population was defined to be localized at the oasis if at least one individual was present at the oasis, . Different criteria are of course possible. Here, we compare the following conditions: “”: the center-of-mass (com) is within a distance of the oasis. “”: the occupancies in demes , (oasis) and exceeds simultaneously. We also checked a combination of the last two conditions. Overall we find that our choice, , is in good agreement with the others, and the predictions from the tuned constraint model (black line).

Appendix B Relaxation spectrum in the case without convection

In the main text, we have analyzed the steady state of the analytically tractable tuned model. Here, we discuss the relaxation towards the steady state for the case without convection. The relaxation dynamics is governed by the eigenmodes and eigenvalues of the linear differential equation for  [1, 2, 20]. The following analysis shows that only a single bound state exists, together with a continuum of extended, delocalized states. These extended modes decay exponentially, with a rate faster than in (rescaled) time.

Figure 12: Spectral properties of the oasis model. (a) Without convection, , the spectrum can be obtained analytically and consists of the single bound state at and a continuous line of eigenvalues, the open interval , for the extended delocalized states. In analogy with [2] we expect the convection to render the eigenvalues complex (in pairs). The real part of these eigenvalues is approaching the bound state with increasing convection until the first extended state reaches the bound state at the the delocalization velocity . (b) Bound state () and one of the delocalized states (). Both states are normalized to fulfill .

The dynamics of the mean density is described by (c.f. (19) and (13))


The mean density can be decomposed into a basis consisting of (right) eigenfunctions of the linear operator . Without convection this operator is hermitian, rendering all eigenvalues and eigenfunctions real valued. Here, we only treat the convection-less case, . The eigenfunctions can then be chosen as an orthogonal basis, satisfying


In addition to (52) we can demand the normalization of , , which is used later in (56).

By the spectral theorem we can write


characterizing the relaxation dynamics of perturbations in the mean density. In (53) the coefficients are the projections of the mean density at initial time onto the eigenstates ,


Note that our constraint fixes the coefficient of the stationary eigenfunction with eigenvalue zero,


In going from (55) to (56) we used the fact that is proportional to with a proportionality constant .

Inserting the operator in (50) into (51) leads to


for and .

The eigenfunction is a linear superposition, , of two functions and , given by


Here, we have introduced the abbreviation . The parameter is fixed by the two boundary conditions


Different cases for have to be distinguished now. and render the term in the eigenfunctions either real or purely imaginary, respectively. First, we treat the case . The boundary condition (61) is fulfilled only for . For the second boundary condition (62) we obtain as a necessary requirement, leading to . This is the (already known) stationary state, . For , the arguments of the hyperbolic functions are (purely) imaginary. Using the identity between hyperbolic functions with complex arguments and trigonometric functions, we obtain oscillating (but bound) functions , which automatically comply with the first boundary condition. Straightforward, but tedious, calculations yield the coefficient , which depends only on and .

These calculations imply that the spectrum of the linear operator is given by


as illustrated in Figure (a)a. is the single bound state, whereas we obtain a continuum of extended, oscillating functions, which correspond to the delocalized states.

Analogous to previous results on the deterministic case [2], we expect that the real part of the eigenvalue of the extended states increases with increasing convection speed. Hence, perturbations relax more slowly than in the convection-less case. At the critical convection speed the relaxation time diverges.



  • [1] K. Dahmen, D. Nelson, and N. Shnerb. Population dynamics and non-hermitian localization. Statistical mechanics of biocomplexity, 527:124–151, 1999.
  • [2] K.A. Dahmen, D.R. Nelson, and N.M. Shnerb. Life and death near a windy oasis. Journal of mathematical biology, 41(1):1–23, 2000.
  • [3] M.M. Desai and D.R. Nelson. A quasispecies on a moving oasis. Theoretical population biology, 67(1):33–45, 2005.
  • [4] B. Drossel. Biological evolution and statistical physics. Advances in Physics, 50(2):209–295, 2001.
  • [5] M. Eigen. Selforganization of matter and the evolution of biological macromolecules. Naturwissenschaften, 58(10):465–523, 1971.
  • [6] M. Eigen, J. McCaskill, and P. Schuster. The molecular quasi-species. Advances in chemical physics, 75:149–263, 1989.
  • [7] R.A. Fisher. The wave of advance of advantageous genes. Annals of Human Genetics, 7(4):355–369, 1937.
  • [8] B. H. Good, I. M. Rouzine, D. J. Balick, O. Hallatschek, and M. M. Desai. Distribution of fixed beneficial mutations and the rate of adaptation in asexual populations. Proceedings of the National Academy of Sciences, 109(13):4950–4955, 2012.
  • [9] O. Hallatschek. The noisy edge of traveling waves. Proceedings of the National Academy of Sciences, 108(5):1783, 2011.
  • [10] J. Hermisson, O. Redner, H. Wagner, and E. Baake. Mutation–selection balance: Ancestry, load, and maximum principle. Theoretical population biology, 62(1):9–46, 2002.
  • [11] V.M. Kenkre. Results from variants of the fisher equation in the study of epidemics and bacteria. Physica A: Statistical Mechanics and its Applications, 342(1):242–248, 2004.
  • [12] V.M. Kenkre and M.N. Kuperman. Applicability of the fisher equation to bacterial population dynamics. Physical Review E, 67(5):051921, 2003.
  • [13] M. Kimura. The neutral theory of molecular evolution. Cambridge Univ Pr, 1985.
  • [14] M. Kimura and G.H. Weiss. The stepping stone model of population structure and the decrease of genetic correlation with distance. Genetics, 49(4):561–576, 1964.
  • [15] A.N. Kolmogorov, I. Petrovsky, and N. Piscounoff. Study of the diffusion equation with growth of the quantity of matter and its application to a biology problem. Bull. Univ. Moscow, Ser. Int. A, 1(6):1–26, 1937.
  • [16] A.L. Lin, B.A. Mann, G. Torres-Oviedo, B. Lincoln, J. Käs, and H.L. Swinney. Localization and extinction of bacterial populations under inhomogeneous growth conditions. Biophysical journal, 87(1):75–80, 2004.
  • [17] M. Lynch, R. Bürger, D. Butcher, and W. Gabriel. The mutational meltdown in asexual populations. Journal of Heredity, 84(5):339–344, 1993.
  • [18] R.A. Neher and B.I. Shraiman. Statistical genetics and evolution of quantitative traits. Reviews of Modern Physics, 83(4):1283, 2011.
  • [19] T. Neicu, A. Pradhan, D.A. Larochelle, and A. Kudrolli. Extinction transition in bacterial colonies under forced convection. Physical Review E, 62(1):1059, 2000.
  • [20] N.M. Shnerb. Extinction of a bacterial colony under forced convection in pie geometry. Physical Review E, 63(1):11906–11906, 2001.
  • [21] V. Tunnicliffe. The biology of hydrothermal vents: ecology and evolution. Oceanography and Marine Biology an Annual Review, 29:319–407, 1991.
  • [22] V. Tunnicliffe, A.G. McArthur, and D. McHugh. A biogeographical perspective of the deep-sea hydrothermal vent fauna. Advances in Marine Biology, 34:353–442, 1998.
  • [23] W. Van Saarloos. Front propagation into unstable states. Physics reports, 386(2-6):29–222, 2003.
  • [24] G. Wächtershäuser. Evolution of the first metabolic cycles. Proceedings of the National Academy of Sciences, 87(1):200, 1990.
  • [25] G.H. Weiss and M. Kimura. A mathematical analysis of the stepping stone model of genetic correlation. Journal of Applied Probability, 2:129–149, 1965.
  • [26] M.C. Whitlock, C.K. Griswold, and A.D. Peters. Compensating for the meltdown: the critical effective size of a population with deleterious and compensatory mutations. In Annales Zoologici Fennici, volume 40, pages 169–183. Helsinki: Suomen Biologian Seura Vanamo, 1964-, 2003.
  • [27] C. Wilke. Quasispecies theory in the context of population genetics. BMC evolutionary biology, 5(1):44, 2005.
  • [28] A. Wolff and J. Krug. Robustness and epistasis in mutation-selection models. Physical Biology, 6:036007, 2009.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description