Survival of the Scarcer
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec Survival of the Scarcer
Alan Gabel, Baruch Meerson, and S. Redner Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215, USA Racah Institute of Physics, Hebrew University, Jerusalem 91904, Israel
We investigate extinction dynamics in the paradigmatic model of two competing species A and B that reproduce ( A → A , B → B ), self-regulate by annihilation (2 A →
0, 2 B → A + B → A , A + B → B ). For a finite system that is in the well-mixed limit, a quasi-stationarystate arises which describes coexistence of the two species. Because of discrete noise, both specieseventually become extinct in time that is exponentially long in the quasi-stationary population size.For a sizable range of asymmetries in the growth and competition rates, the paradoxical situationarises in which the numerically disadvantaged species according to the deterministic rate equationssurvives much longer. PACS numbers: 05.10.Gg, 87.23.Cc, 02.50.-r
In the paradigmatic two-species competition model, apopulation is comprised of distinct species A and B , eachof which reproduce and self regulate by intraspecies com-petitive reactions. In addition, interspecies competitivereactions occur, which are deleterious to both species [1].For large, well-mixed populations, the dynamics can beaccurately described by deterministic rate equations. Forfinite systems, however, fluctuations in the numbers of in-dividuals ultimately lead to extinction, in stark contrastto the rate equation predictions.In this work, we investigate how asymmetric inter-species competition influences the extinction probabilityof each species. In a finite ecosystem, extinction arisesnaturally when multiple species compete for the same re-sources. In such an environment, one species often dom-inates, while the others become extinct [2–6], a featurethat embodies the competitive exclusion principle . A re-lated paradigm appears in the context of competing par-asite strains that exploit the same host population, or inthe fixation of a new mutant allele in a haploid popula-tion whose size is not fixed [7].With asymmetric interspecies competition, we uncoverthe surprising feature that deterministic and stochasticeffects, which originate from the same elemental reac-tions, act oppositely. For sizable asymmetry ranges in thegrowth and competition rates, the situation arises wherethe combined effects of the elemental reactions leads toone species being numerically disadvantaged at the mean-field level, despite its interspecies competitive advantage,but this competitive advantage dominates the other re-action processes at the level of large deviations. Thusthe outcompeted and less abundant species has a higherlong-term survival probability: “survival of the scarcer”. Model:
Asymmetric competition of two species A and Bis defined by the reactions: A −→ A + A B g −→ B + B ,A + A /K −→ B + B /K −→ , (1) A + B ǫ/K −→ B A + B αǫ/K −→ A. The first line accounts for reproduction, the second for in-traspecies competition, and the last for interspecies com-petition. Here K is the environmental carrying capacity,which sets the size of the overall population, ǫ quantifiesthe severity of the competition, while g and α quantifythe asymmetries in the growth and interspecies compe-tition rates, respectively. In our presentation, we focuson the limit K ≫
1. While a general model should alsocontain asymmetry in the intraspecies competition rate,no new phenomena arise by this generalization; for sim-plicity, we study the model defined by Eqs. (1).To probe extinction in two-species competition, we fo-cus on P m,n ( t ), the probability that the population con-sists of m ≥ n ≥ t . In the limitof a perfectly-mixed population, the stochastic reactionprocesses in (1) lead to P m,n ( t ) evolving by the masterequation˙ P m,n ( t ) = ˆ HP m,n = (cid:2)(cid:0) E − − (cid:1) m + g (cid:0) F − − (cid:1) n (cid:3) P m,n + h(cid:0) E − (cid:1) m ( m − K + (cid:0) F − (cid:1) n ( n − K i P m,n + h ǫK ( E −
1) + αǫK ( F − i mnP m,n . (2)Here E and F are the raising and lowering operators [8]for species A and B, respectively; viz. E i P m,n = P m + i,n and F j P m,n = P m,n + j . Deterministic Rate Equations:
First we focus on the av-erage population sizes h m i = P m,n mP m,n and h n i = P m,n nP m,n . From (2), the evolution of these quantitiesis given by ˙ h m i = h m i (cid:18) − h m i K − ǫ h n i K (cid:19) , ˙ h n i = h n i (cid:18) g − h n i K − αǫ h m i K (cid:19) . (3)Here we neglect correction terms of the order of 1 /K and, more importantly, neglect correlations by assumingthat h m i = h m i , h n i = h n i , and h mn i = h m ih n i .We restrict ourselves to the parameter range αǫ < g < /ǫ , which guarantees that the fixed point correspondingto coexistence of both species is stable. The four fixedpoints of the rate equations (3) are then:( m ∗ , n ∗ ) = (0 ,
0) unstable node , = ( K, , (0 , Kg ) saddles , (4)= (cid:16) K − gǫ − αǫ , K g − αǫ − αǫ (cid:17) stable node . If the initial populations of both species are non-zero,they are quickly driven to the stable node (Fig. 1) thatdescribes the steady-state populations in the mean-fieldlimit. The relaxation time toward the stable node, τ r , isindependent of K . These steady-state populations of thetwo species are equal when g ∗ = 1 + αǫ ǫ . (5)For g < g ∗ , the B -population is scarcer. Naively, thescarcer population should be more likely become extinctfirst. However, as we shall show, a proper account of thefluctuations that stem from the underlying elemental re-actions themselves leads to a radically different outcome. FIG. 1. Schematic flow diagram in asymmetric two-speciescompetition for weak competition in the mean field. The un-stable node, the saddles, and the stable node are shown asopen, hatched, and solid, respectively.
Extinction:
The mean-field picture is incomplete becausefluctuations of the population sizes about their fixed-point values are ignored. For large populations (corre-sponding to carrying capacity K ≫ quasi-stationary state where the two species coexist. Thisstate is stable in the mean-field description (heavy dotin Fig. 1). However, an unlikely sequence of deleteriousevents eventually occurs that ultimately leads one popu-lation, and then the other, to extinction. After the firstextinction, the remaining population settles into anotherquasi-stationary state around one of the single-speciesfixed points ( m ∗ , n ∗ ) = ( K,
0) or (0 , Kg ). Eventually alarge fluctuation drives the remaining species to extinc-tion. This second extinction time is typically much longer[by a factor that scales as exp(const. × K )] than the firsttime, because the remaining species does not suffer inter-species competition. Once a species is extinct, there is no possibility of recovery since there is no replenishmentmechanism.The question that we address is: which species typi-cally goes extinct first? The answer is encoded in thedynamics of the two-species probability P m,n ( t ). Duringthe initial relaxation stage, a quasi-stationary probabilitydistribution is quickly reached (Fig. 2). The probabilitydistribution is sharply peaked at the stable fixed pointof the mean-field theory. This probability slowly “leaks”into localized regions near each of the single-species fixedpoints ( m ∗ , n ∗ ) = ( K,
0) and (0 , Kg ). Thus two sharply-peaked single-species distributions start to form. If the( K,
0) peak grows faster then the B species is more likelyto go extinct first. Similarly, a faster growing (0 , Kg )peak means A extinction is more likely. Eventually, theprobability distribution that is localized at one of the twosingle-species fixed points slowly leaks toward the fixedpoint (0 ,
0) that corresponds to complete extinction [9].To determine extinction rates, it is helpful to define P A , P B , and P φ as the respective probabilities thatspecies A is extinct, species B is extinct, or neither isextinct at time t [10]. (Being interested in times muchshorter than the expected extinction time of both species,we can neglect the probability of the latter process.) Bydefinition, these extinction probabilities are P A = X n> P ,n , P B = X m> P m, , P φ = X m,n> P m,n ; (6)these satisfy P A + P B + P φ = 1, up to an exponen-tially small correction that stems from the process whereboth species become extinct simultaneously. In the limit K ≫ τ r , the sums in Eqs. (6) are dominated bycontributions from values of m and n that are close tothe single-species and coexistence fixed points. Moreover,these extinction probabilities evolve according to a set ofeffective coupled equations˙ P A = R A P φ , ˙ P B = R B P φ , ˙ P φ = − ( R A + R B ) P φ , (7)that define R A and R B as the respective extinction ratesfor species A and species B. Solving these equations yieldsthe time dependence of the extinction probabilities P A ( t ) = R A R (cid:0) − e −R t (cid:1) , P B ( t ) = R B R (cid:0) − e −R t (cid:1) , (8)with R = R A + R B . To determine R A and R B , we followthe evolution of the eigenstate of the master equation(2) that determines the leakage of probability from thevicinity of the coexistence point: P m,n ( t ) = Π m,n e −R t , m, n > , (9)where ˆ H Π m,n = −R Π m,n , m, n > , (10)and R is the third -lowest positive non-trivial eigenvalueof the operator ˆ H . The two still-smaller positive non-trivial eigenvalues correspond to the much slower decayof quasi-stationary single-species states and play no rolein the dynamics of the first extinction event. There is alsoa trivial eigenvalue that corresponds to the final state ofcomplete extinction.Combining Eq. (2) with (6)–(9), we obtain the follow-ing expression for the extinction rate of the A species: R A = 1 K X n> (cid:0) ǫ n Π ,n + Π ,n (cid:1) . (11a)As expected, the extinction rate for As involves two pro-cesses: (i) elimination of the last remaining A via compe-tition with Bs and (ii) annihilation of the last remainingpair of As. Similarly, R B = 1 K X m> (cid:0) αǫ m Π m, + Π m, (cid:1) . (11b)To calculate R A and R B , we therefore need to evaluatethe small-population-size tails of Π m,n . This task canbe achieved by applying a variant of Wentzel-Kramers-Brillouin (WKB) approximation, that was pioneered inRefs. [11–14], and was applied more recently to pop-ulation extinction, in particular, for stochastic two-population systems [10, 15–22]. The WKB ansatz forΠ m,n has the form Π m,n = e − KS ( x,y ) , (12)where x = m/K and y = n/K are treated as continuousvariables. Substituting Eq. (12) into (11a) and assuming K ≫
1, gives, to lowest order in 1 /KR A ∼ e − KS A , R B ∼ e − KS B , (13)where S A = S (0 , g ) and S B = S (1 , K ≫
1, the eventual extinction probabili-ties in (8) simply become (up to pre-exponential factorsthat depend on K ), as P A ( t = ∞ ) = 1 −P B ( t = ∞ ) ≃ e − KS A e − KS A + e − KS B . (14)To determine the extinction probabilities explicitly, wetherefore need S A and S B . To this end, we substitute theWKB ansatz (12) in Eq. (10) and Taylor expand S ( x, y )to lowest order in 1 /K . After some algebra, we obtain aneffective Hamilton-Jacobi equation H ( x, y, ∂ x S, ∂ y S ) = −R , with the Hamiltonian H ( x, y, p x , p y ) = x ( e p x −
1) + gy ( e p y − x e − p x − y e − p y − ǫxy ( e − p x − αǫxy ( e − p y − . (15) −8 −6 −4 −2 Population Size P r obab ili t y WKBB speciesA species
FIG. 2. Quasi-stationary probability distributions for speciesA, P m = P n P m,n , and species B, P n = P m P m,n . Parame-ters are K = 100, ǫ = 0 . g = 0 .
45, and α = 0. Symbols aresimulation results, while the solid curve is WKB approxima-tion for the B species distribution. Here p x = ∂ x S and p y = ∂ y S are the canonical momentathat are conjugate to the “coordinates” x and y . Corre-spondingly, S ( x, y ) is the classical action of the system.Since we expect that −R (which now has the meaningof energy in this Hamiltonian system) is expected to beexponentially small, we set it to zero. The Hamiltonianequations of motion ˙ x = ∂H/∂p x , ˙ p x = − ∂H/∂x , etc.,have six finite zero-energy fixed points, and three morefixed points where one or both momenta are at minusinfinity. Only three of the fixed points, however, turnout to be relevant for answering our question about whichspecies typically goes extinct first. These are: F φ = (cid:18) − gǫ − αǫ , g − αǫ − αǫ , , (cid:19) ,F A = (0 , g, ln( gǫ ) , ,F B = (1 , , , ln( αǫ/g )) . (16)A straightforward way to determine S A and S B is bycalculating the action along the activation trajectories .These are zero-energy, but non-zero-momentum trajec-tories of the Hamiltonian system (15) that go from F φ to F A , and from F φ to F B , respectively. These actions are S A = Z F A F φ ( p x dx + p y dy ) , (17)and similarly for S B . In general, these activationtrajectories—separatrices, or instantons—cannot be cal-culated analytically because of the lack of an integral ofmotion that is independent of the energy. However, for ǫ ≪ ǫ = 0, which corresponds totwo uncoupled species. Here the zero-energy activationtrajectories can be easily found. For the F φ → F A sep-aratrix, the B species is unaffected by A extinction so( y, p y ), which correspond to the coordinates of the Bs,remains constant throughout the evolution. As a result,a parametric form of the F φ → F A separatrix is x = 2 e p x e p x +1 , (18a)with y = g and p y = 0 throughout. Similarly, for the F φ → F B separatrix one obtains y = 2 ge p y e p y +1 , (18b)with x = 1 and p x = 0 throughout. Substituting thetrajectories given in (18) into (17) and performing theintegration by parts gives S A = 2(1 − ln 2) and S B =2 g (1 − ln 2) [23, 24].For weak interspecies competition ( ǫ ≪ ǫ . For this purpose, we split the Hamiltonian (15)into unperturbed and perturbed parts, H = H + ǫH ,and similarly expand the action as S = S + ǫS + . . . .Following [15, 18, 25], the correction to the action is S = R ∞−∞ H [ x ( t ) , y ( t ) , p x ( t ) , p y ( t )] dt , where the integralis evaluated along the unperturbed trajectories given byEqs. (18). Performing this integral for S yields the cor-rected actions S A = 2(1 − ln 2) − ǫ (2 g ln 2) ,S B = 2 g (1 − ln 2) − ǫ (2 α ln 2) . (19)Equations (13) and (19) give the analytic expressionfor the extinction probabilities of each species for weakinterspecies competition. Using Eq. (19) and imposingthe condition S A = S B from Eq. (14), we obtain thefollowing condition for equal extinction probability forboth species: g = 1 + α c ǫ c ǫ , (20)where c = (cid:2) (ln 2) − − (cid:3) − . The predictions of Eqs. (14)and (19) are in good agreement with our simulation re-sults (Fig. 3). Phase Diagram:
Comparing Eqs. (5) and (20), one seesthat there is a sizable region in the α - g parameter spacewhere one species has a smaller quasi-stationary popu-lation and yet an (exponentially) smaller probability tofirst become extinct. As an illustration, Fig. 3 shows theprobability for B to become extinct first for fixed α and ǫ . We also produced analogous curves as Fig. 3 at manyvalues of α . From the value of g at which the extinctionprobabilities are equal, we infer the phase diagram shownin Fig. 4. Simulations at larger values of ǫ yield the samequalitative phase diagram. g P B A dominant B dominant
FIG. 3. Probability that species B first becomes extinct asa function of growth rate asymmetry g for α = 0, ǫ = 0 . K = 40. The curve is the prediction of Eq. (14), with actionsgiven by (19), while circles are simulation results. In thehatched region, the quasi-stationary population of B is lessthan that of A, but Bs are less likely to become extinct first. α g B dominant A dominant
FIG. 4. Phase diagram for ǫ = 0 . Conclusion:
In two-species competition, interspeciescompetitive asymmetry leads to the unexpected phe-nomenon of “survival of the scarcer”. The very sameelemental reactions that lead to a disadvantage in thequasi-stationary population size of one species within adeterministic mean-field theory, may also give this speciesa great advantage in its long-term survival when fluctu-ation effects are properly accounted for.AG and SR gratefully acknowledge NSF grant DMR-0906504 and DMR-1205797 for partial financial sup-port. BM was partially supported by the Israel ScienceFoundation (Grant No. 408/08), by the US-Israel Bina-tional Science Foundation (Grant No. 2008075), and bythe Condensed Matter Theory Visitors Program in theBoston University Physics Department. [1] J. D. Murray,
Mathematical Biology I: An Introduction,3rd ed. (Springer, New York, 2001).[2] G. Hardin, Science, , 1292 (1960).[3] J. C. Flores, J. Theor. Biol. , 295 (1988).[4] J. Bengtsson, Nature , 713 (1989).[5] A. E. Noble and W. F. Fagan, arXiv:1102.0052v1.[6] D. Gravel, F. Guichard, and M. E. Hochberg, Ecol. Lett. , 828 (2011).[7] T. L. Parsons and C. Quince, Theor. Popul. Biol. , 121(2007).[8] N. G. Van Kampen, Stochastic Processes in Physics andChemistry, nd ed. (North-Holland, Amsterdam, 2001).[9] See Supplemental Material at http://physics.bu.edu/~redner/pubs/pdf/weaker-sm.pdf ,which contains a movie that shows the time evolution of P m,n ( t ).[10] O. Gottesman and B. Meerson, Phys. Rev. E, , 021140(2012).[11] R. Kubo, K. Matsuo, and K. Kitahara, J. Stat. Phys. ,51 (1973).[12] G. Hu, Phys. Rev. A , 5782 (1987).[13] C. S. Peters, M. Mangel, and R. F. Costantino, Bull. Math. Biol. , 625 (1989).[14] M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J.Chem. Phys. , 5735 (1994).[15] M. I. Dykman, I. B. Schwartz, and A. S. Landsman,Phys. Rev. Lett. , 078101 (2008).[16] A. Kamenev and B. Meerson, Phys. Rev. E , 061107(2008).[17] B. Meerson and P. V. Sasorov, Phys. Rev. E , 041130(2009).[18] M. Khasin and M. I. Dykman, Phys. Rev. Lett. ,068101 (2009).[19] M. Khasin, M. I. Dykman, and B. Meerson, Phys. Rev.E , 051925 (2010).[20] I. Lohmar and B. Meerson, Phys. Rev. E , 051901(2011).[21] A. J. Black and A. J. McKane, J. Stat. Phys. P12006(2011).[22] M. Khasin, B. Meerson, E. Khain, and L. M. Sander,Phys. Rev. Lett. , 138104 (2012).[23] V. Elgart and A. Kamenev, Phys. Rev. E , 041106(2004).[24] D. A. Kessler and N. M. Shnerb, J. Stat. Phys. , 5(2007).[25] M. Assaf, A. Kamenev, and B. Meerson, Phys. Rev. E,78