Table Of ContentMultiple extinction routes in stochastic population models
Omer Gottesman1,2 and Baruch Meerson2
1Weizmann Institute of Science, Rehovot 76100, Israel and
2Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel
Isolated populations ultimately go extinct because of the intrinsic noise of elementary processes.
In multi-population systems extinction of a population may occur via more than one route. We
investigate this generic situation in a simple predator-prey (or infected-susceptible) model. The
predatorandpreypopulationsmaycoexistforalongtimebutultimatelybothgoextinct. Inthefirst
extinctionroutethepredatorsgoextinctfirst,whereasthepreythriveforalongtimeandthenalso
goextinct. Inthesecondroutethepreygoextinctfirstcausingarapidextinctionofthepredators.
Assuminglargesub-populationsizesinthecoexistencestate,wecomparetheprobabilitiesofeachof
2 thetwoextinctionroutesandpredictthemostlikelypathofthesub-populationstoextinction. We
1 alsosuggest aneffectivethree-statemasterequationfortheprobabilities toobservethecoexistence
0 state, thepredator-free state and theempty state.
2
n PACSnumbers: 87.18.Tt,87.23.Cc,02.50.Ga,05.40.Ca
a
J
I. INTRODUCTION II. MODEL AND MEAN-FIELD DYNAMICS
9
2
We assume that predators (F - foxes) and prey (R -
] Isolated populations ultimately go extinct with prob- rabbits)arewellmixedinspacesothatspatialdegreesof
E ability one, even in the absence of adverse environmen-
freedom are irrelevant. The rabbits reproduce naturally,
P tal variations, because of the presence of an absorbing
whereasthefoxesonlyreproducebypredation. Thefoxes
o. state atzeropopulationsize andthe intrinsic noiseofel- and rabbits die or leave with constant per-capita rates,
i ementary processes. Populationextinction risk is an im-
b ingeneraldifferentforeachpopulation. Wealsoaccount,
portant negative factor in viability of small populations
- in a model form, for the competition of rabbits over re-
q [1, 2], whereas extinction of an endemic disease from a
sources, by adding the elementary process 2R R that
[ community [1] is of course a good thing. For decades, →
becomesmoreeffectiveatlargepopulationsizes. Bytak-
quantitativeanalysisofpopulationextinction, causedby
2 ing into account the competition, this model generalizes
v intrinsicandextrinsicnoises,hasbeeninthe focus ofat- the classical Lotka-Volterra model [25]. The elementary
1 tention from bio-mathematicians, see e.g. Ref. [2] for a
processes and their rates are described in Table I. We
3 review. More recently, it has also attracted much atten-
chose the units of time so that the per-capita death rate
3 tion from physicists [3–24], as a highly relevant example
of the foxes is equal to 1. The large parameter N de-
4 of a rare large fluctuation far from thermal equilibrium.
. termines the typical scaling of the sub-population sizes,
2
Inmulti-populationsystemsextinctionofapopulation seebelow,andwewillassumeastronginequalityN 1
1 ≫
1 may occur via more than one route. In this paper we throughout the paper.
1 analyze this generic situation in a simple predator-prey
: model where the predators need the prey for survival.
v Process Type of transition Rate
Thesamemodelcandescribespreadofaninfectious dis-
i Reproduction of rabbits R→2R a
X ease in an isolated community. Although the predator
Predation/reproduction of foxes F +R→2F 1/(ΓN)
r and prey populations (or infected and susceptible pop-
a ulations) may coexist for a long time, they ultimately Death of rabbits R→0 b
bothundergoextinction,andthishappensviaoneoftwo Death of foxes F →0 1
routes. In the first route, that can be called sequential, Competition among rabbits 2R→R 1/N
the predators go extinct first, whereas the prey thrive
for a long time and then also go extinct. In the second
TABLE I:Predator-prey model
route, that can be called (almost) parallel, the prey go
extinct first, causing a rapid extinction of the predators.
The mean-field equations for this model are
Assuming large sub-population sizes in the coexistence
state,weapplyaWKB(afterWentzel,KramersandBril- 1 1
louin) approximation to the pertinent master equation. R˙ = (a b)R RF R2,
− − ΓN − 2N
In this way we evaluate the probabilities of each of the
1
twoextinctionroutes andpredictthe mostlikelypathof F˙ = RF F. (1)
ΓN −
the sub-populations on the way to extinction. We also
suggest an effective three-state master equation for the The same model can be reinterpreted to describe the
evolution of the probabilities to observe the coexistence spreadofaninfectious diseasein anisolated community.
state, the predator-free state and the empty state. Here we re-interpret the rabbits as the susceptibles (S)
2
andthefoxesastheinfected(I).AsintheconventionalSI 1
modelwithpopulationturnover[11,16,26–28],asuscep- HaL
tible individual can become infected upon contact with
another infected, while infected individuals are removed
(recover with immunity, leave or die) with a constant
0.5
per-capita rate. In the conventional SI model the sus- y
ceptibles arrive from outside. In the modified SI model,
M
presented here, there are no arrivals from outside, but 3
the susceptibles can reproduce by giving birth. They M
0 1
are also removed (die or leave), as in the conventional M
2
SI model. Finally, they compete for resources, 2S S, 0 0.5 1 1.5 2 2.5
→
so their population size remains bounded. See the ele-
x
mentary processes and their rates in Table II, where we
1.5
measure time in the units of per-capita removal rate of HbL
the infected. The mean-field equations for the modified
SI model are
1
1 1
S˙ = (a b)S SI S2, y
− − ΓN − 2N M
0.5 3
1
I˙ = SI I. (2)
ΓN −
M
These of course coincide with Eqs. (1) upon the change 0 M1 2
of S to R and I to F. 0 0.5 1 1.5 2 2.5
x
Process Typeof transition Rate
Reproduction of susceptibles S →2S a FIG. 1: (Color online) The mean field phase trajectories of
Infection I+S →2S 1/(ΓN) thesystem for a=2, b=1 andtwo different valuesof Γ. (a)
Γ = 1.8: M3 is a stable node. (b) Γ = 0.4: M3 is a stable
Removalof susceptibles S →0 b
focus.
Removalof infected I →0 1
Competition among susceptibles 2S →S 1/N
a stable node for Γ > Γ = 4(√1+a b 1), and a
0
− −
TABLE II: Susceptible-Infected model for an isolated com- stable focus for Γ < Γ0. Note that Γ0 < Γ∗ for any
a > b. Note also that y¯ is a non-monotonic function
munity 3
of Γ. It vanishes at Γ = 0 and Γ = Γ∗, and reaches
a maximum, Γ2∗/8, at Γ = Γ∗/2 corresponding to the
For concreteness, we will use the predator-prey nota-
optimal predation rate. Figure 1 shows two examples
tioninthefollowing. Introducingtherescaledpopulation
of mean-field trajectories: for Γ0 < Γ < Γ∗ (a) and for
sizes x = R/N and y = F/N, we can rewrite the mean-
Γ < Γ (b). The characteristic relaxation time scale t
0 r
field equations (1) as
of the coexistence state is determined by the realpartof
the eigenvalues of the linear stability matrix at M .
xy x2 xy 3
x˙ =(a b)x , y˙ = y. (3) Let us compare the mean-fielddynamics of this model
− − Γ − 2 Γ −
with those of the classicalLotka-Volterramodel[25]and
The mean-field equations are fully characterized by two of the SI model with population turnover [11, 16, 26].
parameters, a b and Γ. We will be only interested in The competition among the rabbits eliminates one un-
−
the case of a b > 0 and Γ < Γ∗ = 2(a b), when desirable feature of the Lotka-Volterra model: the un-
− −
Eqs. (3) have three fixed points corresponding to non- limited proliferation of rabbits in the absence of foxes.
negative population sizes. The fixed point M (x¯ = Furthermore, there are no neutral cycles in this model,
1 1
0, y¯ = 0) describes an empty system. It is a saddle in contrast to the Lotka-Volterra model. The attracting
1
point: attracting in the y-direction (when there are no fixed point M , describing a unique coexistence state of
3
rabbits),andrepellinginthex-direction. Thefixedpoint the rabbits and foxes, appears instead, with either oscil-
M2 [x¯2 =Γ∗,y¯2 =0]describesanestablishedpopulation latory (underdamped) or non-oscillatory (overdamped)
of rabbits in the absence of foxes. It is also a saddle character of phase trajectories approaching it. With ac-
point: attracting in the x-direction (when there are no count of the discreteness of rabbits andfoxes, and of the
foxes), and repelling in a direction corresponding to the stochastic character of their interaction, this leads to an
introduction of a small number of foxes into the system. exponentially long life-time of the coexistence state, in
The third fixed point M3 [x¯3 = Γ, y¯3 = Γ(Γ∗ Γ)/2] contrast to a much shorter, power-law life-time, charac-
−
is attracting and describes the coexistence state. It is teristic of the classical Lotka-Volterra model [7, 15].
3
On the other hand, the mean-field dynamics (2) is
300
quite similarto thatofthe conventionalSI orSIRmodel
with population turnover [11, 16, 26], except that the 250
fixed point M (0,0), that is absent in the conventional
1
SImodel,nowappears. Asthisfixedpointisrepellingin F 200
d
the x-direction, its presence does not make much differ- an
ence in the mean-field descriptionof the coexistence (or, R 150 (a)
in the context of the SI model, of the endemic state of 100
the disease in the population).
50
0
0 5000 10000 15000
t
120
III. STOCHASTIC ANALYSIS 100
F 80
d
n
A. Two extinction routes a
R 60 (b)
40
Thesituation,however,changesconsiderablywhenone
20
accounts for the stochasticity. Here the empty system
is an absorbing state describing extinction of both sub- 0
populations. Having reached this state via a rare large 0 1000 2000 3000 4000
fluctuation, the system will stay there forever. It is this t
empty state that representsthe only true steady state of FIG. 2: (Color online) Two different stochastic realizations
this system. (This is different from the conventional SI of the model for a = 2, b = 1, Γ = 0.4 and N = 120.
model, where the true steady state describes an estab- Shown are the population sizes of rabbits (solid lines) and
lished infection-free population.) Importantly, the ab- foxes (dashed lines) versus time (measured in Monte Carlo
sorbingstateatzerocantypicallybe reached,atN 1, steps). (a) Thefoxesgo extinctfirst,whereastherabbitsget
via one of two routes. The first route can be called≫se- established around the fixed point M2 and, after a very long
time,also go extinct(notshown). (b)Therabbitsgo extinct
quential: The foxes go extinct first, whereas the rabbits
first, causing a rapid extinction of thefoxes.
thrive for a long time, and then also go extinct. The
presence of this route is caused by the fact that all fox-
free states represent absorbing states for the foxes. The B. Master equation and long-time dynamics
second route is (almost) parallel: the rabbits go extinct
first causing a rapid (almost deterministic) extinction of
Let P (t) be the probability to observe, at time t,
m,n
the foxes.
m rabbits and n foxes, where m,n = 0,1,2.... The
The two different routes to extinction are clearly ob- evolution of Pm,n(t) is described by the master equation
served from Figure 2a and b that shows two different
realizationsof the stochastic processes,listed in Table 1, P˙m,n = HˆPm,n =a[(m−1)Pm−1,n−mPm,n]
for the same set of rate coefficients and for the same ini- + (1/ΓN)[(m+1)(n 1)Pm+1,n−1 mnPm,n]
− −
tial conditions. The parameters are the same as in Fig. + b[(m+1)P mP ]
m+1,n m,n
−
1b. In figure a the foxes go extinct first, whereas the
+ (n+1)P nP
m,n+1 m,n
rabbitpopulationgetsestablished. Thentherabbitswill −
+ (1/2N)[(m+1)mP m(m 1)P ],(4)
also go extinct after a very long time (not shown in the m+1,n− − m,n
figure). In figure b the rabbits go extinct first causing a whereP =0whenanyoftheindicesisnegative. There
i,j
rapid extinction of the foxes. Figure 3a and b shows the arethreequantities ofinterestheredescribingextinction
same trajectories in the R,F plane. One can notice in of the sub-populations involved. The probability of ex-
Figures 2b and 3b that the population of rabbits is very tinction, by time t, of foxes at any number of rabbits is
small close to the time when the foxes go extinct. This equal to ∞ P (t). The evolution of this quantity
m=0 m,0
feature,reproducibleinmanystochasticsimulations,will in time isPdescribed by the equation
obtainanaturalexplanationintheWKBtheory,seesub-
∞ ∞
section C3. More generally,the WKB theory will enable d
P (t)= P (t), (5)
onetocomparethe probabilityofeachofthetwoextinc- dt m,0 m,1
tionroutes,topredictthemostlikelypathoftherabbits mX=0 mX=0
and foxes to extinction, and to evaluate the mean time directly following from Eq. (4). The right-hand side of
to extinction of each sub-population. Eq. (5) is simply the probability of death of the last re-
4
M ofthemean-fieldtheory. Finally,theextinctionprob-
160 3
ability of both sub-populations P (t) corresponds to a
0,0
140 (a) Kronecker-delta probability density. Not only the struc-
120 ture, but the long-time dynamics of these three distribu-
100 tions aredifferent. To clearlysee this point, let us define
F the total “probability contents” of the vicinities of each
80
of the fixed points M , M and M :
1 2 3
60
40 1(t) = P0,0(t), (8)
P
∞
20
(t) = P (t), (9)
0 P2 m,0
0 50 100 150 200 250 300 mX=1
∞ ∞
R
120 3(t) = Pm,n(t). (10)
P
(b) mX=1nX=1
100
At N 1 and t t the sums in Eqs. (9) and (10) are
r
80 ≫ ≫
mostly contributed to by close vicinities of M and M ,
2 3
F respectively. The long-time evolution of , and
60 P1 P2 P3
is described by an effective three-state master equation:
40
˙ (t) = r (t)+r (t),
1 21 2 31 3
20 P P P
˙ (t) = r (t)+r (t),
2 21 2 32 3
P − P P
0 ˙ (t) = (r +r ) (t), (11)
0 20 40 60 80 100 120 P3 − 31 32 P3
R where r is the (yet unknown) transition rate from the
ij
FIG. 3: (Color online) Same as in Fig. 2, but on the R,F vicinity of the fixed point i to the vicinity of the fixed
plane. pointj. Letthe initialconditioncorrespondtothe coex-
istence state around M :
3
maining fox at any number of rabbits. [ (0), (0), (0)]=(0,0,1). (12)
1 2 3
∞ P P P
In its turn, P (t) is the probability of extinc-
n=0 0,n Then the solution of Eqs. (11) is
tion, by time tP, of rabbits at any number of foxes. This
quantity evolves in time according to r e−r21t+(r r )e−(r31+r32)t
32 31 21
(t) = 1+ − ,(13)
∞ ∞ ∞ P1 r r r
d 1 21 31 32
− −
dt P0,n(t)= ΓN nP1,n(t)+b P1,n(t). (6) r32[e−(r31+r32)t e−r21t]
nX=0 nX=1 nX=0 2(t) = − , (14)
P r r r
21 31 32
− −
The first and second terms on the right are the prob- (t) = e−(r31+r32)t. (15)
3
abilities of predation and of natural death of the last P
remaining rabbit, respectively. Finally, the probability Once the transition rates r31, r32 and r21 are known,
P (t) of extinction of both rabbits and foxes by time t Eqs. (13)-(15) provide a valuable “coarse-grained” de-
0,0
is described by the equation scription of the stochastic system in terms of the long-
time evolution of the probabilities to observe the coexis-
dP0,0 =bP (t)+P (t). (7) tencestatearoundM3,thefox-freestatearoundM2 and
dt 1,0 0,1 the extinction state at M1. The coexistence probabil-
ity (t) goes down to zero exponentially in time. The
3
AtN 1,andattimes muchlongerthanthe relaxation P
≫ probability 1(t) of extinction of both sub-populations
time t of the mean-field theory, the quantities that ap- P
r increases monotonically with time, exhibiting two differ-
pear in Eqs. (5)-(7) become sharply peaked around the
ent exponents, and ultimately reaches 1. Finally, the
correspondingfixed points of the mean-field theory. The
probabilityof the fox-freestate (t) first increaseswith
2
structure of these peaks is affected by the presence or P
time,reachesamaximumandthengoesdowntozero. At
absence of absorbing states. At long times, the proba-
intermediate times, t t min[1/r ,1/(r +r )],
r 21 31 32
bility distribution P (t) with m > 0 (which describes ≪ ≪
m,0 Eqs. (13)-(15) read
thedynamicsofrabbitsconditionalonpriorextinctionof
thefoxes)isaone-dimensionaldistributionpeakedatthe (t) r t, (t) r t,
1 31 2 32
P ≃ P ≃
fixed point M ofthe mean-field theory. The probability
2 (t) 1 (r +r )t. (16)
3 31 32
distribution P (t) with m,n > 0 (which describes the P ≃ −
m,n
long-timedynamicsofthecoexistingrabbitsandfoxes)is During this stage of the dynamics, (t) and (t) grow
1 2
P P
a two-dimensionaldistribution peaked at the fixed point linearly with time, whereas the transition 2 1 does
→
5
not show up. What happens at longer times? The tran- 2R, R 0 and 2R R, where the predators are absent
→ →
sition rates are usually widely different in magnitude. [18]. The eigenvalue r +r is also exponentially small
31 32
According to our WKB calculations, presented below, in N, see below. It corresponds to the quasi-stationary
r r + r . Then, at times t 1/(r + r ), distribution π at n > 0, sharply peaked at the fixed
21 31 32 31 32 m,n
≪ ≫
Eqs. (13)-(15) become point M . Our main effort in the following will be to
3
evaluate, with exponential accuracy, the transition rates
(t) 1 r32e−r21t, (17) r31 and r32 that contribute to the eigenvalue r31 +r32.
1
P ≃ − r +r As all the effective transition rates r are exponentially
31 32 ij
r e−r21t small in N 1, we can drop the right-hand side in the
2(t) 32 , (18) eigenvalue p≫roblem
P ≃ r +r
31 32
P3(t) ≃ 0. (19) Hˆπm,n =−(r31+r32)πm,n, n>0 (23)
Nowthereisaprobabilitycurrentfromstate2tostate1, and arrive at the quasi-stationary equation
describing extinction of rabbits in the absence of foxes,
but the transition rates r31 and r32 are still present in Hˆπm,n ≃0, m>0, n>0. (24)
the equations.
It is common, see e.g. Ref. [2], to characterize the
C. WKB approximation
extinction risk of a stochastic population in terms of its
mean time to extinction. From Eqs. (13)-(15), the mean
time to extinction of foxes is 1. General
∞ ∞
τ = dtt[P˙ (t)+P˙ (t)]= dttP˙ (t) For N 1, Eq. (24) can be approximately solved via
F Z0 1 2 −Z0 3 the WKB≫ansatz
1
= , (20)
π =exp[ NS(x,y)], (25)
r +r m,n
31 32 −
whereas the mean time to extinction of both sub- whereSisassumedtobeasmoothfunctionofthecontin-
uousvariablesx=m/N andy =n/N. The useofWKB
populations is
approximation for finding stationary or quasi-stationary
∞
r +r solutions of master equations with a discrete state space
τ = dttP˙ (t)= 21 32 . (21)
RF Z 1 r (r +r ) was pioneered in Refs. [29–32]. By now this approxi-
0 21 31 32
mation, in the space of population sizes or in the mo-
How is this dynamics encoded in the spectral proper- mentum space, has become a standard tool in the anal-
ties of the the linear operator Hˆ, see Eq. (4)? Let λi be ysis of rare large fluctuations of stochastic populations
theeigenvaluesandπ(i) theeigenstatesofHˆ. Theseare [2, 4, 8, 10, 11, 13, 14, 16, 17, 20–24, 33–37]. Plugging
m,n
defined by the relations Eq. (25) into Eq. (24) and Taylor expanding S to first
orderaround(x,y),wearriveatazero-energyHamilton-
Hˆπm(i),n =−λiπm(i),n, i=1,2,... , Jacobi equation H(x,y,∂xS,∂yS) = 0, with the Hamil-
tonian
so that
∞ H(x,y,px,py)=ax(epx −1)+bx(e−px −1)
Pm,n(t)= Ciπm(i),ne−λit, +xy(epy−px 1)+y(e−py 1)+ x2(e−px 1).(26)
Xi=1 Γ − − 2 −
wheretheconstantsC aredeterminedbytheinitialcon- The trajectories are given by the Hamilton’s equations
i
dition P (0). Let us order the eigenvalues so that for the “coordinates”x andy and conjugate “momenta”
m,n
λ1 < λ2 < λ3 < .... The true (empty) steady state px and py:
correspondstothe onlyzeroeigenvalue: λ =0,whereas
1 x2 xy
πm(1,)n = δm0δn0. Equations (13)-(15) and the inequality x˙ = axepx bx+ e−px epy−px,
−(cid:18) 2 (cid:19) − Γ
r <r +r implythat,atN 1,thenexttwoeigen-
v2a1lues a3r1e λ232=r21 andλ3 =r31≫+r32. The quantity r21 y˙ = xyepy−px ye−py,
was found in Refs. [14, 18], and it is exponentially small Γ −
y
in N: p˙ = (b+x)(1 e−px)+ 1 epy−px +a(1 epx),
x
− Γ − −
x (cid:0) (cid:1)
bN (a b)2 b p˙ = 1 e−py + 1 epy−px , (27)
r21 ≃r π −a exp(cid:20)−2N(cid:18)a−b+b lna(cid:19)(cid:21). y − Γ(cid:0) − (cid:1)
(22) where we are only interested in zero-energy trajectories.
The corresponding eigenvector π δ is the quasi- The (zero-energy)invarianthyperplane p =p =0 cor-
m,0 n0 x y
stationarydistributionoftheone-population modelR responds to the mean-field dynamics (3). The invariant
→
6
hyperplanes x = 0 and y = 0 correspond to the rabbit- The one-dimensional quasi-stationary distribution
free and fox-free dynamics, respectively. The Hamilto- π δ , corresponding to the long-lived population of
m,0 n,0
nian problem is characterized by three independent pa- rabbits in the absence of foxes, corresponds to another
rameters a, b and Γ. instanton-like trajectory: the one connecting the fixed
The Hamiltonian flow (27) has five zero-energy fixed points M and F [14, 18]. The accumulated action
2 1
points:
F1 b
M1 = (0,0,0,0), S21 =Z pxdx=2(cid:18)a−b+blna(cid:19) (30)
M2
M2 = (Γ∗,0,0,0),
M3 = [Γ,Γ(Γ∗ Γ)/2,0,0], yields the transition rate r21 which agrees, up to a pre-
−
exponential factor, with Eq. (22). See Refs. [14, 18] for
F = [0,0,ln(b/a),0],
1
detail.
F2 = [Γ∗,0,0, ln(Γ∗/Γ)]. (28)
−
The zero-momentumfixedpoints M , M andM corre-
1 2 3
spond to the three fixed points of the mean field equa- 2. Analytic result for S32 near bifurcation Γ=Γ∗
tions (3), so we keep the same notation for them. The
twoother fixedpoints, F1 andF2, arefluctuationalfixed Because of the lack of an independent integral of mo-
points describing a fox-free state at a non-zero number tion in addition to the Hamiltonian, the canonical equa-
ofrabbits(F2)andanemptysystem(F1). Fluctuational tions(27)arenon-integrableanalyticallyforageneralset
fixed points have a non-zero px or py component and ofparametersa,bandΓ. ClosetothebifurcationΓ=Γ∗,
appear in a broad class of stochastic population mod- however, the fixed points M and F become very close
3 2
els exhibiting extinction in the absence of an Allee ef- to each other, and the action S can be calculated per-
32
fect [2, 10, 11, 14, 17, 18, 20, 22, 24, 27, 38]. They turbatively by exploiting the time separation property
play an important role in the calculations of the quasi- uncoveredinRefs. [10, 11]. Letus define the bifurcation
stationarydistributionsandthemeantimetoextinction. parameter δ = 1 Γ/Γ∗. For δ 1, the fixed points
ItismostlytheirpresencethatmakestheWKBtheoryof M3 and F2 becom−e M3 Γ∗(1≪ δ),δΓ2∗/2,0,0 and
population extinction distinct from the WKB theory of F2 (Γ∗,0,0, δ). Motiv≃ate(cid:2)d by t−hese expression(cid:3)s, we
noise-induced switches between different states that are assu≃methefollo−wingscalings: x=Γ∗+δq1,y =δΓ2∗q2/2,
stable in the mean-field theory [39]. p = δ2p , and p = δp , where q , q , p and p are
x 1 y 2 1 2 1 2
To determine S(x,y) in Eq. (25), one should calcu- (non-canonical)rescaledvariables. The equationsofmo-
late the action accumulated along the (zero-energy) ac- tion become, in the leading order in δ,
tivation trajectory in the phase space of the Hamiltons
equations (27) that exits the fixed point M and ends at
3
(x,y): Γ∗
q˙1 = (q1+Γ∗q2),
− 2
(x,y)
q
S(x,y)=ZM3 pxdx+pydy, q˙2 = δq2(cid:18)1+2p2+ Γ1∗(cid:19),
Γ∗
andthenminimizetheresultwithrespecttopx andpy at p˙1 = 2 (p1−q2p2),
the end point (x,y). To evaluate the transition rates r
31 q p
aannddr(3x2,=weΓn∗e,eydt=oe0v)a,luraetsepeπcmti,vnealyt.poAinstsin(xm=an0y,yo=th0er) p˙2 = δ(cid:18)p1− Γ1∗2 −p2−p22(cid:19). (31)
stochastic population models, there are no phase trajec-
The dynamics of q and p is fast compared to that of
tories that would start at M and end at fixed points 1 1
3
q and p . As a result, q and p adiabatically follow
M or M . There are, however,instanton-like activation 2 2 1 1
1 2
trajectories that start, at t = , at M and enter, at the slowly varying q2 and p2: q1 = Γ∗q2, p1 = q2p2.
3 −
−∞ Plugging these expressions in the equations for q and
t= , the fixed point F or F , respectively. The accu- 2
1 2
∞ p , we obtain equations
mulated actions S and S along these instantons, 2
31 32
F1 F2 q˙2 =δq2(1+2p2 q2), p˙2 =δp2(2q2 1 p2), (32)
S = p dx+p dy and S = p dx+p dy, − − −
31 x y 32 x y
Z Z
M3 M3 derivable from the reduced “universal” Hamiltonian
H (q ,p ) = δq p (p q + 1) [10, 11, 14, 38]. This
yield,withexponentialaccuracy,thetransitionsratesr r 2 2 2 2 2 2
31 −
and r : problem is easily soluble, and one obtains
32
r31 ∼ exp(−NS31) S32 δ2Γ2∗ 0(q2 1)dq2 = (Γ∗−Γ)2. (33)
r32 exp( NS32). (29) ≃ 2 Z1 − 4
∼ −
7
0.75 HaL HaL
0.75
0.5
y
0.25 M
3 F
F 2 0.5
0 1
y
0 0.5 1 1.5 2
M
3
0.25
x
0.2
HbL F1
0.1
y M 0 0.5 1
p 0 3
F x
-0.1 1 F2
-0.6 -0.4 -0.2 0 0.07
0
px py -0.07 F1 M3
-0.14 HbL
FIG. 4: The numerically found instantons M3F1 and M3F2 -0.6 -0.4 -0.2 0
for a = 2, b = 1 and Γ = 1.8, when the fixed point M3
is a node. Shown are the xy-projections (a) and the pxpy px
projections (b).
FIG. 5: The numerically found instanton M3F1 for a =
2, b = 1 and Γ = 0.4, when the fixed point M3 is a fo-
3. Numerical calculations of S31 and S32 cus. Shownarethexy-projection(a)andthepxpy projection
(b).
For a general set of parameters a, b and Γ, one can
find the instantons M F and M F numerically: either
3 1 3 2
by shooting [10, 11], or by iterating the equations for x˙ HaL
0.6
andy˙ forwardintime,andequationsforp˙ andp˙ back-
x y
0.4
ward in time [4, 24, 40]. Here we employed the shooting
y M
method and explored different parts of the parameter 0.2 3
F
space a, b and Γ. Figure 4 shows typical examples of 0 2
the numericallyfoundinstantonsM F andM F inthe
3 1 3 2 0 0.5 1 1.5 2
case when M is a node. Figures 5 and 6 refer to the
3
case when M3 is a focus. As one can see, the instanton x
M F is qualitatively similar, in both regimes, to that
3 2 0
found in the conventional SI model [11]. The instanton
M
M3F1 isdifferent, asit correspondsto extinctionofboth 3
sub-populations, absent in the conventionalSI model.
-0.5
NotethatFig.6aresemblesFig.3a,andFig.5aresem-
bles Fig. 3b. In particular, the instanton M F passes
3 2 y
close to the point (0,0). This explains the feature ob- p
served, for the same set of parameters a, b and Γ, in -1
stochasticsimulations,seeFig. 2a. Asinotherstochastic
systemsexhibitingrarelargefluctuations,oneshouldex-
ptreacjtecthtoartiebsy,acvoenrdaigtiionngaolvoenraaagliavregneenxutminbcteiroonfrsotoucthe,aostniec -1.5 HbL F2
will obtain the corresponding instanton up to an error -0.18-0.12-0.06 0 0.04
that vanishes as N .
→∞
p
Figure 7a shows the accumulated actions S31 and S32 x
forfixeda=2andb=1anddifferentvaluesofΓ(param-
eterized by δ = 1 Γ/Γ∗). The first observation is that FIG. 6: The numerically found instanton M3F2 for a =
−
S31 >S32 for all δ, implying that the effective transition 2, b = 1 and Γ = 0.4, when the fixed point M3 is a fo-
rate r32 is exponentially greater than r31. As expected, cus. Shownarethexy-projection(a)andthepxpy projection
S vanishes at δ 0 and δ 1 and has a maximum (b).
32
→ →
closetoδ =1/2,orΓ=Γ∗/2. Thisbehaviorfollowsthat
8
0.6ìì HaL toholdwhena isonly slightlygreaterthanb(notshown
0.5 here).
s 0.4 When amplified by the very large factor N 1 in
ion 0.3 the exponent, see Eq. (29), the double inequality≫S32 <
t
c S <S leads to
a 0.2 31 21
0.1
r r r , (34)
0 21 ≪ 31 ≪ 32
0 0.2 0.4 0.6 0.8 1
implyingthatthesequentialextinctionroute(foxesfirst,
∆
rabbits second) is much more likely than the (almost)
ìì HbL parallelextinctionroute. Usingtheinequalityr r ,
32 31
1.5 we canfurther simplify Eqs.(17)-(21). Inparticula≫r,the
s mean time to extinction of foxes becomes simply τ
n F
o 1 1/r , whereas the mean time to extinction of both sub≃-
ti 32
ac 0.5 population is τRF ≃ 1/r21. Here the quantities τF and
τ are determined by the “kinetic bottlenecks” of the
RF
effective transitions, as expected.
0
0 0.2 0.4 0.6 0.8 1 However,forasufficientlyhighpredationrate,andnot
∆ toolargeN,thetransitionratesr31 andr32arecompara-
ble, and Eqs. (17)-(21) should be used in their complete
form.
FIG. 7: (Color online) Accumulated actions S31 (solid line)
andS32 (dashedline)asfunctionsofδ=1−Γ/Γ∗ forafixed
b = 1 and two different values of a: 2 (a) and 3 (b). The
dotted lines show the analytical asymptote from Eq. (33). 4. Proximity of instantons M3F1 and M3F2 at small Γ
The diamonds denotethe valuesof S21 from Eq. (30).
As is evident from Fig. 7, the actions S and S ap-
31 32
proach each other closely for sufficiently large δ, that is
of the mean-field population size of foxes in the coexis- forhigh predationrates. The reasonfor itbecomes clear
tence state. S31 behaves differently: it has a maximum upon comparison of the instantons M3F1 and M3F2 for
at δ 0 (that is, at Γ Γ∗), goes down monotonically sufficiently small Γ, see Fig. 8. One can see that the in-
with→an increase of δ, a→nd tends to zero as δ approaches stantonM3F2 almostcoincideswith the instantonM3F1
1 (that is, Γ approaches 0). This behavior can be qual- until, close to the fixed point F1, it abruptly changes its
itatively understood from the Γ-dependence of the fixed directionandgoestowardthefixedpointF2. Onthelat-
points, see Eq. (28). There is, however, one surprising ter segment of the trajectory the values of y and px stay
feature here. As δ approaches zero, the fixed point M close to zero, so the contribution of this segment to the
3
becomes closer and closer to the fixed point M2. One action S32 is negligible. In other words, the most likely
might expect, therefore, that the instanton M F and routeto extinctionofthe foxesinthis regimeofparame-
3 1
the accumulated action S would approach the instan- ters involvesa drastic decreasein the number of rabbits:
31
ton M F and the action S , respectively. This is not almost all the way to their extinction! Then, a few re-
2 1 21
whathappens. Extrapolatingthe numericallycalculated maining rabbits start reproducing, at a very small num-
values of S towardδ =0, we obtain S (δ =0) 0.48. beroffoxes,andtherabbitpopulationgetsreestablished.
31 31
This is considerably less than S (δ = 0) = 0.6≃137..., Note that, in the WKB formalism, the reestablishment
21
obtainedfromEq.(30),seeFig.7a. Furthermore,thenu- of the rabbits, accompanied by extinction of the foxes,
mericallyfoundinstantonsM3F1 atsmallδ aremarkedly occursviathe“fluctuational”fixedpointF2,ratherthan
different from the instanton M2F1 for which y py 0. via directly approaching the mean-field point M2.
≡ ≡
Notonlythe instantonsM F exhibitrelativelylargein-
3 1
termediate values of y and p , but the maximum value
y
of p along the instanton actually increases as δ goes 5. Two modifications of the model
y
down. That is, as δ 0, or Γ Γ∗, the fluctuations in
→ →
the number of foxes play an important role in the joint One central result of this work is the strong inequal-
extinction of the rabbits and foxes. ity r r observed in most of the parameter space
31 32
≪
We also observedthese features for other values of pa- thatweexplored. Thatis,thesequentialextinctionroute
rameters a, b and Γ from the coexistence region 0 < (first the predators then, after a long time, the prey) is
Γ < Γ∗. For example, Fig. 7b shows the δ-dependence usuallymuchmorelikelythanthe(almost)parallelroute.
of the accumulated actions S and S for a = 3 and Inotherwords,thepredatorsareusuallymuchmorevul-
31 32
b = 1. Here the two actions are very close to each nerable extinction-wise than their prey. In retrospect,
other in a broader region of δ, but the double inequality this feature is hardly surprising: As the predators are
S < S < S still holds. Furthermore, it continues dependent on their prey for survival, they can be at best
32 31 21
9
3 modified ones is that, as we improve the conditions for
HaL
thepredators,thisconvergenceoccursatlowerpredation
rates.
2 M3
y
IV. CONCLUSIONS
1
Weconsideredasimplestochasticpredator-preymodel
F F initscoexistenceregionandobservedthatthismodelex-
0 1 2
hibitstwodifferentroutesofextinctionofeachofthepop-
0 1 2 3 4
ulations: the sequential and the (almost) parallel. This
x multiplicity of the extinction routes can be conveniently
accountedforbyaneffectivethree-statemasterequation
M forthe probabilitiesto observethe coexistencestate, the
0 3
F predator-free state and the empty state. The WKB ap-
1
proximation yields the effective transition rates between
-0.3 these three states, as well as the most likely paths of
the sub-populations to extinction: the instantons. We
y
p showed numerically that the parallel extinction route is
-0.6 usuallymuchlesslikelythanthesequentialone,implying
a great robustness of the prey against predation. For a
sufficiently high predation rate, however, the two routes
-0.9 HbL F may have comparable probabilities.
2
Therestofparametersbeing fixed,ourmodelpredicts
-0.9 -0.6 -0.3 0
an optimal predation rate so that the mean time to ex-
p tinction of the predators is maximum. Not surprisingly,
x
this optimal predationrate is close to that for which the
quasi-stationarypopulationsizeofthepredatorsismaxi-
FIG.8: (Color online) xy (a) andpxpy (b)projections of the
mum. Asurprisingresultisthat,forsufficientlyhighpre-
instantons M3F1 (solid lines) and M3F2 (dashed lines) for
a = 3, b = 1 and Γ = 1.5. For these parameters δ = 0.625, dationrates,theoptimalpathtoextinctionofthepreda-
and the actions S31 and S32 are very close, see Fig. 7b. tors almost coincides with the optimal path to the joint
extinctionofthe predatorsandpreyuntil apoint(corre-
sponding to an almost zero size of each sub-population)
as resilient extinction-wise as the prey. (This also gives isreachedwherethetwooptimalpathsdepartfromeach
an intuitive explanation to the proximity of the instan- other: the predatorpopulation continues moving toward
tonsM F andM F athighpredationrates,reportedin extinction whereas the prey population reestablishes it-
3 2 3 1
the previoussubsection.) Willthis featureholdforother self. Thatis, fora highpredationrate,the predatorsare
predator-preymodels? Togetinsight,weintroducedtwo more likely to reach extinction by consuming nearly all
minormodificationsofourmodel,eachofthemendowing the prey. Thisphenomenonappearstobe robustandin-
the predatorswith moreresiliencebut stillkeepingthem dependent of details of the predator-prey model, as long
dependent on the prey for survival. In one modification as predators still need prey for their survival.
we added a direct reproduction process F 2F with Finally, all our results can be reformulated in terms
→
rate d<1. This can model a large amount of additional of the SI epidemic model for an isolated community, see
small prey for the foxes, say mice. The direct reproduc- Table 2, where both the infected, and the susceptible
tion reduces the effective death rate of the foxes. Still, populations are prone to extinction.
as long as d < 1, the foxes go extinct deterministically
in the absence of rabbits. In another modification of the
model we replaced the predation process F +R 2F
→ ACKNOWLEDGMENTS
by a more efficient one F + R 3F. Our numerics
→
showed that, for both of these modified models, the in-
equality S > S continues to hold. Furthermore, in We acknowledge useful discussions with Michael
3,1 3,2
both models we observed, for sufficiently high predation Khasin and Pavel Sasorov. This work was supported by
rates,theconvergenceoftheinstantonsM F andM F the Israel Science Foundation (Grant No. 408/08) and
3 2 3 1
toeachotheruntilaclosevicinityofthefixedpointF is by the US-Israel Binational Science Foundation (Grant
1
reached. Adifferencebetweentheoriginalmodelandthe No. 2008075).
10
[1] M. S. Bartlett, Stochastic Population Models in Ecology P07018 (2010).
and Epidemiology (Wiley, NewYork, 1961). [22] B. Meerson and P.V. Sasorov, Phys. Rev. E 83, 011129
[2] O. Ovaskainen and B. Meerson, Trends in Ecology and (2011).
Evolution 25, 643 (2010). [23] M. Khasin and M. I. Dykman,Phys. Rev. E 83, 031917
[3] J.W.TurnerandM.Malek-Mansour,PhysicaA93,517 (2011).
(1978). [24] I. Lohmar and B. Meerson, Phys. Rev. E 84, 051901
[4] V. Elgart and A. Kamenev, Phys. Rev. E 70, 041106 (2011).
(2004). [25] A.J. Lotka, Proc. Natl. Acad. Sci. U.S.A. 6, 410 (1920);
[5] C.R. Doering, K.V. Sargsyan, and L.M. Sander, Multi- V. Volterra, Le¸cons sur la Th´eorie Math´ematique de la
scale Model. and Simul. 3, 283 (2005). Lutte Pour la Vie (Gauthier-Villars, Paris, 1931); J.
[6] M. Assaf and B. Meerson, Phys. Rev. Lett. 97, 200602 D. Murray, Mathematical Biology: I. An Introduction
(2006); Phys.Rev. E75, 031122 (2007). (Springer, Berlin, 1989).
[7] T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E [26] H. Hethcote, Math. Biosci. 28, 335 (1976); SIAM Rev.
74, 051907 (2006). 42, 599 (2000); I.N˚asell, J. R. Stat. Soc. Ser. B 61, 309
[8] D.A. Kessler and N. M. Shnerb,J. Stat. Phys.127, 861 (1999).
(2007). [27] O. A. van Herwaarden and J. Grasman, J. Math. Biol.
[9] M. Assaf and B. Meerson, Phys. Rev.Lett. 100, 058105 33, 581 (1995).
(2008). [28] O. A. van Herwaarden, J. Math. Biol. 35, 793 (1997).
[10] M. I. Dykman, I. B. Schwartz, and A. S. Landsman, [29] R. Kubo,K. Matsuo, and K. Kitahara, J. Stat.Phys. 9,
Phys. Rev. Lett. 101, 078101 (2008); I. B. Schwartz, L. 51 (1973).
Billings, M. Dykman, and A. Landsman, J. Stat. Mech. [30] G. Hu,Phys.Rev. A 36, 5782 (1987).
P01005 (2009). [31] C.S. Peters, M. Mangel, and R. F. Costantino, Bull.
[11] A. Kamenev and B. Meerson, Phys. Rev. E 77, 061107 Math. Biol. 51, 625 (1989).
(2008). [32] M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J.
[12] A.Kamenev,B. Meerson, and B. Shklovskii, Phys.Rev. Chem. Phys.100, 5735 (1994).
Lett.101, 268103 (2008). [33] B. Meerson and P. V. Sasorov, Phys. Rev. E 78,
[13] M. Assaf, A. Kamenev, and B. Meerson, Phys. Rev. E 060103(R) (2008).
78, 041123 (2008). [34] C. Escudero and A. Kamenev, Phys. Rev. E 79, 041149
[14] M. Assaf, A. Kamenev, and B. Meerson, Phys. Rev. E (2009).
79, 011127 (2009). [35] B.Meerson,P.V.Sasorov,andY.Kaplan,Phys.Rev.E
[15] M. Parker and A. Kamenev, Phys. Rev. E 80, 021129 84, 011147 (2011).
(2009); J. Stat. Phys. 141, 201 (2010). [36] B. Meerson and P.V. Sasorov, Phys. Rev. E 84,
[16] B. Meerson and P.V. Sasorov, Phys. Rev. E 80, 041130 030101(R) (2011).
(2009). [37] M. Assaf, E. Roberts, and Z. Luthey-Schulten, Phys.
[17] M. Khasin and M. I. Dykman, Phys. Rev. Lett. 103, Rev. Lett.106, 248102 (2011).
068101 (2009). [38] V. Elgart and A. Kamenev, Phys. Rev. E 74, 041101
[18] M. Assaf and B. Meerson, Phys. Rev. E 81, 021116 (2006).
(2010). [39] M. I. Freidlin and A. D. Wentzell, Random Perturba-
[19] M.Khasin,B.Meerson,andP.V.Sasorov,Phys.Rev.E tionsofDynamicalSystems (Springer-Verlag,NewYork,
81, 031126 (2010). 1998), 2nd ed.
[20] M. Khasin, M. I. Dykman, and B. Meerson, Phys. Rev. [40] A. I. Chernykh and M. G. Stepanov, Phys. Rev. E 64,
E 81, 051925 (2010). 026306 (2001).
[21] M.Assaf, B. Meerson, and P.V.Sasorov, J. Stat.Mech.
