Survival of the scarcer in space

Survival of the Scarcer in Space Renato Vieira dos Santos Departamento de Fı́sica, Instituto de Ciências Exatas, Universidade Federal de Minas Gerais, CP 702, CEP 30161-970, Belo Horizonte, Minas Gerais, Brasil. arXiv:1304.5956v2 [q-bio.PE] 3 Jul 2013 e-mail: econofisico@gmail.com Ronald Dickman Departamento de Fı́sica, Instituto de Ciências Exatas, and National Institute of Science and Technology for Complex Systems, Universidade Federal de Minas Gerais, CP 702, CEP 30161-970, Belo Horizonte, Minas Gerais, Brasil. November 13, 2018 Abstract The dynamics leading to extinction or coexistence of competing species is of great interest in ecology and related fields. Recently a model of intra- and interspecific competition between two species was proposed by Gabel et al. [Phys. Rev. E 87 (2013) 010101], in which the scarcer species (i.e., with smaller stationary population size) can be more resistant to extinction when it holds a competitive advantage; the latter study considered populations without spatial variation. Here we verify this phenomenon in populations distributed in space. We extend the model of Gabel et al. to a d-dimensional lattice, and study its population dynamics both analytically and numerically. Survival of the scarcer in space is verified for situations in which the more competitive species is closer to the threshold for extinction than is the less competitive species, when considered in isolation. The conditions for survival of the scarcer species, as obtained applying renormalization group analysis and Monte Carlo simulation, differ in detail from those found in the spatially homogeneous case. Simulations highlight the speed of invasion waves in determining the survival times of the competing species. Keywords: Demographic stochasticity, Doi-Peliti mapping, Renormalization group, Monte Carlo simulation. 1 Introduction Coexistence of species and maintenance of species diversity are key issues in ecology as well as in conservation and restoration. A key idea in this context is the competitive exclusion principle, which asserts that similar species competing for a limited resource cannot coexist [1]. Individual-based stochastic models, both with and without 1 spatial structure, are useful for analyzing these questions [2]. There is some evidence [3] to suggest that spatial structure should facilitate coexistence of similar competitors. This is because when dispersal and interactions are localized, individuals tend to interact with conspecific neighbors more frequently than would be suggested by the overall densities at the landscape level; this is a common pattern in natural communities [4]. Computer simulations suggest that such spatial structures could lead to long-term coexistence of similar competitors and hence to the maintenance of high levels of biodiversity [5]. This mechanism has been encapsulated in the segregation hypothesis [6], which states that intraspecific spatial aggregation promotes stable coexistence by reducing interspecific competition. In the context of competition and extinction, the question of population size plays a fundamental role. In the simplest birth-and-death models, for example, the extinction probability of a viable, but initially small population is high, but decreases rapidly with population size.1 Recently, Gabel, Meerson, and Redner [7] proposed a stochastic model of two-species competition in which the species with smaller population size is, under certain conditions, less susceptible to extinction than the more populous species. This phenomenon, dubbed survival of the scarcer (SS) in [7], is somewhat surprising since conventional wisdom on population dynamics suggests that a smaller population is more susceptible to extinction due to demographic fluctuations. The authors of [7] use a variant of the Wentzel-Kramers-Brillouin-Jeffreys (WKBJ) approximation to obtain the tails of the quasi-stationary (QS) probability distribution [8, 9]. (Here and below, the quasi-stationary regime characterizes the behavior at long times, conditioned on survival.) They show that asymmetry in interspecific competition can induce survival of the scarcer species. While the analysis of [7] holds for well mixed (spatially uniform) populations, in view of the above mentioned segregation hypothesis it is natural to inquire whether the same conclusions apply for populations distributed in space. In this paper we study a stochastic model similar to that of [7], but with organisms located on a d−dimensional lattice, and investigate the conditions required for survival of the scarcer in space (SSS). On reproduction, a daughter organism appears at the same site as the mother; competition only occurs between organisms at the same site. In addition to these reactions, organisms also diffuse on the lattice. We study the model using two complementary approaches: field-theoretic analysis and direct numerical (Monte Carlo) simulation. Since the field-theoretic study requires specialized techniques, we summarize the method for readers less familiar with this approach. Starting from the master equation, we construct a representation involving creation and annihilation operators. This allows us to obtain a functional integral formulation, that is, a mapping of the stochastic process to a field theory. This now standard procedure is known as the Doi-Peliti mapping [10, 11, 12]. We consider the population dynamics in a critical situation, i.e., competing populations close to extinction. The dynamic renormalization group (DRG) is used to study the nonequilibrium critical dynamics of the populations, specifically, to determine how the model parameters transform under rescaling of length and time. The analysis furnishes a set of reduced stochastic evolution equations that describe the population dynamics for parameters close to those of a fixed point of the DRG. Analysis of these equations lets us map out regions of coexistence and of single-species survival in the space of reproduction and competition rates. In particular, we find regions in which the species predicted by mean-field theory to have the larger population in 1 See Ch. 2 of [2]. 2 fact goes extinct. We use Monte Carlo simulations to probe the QS population densities and lifetimes of the two species. We find SSS in a small but significant region of parameter space, in which the less populous species is more competitive. The simulations yield insights into the spatiotemporal pattern of population growth and decay, the speed of species invasion, and the QS probability distribution. The remainder of this paper is organized as follows. In section 2 we describe the model and write the effective action associated with the stochastic process. A brief analysis of mean-field theory, valid in dimensions greater than the critical dimension is performed. In section 3 we use some known results on multispecies directed percolation to obtain the renormalization group flow in parameter space. This analysis, combined with numerical simulations of the associated stochastic partial differential equations, allows us to determine the regions in parameter space in which SSS occurs. Results of Monte Carlo simulations are presented in Section 4. Section 5 closes the paper with our conclusions. 2 Model The model proposed in [7] may be described, with some minor modifications in the rates as specified below, by the following set of reactions:          β α A ⇀ A+A B ⇀ B+B β′ α′ A + A ⇀ ∅ (A) B + B ⇀ ∅ (B) ζ (1) ξ A+B ⇀ B A+B ⇀ A The reactions in the first line correspond to reproduction of species A and B. Those in the second line describe intraspecific competition resulting in the death of one individual (A/B) or of both (∅), while the third line represents interspecific competition. α (β), α′ (β ′ ) and ζ (ξ) are the rates of the reactions associated with species A (B). The relation between the rates in equation (1) and those of the original model [7] is (in the case of mutual annihilation): α ↔ 1, β ↔ g, α′ ↔ 1/K, β ′ ↔ 1/K, ζ ↔ ǫ/K, and ξ ↔ δǫ/K. δ < 1 implies an asymmetry in competition between species that favors B at the expense of A. In this case species B is more competitive than A; the opposite occurs if δ > 1. The authors of Ref. [7] constructed a phase diagram in the g-δ plane (shown schematically in Figure 1). The phase diagram includes results of both mean-field theory (MFT) and analysis of the master equation for a stochastic population model in a well mixed system, i.e., without spatial structure. In MFT, the populations of the two species are equal when gmf (δ) = (1+δǫ)/(1+ǫ). Analysis of the master equation in the QS regime yields  −1 equal extinction probabilities when gss (δ) = (1 + mδǫ)/(1 + mǫ) with m ≡ (ln 2)−1 − 1 . These conditions are plotted in figure 1; they intersect at the point (g, δ) = (1, 1). The phase diagram includes two “normal” regions, in which the more populous species is less likely to go extinct, and two “anomalous” regions, in which the scarcer species is more likely to survive. Thus in region I (where δ < 1 and B is more competitive in the sense defined above), species A is more numerous in the quasi-stationary state and simultaneously more susceptible to extinction, while in region II (where δ > 1 and A is more competitive), species B is more numerous, yet more likely to go extinct first. Regions I and II are anomalous since they correspond to an 3 g B dominant I II A dominant δ Figure 1: Phase diagram of well mixed system in the g-δ plane, as furnished by the analysis of Ref. [7]. Dashed blue curve: gmf = (1 + δǫ)/(1 + ǫ). Full red curve: gss = (1 + mδǫ)/(1 + mǫ). inversion of the usual dictum that susceptibility to extinction grows with diminishing population size. Our goal in this work is to determine whether a stochastic model of two-species competition with spatial structure, in which organisms diffuse in d−dimensional space, exhibits a similar phenomenon. In the following subsection we develop a continuum description for this spatial stochastic process. 2.1 Effective action We consider a stochastic implementation of the reactions of equation (1) on a d-dimensional lattice, adding diffusion (nearest-neighbor hopping) of individuals of species A and B at rates DA and DB respectively. Following the standard Doi-Peliti procedure mentioned in the Introduction, ignoring irrelevant terms (in the renormalization group sense) and with the so-called Doi shift [12] already performed, we obtain the following effective action: S(φ̄, φ, ψ̄, ψ) = Z dd x Z + Z dd x Z  dt φ̄[∂t + DA (σA − ∇2 )]φ − αφ̄2 φ + jα′ φ̄φ2 + α′ φ̄2 φ2 + ζψφφ̄  dt ψ̄[∂t + DB (σB − ∇2 )]ψ − β ψ̄ 2 ψ + jβ ′ ψ̄ψ 2 + β ′ ψ̄ 2 ψ 2 + ξφψ ψ̄ , (2) where σA ≡ −α/DA , σB ≡ −β/DB and j ≡ 1 (2) for individual death (mutual annihilation). Action S has four fields. The expected values of the fields φ and ψ (evaluated by performing functional integrals over the four fields, using the weight exp[−S]), represent the mean population densities of species A and B, respectively. The fields with overbars have no immediate physical interpretation, but are related to intrinsic fluctuations of the stochastic dynamics. In particular, terms quadratic in φ and ψ correspond to noise in the associated stochastic 2 2 evolution equations [13]. The terms ∝ φφ and ψψ correspond, respectively, to intrinsic fluctuations of species A and B. In the expansion of the action in powers of the fields, the lowest order term involving the product φψ would be φψφψ. Since it is of fourth order in the fields, this term is irrelevant, that is, the coefficient of this term flows to zero under repeated renormalization group transformations; it has therefore been excluded from the action. Fluctuations of one species nevertheless affect the other species via the coupling terms ∝ ζ and ξ in Eq. (2). −1 −1 ψ, φ, ψ̄ → θB ψ̄, and ψ → θB Using the definitions of σA and σB and the change of variable φ̄ → θA φ̄, φ → θA q q jα′ jβ ′ where θA ≡ α , θB ≡ β , we can write the effective action in the form: 4 S(φ̄, φ, ψ̄, ψ) = Z dd x Z n dt φ̄[∂t + DA (σA − ∇2 )]φ + ψ̄[∂t + DB (σB − ∇2 )]ψ + gA φφ̄(φ − φ̄) + gB ψ ψ̄(ψ − ψ̄) o hB φψ ψ̄ + hA φψ φ̄ , + with   gA       hA   gB      hB Dimensional analysis yields, (3) √ jαα′ ≡ −1 θA ζ √ jββ ′ ≡ ≡ (4) −1 θB ξ. ≡ d [hA ] = [hB ] = p2− 2 = [gA ] = [gB ], where [X] denotes the dimensions of X, and p has dimensions of momentum [12]. The upper critical dimension is therefore dc = 4, as for single-species processes such as the contact process and directed percolation. 2.2 Mean-field approximation For d ≥ dc , mean-field analysis, which ignores fluctuations in φ and ψ, yields the correct critical behavior. Imposing the conditions δS δφ = 0, δS δ φ̄ = 0, δS δψ   ∂φ ∂t = ∂ψ ∂t DA ∇2 φ + αφ − gA φ2 − hB φψ = DB ∇2 ψ + βψ − gB ψ 2 − hA φψ.  = 0 and δS δ ψ̄ = 0, we obtain the mean-field equations (5) Let L be a typical length scale and define X = gA φ/α, Y = gB ψ/β, s = DA t/L2 , and define a rescaled coordinate x′ = x/L. Then equation (5) can be written in dimensionless form as    ∂X ∂s = ∂Y ∂s =
∇2 X + γ X − X 2 − aXY ≡ ∇2 X + f (X, Y )
D∇2 Y + γ cY − cY 2 − bXY ≡ D∇2 Y + g(X, Y ) (6) with γ ≡ αL2 /DA , a ≡ βhB /(αgB ), b ≡ βhA /(αgB ), c ≡ gA β 2 /(gB α2 ), D ≡ DB βgA /(DA αgB ), f (X, Y ) ≡

