The effect of spatial heterogeneity on the extinction transition in stochastic population dynamics

Stochastic logistic-type growth on a static heterogeneous substrate is studied both above and below the drift-induced delocalization transition. Using agent-based simulations, the delocalization of the highest eigenfunction of the deterministic operator is connected with the large N limit of the stochastic theory. It is seen that the localization length of the deterministic theory controls the divergence of the spatial correlation length with N at the transition. It is argued that, in the presence of a strong wind, the extinction transition belongs to the directed percolation universality class, as any finite colony made of discrete agents is washed away from a heterogeneity with compact support. Some of the difficulties in the analysis of the extinction transition in the presence of a weak wind, where there is a localized active state, are discussed.


2
In general, the use of deterministic rate equations and the neglect of the demographic (shot) noise may be justified if the bacterial density is large [12,13]. In a dilute system, e.g. close to the extinction transition, stochasticity dominates the dynamics and leads to substantial deviations from the predictions of the deterministic theory [14]. Studying the density dependence and the approach to the large density limit of the stochastic system enables one to track how the deterministic dynamics is modified in the presence of fluctuations.
The issues to be addressed here have two aspects: practical and conceptual. Experimentally, both the previous works [10,11] and a recently suggested setup [9] have concentrated on the identification and characterization of the extinction transition, where the effect of stochasticity is large. It is thus crucial to explore this regime theoretically. The fundamental aspects are also intriguing. For example, what defines the large density limit on a fluctuating substrate? Is it necessary that the density be large everywhere in order to allow the use of deterministic equations? These questions are analyzed and answered below for the single oasis case; before presenting it, we begin with a brief review of the deterministic theory [2,6,7].
In the context of population dynamics, the deterministic description of logistic growth of a motile population (with space-time-dependent density c(x, t)) on a 1d static spatially heterogeneous substrate is given by Here, D is the diffusion rate, v is the drift, a(x) is the local net growth/death rate and the nonlinear (Verhulst) term controls the saturation. Below the extinction transition one may neglect the nonlinear term, and the time evolution of the system,ċ = Lc, is determined by the linearized operator The real part of an eigenvalue of L reflects the overall growth/death rate of a colony that takes the spatial profile of the corresponding eigenvector. The picture is simplest for a homogeneous environment (a = a 0 ≡ σ − µ, where a 0 , the difference between the birth rate σ and the death rate µ, is independent of spatial location 2 ). Here L admits a full set of delocalized harmonic functions and the extinction/proliferation transition takes place when the highest eigenvalue, a 0 , crosses zero. Above the transition one gets the celebrated Fisher-Kolomogorov-Petrovsky-Piscounov (FKPP) equation, a generic description of an invasion of a stable state (c * = a 0 /b) into an unstable one (c * = 0). In the asymptotic long-time limit, this system supports a front that travels with constant speed v F = 2 √ Da 0 . The drift corresponds to a simple Galilean transformation. Things change when the translational invariance is broken, i.e. in the presence of spatial inhomogeneity. In the absence of drift (v = 0), L is a real symmetric operator that may support both extended and localized eigenstates. The eigenstates of a quantum particle in a single potential well, for example, are localized inside the well for low energies (small negative eigenvalues), and are extended above some threshold energy (large negative value).
In the presence of drift, or other non-Hermitian perturbation [1], the system undergoes another phase transition beyond which all the localized wavefunctions have become extended, and the corresponding eigenvalues have migrated from the real axis to the complex plane. This 3 transition was first analyzed by Hatano and Nelson [2] in the context of flux lines in high T c superconductors with columnar defects subjected to a tilted external magnetic field. Since then, many authors have considered this transition in different fields, e.g. hydrodynamics [3], random lasers [4] and quantum dots [5] among many others.
Two main types of heterogeneous growth are considered in the literature [6]- [8], [10,11]: a 'single oasis' case, where the growth rate is larger on a single domain, and the disordered case, where a(x) = a 0 + δa(x), δa being taken from some random distribution with zero mean. In both cases, the spectrum of the linear evolution operator L(v = 0) admits localized wavefunctions; if {φ 0 n (x), n } is the set of eigenfunctions and the corresponding eigenvalues, at least some eigenstates in the tail of the spectrum (or all the states for a disordered potential below 2d) are exponentially localized. The effect of small drift on the localized eigenstates is trivial: This 'gauge invariance' breaks down at v c n = 2D/ξ n , where ξ n is the localization length of the nth eigenstate. Above v c n this eigenstate delocalizes and the boundary conditions begin to play an important role: e.g. for periodic boundary conditions, the eigenvalues that correspond to delocalized eigenstates become complex [15]. The spectrum then takes the form of a 'bubble' in the complex plane, where the localized eigenstates correspond to the spectral points in the tail, since the localization length in the center of the band is larger. A non-Hermitian 'mobility edge' appears between the two regimes. Increasing v even more, the bubble spreads and captures more and more spectral points, and at the end the ground state also delocalizes. The Perron-Frobenius theorem [16] ensures that the highest eigenstate stays on the real line; the delocalization transition is identified by the breakdown of the trivial gauge (equation (3)) and the vanishing of the spectral gap [6,15]. Figure 1 shows some examples of the spectrum of L, together with a sketch of the phase diagram, for the single oasis scenario. Here, a discrete space model (1000 sites with periodic boundary conditions) has been used with a local two-particle annihilation interaction.
On a continuous substrate one should define an interaction length scale within which particles interact. This deterministic problem and the spectrum of L has been discussed in detail previously [2,6,7,15]. In the absence of drift, there is a single localized state at the right edge of the spectrum (if the oasis is large, a few localized states exist), followed by a continuum of states that correspond to extended eigenfunctions. Even a small drift is enough to push the delocalized eigenvalues to the complex plane, but the localized state only develops a slight asymmetry with almost no effect on the eigenvalue. Only for high enough drift does the highest eigenstate delocalize and the gap disappear. A change of a 0 corresponds to a rigid shift of the whole spectrum along the real line. Thus, three regimes exist in the drift-proliferation parameter space: the extinction region, where the real part of all the eigenvalues is negative; the localized region, where only the localized states admit positive growth rate; and the proliferation regime, where both localized and extended states may grow. Above v c 0 only the extinction and the proliferation regions exist.
How much of this picture, however, survives in a fully stochastic model? As we mentioned above the deterministic dynamics is an approximate description of the system. It becomes exact when the effect of stochasticity vanishes, i.e. when the density of agents is large, since the relative fluctuations scale with 1/ √ N . Technically, the exact stochastic master equation can then be replaced by a deterministic description using the Kramers-Moyal expansion, or more rigorously by van Kampen's expansion and related methods [12,13]. Joo and Lebowitz [17] (1) is the extinction region where the highest (in the real part sense) eigenvalue is negative. In region (2), only the localized state admits an eigenvalue with positive real part, and in regions (3a) and (3b) extended eigenstates become 'active'. In the right side of the figure, the highest eigenstate (right panels, φ 0 (x) is plotted versus x where the oasis is in the center) and the corresponding spectrum in the complex plane (left panels, Im( n ) versus Re( n ), the spectrum consists of a collection of points, each point corresponds to a single eigenvalue, the number of eigenvalues being equal to the number of points in the lattice) are plotted for the cases of no drift (v = 0), small drift (v = 0.1) and large drift (v = 0.5) for a single oasis. The deterministic phase transition happens when the rightmost spectral point touches zero; a change of a 0 shifts the spectrum rigidly, thus increasing a 0 (e.g. by decreasing µ) and keeping v fixed the system may undergo the deterministic transition either by first exciting a localized state (as exemplified by path 1, corresponding to the middle inset) or an extended state (path 2, corresponding to the bottom inset).
have already pointed out that in the limit of large N one should expect a population density distribution that follows the spatial features of the active eigenstates, i.e. the eigenstates for which Re( ) > 0. On the other hand, for a dilute system, i.e. close to the extinction transition the deterministic dynamics cannot be expected to be a reliable guide. In particular, the question is the fate of the sharp delocalization transition.
It should be noted that Brunet and Derrida [18] have already discussed the modification of the Fisher velocity as a result of demographic stochasticity, but they considered a different regime. These authors assumed that the distance from the extinction transition is fixed and positive, then took N to infinity and found that the velocity converges slowly to its deterministic value. Here we consider the N → ∞ limit at, or close to, the extinction transition. In particular, at the extinction point the Fisher velocity is zero, independent of how large N is.
Let us first present some general considerations. Janssen and Grassberger [19] suggested long ago that the extinction transition to a single absorbing state on a homogeneous substrate belongs (in the absence of special additional symmetries) to the directed percolation (DP) equivalence class, independent of the microscopic details of the stochastic process. DP is a continuous transition and the correlation length and correlation times diverge at the transition point with their characteristic exponents (see [14] for a general review). On a homogeneous substrate the correlation length is the only length scale of the problem. On a static heterogeneous substrate, on the other hand, another length scale appears-the localization length. How do these two quantities relate to each other? What are the properties of the stochastic extinction transition below and above the deterministic delocalization transition? In what sense is equation (1) a deterministic limit of a stochastic process when the effect of stochasticity is important, i.e. close to the extinction transition? Recently, we have addressed this last question numerically in the homogeneous case [20]; here we extend our numerics to the case of a nonuniform environment.
In order to simulate the heterogeneous system in the large N limit, we use an individualbased model that allows for an accurate determination of the transition point even for large N . We consider a logistic growth process on a one-dimensional lattice with periodic boundary conditions; the state of the system at t + t is determined by its state at t using the following procedure [21]. The number of agents at the ith lattice site is an integer n i , and each cycle of the Monte Carlo simulation involves two consecutive steps. The first step is the reaction: each of the agents at the site produces an offspring with probability (σ 0 + δσ i )(1 − n i /N ) t, and dies with probability µ t, so that the number of birth and death events at the site are given by binomial deviates. (While the deterministic results depend only on the difference σ − µ, to get generic results for the stochastic system one must include µ [20,22]. Without the death process, the deterministic limit of the extinction transition takes place when the birth rate vanishes, and does not belong to the DP equivalence class.) Here σ 0 is the bulk reproduction rate, and δσ = 0.2δ i,L/2 is the localized additional growth rate on the 'oasis', which we take to be a single lattice site. N represents the carrying capacity of a single site, beyond which the net growth vanishes. In the second step of the cycle, the diffusion step, every agent migrates with probability 2D t. Each migrating agent chooses its destination-to the left with probability q L = (1 + ν)/2 or to the right with probability 1 − q L . The migration takes place via a parallel update scheme; the n i are updated only after the diffusion step is completed. When t is small this Monte Carlo procedure converges, in the N → ∞ limit, to the Euler method of integrating ordinary differential equations.
In the linearized deterministic limit this model corresponds to an L-dimensional map, where L is the number of sites. This map is given by the multiplication of the reaction matrix, R i, j = δ i, j [1 + t (σ 0 + δσ i − µ)], by the diffusion matrix that takes the form (up to the boundary conditions) Diagonalizing the product D R one finds the highest eigenvalue˜ 0 and the corresponding eigenvector, φ i 0 ; adding another death process, where each particle in the MC simulation is selected to die after any cycle with probability 1/˜ 0 , ensures that the system is exactly at the transition point for N → ∞. In different words, the agent-based system is simulated with a parameter set that ensures 0 = 0 in the deterministic limit. Clearly, a system with finite carrying capacity N is always closer to extinction than the deterministic system when all other parameters are equal. This implies that, scaling the parameters as described above and increasing N , the system is always in the extinction phase and reaches the transition exactly at N = ∞. Below the extinction transition one can ignite the system many times, starting from a single particle located on the oasis, and monitor the overall number of birth events that took place at any spatial location, integrated over time. It turns out [20] that the particle density obtained in that way yields a smooth exponential profile, and the spatial correlation scale, ξ ⊥ , is extracted from the exponential falloff. In figure 2, we examine ξ ⊥ . Plotting it versus N on a log-log scale reveals the real meaning of the deterministic delocalization transition: below v c , i.e. when φ 0 is localized, the correlation length associated with the stochastic process first grows and then saturates to the deterministic localization length ξ 0 /[1 − vξ 0 /(2D)] (here ξ 0 is the localization length associated with the ground state of the deterministic evolution operator, see equation (3)). This demonstrates the fact that the state that becomes active at the transition is localized and the correlation length of the stochastic process cannot grow beyond this deterministic length.
On the other hand, above v c the correlation length grows unboundedly with N . The data are consistent with ξ ⊥ ∼ N 0.5 , which is (up to logarithmic corrections), the lifetime of a wellmixed system at the deterministic transition point [23]. This reflects the fact that the delocalized system is not really one-dimensional but rather '0 + 1' dimensional, with the spatial direction playing the role of time. Life in that system is a result of a drift from a source, not of uniform growth, and the lifetime of the colony at the oasis determines the spatial extent reached by its descendants.
The signature divergence of the correlation length at the delocalization transition thus does not survive when the transition is approached by increasing N in the localized regime. Nonetheless, there is a qualitative difference between the dynamics above and below the delocalization velocity even at finite N . For finite N and with constant drift velocity v, the system undergoes an extinction transition as a 0 = σ 0 − µ decreases. This may happen either via the localized phase (e.g. along path 1 shown by the arrow in figure 1), or directly to the delocalized phase (path 2 in figure 1). While for N → ∞ the transition happens when Re( 0 ) touches zero, for any finite N the transition takes place when a finite region of the upper part of the spectrum is above zero (inset of figure 3). For a single oasis (or otherwise when the number of oases is finite) all the localized states decay in the long run as a result of demographic stochasticity; only when the extended eigenstates are 'excited' (their eigenvalues cross to the positive real part of the spectrum) will the system be in its active phase. As a result, the scenarios 1 and 2 can be seen to differ significantly. Let us first consider path 2. Intuitively, above v c the colony is carried off the oasis by the wind, thus the large-scale properties of the system are identical with a homogeneous substrate with drift in the thermodynamic limit. More precisely, for finite N the transition occurs when a finite part of the spectrum, made of delocalized states, is already 'excited'  Close to the transition point for finite N , on the other hand, many linear states are already excited and the growth of the colony is unaffected by the nonlinear competition at short times. Only after the characteristic time ξ || does nonlinearity suppress the growth, leading to extinction. Within this growth period the system behaves deterministically and a 'Fisher front' starts to invade the empty region. As the wind velocity is larger than the Fisher velocity above v c [6], the maximum of C T (x) is shifted in the direction of the wind, as demonstrated in figure 3. This second regime vanishes at the deterministic limit; accordingly, the detachment of the peak from the nucleation point disappears upon increasing N .
The situation is completely different along path 1. The highest state is now localized, and its nonlinear interaction differs substantially from the interaction between extended states. If the localization length is finite the oasis region decouples from the rest of the system in the thermodynamic limit and the DP dynamics happens in parallel with the zero-dimensional stochastic process on the oasis. This decoupling, however, is impossible at the bulk DP transition, where ξ ⊥ diverges (see in this regard, [24]), and the oasis influences the entire system. A related issue is the transition in the presence of a finite density of randomly distributed oases: below v c a nonuniversal Griffiths phase appears between the active and the inactive parameter regions [25]. In the deterministic limit only the highest localized state becomes active at the transition, thus the Griffith phase admits no deterministic limit, and its width shrinks to zero. These last two observations suggest that the deterministic description of the system by means of excited localized eigenstates is insufficient, as the convergence of a finite N system to the deterministic limit is a very subtle issue, to be addressed in subsequent work.