γ X − X 2 − aXY and g(X, Y ) ≡ γ cY − cY 2 − bXY . Coexistence is possible in the stationary state if the conditions c > b and a < 1 hold, implying hA < βgA /α and hB < αgB /β. These conditions depend monotonically on parameters α, β, gA and gB , analogous to the coexistence conditions found in [7]. The same is not true in spatial stochastic theory, as we shall see in the next section. 3 Renormalization group (RG) flow In this section we apply a renormalization group analysis to determine regions of coexistence and of singlespecies survival in the space of reproduction and competition rates. Some years ago, Janssen analyzed a class 5 of reactions of the kind defined in equation (1) [14]. This work considered multi-species reactions of the form Xi ↔ 2Xi Xi → ∅ Xi + Xj → kXi + lXj , (7) where i and j are species indices and k, l are either zero or unity; this process is called multicolored directed percolation (MDP), with different colors referring to different species. The RG analysis of [14] shows that knowledge of the directed percolation fixed points is sufficient to determine the fixed points of the full MDP. As shown in that work, the parameter combinations gA /DB ≡ uA and gB /DA ≡ uB , related to intraspecific competition, flow under renormalization group transformations to the stable DP fixed point u∗A = u∗B = u∗ = 2ǫ/3 with ǫ = 4 − d > 0.2 Therefore, regarding the renormalization of intraspecific competition parameters, inclusion of other species behaving as in DP does not alter the fixed point. Janssen’s analysis shows that there are four fixed points for interspecific competition parameters hA /DA ≡ ∗ ∗ vA and hB /DB ≡ vB , which from reactions (1) are related to ξ and ζ. The first two are (vA , vB ) = (0, 0), which is unstable, and ∗ ∗ (vA , vB ) =  DA + DB 2 DA + DB 2 ǫ, ǫ DA 3 DB 3  (8) which is a hyperbolic fixed point if DA = DB .3 For DA 6= DB , it was conjectured [14] that the stability of the fixed points remains the same despite small changes in the renormalization group flow diagram topology. The ∗ ∗ ∗ ∗ other two fixed points for DA = DB are stable; they are given by (vA , vB ) = (0, 2u∗ ) and (vA , vB ) = (2u∗ , 0). To summarize, for DA = DB ≡ D0 , the interspecies competition parameters flow to one of the following four d−dependent fixed points (denoted as O,F,G, and H below), (vA , vB ) −1 −1 ξ θA ζ θB , D0 D0  =  → n O : (0, 0), F : (u∗ , u∗ ), o G : (2u∗ , 0), H : (0, 2u∗ ) (9) Figure 2 shows the renormalization group flow diagram. The unstable fixed point O corresponds to complete decoupling between the two species, and the hyperbolic point F to symmetric coupling. In the symmetric subspace, any nonzero initial value of v ≡ vA = vB , flows to F. For vB > vA (δ > 1), (vA , vB ) flows to G, so that species A is effectively unable to compete with B. Similarly, if vA > vB (δ < 1), the flow attains point H. In the following subsection we discuss the stationary states associated with these fixed points. 3.1 Steady state close to critical points Population dynamics in the critical regime is governed by a pair of coupled stochastic partial differential equations (SPDE), which are readily deduced from the action of equation (3)[12]: 2 Not 3 In to be confused with parameter ǫ used in equation 1. the nomenclature of [14], if species are of the same flavor. 6 vB IV G I III F II H O vA Figure 2: RG flow of the interspecific competition parameters vB and vA . Only the first quadrant is relevant due to the requirement of positive rates. There is an unstable fixed point O at (vA , vB ) = (0, 0) (black dot); a hyperbolic fixed point F at (vA , vB ) = (u∗ , u∗ ) (red dot), and two stable fixed points G and H, located at (vA , vB ) = (0, 2u∗ ) and (vA , vB ) = (2u∗ , 0) respectively (blue dots). The separatrix vA = vB is the boundary between the basins of attraction of G and H. ∂φ = D0 ∇2 φ + αφ − gA φ2 − hB φψ + η1 ∂t (10) ∂ψ = D0 ∇2 ψ + βψ − gB ψ 2 − hA ψφ + η2 ∂t (11) and where the noise terms η1 (x, t) and η2 (x, t) satisfy hη1 (x, t)i = hη2 (x, t)i = 0, and hη1 (x, t)η1 (x, t)i = 2gA φ(x, t)δ d (x − x′ )δ(t − t′ ), (12) hη2 (x, t)η2 (x, t)i = 2gB ψ(x, t)δ d (x − x′ )δ(t − t′ ). (13) Note that these multiplicative noise terms depend on the square root of their respective fields, and were obtained directly from the action, without any additional hypotheses. 3.1.1 SPDE in the vicinity of point F In the symmetric subspace, δ = 1, the RG flow is to the fixed point F, and parameters gA /D0 , gB /D0 , hA /D0 , and hB /D0 take the associated values. Therefore the SPDEs in (10) and (11) can be written as [15]: ∂φ = D0 ∇2 φ + φ − D0 u∗ φ2 − D0 u∗ φψ + η1 , ∂t (14) ∂ψ = D0 ∇2 ψ + gψ − D0 u∗ ψ 2 − D0 u∗ ψφ + η2 ∂t (15) where we set α = 1 and β = g. With the rescaling φ → φ/(D0 u∗ ), ψ → ψ/(D0 u∗ ) and x → (1/D0 )1/2 x, we can write equations (14) and (15) in the form: 7 1.0 1.0 0.8 Ψ Population Density Population Density Φ 0.6 0.4 Ψ 0.6 0.4 0.2 0.2 0.0 Φ 0.8 0 100 200 300 0.0 400 0 100 200 300 400 Time Time (a) Symmetric case with σ = 0.0 and g = 1. Initial populations are φ(x, 0) = 0.4 and ψ(x, 0) = 0.6. (b) Symmetric case with σ = 0.03 and g = 1. Initial populations are φ(x, 0) = 0.4 and ψ(x, 0) = 0.6. Figure 3: Numerical simulations symmetric SPDE. ∂φ = ∇2 φ + φ − φ2 − φψ + η1 , ∂t (16) ∂ψ = ∇2 ψ + gψ − ψ 2 − ψφ + η2 ∂t (17) with rescaled noise 2− d 2 (u∗ )2 φ(x, t)δ d (x − x′ )δ(t − t′ ) (18) 2− d 2 (u∗ )2 ψ(x, t)δ d (x − x′ )δ(t − t′ ). (19) hη1 (x, t)η1 (x, t)i = 2 (D0 ) hη2 (x, t)η2 (x, t)i = 2 (D0 ) 3 For d < 4, the rescaled noise intensity σ grows with diffusion rate D0 ; for d = 1 we have σ 2 = 2D02 u∗2 with u∗2 ≈ 2. Figures 3(a) and 3(b) show results of numerical simulations of equations (16) and (17) in d = 1 with g = 1. They show, for a one-dimensional lattice of L = 128 sites in the interval (−1, 1), the population densities averaged over 1000 Monte Carlo realizations for two different values of the noise intensity. The time interval T is T = 400, partitioned into N = 40000 steps so that ∆t ≡ T /N = 400/40000 = 0.01. Initial conditions are φ(x, 0) = 0.4 and ψ(x, 0) = 0.6. For σ = 0, population densities are almost time-independent as shown in Figure 3(a). This is not the case for larger noise values. In this case, the two populations tend to the same value (see figure 3(b)). Although Figs. 3(a) and 3(b) represent averages over many realizations, we have verified coexistence in individual runs. These simulations were performed using standard integration techniques for the Langevin equation with the XMDS2 software [16]. Since we are interested in showing only the initial temporal trends of the two population densities in the vicinity of fixed points, the use of standard techniques of integrating the Langevin equations is sufficient to reveal the qualitative nature of the solutions. Near a phase transition to an absorbing state, one of the population densities tends to zero and the standard numerical integration scheme fails. This standard algorithm is not suitable for extracting the more accurate results associated with critical exponents or asymptotic decays. In this case we would have to use more sophisticated algorithms such as those proposed in [17, 18, 19]. 3.1.2 SPDE in the vicinity of a stable fixed point Now consider the case δ 6= 0. With δ > 1, the parameters flow to point G and Eqs. (10) and (11) become 8 ∂φ = D0 ∇2 φ + φ − u∗ φ2 + η1 , ∂t (20) ∂ψ = D0 ∇2 ψ + gψ − u∗ ψ 2 − 2u∗ ψφ + η2 , ∂t (21) where we put α = 1 and β = g. With a rescaling similar to that used above, we have ∂φ = ∇2 φ + φ − φ2 + η1 , ∂t (22) ∂ψ = ∇2 ψ + gψ − ψ 2 − 2φψ + η2 ∂t (23) with η1 and η2 satisfying (18) and (19), respectively. For δ < 1 the situation is completely analogous to the case δ > 1, and the resulting equations are as above with the exchange of φ and ψ and changing the position of the factor g accordingly. If we neglect diffusion and noise terms, we have a system of ordinary differential equations having a fixed point associated with the stable coexistence state given by (φ∗ , ψ ∗ ) = (1, g − 2).4 Therefore, the condition for coexistence is g > 2. If g < 2, only species A survives for δ > 1. Figure (4) shows the result of a simulation with g = 2.2 in which the equations are integrated including both diffusion and noise. Numerical experiments indicate that the effect of these terms is to increase the value of g ' 2 for coexistence. Higher noise values imply higher thresholds g for coexistence. Φ Population Density 1.5 Ψ 1.0 0.5 0.0 0 100 200 300 400 500 Time Figure 4: Asymmetric case (δ = 1.5) with σ = 0.03 and g = 2.2. Initial populations for A and B species are φ(x, 0) = 0.1 and ψ(x, 0) = 0.9, respectively. For δ < 1 similar reasoning applies. The results are summarized in the phase diagram of Fig. 5. There are four stationary phases: a pure-A phase for δ > 1 and g < 2, a symmetric pure-B phase for δ < 1 and g < 1/2, and two (disjoint) coexistence phases. The latter arise when the less competitive species proliferates sufficiently faster than the more competitive one. 3.1.3 Nature of phase transitions From the equations that omit the diffusion and noise terms, one may infer the nature of the phase transitions in the diagram of Fig. 5. From Table 1, which shows the fixed points of the equations for different values of δ, we see that in general the values in the final column do not match for δ 6= 1. This fact indicates that the vertical blue line in the diagram (5), i.e., the line δ = 1, is a line of discontinuous phase transitions. Similarly, examining the 4 For the case δ < 1, fixed points are (φ∗ , ψ∗ ) = (1 − 2g, g), and coexistence condition is 0 < g < 1/2. 9 g AB B r II A I p AB δ b Figure 5: Stationary phase diagram in the g-δ plane for stochastic model with spatial structure. The reference points have coordinates: b = (1, 0), p = (0, 1/2), r = (3, 2). For δ > 1, when g < 2 only species A is present, while for g > 2 there is coexistence. Similarly, for δ < 1, when g > 1/2 only species B is present, and for g < 1/2 there is coexistence. The diagonal line connecting points p and r is the equal-population criterion furnished by mean-field theory, gmf (δ) = (1 + ǫδ)/(1 + ǫ), with ǫ = 1. In regions I and II, the more populous species (A or B, respectively), as predicted by mean-field theory, in fact goes extinct. The thick, dashed, blue √ and black curves represent gmf and gsw for the value of ǫ that maximizes the area between them, i.e., ǫ = 1/ m ≈ 0.66. Table 1: Simplified equations and their fixed points for different values of δ δ =1 >1 <1 Equation Fixed Point 2 2 φ̇ = φ − φ − φψ and ψ̇ = gψ − ψ − φψ φ̇ = φ − φ2 and ψ̇ = gψ − ψ 2 − 2φψ φ̇ = gφ − φ2 and ψ̇ = ψ − ψ 2 − 2φψ (φ∗ , ψ ∗ ) = (0, g) (φ∗ , ψ ∗ ) = (1, 2 − g) (φ∗ , ψ ∗ ) = (g, 1 − 2g) expressions in the final column of Table 1, we infer the nature of phase transitions along the horizontal lines at g = 1/2 and g = 2. In regions of coexistence, the densities of the minority species grow continuously from zero; thus these transitions are continuous. One should of course recognize that the predictions for the phase diagram are essentially qualitative; quantitative predictions for nonuniversal properties such as phase boundaries are not possible once irrelevant terms have been discarded. Regarding the nature of the transitions, we expect the continuous ones, as furnished by this mean-field-like analysis, to in fact belong to the DP universality class [14]. The line of discontinuous transitions contains a portion (between levels p and r) that is between absorbing subspaces (only on species present in each phase), and other portions that separate (from the viewpoint of the species undergoing extinction) an active and an absorbing phase. Discontinuous transitions of this nature are not expected to occur in one dimension [20]. According to mean-field theory, a point in region II corresponds to ρB > ρA . If we increase the competitiveness of species A, increasing δ while maintaining g fixed, there will be a point beyond which species B no longer has the greater population density. Improvement of the competitiveness of species A will make it the more populous species if organisms are well mixed. In the case of a spatially structured population, by contrast, any competitive advantage of species A, no matter how small, is sufficient to dominate the B population (if the latter does not proliferate very rapidly) when the population densities are very low, as is the case near criticality. Thus in a situation in which mean-field theory predicts species B to be the majority, we can have species A exclusively. This can be seen as a strong version of the survival of the scarcer phenomenon. 10 To compare the predictions of the mean-field and spatially structured descriptions, a pair of dashed lines are plotted in Fig. 5, representing equations for gmf and gsw as defined previously, for parameters such that the area between them is maximum, thereby maximizing the region in which SS occurs. [The maximum area √ is obtained using ǫ = 1/ m with m = [(ln 2)−1 − 1]−1 , in Eq. (1)]. Given the much larger area of the region exhibiting SS in the spatial model, compared with that in the mean-field analysis (i.e., the area between the two dashed lines for 1 < δ < 3,) we see that the SS phenomenon can be greatly intensified when spatial structure is included. Analogous reasoning applies for region I (green). (For the parameter values used, however, this reasoning does not apply for δ > 3, since in this case the area between the dashed lines can be arbitrarily large. Also we no longer have a region analogous to region II with species A only.) In this way we see that the phase diagram for the spatial stochastic model is rather different from that found in [7]. Moreover, under certain conditions, the survival of the scarcer phenomenon is strengthened. 4 Simulations Since the preceding analysis involves approximations whose reliability is difficult to assess, we perform simulations of a simple lattice model to verify SSS. This approach permits us to access the QS regime of the spatial model, in finite systems. We consider the following spatial stochastic process, defined on a lattice of Ld sites. Each site i is characterized by nonnegative occupation numbers ai and bi for the two species. The organisms (hereafter “particles”) of the two species evolve according to the following rates: • Particles of either species hop to a neighboring site at rate D. • Particles of species A(B) reproduce at rate λA (λB ). • At sites with two or more A particles, mutual annihilation occurs at rate αA ai (ai − 1), and similarly for B particles, at rate αB bi (bi − 1). • At sites having both A and B particles, the competitive reaction B → 0 occurs at rate ζA ai bi and the reaction A → 0 occurs at rate ζB ai bi . The model is implemented in the following manner. In each time step, of duration ∆t, each process is realized, at all sites, in the sequence: hopping, reproduction, annihilation, competition. In the hopping substep, each site is visited, and the probability of an A particle hopping to the right (under periodic boundaries) is set to ph = Dai ∆t/(2d). A random number z, uniform on [0, 1) is generated, and if z < ph the site is marked to transfer a particle to the right. Once all sites have been visited, the transfers are realized. The same procedure is applied, in parallel, for B particles hopping to the right. Subsequently, hopping in the other 2d − 1 directions is realized in the same manner. In the reproduction substep, the A particle reproduction probability at each site i is taken as pc = λA ai ∆t. A random number z is generated and a new A particle created at this site if z < pc . An analogous procedure is applied for reproduction of B particles. In the annihilation substep, a probability pa = αA ai (ai −1)∆t is defined at each site, and mutual annihilation (ai → ai − 2) occurs if a random number is < pa . Again, an analogous procedure is applied for annihilation of B particles. 11 Finally, at sites harboring both A and B particles, probabilities pA = ζA ai bi ∆t and pB = ζB ai bi ∆t are defined. The process B → 0 occurs if a random number z < pA , while the complementary process A → 0 occurs if pA ≤ z < pA + pB . Note that at most one of the processes (A → 0 and B → 0) can occur at a given site in this substep. The time step ∆t is chosen to render the reaction probabilities relatively small. To do this, before each step we scan the lattice and determine the maximum, over all sites, of ai , bi , and ai bi . With this information we can determine the maximum reaction rate over all sites and all reactions. Then ∆t is taken such that the maximum reaction probability (again, over sites and reactions) be 1/5. In this way, the probability of multiple reactions (e.g., two A particles hopping from i to i + 1, etc.) is at most 1/25, and can be neglected to good approximation, particularly in the regime in which occupation numbers are typically small. We simulate of the process on a ring (d = 1) using quasistationary (QS) simulations, intended to sample the QS probability distribution [21, 22]. In these studies the system is initialized with one A and one B particle at each site. When either species goes extinct (an absorbing subspace for the process), the simulation is reinitialized with one of the active configurations (having nonzero populations for both species) saved during the run. Following a brief transient, the particle densities fluctuate around steady values. We search for parameter values such that species A is more numerous, despite being much less competitive than species B (i.e., ζB >> ζA ). Given the large parameter space, certain rates are kept fixed in the study: We set D = αA = αB = 1/4. To understand the competitive dynamics we first need to determine the conditions for single-species survival. We use QS simulations as well as spreading simulations [23] to estimate the critical value of λ for survival of single species, given α = 1/4. In spreading simulations the initial configuration is that of a single site with one particle, and all other sites empty. We search for the value of λ associated with a power-law behavior of the survival probability P (t) and the mean population size n(t). These studies yield λc = 1.267(1) as the critical point for survival of a single species, i.e., without interspecies competition. (Here and in the following, figures in parentheses denote statistical uncertainties in the final digit.) The phase transition is clearly continuous; details on scaling behavior will be reported elsewhere. In the two-species studies we set λA = 1.6 and ζA = 0.005, so that species A is well above criticality but weakly competitive. We then vary λB and ζB , monitoring the QS population densities ρA and ρB of the two species, as well as their QS lifetimes, τA and τB . The latter are estimated by counting the number of times, in a long QS simulation, that one or the other species goes extinct. Survival of the scarcer is then characterized by the conditions ρA > ρB and τA < τB . For the parameter values studied, we observe SSS for λB ≈ λc and ζB ≫ ζA . Figure 6 illustrates the variation of the population densities, and of the ratio τA /τB , as ζB is varied for fixed λB = 1.25, on a ring of 200 sites. In this case, ρA > ρB throughout the range of interest; the scarcer species, B, survives longer than A for ζB ≥ 0.25. Figure 7 shows a typical evolution, for ζB = 2 and other parameters as in Fig. 6. (The spatiotemporal pattern observed for other parameter sets exhibiting SSS is qualitatively similar.) The system is divided into A-occupied regions, B-occupied regions, and voids. The voids arising in B-occupied regions are larger, as this species is nearer criticality. Species A, which is well above criticality, rapidly invades empty regions, but is in turn subject to invasion by the more competitive species B. Sites bearing both species are quite rare, occurring only at the frontiers between regions occupied by a single 12 species. Figure 6: Quasistationary population densities ρA (blue) and ρB (red), and lifetime ratio τA /τB (black) versus ζB , for parameters as specified in text. Figure 8 shows the result of a systematic search in the ζB -λB plane on a ring of L = 200 sites. These data represent averages of 200 realizations, each lasting 105 time units. SSS is observed in the region bounded by the two curves; the lower curve corresponds to τA = τB while on the upper we have ρA = ρB . There is an interval of λB values, including λc , on which SSS is observed for any value of ζB greater than a certain minimum. For larger values of λB , SSS occurs only within a narrow range of ζB . In the latter region we have α > ζB > ζA , so that intraspecies competition is stronger than competition between species. Although we use a rather small system to facilitate the search, SSS is not restricted to this system size. For L = 800 for example, the range of λB values admitting SSS is more restricted (from about 1.24 to 1.27) but at the same time the effect can be more dramatic. For λB = 1.27 and ζB = 0.2, for example, we find ρA = 0.74, ρB = 0.49, and τB ≈ 100τA . Studies using λA = 1.35 and ζA = 0.01 yield a qualitatively similar diagram, leading us to conjecture that the form of the region exhibiting SSS in one dimension is generically that shown in Fig. 8. To summarize, we verify SSS in a limited but significant region of parameter space. 13 Figure 7: A sample evolution on a ring of 200 sites, of duration 4000 time units (with time increasing downward). Blue: sites with ai > 0 and bi = 0; red: sites with ai = 0 and bi > 0; black: sites with both species present. In the inset of Fig. 8 the data are plotted in the δ-λB plane, for comparison with the theoretical predictions shown in Fig. 5. (Recall that parameter g corresponds to the ratio λB /λA , which is smaller than unity here. The parameter δ = ζA /ζB , i.e., the ratio of the interspecies competition rates.) The region exhibiting SSS is rather different from the predictions of both MFT and stochastic field theory. In particular the lower boundary in the δ-λB plane is not horizontal and it is unclear whether the SSS region extends to δ = 1. (As λB is increased, SSS is observed in an ever more limited range of ζB values, making numerical work quite time-consuming.) Here it is important to recall that irrelevant terms (in the renormalization group sense) are discarded in the theoretical analysis, so that predictions for the phase diagram are only qualitative. The evolution shown in Fig. 7 suggests that survival of a given species depends on the speed of invasion into territory occupied by the competing species; we study this issue via the following procedure. We first use simulations to prepare collections of configurations drawn from the QS distribution for each species in isolation. (We saved 1000 such configurations, on a ring of L = 200 sites, for each λ value of interest.) To determine the speed of invasion we prepare initial configurations for the two-species system by placing a randomly chosen, saved configuration with species A only at sites 1,...,L, beside a similar configuration, with species B only, occupying sites L + 1, ..., 2L. In the subsequent dynamics, the two species interact at the boundary, leading to invasion by one species or the other in different realizations. Fig. 9 shows the evolution of the A and B density profiles, averaged over many realizations (3000 or more). Given the average density profiles ρA (x, t) and ρB (x, t) we identify the interface position Xi (t) of species i via the condition ρi [Xi (t), t] = ρi /2, with ρi the QS density of species i in isolation. Plotting the interface positions versus time, we find that they attain steady velocities 14 zz z zz zz zz zzzz zz z z z z z z z z z z z z z z Figure 8: Simulations in one dimension: SSS is observed in the region bounded by the two curves; see text for parameters. Inset: the same data plotted using δ = ζA /ζB instead of ζB on the horizontal axis. after a certain initial transient; the velocities are plotted versus ζB in Fig. 10. (For the system size used here, steady velocities are attained after about 500 time units.) The interface velocity vB of species B increases slowly with ζB , but that of species A (vA ) falls rapidly, becoming negative for ζB greater than about 0.13. For the parameters considered here, the mean lifetime of the scarcer species (B) is longer for ζB ≥ 0.113; very near this value, vA becomes smaller than vB . This suggests that, as one would expect, the preferential extinction of the more populous species is associated with a smaller rate of spreading, so that A domains eventually die out due to invasion by B. To close this discussion of simulation results, we present in Fig. 11 a portrait of the QS probability distribution in the SSS region. The basic features identified in the spatially uniform case (i.e., a well defined maximum corresponding to coexistence with NA > NB , as well as single-species QS distributions), persist in the presence of spatial structure. We may speculate that the higher extinction probability for species A is associated with the more diffuse barrier in the small-NA region (near the y-axis) as compared with the corresponding barrier near the x axis. More detailed simulations, including other parameter sets, larger systems, and simulations on the square lattice, are planned for future work. 15 Figure 9: Mean density profiles ρA (x, t) (maximum at right) and ρB (x, t) (maximum at left) starting from neighboring single-species domains. Parameters: λA = 1.6, λB = 1.29, ζA = 0.005, ζB = 0.16. 5 Conclusion We study a spatial stochastic model of two-species competition, inspired by the spatially uniform system analyzed in Ref. [7]. Based on Janssen’s results of for multicolored directed percolation [14], we find that the flow of parameters under the renormalization group is such that a small competitive advantage may be amplified to the point of excluding the more numerous, but less competitive species, unless the latter reproduces rapidly. Thus the less competitive population can go extinct, regardless of its size as given by mean-field theory. These effects tend to be stronger, the smaller the dimensionality d of the system. Our result shows that the survival of the scarcer phenomenon observed in [7] persists, in some cases in a stronger form, in the presence of spatial structure. These differences are evident when comparing Figs. 1 and 5. Monte Carlo simulations of a one-dimensional lattice model verify the existence of SSS when the more competitive species has a reproduction rate near the critical value, while that of the less competitive one is above criticality. The latter species, although more numerous in the quasi-stationary regime, loses out due to invasion by the more more competitive species. Our results suggest that, in an ecological setting, a species with smaller population density can in fact outcompete one with a higher density, even when the populations are not well mixed, so that the spatiotemporal pattern is one of shifting domains. Under conditions that favor survival of the scarcer, when both species are 16 Figure 10: Steady-state interface velocities for species A (blue) and B (red) versus ζB ; other parameters as in Fig. 9. Figure 11: Quasi-stationary probability distribution for parameters λA = 1.6, λB = 1.26, ζA = 0.005, and ζB = 0.3, system size L = 200. Red corresponds to maximum and violet to zero probability. The maxima along the x and y axes correspond to the single-species QS distributions. present in the same niche, the more competitive species (B) tends to survive longer, even if its density is smaller than the that of the less competitive species (A). In isolation however, species A is less susceptible to extinction than is species B, again considered in isolation. This suggests that, as long as both species are present globally, a cyclical dynamics of the form A → B → 0 → A,..., may be observed in localized regions. That is, a region occupied by species A is susceptible to invasion by B, which in turn is susceptible to extinction, permitting return of species A, and so on. Sequences of this kind are indeed seen in Fig. 7. Given the simplicity of the 17 model and the limitations of the present theoretical and numerical analysis, the above conclusions should be seen as preliminary at best. They may nevertheless suggest directions for further research in ecological modelling. Acknowledgements This work is supported by CNPq, Brazil. References [1] Hardin, G. (1960) Science 131(3409), 1292. [2] Renshaw, E. (1991) Modelling biological populations in space and time, Cambridge University Press, . [3] Murrell, D. J., Purves, D. W., and Law, R. (2001) Trends Ecol Evol 16(10), 529. [4] Condit, R., Ashton, P. S., Baker, P., Bunyavejchewin, S., Gunatilleke, S., Gunatilleke, N., Hubbell, S. P., Foster, R. B., Itoh, A., LaFrankie, J. V., et al. (2000) Science 288(5470), 1414. [5] Silander Jr, J. A. and Pacala, S. W. (1985) Oecologia 66(2), 256. [6] Pacala, S. W. and Deutschman, D. H. (1995) Oikos p. 357. [7] Gabel, A., Meerson, B., and Redner, S. (2013) Phys. Rev. E 87, 010101. [8] Assaf, M. and Meerson, B. (2006) Phys. Rev. E 74(4), 041115. [9] Assaf, M. and Meerson, B. (2010) Phys. Rev. E 81(2), 021116. [10] Doi, M. (1976) J Phys A-Math Gen 9(9), 1479. [11] Peliti, L. (1985) J. Physique 46, 1469. [12] Tauber, U. C., Howard, M., and Vollmayr-Lee, B. P. (2005) J Phys A-Math Gen 38(17), 79. [13] Cardy, J. (1996) Scaling and Renormalization in Statistical Physics, Cambridge University Press, . [14] Janssen, H. (2001) J Stat Phys 103(5), 801. [15] Barabasi, A. L. and Stanley, E. (1995) Fractal Concepts in Surface Growth, Cambridge University Press, . [16] Dennis, G. R., Hope, J. J., and Johnsson, M. T. (2013) Comput Phys Commun 184(1), 201. [17] Dickman, R. (1994) Phys. Rev. E 50(6), 4404. [18] Moro, E. (2004) Phys. Rev. E 70(4), 045102. [19] Dornic, I., Chaté, H., and Muñoz, M. A. (2005) Phys. Rev. Lett. 94, 100601. [20] Hinrichsen, H. (2000) Adv Phys 49, 815. [21] De Oliveira, M. M. and Dickman, R. (2005) Phys. Rev. E 71(1), 016129. 18 [22] Dickman, R. and De Oliveira, M. M. (2005) Physica A 357(1), 134. [23] Grassberger, P. and De La Torre, A. (1979) Ann Phys (NY) 122(2), 373. 19 This figure "f2.png" is available in "png" format from: http://arxiv.org/ps/1304.5956v2