- Email: [email protected]
Persistence in density dependent stochastic populations
Persistence in Density Dependent
Stochastic Populations*
CHARLES TIER ANDFLOYD B. HANSON Departmeni of Mathematics, Unioersi@ of Illinois at Chicago Circle, Chicago, Iilinois 60680 Receiwd
26 Februaty 1980; revised 28 Ju!v 1980
ABSTRACT Persistence, as measured by time to extinction, is studied in a density dependent population that is subject to small environmental and demographic randomness. Diffusion processes are formally derived from branching processes in constant and random environments. The moment generating function of the extinction time, which satisfies a second order ordinary differential equation, is found asymptotically in the limit of small diffusion and is related to the diffusion limit of the Galton-Watson process and the OmsteinUhlenbeck process. Extinction occurs with probability one, though the mean and variance of the extinction time are found to be exponentially large and suggest the extinction time is exponentially distributed. The notion of persistence is compared with other qualitative measures of stability. Four examples are studied and compared.
1.
INTRODUCTION
Persistence, as measured by time to extinction, has been suggested by Ludwig [19] as a qualitative measure of the effects of random fluctuations on a species. The inclusion of small random fluctuations or stochastic effects may lead to the eventual extinction of a species even when its population size is near a nontrivial stable deterministic equilibrium. The extinction time r is then a random variable which is the first passage time to the zero population size. If extinction occurs with probability one, persistence can be measured by the moments of the extinction time, such as the mean and the variance. We shall study persistence using diffusion theory with population models that are deterministically of the logistic type but include small random fluctuations due to demographic and environmental stochasticity. Hence we shall study the diffusion process X(t) with equations involving the back*Research supported by the National Science Foundation under Grants MCS 78-01167 and MCS 79-01718. MATHEMATICAL
BIOSCIENCES
53:89- 117 (1981)
OElsevier North Holland, Inc., 1981 52 Vanderbilt Ave., New York, NY 10017
89 0025-5564/81/010089+29$02.50
90
CHARLES
TIER AND
FLOYD
B. HANSON
ward operator
fqx)uxx +wx)ux with the proper auxiliary conditions.
Here M(x)
(1.1) and V(x)
are the infinites-
imal mean and variance, respectively, of the diffusion process X(r), and X(0)=x is the initial population size. Demographic stochasticity or sampling fluctuation is due to randomness in the survival ability or fertility of individuals within the population. Typically, models with demographic stochasticity arise from discrete variable stochastic processes such as birth and death processes or branching processes. A diffusion process can arise as a limit approximation of the discrete process. For example, the diffusion limit of the Galton-Watson process has infinitesimal mean and variance (Feller [8]) M(x)=rx,
(1.2)
v(x)=px. Here there are no density dependent
effects, since the per capita birth rate
M(x)/x and variance V(x)/x are constants. Density dependent diffusion models arise when offspring distributions
or
per capita birth and death rates depend upon the current population size. Goel and Richter-Dyn [ 1 l] obtained density dependent diffusion models as formal limits of general birth and death processes. A rigorous proof of the convergence of a continuous time, density dependent branching process to a diffusion process was given by Lipow [ 171 and by Kurtz [ 161 within the context of general continuous time Markov infinitesimal mean and variance of the form
M(x)=-v(x), V(x)=j3x.
chains.
Lipow
derived
an
(l-3)
Diffusion models with (1.3) in general lead to extinction with probability one and a finite mean extinction time if the per capita growth rate p(x) is integrable near x=0. Thus, in this case, the mean extinction time can be used as a measure of persistence. Environmental stochasticity results from random fluctuations in the environment affecting the population as a whole, beyond the effects the environment has on the variance of the offspring distribution. Environmental stochasticity can be introduced into a standard deterministic model by replacing one or more of the parameters with random processes. This type of problem has been analyzed extensively by May [23], Feldman and Roughgarden [7], and Turelli [26]. For example, if Gaussian white noise is
PERSISTENCE
IN STOCHASTIC
91
POPULATIONS
introduced into parameters in the logistic equation a!x/dt =rx( 1 -x/k), a diffusion process is obtained using Ito calculus with infinitesimal mean and variance M(x)=rx(l
k, r>O,
-x/k),
V(x)=x’C(x),
Z(x) > 0,
(1.4)
where t(x) is related to the variance of the noise. For the diffusion model with (1.4), the zero population size is not attainable. Hence either the probability of extinction is equal to one and the mean extinction time is infinite, or the probability of exit is equal to zero. In this case, persistence is not useful, and other qualitative measures of stability have been used, some involving the stationary distribution by May [23] and Turelli [26]. Models with both environmental and demographic stochasticity have been investigated by Athreya and Karlin [2], Keiding [14,15], Barbour [4], and Ludwig [20]. Keiding [IS] and Ludwig [20] formally derive diffusion approximations to discrete time, density dependent branching processes with a random environmental parameter in the mean of the offspring distribution. Their diffusion approximations have an infinitesimal mean and variance M(x)=rx(l
-x/k),
k, r>O,
v(x)=xd(x)+x2c?(x),
ci(x)>O,
5(x)>O,
(1.5)
where x is a scaled population size. In V(x), the term xd(x) arises from the variance of the demographic stochasticity, and the term x’c?(x) arises from the variance of the environmental stochasticity. The deterministic approximation, 2, satisfies the logistic equation di/dt = M(P) = r?( I- 2/k). In Sec. 2, we formally derive a diffusion equation with coefficients (1.5) following Ludwig [20]. Ludwig investigated, asymptotically, the mean extinction time when the environmental stochasticity dominated the demographic stochasticity, i.e., C?(X)= 1/N, N > 1, and g(x) = (Yin (1.5). Analogous models arise in population genetics with the demographic stochasticity, called random genetic drift, due to sampling error and the environmental stochasticity due to random selection (see Gillespie [lo]; Karlin and Levikson [ 131). The extinction time problem is the first passage problem in probability theory. As we derive in Appendix 1, the moment generating function of the extinction time 7, u(x, s), satisfies
~V(x)u,, +M(x)u,
=su,
o
(1.6)
92
CHARLES
TIER AND
FLOYD
B. HANSON
FIG. 1. The interval x > 0 with the four regions B, L, T, and R indicated. solutions us, UL, UT, and uR are given in their respective regions of validity.
The
If M(x) and V(x) are given by (1.5) with d(0) > 0, then x = 0 is an exit boundary and x = cc is an entrance. The analysis of (1.6) subject to the proper boundary conditions is difficult except for certain special cases (Darling and Siegert [5]; Ricciardi [25]; Goel and Richter-Dyn [ 111). In Sec. 3, we present our main results, where we find asymptotic solutions of (1.6) with (1.5) but with d(x)=d(x)/N and g(x)=e(x)/N, N ~1. Here N is some reference population size, and its reciprocal measures the smallness of the variance. The WKB approximation technique and boundary layer analysis of Hanson and Tier 1121, extended to the infinite domain, are used for large N. We find that the solution has different approximations depending upon the value of x (see Fig. 1). For x near the stable deterministic equilibrium k, the process is approximated by the Omstein-Uhlenbeck process to leading order. This is analogous to the approximation of discrete processes about an equilibrium using the Grnstein-Uhlenbeck process by McNeil and Schach [21] and Barbour [3]. Similar results have been derived by Feller [8], by Ludwig [18], and in nonbiological contexts by Matkowsky and Schuss [22]. Near x=0, we find that the process is approximated by the diffusion limit of the Galton-Watson branching process. This is a result of the demographic stochasticity. This approximation and the one near x = k are then connected together using the WKB solution and asymptotic matching. The usual boundary conditions on the diffusion limit of the Galton-Watson process and the Ornstein-Uhlenbeck process are thus replaced by matching conditions. The matching process is described in Appendix 2. Using the asymptotic solution for N>l we find in Sec. 4 that extinction occurs with probability asymptotic to 1. The mean and the variance of the extinction time for x away from 0 are found to be E[7]--h,, Var[ T]-A~, where A,, given by (4.7) is exponentially large in N. This indicates that the population persists for long periods of time before extinction. In addition this is evidence that the extinction time 7 is exponentially distributed for N>l. For x near 0, we find that E[T] and Var[~] change rapidly with x,
PERSISTENCE IN STOCHASTIC POPULATIONS
93
which is expected and is typical of boundary layer behavior. We also compare our results for persistence with the notion of stability used by other authors [23,26] and find that they are qualitatively similar. In Sec. 5, we consider four density dependent models, one with a constant environment and three others that differ by the manner in which the environmental stochasticity is introduced. We are able to conclude that the introduction of environmental stochasticity in models with demographic stochasticity decreases the persistence over models with a constant environment since it decreases the mean extinction time. This is the result of the environmental stochasticity increasing the probability of the population size becoming small where the demographic stochasticity can cause extinction. 2.
THE DIFFUSION
APPROXIMATION
We shall now formally derive diffusion approximations for discrete time, density dependent branching processes in a constant environment and in a random environment following Keiding [15] and Ludwig [20]. Let Z, be a discrete random variable equal to the number of individuals in the population in generation n. We assume that the generations are nonoverlapping and the individuals within each generation reproduce independently. Let the number of offspring born to individual i in generation n be given by the random variable Bi, all of which are independent and identically distributed. CONSTANT
ENVIRONMENT
For the density dependent branching process in a constant environment, the distributions of the Bi are defined by the probabilities pi(z) = Pr( Bi = A-C =z) with mean and variance h(z)=E[B,IZ,=z],
(2.1) v(z)=Var[BilZn
=z],
respectively. Within each generation, the total population Z, is assumed to be fixed and hence the distribution and moments of Bi are conditioned on Z” =z. The number of individuals in generation IZ+ 1 is the sum of the offspring of all the Z,, individuals, and hence
(2.2) Using (2.2) with (2.1), we can compute
the first and second moments
of the
CHARLES
94 increment
TIER AND FLOYD
B. HANSON
AZ, = Z,,, , - Z, as E[AZ,IZ,,=z]=z[h(z)-11, (2.3) E[(AZ,)2jZ,
=z] =zv(z)+z~[~(z)-
112.
To formally obtain a diffusion approximation to the branching process (2.2), we must transform the discrete time, discrete size branching process into a new continuous time, continuous size process. Thus, we define the continuous time variable t=n At, where the increment At<1 is the generation time. A continuous population size variable, as a function of t, is defined by X(t) = Z,, /L, where L is some large reference population size. In addition, we require that small changes in X(t) occur in small time increments At. This requires that the offspring mean h be approximately equal to 1, i.e., near replacement. Thus we assume h(z)=l+B/k(z/L),
(24
where 6 measures the deviation from replacement. We also assume that v(z)=d(z/L). The diffusion approximation X(t) is obtained in the limit as iI-+0 when we choose At = 6 and L = 1/a. Lipow [ 171 proved rigorously the convergence of a continuous time, density dependent process in a constant environment to a diffusion process with a(x) a constant [see (1.3)]. The diffusion approximation can also be made rigorous by using the technique of Ethier and Norman [6]. If we introduce the new continuous variables into (2.3), the infinitesimal mean and variance of the diffusion approximation X(t) are given by
M(x)=
lim Al-+0
E[AX(t)lX(t)=xl =xpFL(x) , At (2.5)
,
At-0
where x= 8z represents a scaled population size. If the higher central moments of Bi are suitably bounded, then the higher moments of AX(t) satisfy E[(AX(t))‘IX(t)=x] At as At+O+.
=O[(At)“‘-‘I,
j>2,
(2.6)
PERSISTENCE
IN STOCHASTIC
95
POPULATIONS
The deterministic approximation to X(t), P, satisfies dz?/dt=fp(2), which is characteristic of deterministic population models. If we assume that the per capita growth rate is p(x) =r( 1 -x/k), then the infinitesimal mean and variance are M(x)=rx(l
-x/k), (2.7)
V(x)=xci(x),
and the deterministic approximation is the logistic equation with scaled carrying capacity k. The diffusion approximation with (2.7) is a stochastic generalization of the logistic equation. This problem with &x)=/3/N is discussed in the examples in Section 5. RANDOM
ENVIRONMENT
We now modify the density dependent branching process described above by allowing random fluctuations in the environment from generation to generation. These fluctuations are accounted for by introducing a random environmental parameter, W,, into the offspring distribution. W, are independent with mean zero, variance 1, and bounded higher moments, and are independent of Z,,, , m < n . The offspring distribution is now defined by pi(z, w)=Pr(B, =jl Z, = Z,W” =w), which are conditioned on both the population size z and the environment w in generation n. Let h(z, w) and Y(Z, w) be the conditional mean and variance, respectively, of the offspring distribution. We also assume the average of Y(Z, w) over the environmental variable W, is
E[vlZ,
=z]=d(z/L),
(2.8)
where a is now the expected value of v. To obtain the diffusion approximation, we assume t = n At, X(t) = Z, /L, the near replacement condition
h(z, w)=
1+ Vl%(z/L)
w+Gp(z/L),
(2.9)
and At = 6= l/L. The fluctuations due to the random environment are assumed to be of order fl . This may be motivated by assuming that these fluctuations are due to the sum of a large number of independent, identically distributed random variables. Then fluctuations about the mean are of the order of one over the square root of the total number in the sum.
96
CHARLES
The infinitesimal
FLOYD
B. HANSON
mean and variance of the diffusion approximation
are computed by iterated expectation M(x)=
TIER AND
lim
At-+0
using (2.3), (2.8)
EIAx(t)lx(t)=xl At
X(t)
and (2.9) as
’
=&moE{EIAZ~]Zn,Wn]]Z,,=x/At,n=t/At}, ‘V(X), v(x)=
(2.10)
lim E[(AX)*IX(t)=x] At At-t0
’
= lim SE{E[(AZ,)*IZ,,W,]IZ,,=x/At,n=t/At}, At-+0
If the higher central moments of W, and Bi are suitably bounded, then the higher moments of AX satisfy a condition of the form (2.6). The diffusion approximation for the constant environment model is a special case of the random environment model with 5(x)=0. Both models have the same mean, since the environmental parameter has mean zero. In V(x) in (2.10), the term x&x) arises from the variance of the offspring distribution, while the term x*P(x) arises from the fluctuations in the environment. The models of May [23], Feldman and Roughgarden [7], and Turelli [26] with random environments but no demographic stochasticity can be obtained by setting &x)=0. Their models were obtained, for example, by replacing with r/k
r in the logistic equation dx/dt=rx(l fixed and y as normalized Brownian
stochastic
differential
-x/k) by r+ ~6 dy/dt, motion. This leads to the
equation dx=rx(l-x/k)dt+&xdy
and a diffusion,
(2.11)
using Ito calculus, with mean and variance M(x)=rx(l V(x)=ax2.
-x/k),
(2.12)
This corresponds to setting d(x) = 0 and t?(x) = a in (2.10). This completes the formal derivation of the diffusion approximation. Our principal interest is the analysis of diffusion equations where V(x) is small relative to M(x). Therefore, we replace d(x) by d( x)/N and C?(X) by .$x)/N in (2.10), which gives (3.1). This scaling assumes that the stochastic effects are small compared to deterministic effects but the environmental fluctuations are the same order as the demographic effects.
PERSISTENCE
3.
IN STOCHASTIC
ASYMPTOTIC
97
POPULATIONS
SOLUTION
OF EXTINCTION
PROBLEM
We shall now construct an asymptotic solution for the moment generating function of the extinction time r, u(x, s; N). As we show in Appendix 1, u satisfies
&X)%x+m(x)u,
o
=su,
(3.1)
where m(x)=rx(l-x/k),
k, r>O,
u(x)=xd(x)+x2e(x),
d(x)>O,
e(x)>O,
(3.2)
and d(x) and e(x) are bounded for x< co. The equation (3.1) is (1.6) where the coefficients (3.2) are the same as (1.5) but with a(x), g(x), and V(x) replaced by d(x)/N,e(x)/N, and u(x)/N, respectively. The boundary conditions for (3.1) with (3.2) as discussed in Appendix 1 are ~(0, s; N)=
1,
(3.3a)
lim u(x, s; N) < 00, x-m
(3.3b)
since x=0 is an exit boundary and x= cc is an entrance boundary. We seek an asymptotic solution valid as N-+oc using the WKB method (O’Malley [24]; Hanson and Tier [12]). To do so, we assume that the solution is of the form X’p(~,s;N)dz
.
(3.4)
Substituting (3.4) into (3.1) yields the Ricatti equation (p2+2_ m(x) --- 2s u(x) ‘- NV(X)
(3.5)
We assume thatcp has an asymptotic expansion of the form p(x,s;
N)=cp,(x,s)+
w
+o(-$),
A’--,
(3.6)
substitute (3.6) into (3.5), and equate coefficients of powers of l/N to zero; then we find to leading order the two solutions cp; =0
and
‘p:=2@‘(x),
(3.7)
98
CHARLES
TIER AND FLOYD
B. HANSON
where (3.8) The next order terms give
cpz=S/m(x>, cpz’= -rp;
(3.9)
-w/w.
Since cp; are not integrable at either x=0 or x= k, we see that the asymptotic validity of the WKB approximation will fail in some small neighborhood of these points, i.e., regions B and Tin Fig. 1. The point x=0 is a singular point, and the point x= k is a turning or transition point. The WKB approximation (3.4) is valid in region L in Fig. 1, and a similar approximation is valid in region R in Fig. 1 but subject to the boundary condition (3.3b). The different solutions in Fig. 1 are defined by
I 1
uB,
u(x, s)w
UL2
UT? uR>
xin xin xin xin
B, L, T, R.
The approximate solutions, uR in R and uL in L, are given by a linear combination of the two solutions of the form (3.4) with (3.7)-(3.9), which we give as
(3.10)
lqx,
N) =
( -&J
exp;;fl;;,x” ,
(3.11)
with m(x) and u(x) given by (3.2). We note that A(x; N) grows exponentially with N, since cP(x)>O. The constants aR,aL, b,, and b, are to be
PERSISTENCE
IN STOCHASTIC
99
POPULATIONS
determined, and
“‘.
(3.12)
Applying the boundary condition (3.3b) to uR given by (3.10), which is valid for x > k, we find that (3.13)
b, -0,
since A(x; N) is always unbounded as x+00 for fixed N. The approximate solution in the transition region T, u&, s; NT), is found by introducing the stretched variable {= Nr(x- k) into (3.1) with NT given by (3.12). Then we obtain to leading order the transition region equation (3.14)
UT,55
The left side of (3.14) is the backward operator of the Ornstein-Uhlenbeck equation, which has been derived in [8], [21], and [3]. The solution of (3.14) is a linear combination of the standardized parabolic cylinder functions [ 11,U(s/r - 0.5, {) and V(s/r - 0.5,{), in the form UT-[
U,U(f
-0.5,
{)+b,V(;
-0.5,
I)] exp( $),
(3.15)
where aT and b, are to be determined. From the theory of asymptotic analysis of differential equations (O’Malley [24]), we see that the transition region approximation UT is valid for at least Ix- kl= O(N -‘/‘). Clearly, UT is not a good approximation to the solution at points that are order one from x = k, such as the boundary point. To obtain consistent approximations to the solution we require that UT be consistent or match with us and uL as N+o3. The procedure of matching us and UT consists of expanding UT for {-+oc and uR as x+k +, and then choosing aR, b,, a,, and bT so that the leading terms of uR and uT are identical. We sketch the matching in Appendix 2; complete details are in Hanson and Tier [12]. From Appendix 2, we find 6, =O, aR =aT.
(3.16)
The constant ar will be determined by the boundary condition at x =O. Similarly, uT as I+ - 00 is matched to uL as x+k -. From Appendix 2, we
100
CHARLES
TIER AND FLOYD
B. HANSON
find aL =a+x~,
b,
r
(3.17)
=T&
In order to extend the results to region B we must construct a new solution, since uL is not valid near x=0. The solution in the boundary region B, u,(.$, s; N,), is found by introducing the stretched variable (3.18)
(=NBx
where (3.19)
Ns -2Nr/u’(O). Substituting .$ for x in (3.1) and expanding as N+cc leading term the boundary layer equation
‘tu B. Et +,I$,,< =
with 5 fixed gives as a
fUB.
(3.20)
We note that the left side of (3.20) is the backward operator for the diffusion limit of the Galton-Watson process [see (1.2)]. The solution of (3.20) satisfies the boundary condition at x = 0 with a matching condition with uL at [= cc instead of the boundary condition associated with (1.2). The general solution of (3.20) is
UB-[UBw~(-~,~)+bBw2(-f,E)lexp(-~), (3.21) where we define
W,(k,x)=Wk,,/z(x) w,(k,x)=cos(ka)Re[
W- k,1,2(-X)]+r(l--)s~2(~~)~~.l/2(x)~ (3.22)
The functions W,,,,,(x) and M h,,,Z(x) are Whittaker functions [l]. The constants aB and bB are to be determined. The boundary region solution uB is defined for 0 < 4 < co, and we require that it satisfy the boundary condition at x= 0 and that uB as &cc match uL The size of region B is at least x= 0(1/N). The matching, as x+0+.
PERSISTENCE
IN STOCHASTIC
sketched in Appendix
POPULATIONS
101
2, gives a,rA,(k2NBNT)S’r
a,=
’
Vs/r) b, =a+os
(
7
>
(3.23)
(k2NBNT)-S’r,
where A,=
N)=d(O)(
liiA(x;
r~~(k))“2ex~~2Nm(0)].
Using (3.23) in (3.21) at x = 0 with the results that W,( k, 0) = l/I’( 1 -k) and W2( k, 0) = (cos ka)/r( 1 + k), the absorbing boundary condition u,(O, s) = 1 gives 1
a T=
.
r( k2NBNT)s’rh0
(3.24)
+ cos2(m/r)(k2NBNT)-S”
r(s/r)r(i +s/T)
qi -+I
We summarize the solution with aT given by (3.24) and A(x; N) by (3.11) in the different regions in Fig. 1 as: Region R: (3.25) Region T: u,-a,(l(S-0.5,l)exp($),
(3.26)
where
{=&(x-k),
NT
= [2rN/o( k)]“‘;
Region L:
s/r N,k;x-k, +qx.
,
N)
(NTklx-kl/x)“” ; r(s/r)
1
(3.27)
102
CHARLES
TIER AND FLOYD
B. HANSON
Region B:
r(
( )i 6
- 2
u,-a,exp
k2NBNJ”AOW,( -s/r,
5)
Us/r)
(7 >(k2NBNp’w*(-s/r,
+cos
5)
1 )
(3.28)
where W, and W, are given in (3.22) and N, =2rN/o’(O).
5=N*x, 4.
MOMENTS
OF EXTINCTION
TIME AND PERSISTENCE
Since u(x, s) is the moment generating function of r, its moments are given by E[ 7”]
qYf(
x,r)dt=(-l)~~u(x,s)~,_,’
(4.1)
where f(x, t) is the probaslity density function of r discussed in Appendix f (x, t ) dt = u( x, 0) represents the ultimate proba1. The zeroth moment / bility of extinction and ii found by evaluating the limit as s-0 in us, uL, tlr, and uR, given by (3.25)-(3.28). Using 1 ~ I(s/r)
=
s+o(s*), r
(4.2)
as s+O, we find that lima,=I(l)= s+Q
1
(4.3)
uL, f+, u,-1.
(4.4)
and hence limu,,
s-0
Hence r is a proper random variable, and the probability of extinction is asymptotically equal to 1 for all x.
PERSISTENCE
IN STOCHASTIC
103
POPULATIONS
If we denote the mean and variance of r for x in region T by Er[ r] and Var,[ r], respectively, then from (4.1)
Var,[r]=
a2
lim
_?ur(x,s)-
s-0+
[
(
as
11 2
&ur(x,s)
Differentiating (3.26) with respect to S, letting s+O+,
(4.5)
.
and using (4.3) and
aa, Jl-i+ T”-Ao, we find after a short calculation that
(4.6) Here A0 is exponentially large and is given by
A0
exp(2NQ(O)). =dP)( r1&))“2
(4.7)
A similar calculation for the mean and variance of r in regions 8, L, and R gives
EB[~]-Ao(l-e-E),
N+Co,
,
(4.8)
and Var,[r]-Var,[r]-At, Var,[r]-A:(1
-eezE),
I
N-XCJ,
(4.9)
where A0 is given by (4.7) and .$= 2rNx/u’(O). In region L we have neglected A(x; N) as exponentially small compared to A,. The moments of r can also be computed by a modification of the method of Matkowsky and Schuss [22] as described in Hanson and Tier [12]. It is interesting to note that the means and variances of r for x in L, T, and R are identical and independent of x. Thus, except for x in B, E[ r]-A,, and Var[r]-11; are characteristic of an exponential distribution with parameter l/A,. Here the mean and variance are exponentially large as N+m, since A,, is exponentially large. This can be interpreted by the fact
104
CHARLES
TIER AND
FLOYD
B. HANSON
that the population is strongly attracted to the stable equilibrium at x = k in the presence of small random fluctuations and hence lingers near k. In addition the mean and the variance in regions L, T, and R are independent of x, and we are unable to distinguish populations that start there. Eventually, the random fluctuations cause the population size to become small, so that the demographic effects are important and can lead to extinction. For initial values near x=0 (i.e., in region B), the mean and the variance of r, given by (4.8) and (4.9) are very sensitive to changes in x. The mean extinction time increases rapidly with x from the value EB[7] =0 at x=0, where absorption is certain. Thus, the larger the initial population size, the longer the mean extinction time. From (4.8) with ,$=2rNx/o’(O) we see that a key factor in measuring the persistence of a population with a small initial size is the ratio rN/u’(O). This quantity represents the ratio of the growth rate r to the rate of change of the demographic variance at x=0. The larger the value of r or the smaller the value of u’(O), the greater the persistence of the population for large fixed N. Other authors (May [23], Feldman and Roughgarden [7], Turelli [26]) have examined the “stability” of populations in a random environment in the absence of demographic effects which corresponds to d(x)=0 in (3.2). With d( x)=0 in (3.2), extinction is not possible, since the zero population size is not attainable (x = 0 is an entrance, or a natural boundary in the terminology of Feller [9]). A qualitative measure of the effects of the random environment is whether a stationary distribution exists which depends upon the behavior near x= 0. In the case when no stationary distribution exists, the mean extinction time is infinite and hence it is not a good qualitative measure. The stationary distribution li, if it exists, is found using the forward equation and is
where C is chosen to normalize the distribution, u(x) =x*e(x) > 0, and @J(X) is given by (3.8). To distinguish populations in a random environment that have stationary distributions, May [23] and Turelli [26] have measured the spread or dispersion of the stationary distribution about the deterministic equilibrium k, which is the mean of the stationary distribution if the noise is small. Turelli’s account [26, pp. 250,264] of May’s “dispersion”criterion [23] is that the square of the coefficient of variation, variance/(mean)2, of the stationary distribution (4.10) locally about x = k must be bounded by some constant (e.g. 1). By expanding O(x)= -m’(k)(x-k)*/2u(k) in (4.10), we find that their criterion gives -m’(k)>o(k)/2Nk*.
(4.11)
PERSISTENCE
IN STOCHASTIC
105
POPULATIONS
It is important in their criterion that the noise be small; otherwise the stationary distribution may not be peaked about k. In our problem we have explicitly given the size of the noise as l/N, which is small when N> 1. Then the criterion (4.11) is always satisfied and is not particularly useful in distinguishing populations. May [23] and Turelli [26] developed an approximation method that could be applied directly to the forward equation, which gives the criterion (4.11). Their method is related to the approximation (3.14) in the transition region T and the Omstein-Uhlenbeck process. The problem we study, (3.1) with (3.2), contains the added aspect of demographic stochasticity, which is more realistic. With the demographic stochasticity included, extinction is possible (x=0 is an exit boundary) and the stationary distribution is identically 0. As suggested by Ludwig [19], persistence (using the extinction time) provides an absolute criterion for measuring the qualitative effects of the randomness. The importance of the asymptotic result for the mean extinction time (4.6) derived in this section is that it provides a measure of persistence and is valid for all x. It was derived with the assumption of small noise, N>>l, consistent with May [23] and Turelli [26]. The quantity A, given by (4.7) distinguishes different mean extinction times, and for N> 1 it is dominated bykthe exponent 2N@(O) in the exponential function. The integral Q(O) = s [ m(x)/u(x)] dx depends upon 0
the behavior of m(x)/u(x) over the entire interval [0, k] and not just near x- k. While it is better to calculate the integral a(O), a simplified approximation of Q(O) can be found using tangent line approximations at x =0 for XE(O, k/2) and at x=k for xE(k/2, k) to give m(x)k ___ 20(x) >>1.
+ (m/u)‘(x)@ 8
]-::;:y
) (4.12)
This formula gives an asymptotic analogy of the condition (4.11) derived by the criterion of May, since 2NcP(O)> 1 when N>> 1. It represents an improvement over (4.1 l), since it uses behavior not only near x =k but also near x=0. We shall use (4.12) in the examples in Sec. 5. 5.
EXAMPLES
We now apply the results in Section 4 to study the persistence of four density dependent stochastic population models. All four models have the logistic equation as an infinitesimal mean, but they differ in the manner in which the environment is taken into account.
106
CHARLES
TIER AND FLOYD
B. HANSON
Model I (Constant Environment). The diffusion approximation in a constant environment is characterized by (3.1) with m(x)=rx(l-x/k),
(5.1)
u(x)=Px,
i.e., e(x)=0 in (3.2). This type of model is described by Goel and RichterDyn [ 111, Lipow [ 171, and Hanson and Tier [ 121. Using (4.6)-(4.9), we find that
E[r] -(
$$)“‘exp(
$k)=Ai,\
Var[ r] -A:,
x inL, T, R,
(5.2)
J
and E,[T] -A,[
I-exp(
?)I, x inB,
Var,[T] -A:[
1 - exp( ?)I.
(5.3)
1 Model 2
The diffusion approximation is characterized by m(x)=rx(l-x/k), (5.4)
u(x)=px+cXx*,
i.e., e(x)=cx, a constant, in (3.2). Models with environmental randomness as reflected in the term ox* were discussed in [23], [7], and [26] without demographic randomness and can be interpreted as the result of a Gaussian white noise perturbation of r with r/k fixed using Ito calculus. This model including demographic randomness was studied in [20] and [15]. Using (4.6)-(4.9), we find that for x in L, T, or R
E[rl
-
Var[ ~1 -A;, where
TP Nkr3(l +ak//3)
1 I/*
exp[2N@2(0)]=A2,
(5.5)
PERSISTENCE IN STOCHASTICPOPULATIONS
107
while for x in B E,[T]
Var,[7] Model 3
-A2[1-exp(-2A+x/P)19
(5.6)
-AQl-exp(-4NrxlP)l.
The diffusion approximation is characterized by m(x)=rx(l-x/k),
(5.7)
o(x)=Px+(yx2(1-x/k)‘.
Here e(x) = o( 1 - x/k)2 in (3.2) can be interpreted as the result of introducing Gaussian white noise in r with k fixed. This model was discussed by Keiding [ 151. Using (4.6)-(4.9) we find that for x in L, T, or R E[ 71 -[ ?rp/Nkr3]“2exp[2N@3(0)]
=A,,
Var[ 7 ] -A:,
(5.8)
where
For x in B,
-A3[l-exp(-2~rx/P)1,
EB[7]
Var,[r] Model 4
-A’,[l-q$--4Nrx/kQl.
(5.10)
The diffusion approximation is characterized by
m(x)=rx(l
-x/k),
v(x)=px+ax4/k2.
(5.11)
This model was mentioned by Keiding [15]. Without demographic stochasticity (/3= 0), this diffusion was studied by Feldman and Roughgarden [7] and Turelli [26], and can be interpreted as the result of introducing Gaussian white noise in r/k with r fixed. Using (4.6)-(4.9), we find that for x in L, T, or R,
E[rl-
4 Nr3k(l
Var[ 71 -A’4
+ak/P)
1 l/2
exp[2NQ4(0)]=A4,
(5.12)
CHARLES
108
TIER AND FLOYD
B. HANSON
where (5.13) For x in B, E8[7] Var,[r]
-A4[l-exp(-2A+x/P)19
(5.14)
-A24[1-exp(-4Nrx/P)le
We can now compare the persistence of the four models described above as measured by the mean extinction time E [ T]. First we note that if cx= 0, models 2,3, and 4 reduce to model 1 since the environmental randomness enters the models only through the term containing CLThe parameter cx is the same in models 2,3, and 4 by our scaling and measures the variance of the environmental noise as in (2.12). In Fig. 2, we sketch E[T] versus x for N>>l. The graphs were sketched qualitatively using EB[7] in (5.3), (5.6), (5. lo), and (5.14), since these results are uniformly valid for x in all regions. The rapid change in the mean as x varies through region B is illustrated by enlarging region B as in Fig. 2. In addition, we obtain an ordering of the models in terms of persistence. This ordering is the same as the ordering of
i =l
@I
I I I i=4
p----_ A---l ’
0
i=2
X
k
FIG. 2. The mean extinction time I?[ T] for models 1 to 4 sketched qualitatively from the asymptotic results (5.3), (5.6), (5.10) and (5.14) for N>l. The regions B, L, T, and R are given.
PERSISTENCE
IN STOCHASTIC
109
POPULATIONS
5O.Or
0.04 0.0
5.0
10.0
Q Fro. 3. log Ai is calculated from the asymptotic a function of (I with r=k= 1, N= 100, and fi= 1.
results (5.2),(5.5),(5.8),
and (5.12) aa
the terms Ai, i= 1, . . . ,4. For the case in Fig. 2, we have A,>A,>A,>A,.
(5.15)
This ordering can be shown to be true if N> 1, or more precisely, N//3>>1, since @i(O) = rk/2/3 > ad(O) > aa, > Qz(0). In Figs. 3-5, we graph Ai versus a for different values of p. Here hi, i= 1, . . . ,4, is used to distinguish the mean extinction times of the four models using (5.2), (5.5), (5.9, and (5.12). For these figures, we have chosen r= k= 1 and N= 100. We see from Figs. 3-5 that Ai decreases, and hence the mean extinction time decreases as a increases. As a increases, the size of the random fluctuations increases, which in turn increases the chance that the sample paths of the population reach extinction. As B increases, but with the asymptotic results still valid, the ordering A,>A,>A,>A,
(5.16)
is possible for certain ranges of a. Figure 3 with fi= 1 has the ordering (5.15). In Fig. 4 with /3=5, the ordering changes from (5.15) to (5.16) at about a= 6.7. When p= 10 as in Fig. 5, the ordering is given by (5.16) for all a. We note the scale of the ordinate is different in Figs. 3-5.
110
CHARLES
TIER AND FLQYD
B. HANSON
i=l
6.0-
FIG. 4. log Ai is calculated from the asymptotic results (5.2), (5.5), (5.8), and (5.12) as afunctionofawithr=k=l,N=lOO,and~=5.
0.01 0.0
5.0
10.0
Q FIG. 5. log Ai is calculated from the asymptotic results (5.2), (5.5), (5.8), and (5.12) as afunctionofawith/Pk=l,N=100,andB=10.
PERSISTENCE
IN STOCHASTIC
POPULATIONS
111
50.0 r i =l
Fro. 6. log Ai is calculated from the asymptotic results (5.2), (5.5), (5.8), and (5.12) as afunctionof/3withr=k=l,N=lCQ,andn=5.
For N sufficiently large,- our asymptotic approximations involving the Ai given by (5.2), (5.5) (5.8), and (5.12) imply the ordering (5.15) should be valid. However, for fixed large N with specific values of a,j3, r, and k, the numerical evaluation of our asymptotic approximation can also yield the ordering (5.16) when fl is sufficiently large. One reason is that as /3becomes large the mean extinction time decreases significantly and the difference between the four models becomes small. This is shown in Fig. 6, where Ai is plotted versus 8. Further, N always appears in the ratio N//I, as is clearly seen in (5.2), (5.5), (5.9, and (5.12). Thus an increase in p is equivalent to a decrease in N, which reduces the importance of the exponent 2N@(O) and increases the importance of the coefficient multiplying the exponential function in Ai. This can cause a change in the ordering of the models. The ordering of the models in (5.15) for N//3>>1 coincides precisely with the ordering of the respective values of u(x) in the four models in the neighborhood of x=0. Hence for N/P>> 1, the size of the fluctuations near extinction (X =0 ‘) is crucial for distinguishing the four models. As B increases, the asymptotic inequality N/P>> 1 is weakened and the size of the random fluctuations near x = k becomes more important. Figures 3- 5 show
CHARLES
112
TIER AND FLOYD
B. HANSON
that the values of A, and A4 are close and an increase in /I can cause the ordering between them to change. The ordering (5.15) can easily be seen from the tangent line approximation (4.12) of the four models, giving
2iWi(0)+[
l_$Qi].
i=l,...,
4.
(5.17)
Here the factors Qi are defined as Qi = 0 for model 1, Q2 = 1 + l/(1 +ak/P) for model 2, Qs = 1 for model 3, and Q4 = l/(1 +cuk/P) for model 4. Clearly,
which gives the ordering (5.15) when Ai, for N>> 1, is dominated by the exponent 2 N@(O). In general, using our results, we can conclude that under similar demographic circumstances, populations in a constant environment have the greatest persistence, as in model 1. The persistence of populations subject to environmental randomness depends upon the manner in which this noise is introduced. Model 2 was the least persistent, which is reasonable, since the noise as measured by u(x) is largest of all the models for each x. The difference in persistence of models 3 and 4 is less clear, since near x = 0 the noise in model 3 is greater than in model 4, but the ordering changes near x= k. The asymptotic results help to distinguish these models in order of increasing persistence where intuition fails.
APPENDIX
1.
THE EXTINCTION
TIME PROBLEM
The extinction problem for a stochastic population is the first passage problem in probability theory (Darling and Siegert [5]; Ricciardi [25]). Let r be the first passage time to the zero population size or the extinction of the population defined on (0, co); then 7(X)=SUp{t~0
(Al.l)
of r is given by
F(x,t)=Pr{7
where p(x, t, y) is the probability
density
function
(Al ,2) for X(t)=y
given that
PERSISTENCE
X(0)=x
IN STOCHASTIC
113
POPULATIONS
and X(t*)
From diffusion theory,p(x,t,y)
P,=;~wPx,+wx)Px~
o
jy=qx-Y),
satisfies (Al .3)
with the proper boundary conditions. The probability density function of r, f(x, t), is then given by (Al .4)
f(x,I)=F;(x,t)=-lga~,(x,t,y)dy with a moment generating function t((x,s)=
s0
mf(x,t)e-“dr.
(Al .5)
Using (A1.4) with (A1.3) we see that u(x, s) satisfies
tfqx)%,
+M(x)u,
=su,
O
(Al .6)
The proper boundary conditions depend upon whether or not the process can exit from the boundary points x= 0 and x = cc. For M(x) and V(x) of the form (1.5), the boundary points x= 0 and x= cc are singular points. These points must be classified using Feller’s theory of singular boundary points [9]. For d(x) > 0 and e(x) bounded, x = 0 is an exit boundary and hence the proper boundary condition is u(O,s)=l.
(Al .7)
We find the boundary point x = cc to be an entrance boundary, and the proper boundary condition is lim u(x,s)
APPENDIX
2.
ASYMPTOTIC
(Al .8)
MATCHING
We now give a sketch of the matching of the asymptotic approximations UB, UL, UT9 and uR. Complete details for the more general problem in a finite domain are in Hanson and Tier [12]. In order for u to be correctly represented by the asymptotic approximations us, uL, uT, and uR, there must be some overlap in the asymptotic validity of each neighboring pair of approximations. The matching consists of equating neighboring approximations in the overlap regions. This not only insures that we have a valid
114
CHARLES
TIER AND FLOYD
B. HANSON
approximation to the solution for all X, but allows us to write all constants in terms of a,, since b, =O. THE MATCHING
OF uR AND uT
In order to match uR and uT we require that uR as x-+k + asymptotically match ur as {-+ + co. Hence as x-+k + , (3.10) with (3.11) and (3.13) gives uR(x)-aR[NT(x-k)]-S’r.
W.1)
For ur. we use the asymptotic form for the standardized parabolic cylinder functions U and I/ from [ 11, for I++ co, U(a,S)-_S-“-0.5exp(-~2/4), V(a, S)-
Ir-
fl
(42.2)
+“-0.5exp( +12/4),
(A2.3)
and we obtain u*({)-I&
“/‘exp({‘/2).
--s’r
The comparison of (A2.4) with (A2.1) and I= A$.( x -k) and aT =uR as given in (3.16). THE MATCHING
(A2.4)
requires that b, = 0
OF uT AND uL
We require that ur as {+ - co match uL as x-+k -. As x+k -, we must retain the most significant parts of the expansion of @ and @’ which appear in the function A defined in (3.11), namely Q(x)-@“(k)(x-k)‘/2=r(xk)‘/2u(k) and a’(x)-@“(k)(x-k)= +r(k-x)/o(k). Therefore (3.10) becomes uL(x)-aJNr(k-x)1-“”
+~~[%(‘=)I””
[
&
1
r/Zexp[ rN(x-k)2/u(k)] r(k-x)/u(k)
* (A2.5)
Since the expansions for U and V given in (A2.2) and (A2.3) are not valid as 14 - co, we must use the connection formulas (see [l]),
U(a,S)=
nV(a,-l) l?(a+0.5)
-sin(aa)U(a,-1)
w4
PERSISTENCE
IN STOCHASTIC
115
POPULATIONS
and V(a,~)=~(a+0.5)cos2(sa)U(a,-S)/a+sin(au)V(u,-S), (fQ.7) so that when I+ - co, (A2.2) and (A2.3) may be used with 1 replaced by -l for -I-++ co. The result, using the fact that b, -0, is
m
uT(.+u,
ITi+;“~:““”
+cos(““) r
s r
The comparison NT = d2rN/o( k) in (3.17).
,S, +].
(AM)
of (A2.8) and (A2.5) with { = Nr(x - k) and yields al.-u~cos(as/r) and b,-u,r/T(s/r) as given
THE MATCHING OF uL AND us
Finally, we require that uL as x-+0 + match uB as .$+ + co. As x+0 + , it follows from (3.10) with (3.11) that
udx)-uL[
q]“‘&exp[
&]“‘“,,[
-s
],
(A2.9)
where we have used a( x)-@(O) - rx/u’(O), W(x)- r/u’(O) as x+0 + and A, = lim,, h(x; N). For the leading term asymptotic expansion of uB, we need the well-known formula for the Whittaker function, as
Wr(k, 5) = W,, 0.5(5)-5ke-f/2
.++cc,
(A2.10)
and the “well-defined” complete asymptotic expansion, w,( k, O-5
-ke +‘I2
as
.++oo,
(A2.11)
which was selected for only having exponentially large terms in its complete expansion. Hence from (3.21) and (3.22), we have (A2.12)
~,(~)--a,E-“/‘e-6+b,~+~/~.
The comparison of (A2.12) and (A2.9) along with (3.17), E=Nex, and (3.19) for Nsr yields u,rA,[
k2NBNT] ‘/I
a,-
b,-a,cos
r(s/r) (
y
which are identical to those in (3.23).
>
’ [ k2NBNT]-s’r
116
CHARLES
TIER AND
FLOYD
B. HANSON
REFERENCES 1 2
M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Nat. Bur. Stand. Math. Series No. 55, Washington, 1964. K. B. Athreya and S. Karlin, On branching processes with random environments I, Ann. Math. Stat. 42:1499-1520 (1971).
3
A. D. Barbour, Quasi-stationary distributions ounces in AppI. Probability 8:296-314 (1976).
4
A. D. Barbour, Second order limit theorem for the Markov branching process in random environments, Stochastic Processes Appl. 4~33-40 (1976). D. A. Darling and A. J. E. Siegert, The fit passage problem for a continuous Markov process, Ann. Math. Stat. 24~624-639 (1953).
5 6 7
in Markov
population
processes,
Ad-
S. N. Ethier and M. F. Norman, An error estimate for the diffusion approximation to the Wright-Fisher model, Proc. Nat. Acad. Sri. U.S.A. 74:50965098 (1977). M. W. Feldman and J. Roughgarden, A population’s stationary distribution and chance of extinction in a stochastic environment with remarks on the theory of species packing, Theoret. Population BioIogv 7: 197-207 (1975).
8
W. Feller, Diffusion processes in genetics, in Proceedings of the Second Berkeley Symposium, Univ. of Calif. Press, Berkeley, 1951, pp. 227-246.
9
W. Feller, The parabolic differential equations and transformations, Ann. Math. 55:468-519 (1952).
10 11 12 13 14 15
the associated
semi-groups
of
J. H. Gillespie, The effects of stochastic environments on allele frequencies in natural populations, Theoret. Population BioIogv 3:241-248 (1972). N. S. Goel and N. Richter-Dyn, Stochastic Models in BioIogv, Academic, New York, 1974. F. B. Hanson and C. Tier, An asymptotic solution of the first passage time problem for singular diffusion in population biology, SIAM J. AppI. Math., 40 (1981). S. Karlin and B. Levikson, Temporal fluctuations in selection intensities: case of small population size, Theoret. Population BioIogv 6:383-412 (1974). N. Keiding, Extinction and exponential growth in random environments, Theoret. Population BioIogv 8:49-63 (1975). N. Keiding, Population growth and branching processes in a random environment, in Proceedings of the 9th Biometric IntemationaI Conference, Invited Pqoers, Vol. II, 1976, pp. 149- 165.
16
T. G. Kurtz, Strong approximation theorems Stochastic Processes AppI. 6:223-240 (1978).
17
C. Lipow, Limiting diffusions for population size dependent branching processes, J. AppI. Probabiliry 14: 14-24 (1977). D. Ludwig, Stochastic PopuIation Theories, Springer, New York, 1974. D. Ludwig, Persistence of dynamical systems under random perturbations, SIAM Reu. 17:605-640 (1975). D. Ludwig, A singular perturbation problem in the theory of population extinction, SIAM-AMS Proc. 10:87- 104 (1976). D. R. McNeil and S. Schach, Central limit analogies for Markov population processes, Proc. Roy. Sot. Stat. Ser. B 35: l-23 (1973). B. J. Matkowsky and 2. Schuss, The exit problem for randomly perturbed dynamical systems, SIAM J. AppI. Math. 33:365-382 (1977). R. M. May, Stabi& and Cornplexiry in Model Ecosystems, Princeton U. P., Princeton, 1973.
18 19 20 21 22 23
for density
dependent
Markov
chains,
PERSISTENCE
IN STOCHASTIC
POPULATIONS
117
24
R. E. O’Malley, Jr., Introduction 1974.
25
L. Ricciardi, Di@sion Processes and Related Topics in Biology, Springer, New York, 1977. M. Turelli, A re-examination of stability in random varying versus deterministic
26
to Singular Perturbations,
Academic,
New York,
environments with comments on the stochastic theory of limiting similarity, Theoret. Population Biology 13:244-267 (1978).
Stochastic Populations*
CHARLES TIER ANDFLOYD B. HANSON Departmeni of Mathematics, Unioersi@ of Illinois at Chicago Circle, Chicago, Iilinois 60680 Receiwd
26 Februaty 1980; revised 28 Ju!v 1980
ABSTRACT Persistence, as measured by time to extinction, is studied in a density dependent population that is subject to small environmental and demographic randomness. Diffusion processes are formally derived from branching processes in constant and random environments. The moment generating function of the extinction time, which satisfies a second order ordinary differential equation, is found asymptotically in the limit of small diffusion and is related to the diffusion limit of the Galton-Watson process and the OmsteinUhlenbeck process. Extinction occurs with probability one, though the mean and variance of the extinction time are found to be exponentially large and suggest the extinction time is exponentially distributed. The notion of persistence is compared with other qualitative measures of stability. Four examples are studied and compared.
1.
INTRODUCTION
Persistence, as measured by time to extinction, has been suggested by Ludwig [19] as a qualitative measure of the effects of random fluctuations on a species. The inclusion of small random fluctuations or stochastic effects may lead to the eventual extinction of a species even when its population size is near a nontrivial stable deterministic equilibrium. The extinction time r is then a random variable which is the first passage time to the zero population size. If extinction occurs with probability one, persistence can be measured by the moments of the extinction time, such as the mean and the variance. We shall study persistence using diffusion theory with population models that are deterministically of the logistic type but include small random fluctuations due to demographic and environmental stochasticity. Hence we shall study the diffusion process X(t) with equations involving the back*Research supported by the National Science Foundation under Grants MCS 78-01167 and MCS 79-01718. MATHEMATICAL
BIOSCIENCES
53:89- 117 (1981)
OElsevier North Holland, Inc., 1981 52 Vanderbilt Ave., New York, NY 10017
89 0025-5564/81/010089+29$02.50
90
CHARLES
TIER AND
FLOYD
B. HANSON
ward operator
fqx)uxx +wx)ux with the proper auxiliary conditions.
Here M(x)
(1.1) and V(x)
are the infinites-
imal mean and variance, respectively, of the diffusion process X(r), and X(0)=x is the initial population size. Demographic stochasticity or sampling fluctuation is due to randomness in the survival ability or fertility of individuals within the population. Typically, models with demographic stochasticity arise from discrete variable stochastic processes such as birth and death processes or branching processes. A diffusion process can arise as a limit approximation of the discrete process. For example, the diffusion limit of the Galton-Watson process has infinitesimal mean and variance (Feller [8]) M(x)=rx,
(1.2)
v(x)=px. Here there are no density dependent
effects, since the per capita birth rate
M(x)/x and variance V(x)/x are constants. Density dependent diffusion models arise when offspring distributions
or
per capita birth and death rates depend upon the current population size. Goel and Richter-Dyn [ 1 l] obtained density dependent diffusion models as formal limits of general birth and death processes. A rigorous proof of the convergence of a continuous time, density dependent branching process to a diffusion process was given by Lipow [ 171 and by Kurtz [ 161 within the context of general continuous time Markov infinitesimal mean and variance of the form
M(x)=-v(x), V(x)=j3x.
chains.
Lipow
derived
an
(l-3)
Diffusion models with (1.3) in general lead to extinction with probability one and a finite mean extinction time if the per capita growth rate p(x) is integrable near x=0. Thus, in this case, the mean extinction time can be used as a measure of persistence. Environmental stochasticity results from random fluctuations in the environment affecting the population as a whole, beyond the effects the environment has on the variance of the offspring distribution. Environmental stochasticity can be introduced into a standard deterministic model by replacing one or more of the parameters with random processes. This type of problem has been analyzed extensively by May [23], Feldman and Roughgarden [7], and Turelli [26]. For example, if Gaussian white noise is
PERSISTENCE
IN STOCHASTIC
91
POPULATIONS
introduced into parameters in the logistic equation a!x/dt =rx( 1 -x/k), a diffusion process is obtained using Ito calculus with infinitesimal mean and variance M(x)=rx(l
k, r>O,
-x/k),
V(x)=x’C(x),
Z(x) > 0,
(1.4)
where t(x) is related to the variance of the noise. For the diffusion model with (1.4), the zero population size is not attainable. Hence either the probability of extinction is equal to one and the mean extinction time is infinite, or the probability of exit is equal to zero. In this case, persistence is not useful, and other qualitative measures of stability have been used, some involving the stationary distribution by May [23] and Turelli [26]. Models with both environmental and demographic stochasticity have been investigated by Athreya and Karlin [2], Keiding [14,15], Barbour [4], and Ludwig [20]. Keiding [IS] and Ludwig [20] formally derive diffusion approximations to discrete time, density dependent branching processes with a random environmental parameter in the mean of the offspring distribution. Their diffusion approximations have an infinitesimal mean and variance M(x)=rx(l
-x/k),
k, r>O,
v(x)=xd(x)+x2c?(x),
ci(x)>O,
5(x)>O,
(1.5)
where x is a scaled population size. In V(x), the term xd(x) arises from the variance of the demographic stochasticity, and the term x’c?(x) arises from the variance of the environmental stochasticity. The deterministic approximation, 2, satisfies the logistic equation di/dt = M(P) = r?( I- 2/k). In Sec. 2, we formally derive a diffusion equation with coefficients (1.5) following Ludwig [20]. Ludwig investigated, asymptotically, the mean extinction time when the environmental stochasticity dominated the demographic stochasticity, i.e., C?(X)= 1/N, N > 1, and g(x) = (Yin (1.5). Analogous models arise in population genetics with the demographic stochasticity, called random genetic drift, due to sampling error and the environmental stochasticity due to random selection (see Gillespie [lo]; Karlin and Levikson [ 131). The extinction time problem is the first passage problem in probability theory. As we derive in Appendix 1, the moment generating function of the extinction time 7, u(x, s), satisfies
~V(x)u,, +M(x)u,
=su,
o
(1.6)
92
CHARLES
TIER AND
FLOYD
B. HANSON
FIG. 1. The interval x > 0 with the four regions B, L, T, and R indicated. solutions us, UL, UT, and uR are given in their respective regions of validity.
The
If M(x) and V(x) are given by (1.5) with d(0) > 0, then x = 0 is an exit boundary and x = cc is an entrance. The analysis of (1.6) subject to the proper boundary conditions is difficult except for certain special cases (Darling and Siegert [5]; Ricciardi [25]; Goel and Richter-Dyn [ 111). In Sec. 3, we present our main results, where we find asymptotic solutions of (1.6) with (1.5) but with d(x)=d(x)/N and g(x)=e(x)/N, N ~1. Here N is some reference population size, and its reciprocal measures the smallness of the variance. The WKB approximation technique and boundary layer analysis of Hanson and Tier 1121, extended to the infinite domain, are used for large N. We find that the solution has different approximations depending upon the value of x (see Fig. 1). For x near the stable deterministic equilibrium k, the process is approximated by the Omstein-Uhlenbeck process to leading order. This is analogous to the approximation of discrete processes about an equilibrium using the Grnstein-Uhlenbeck process by McNeil and Schach [21] and Barbour [3]. Similar results have been derived by Feller [8], by Ludwig [18], and in nonbiological contexts by Matkowsky and Schuss [22]. Near x=0, we find that the process is approximated by the diffusion limit of the Galton-Watson branching process. This is a result of the demographic stochasticity. This approximation and the one near x = k are then connected together using the WKB solution and asymptotic matching. The usual boundary conditions on the diffusion limit of the Galton-Watson process and the Ornstein-Uhlenbeck process are thus replaced by matching conditions. The matching process is described in Appendix 2. Using the asymptotic solution for N>l we find in Sec. 4 that extinction occurs with probability asymptotic to 1. The mean and the variance of the extinction time for x away from 0 are found to be E[7]--h,, Var[ T]-A~, where A,, given by (4.7) is exponentially large in N. This indicates that the population persists for long periods of time before extinction. In addition this is evidence that the extinction time 7 is exponentially distributed for N>l. For x near 0, we find that E[T] and Var[~] change rapidly with x,
PERSISTENCE IN STOCHASTIC POPULATIONS
93
which is expected and is typical of boundary layer behavior. We also compare our results for persistence with the notion of stability used by other authors [23,26] and find that they are qualitatively similar. In Sec. 5, we consider four density dependent models, one with a constant environment and three others that differ by the manner in which the environmental stochasticity is introduced. We are able to conclude that the introduction of environmental stochasticity in models with demographic stochasticity decreases the persistence over models with a constant environment since it decreases the mean extinction time. This is the result of the environmental stochasticity increasing the probability of the population size becoming small where the demographic stochasticity can cause extinction. 2.
THE DIFFUSION
APPROXIMATION
We shall now formally derive diffusion approximations for discrete time, density dependent branching processes in a constant environment and in a random environment following Keiding [15] and Ludwig [20]. Let Z, be a discrete random variable equal to the number of individuals in the population in generation n. We assume that the generations are nonoverlapping and the individuals within each generation reproduce independently. Let the number of offspring born to individual i in generation n be given by the random variable Bi, all of which are independent and identically distributed. CONSTANT
ENVIRONMENT
For the density dependent branching process in a constant environment, the distributions of the Bi are defined by the probabilities pi(z) = Pr( Bi = A-C =z) with mean and variance h(z)=E[B,IZ,=z],
(2.1) v(z)=Var[BilZn
=z],
respectively. Within each generation, the total population Z, is assumed to be fixed and hence the distribution and moments of Bi are conditioned on Z” =z. The number of individuals in generation IZ+ 1 is the sum of the offspring of all the Z,, individuals, and hence
(2.2) Using (2.2) with (2.1), we can compute
the first and second moments
of the
CHARLES
94 increment
TIER AND FLOYD
B. HANSON
AZ, = Z,,, , - Z, as E[AZ,IZ,,=z]=z[h(z)-11, (2.3) E[(AZ,)2jZ,
=z] =zv(z)+z~[~(z)-
112.
To formally obtain a diffusion approximation to the branching process (2.2), we must transform the discrete time, discrete size branching process into a new continuous time, continuous size process. Thus, we define the continuous time variable t=n At, where the increment At<1 is the generation time. A continuous population size variable, as a function of t, is defined by X(t) = Z,, /L, where L is some large reference population size. In addition, we require that small changes in X(t) occur in small time increments At. This requires that the offspring mean h be approximately equal to 1, i.e., near replacement. Thus we assume h(z)=l+B/k(z/L),
(24
where 6 measures the deviation from replacement. We also assume that v(z)=d(z/L). The diffusion approximation X(t) is obtained in the limit as iI-+0 when we choose At = 6 and L = 1/a. Lipow [ 171 proved rigorously the convergence of a continuous time, density dependent process in a constant environment to a diffusion process with a(x) a constant [see (1.3)]. The diffusion approximation can also be made rigorous by using the technique of Ethier and Norman [6]. If we introduce the new continuous variables into (2.3), the infinitesimal mean and variance of the diffusion approximation X(t) are given by
M(x)=
lim Al-+0
E[AX(t)lX(t)=xl =xpFL(x) , At (2.5)
,
At-0
where x= 8z represents a scaled population size. If the higher central moments of Bi are suitably bounded, then the higher moments of AX(t) satisfy E[(AX(t))‘IX(t)=x] At as At+O+.
=O[(At)“‘-‘I,
j>2,
(2.6)
PERSISTENCE
IN STOCHASTIC
95
POPULATIONS
The deterministic approximation to X(t), P, satisfies dz?/dt=fp(2), which is characteristic of deterministic population models. If we assume that the per capita growth rate is p(x) =r( 1 -x/k), then the infinitesimal mean and variance are M(x)=rx(l
-x/k), (2.7)
V(x)=xci(x),
and the deterministic approximation is the logistic equation with scaled carrying capacity k. The diffusion approximation with (2.7) is a stochastic generalization of the logistic equation. This problem with &x)=/3/N is discussed in the examples in Section 5. RANDOM
ENVIRONMENT
We now modify the density dependent branching process described above by allowing random fluctuations in the environment from generation to generation. These fluctuations are accounted for by introducing a random environmental parameter, W,, into the offspring distribution. W, are independent with mean zero, variance 1, and bounded higher moments, and are independent of Z,,, , m < n . The offspring distribution is now defined by pi(z, w)=Pr(B, =jl Z, = Z,W” =w), which are conditioned on both the population size z and the environment w in generation n. Let h(z, w) and Y(Z, w) be the conditional mean and variance, respectively, of the offspring distribution. We also assume the average of Y(Z, w) over the environmental variable W, is
E[vlZ,
=z]=d(z/L),
(2.8)
where a is now the expected value of v. To obtain the diffusion approximation, we assume t = n At, X(t) = Z, /L, the near replacement condition
h(z, w)=
1+ Vl%(z/L)
w+Gp(z/L),
(2.9)
and At = 6= l/L. The fluctuations due to the random environment are assumed to be of order fl . This may be motivated by assuming that these fluctuations are due to the sum of a large number of independent, identically distributed random variables. Then fluctuations about the mean are of the order of one over the square root of the total number in the sum.
96
CHARLES
The infinitesimal
FLOYD
B. HANSON
mean and variance of the diffusion approximation
are computed by iterated expectation M(x)=
TIER AND
lim
At-+0
using (2.3), (2.8)
EIAx(t)lx(t)=xl At
X(t)
and (2.9) as
’
=&moE{EIAZ~]Zn,Wn]]Z,,=x/At,n=t/At}, ‘V(X), v(x)=
(2.10)
lim E[(AX)*IX(t)=x] At At-t0
’
= lim SE{E[(AZ,)*IZ,,W,]IZ,,=x/At,n=t/At}, At-+0
If the higher central moments of W, and Bi are suitably bounded, then the higher moments of AX satisfy a condition of the form (2.6). The diffusion approximation for the constant environment model is a special case of the random environment model with 5(x)=0. Both models have the same mean, since the environmental parameter has mean zero. In V(x) in (2.10), the term x&x) arises from the variance of the offspring distribution, while the term x*P(x) arises from the fluctuations in the environment. The models of May [23], Feldman and Roughgarden [7], and Turelli [26] with random environments but no demographic stochasticity can be obtained by setting &x)=0. Their models were obtained, for example, by replacing with r/k
r in the logistic equation dx/dt=rx(l fixed and y as normalized Brownian
stochastic
differential
-x/k) by r+ ~6 dy/dt, motion. This leads to the
equation dx=rx(l-x/k)dt+&xdy
and a diffusion,
(2.11)
using Ito calculus, with mean and variance M(x)=rx(l V(x)=ax2.
-x/k),
(2.12)
This corresponds to setting d(x) = 0 and t?(x) = a in (2.10). This completes the formal derivation of the diffusion approximation. Our principal interest is the analysis of diffusion equations where V(x) is small relative to M(x). Therefore, we replace d(x) by d( x)/N and C?(X) by .$x)/N in (2.10), which gives (3.1). This scaling assumes that the stochastic effects are small compared to deterministic effects but the environmental fluctuations are the same order as the demographic effects.
PERSISTENCE
3.
IN STOCHASTIC
ASYMPTOTIC
97
POPULATIONS
SOLUTION
OF EXTINCTION
PROBLEM
We shall now construct an asymptotic solution for the moment generating function of the extinction time r, u(x, s; N). As we show in Appendix 1, u satisfies
&X)%x+m(x)u,
o
=su,
(3.1)
where m(x)=rx(l-x/k),
k, r>O,
u(x)=xd(x)+x2e(x),
d(x)>O,
e(x)>O,
(3.2)
and d(x) and e(x) are bounded for x< co. The equation (3.1) is (1.6) where the coefficients (3.2) are the same as (1.5) but with a(x), g(x), and V(x) replaced by d(x)/N,e(x)/N, and u(x)/N, respectively. The boundary conditions for (3.1) with (3.2) as discussed in Appendix 1 are ~(0, s; N)=
1,
(3.3a)
lim u(x, s; N) < 00, x-m
(3.3b)
since x=0 is an exit boundary and x= cc is an entrance boundary. We seek an asymptotic solution valid as N-+oc using the WKB method (O’Malley [24]; Hanson and Tier [12]). To do so, we assume that the solution is of the form X’p(~,s;N)dz
.
(3.4)
Substituting (3.4) into (3.1) yields the Ricatti equation (p2+2_ m(x) --- 2s u(x) ‘- NV(X)
(3.5)
We assume thatcp has an asymptotic expansion of the form p(x,s;
N)=cp,(x,s)+
w
+o(-$),
A’--,
(3.6)
substitute (3.6) into (3.5), and equate coefficients of powers of l/N to zero; then we find to leading order the two solutions cp; =0
and
‘p:=2@‘(x),
(3.7)
98
CHARLES
TIER AND FLOYD
B. HANSON
where (3.8) The next order terms give
cpz=S/m(x>, cpz’= -rp;
(3.9)
-w/w.
Since cp; are not integrable at either x=0 or x= k, we see that the asymptotic validity of the WKB approximation will fail in some small neighborhood of these points, i.e., regions B and Tin Fig. 1. The point x=0 is a singular point, and the point x= k is a turning or transition point. The WKB approximation (3.4) is valid in region L in Fig. 1, and a similar approximation is valid in region R in Fig. 1 but subject to the boundary condition (3.3b). The different solutions in Fig. 1 are defined by
I 1
uB,
u(x, s)w
UL2
UT? uR>
xin xin xin xin
B, L, T, R.
The approximate solutions, uR in R and uL in L, are given by a linear combination of the two solutions of the form (3.4) with (3.7)-(3.9), which we give as
(3.10)
lqx,
N) =
( -&J
exp;;fl;;,x” ,
(3.11)
with m(x) and u(x) given by (3.2). We note that A(x; N) grows exponentially with N, since cP(x)>O. The constants aR,aL, b,, and b, are to be
PERSISTENCE
IN STOCHASTIC
99
POPULATIONS
determined, and
“‘.
(3.12)
Applying the boundary condition (3.3b) to uR given by (3.10), which is valid for x > k, we find that (3.13)
b, -0,
since A(x; N) is always unbounded as x+00 for fixed N. The approximate solution in the transition region T, u&, s; NT), is found by introducing the stretched variable {= Nr(x- k) into (3.1) with NT given by (3.12). Then we obtain to leading order the transition region equation (3.14)
UT,55
The left side of (3.14) is the backward operator of the Ornstein-Uhlenbeck equation, which has been derived in [8], [21], and [3]. The solution of (3.14) is a linear combination of the standardized parabolic cylinder functions [ 11,U(s/r - 0.5, {) and V(s/r - 0.5,{), in the form UT-[
U,U(f
-0.5,
{)+b,V(;
-0.5,
I)] exp( $),
(3.15)
where aT and b, are to be determined. From the theory of asymptotic analysis of differential equations (O’Malley [24]), we see that the transition region approximation UT is valid for at least Ix- kl= O(N -‘/‘). Clearly, UT is not a good approximation to the solution at points that are order one from x = k, such as the boundary point. To obtain consistent approximations to the solution we require that UT be consistent or match with us and uL as N+o3. The procedure of matching us and UT consists of expanding UT for {-+oc and uR as x+k +, and then choosing aR, b,, a,, and bT so that the leading terms of uR and uT are identical. We sketch the matching in Appendix 2; complete details are in Hanson and Tier [12]. From Appendix 2, we find 6, =O, aR =aT.
(3.16)
The constant ar will be determined by the boundary condition at x =O. Similarly, uT as I+ - 00 is matched to uL as x+k -. From Appendix 2, we
100
CHARLES
TIER AND FLOYD
B. HANSON
find aL =a+x~,
b,
r
(3.17)
=T&
In order to extend the results to region B we must construct a new solution, since uL is not valid near x=0. The solution in the boundary region B, u,(.$, s; N,), is found by introducing the stretched variable (3.18)
(=NBx
where (3.19)
Ns -2Nr/u’(O). Substituting .$ for x in (3.1) and expanding as N+cc leading term the boundary layer equation
‘tu B. Et +,I$,,< =
with 5 fixed gives as a
fUB.
(3.20)
We note that the left side of (3.20) is the backward operator for the diffusion limit of the Galton-Watson process [see (1.2)]. The solution of (3.20) satisfies the boundary condition at x = 0 with a matching condition with uL at [= cc instead of the boundary condition associated with (1.2). The general solution of (3.20) is
UB-[UBw~(-~,~)+bBw2(-f,E)lexp(-~), (3.21) where we define
W,(k,x)=Wk,,/z(x) w,(k,x)=cos(ka)Re[
W- k,1,2(-X)]+r(l--)s~2(~~)~~.l/2(x)~ (3.22)
The functions W,,,,,(x) and M h,,,Z(x) are Whittaker functions [l]. The constants aB and bB are to be determined. The boundary region solution uB is defined for 0 < 4 < co, and we require that it satisfy the boundary condition at x= 0 and that uB as &cc match uL The size of region B is at least x= 0(1/N). The matching, as x+0+.
PERSISTENCE
IN STOCHASTIC
sketched in Appendix
POPULATIONS
101
2, gives a,rA,(k2NBNT)S’r
a,=
’
Vs/r) b, =a+os
(
7
>
(3.23)
(k2NBNT)-S’r,
where A,=
N)=d(O)(
liiA(x;
r~~(k))“2ex~~2Nm(0)].
Using (3.23) in (3.21) at x = 0 with the results that W,( k, 0) = l/I’( 1 -k) and W2( k, 0) = (cos ka)/r( 1 + k), the absorbing boundary condition u,(O, s) = 1 gives 1
a T=
.
r( k2NBNT)s’rh0
(3.24)
+ cos2(m/r)(k2NBNT)-S”
r(s/r)r(i +s/T)
qi -+I
We summarize the solution with aT given by (3.24) and A(x; N) by (3.11) in the different regions in Fig. 1 as: Region R: (3.25) Region T: u,-a,(l(S-0.5,l)exp($),
(3.26)
where
{=&(x-k),
NT
= [2rN/o( k)]“‘;
Region L:
s/r N,k;x-k, +qx.
,
N)
(NTklx-kl/x)“” ; r(s/r)
1
(3.27)
102
CHARLES
TIER AND FLOYD
B. HANSON
Region B:
r(
( )i 6
- 2
u,-a,exp
k2NBNJ”AOW,( -s/r,
5)
Us/r)
(7 >(k2NBNp’w*(-s/r,
+cos
5)
1 )
(3.28)
where W, and W, are given in (3.22) and N, =2rN/o’(O).
5=N*x, 4.
MOMENTS
OF EXTINCTION
TIME AND PERSISTENCE
Since u(x, s) is the moment generating function of r, its moments are given by E[ 7”]
qYf(
x,r)dt=(-l)~~u(x,s)~,_,’
(4.1)
where f(x, t) is the probaslity density function of r discussed in Appendix f (x, t ) dt = u( x, 0) represents the ultimate proba1. The zeroth moment / bility of extinction and ii found by evaluating the limit as s-0 in us, uL, tlr, and uR, given by (3.25)-(3.28). Using 1 ~ I(s/r)
=
s+o(s*), r
(4.2)
as s+O, we find that lima,=I(l)= s+Q
1
(4.3)
uL, f+, u,-1.
(4.4)
and hence limu,,
s-0
Hence r is a proper random variable, and the probability of extinction is asymptotically equal to 1 for all x.
PERSISTENCE
IN STOCHASTIC
103
POPULATIONS
If we denote the mean and variance of r for x in region T by Er[ r] and Var,[ r], respectively, then from (4.1)
Var,[r]=
a2
lim
_?ur(x,s)-
s-0+
[
(
as
11 2
&ur(x,s)
Differentiating (3.26) with respect to S, letting s+O+,
(4.5)
.
and using (4.3) and
aa, Jl-i+ T”-Ao, we find after a short calculation that
(4.6) Here A0 is exponentially large and is given by
A0
exp(2NQ(O)). =dP)( r1&))“2
(4.7)
A similar calculation for the mean and variance of r in regions 8, L, and R gives
EB[~]-Ao(l-e-E),
N+Co,
,
(4.8)
and Var,[r]-Var,[r]-At, Var,[r]-A:(1
-eezE),
I
N-XCJ,
(4.9)
where A0 is given by (4.7) and .$= 2rNx/u’(O). In region L we have neglected A(x; N) as exponentially small compared to A,. The moments of r can also be computed by a modification of the method of Matkowsky and Schuss [22] as described in Hanson and Tier [12]. It is interesting to note that the means and variances of r for x in L, T, and R are identical and independent of x. Thus, except for x in B, E[ r]-A,, and Var[r]-11; are characteristic of an exponential distribution with parameter l/A,. Here the mean and variance are exponentially large as N+m, since A,, is exponentially large. This can be interpreted by the fact
104
CHARLES
TIER AND
FLOYD
B. HANSON
that the population is strongly attracted to the stable equilibrium at x = k in the presence of small random fluctuations and hence lingers near k. In addition the mean and the variance in regions L, T, and R are independent of x, and we are unable to distinguish populations that start there. Eventually, the random fluctuations cause the population size to become small, so that the demographic effects are important and can lead to extinction. For initial values near x=0 (i.e., in region B), the mean and the variance of r, given by (4.8) and (4.9) are very sensitive to changes in x. The mean extinction time increases rapidly with x from the value EB[7] =0 at x=0, where absorption is certain. Thus, the larger the initial population size, the longer the mean extinction time. From (4.8) with ,$=2rNx/o’(O) we see that a key factor in measuring the persistence of a population with a small initial size is the ratio rN/u’(O). This quantity represents the ratio of the growth rate r to the rate of change of the demographic variance at x=0. The larger the value of r or the smaller the value of u’(O), the greater the persistence of the population for large fixed N. Other authors (May [23], Feldman and Roughgarden [7], Turelli [26]) have examined the “stability” of populations in a random environment in the absence of demographic effects which corresponds to d(x)=0 in (3.2). With d( x)=0 in (3.2), extinction is not possible, since the zero population size is not attainable (x = 0 is an entrance, or a natural boundary in the terminology of Feller [9]). A qualitative measure of the effects of the random environment is whether a stationary distribution exists which depends upon the behavior near x= 0. In the case when no stationary distribution exists, the mean extinction time is infinite and hence it is not a good qualitative measure. The stationary distribution li, if it exists, is found using the forward equation and is
where C is chosen to normalize the distribution, u(x) =x*e(x) > 0, and @J(X) is given by (3.8). To distinguish populations in a random environment that have stationary distributions, May [23] and Turelli [26] have measured the spread or dispersion of the stationary distribution about the deterministic equilibrium k, which is the mean of the stationary distribution if the noise is small. Turelli’s account [26, pp. 250,264] of May’s “dispersion”criterion [23] is that the square of the coefficient of variation, variance/(mean)2, of the stationary distribution (4.10) locally about x = k must be bounded by some constant (e.g. 1). By expanding O(x)= -m’(k)(x-k)*/2u(k) in (4.10), we find that their criterion gives -m’(k)>o(k)/2Nk*.
(4.11)
PERSISTENCE
IN STOCHASTIC
105
POPULATIONS
It is important in their criterion that the noise be small; otherwise the stationary distribution may not be peaked about k. In our problem we have explicitly given the size of the noise as l/N, which is small when N> 1. Then the criterion (4.11) is always satisfied and is not particularly useful in distinguishing populations. May [23] and Turelli [26] developed an approximation method that could be applied directly to the forward equation, which gives the criterion (4.11). Their method is related to the approximation (3.14) in the transition region T and the Omstein-Uhlenbeck process. The problem we study, (3.1) with (3.2), contains the added aspect of demographic stochasticity, which is more realistic. With the demographic stochasticity included, extinction is possible (x=0 is an exit boundary) and the stationary distribution is identically 0. As suggested by Ludwig [19], persistence (using the extinction time) provides an absolute criterion for measuring the qualitative effects of the randomness. The importance of the asymptotic result for the mean extinction time (4.6) derived in this section is that it provides a measure of persistence and is valid for all x. It was derived with the assumption of small noise, N>>l, consistent with May [23] and Turelli [26]. The quantity A, given by (4.7) distinguishes different mean extinction times, and for N> 1 it is dominated bykthe exponent 2N@(O) in the exponential function. The integral Q(O) = s [ m(x)/u(x)] dx depends upon 0
the behavior of m(x)/u(x) over the entire interval [0, k] and not just near x- k. While it is better to calculate the integral a(O), a simplified approximation of Q(O) can be found using tangent line approximations at x =0 for XE(O, k/2) and at x=k for xE(k/2, k) to give m(x)k ___ 20(x) >>1.
+ (m/u)‘(x)@ 8
]-::;:y
) (4.12)
This formula gives an asymptotic analogy of the condition (4.11) derived by the criterion of May, since 2NcP(O)> 1 when N>> 1. It represents an improvement over (4.1 l), since it uses behavior not only near x =k but also near x=0. We shall use (4.12) in the examples in Sec. 5. 5.
EXAMPLES
We now apply the results in Section 4 to study the persistence of four density dependent stochastic population models. All four models have the logistic equation as an infinitesimal mean, but they differ in the manner in which the environment is taken into account.
106
CHARLES
TIER AND FLOYD
B. HANSON
Model I (Constant Environment). The diffusion approximation in a constant environment is characterized by (3.1) with m(x)=rx(l-x/k),
(5.1)
u(x)=Px,
i.e., e(x)=0 in (3.2). This type of model is described by Goel and RichterDyn [ 111, Lipow [ 171, and Hanson and Tier [ 121. Using (4.6)-(4.9), we find that
E[r] -(
$$)“‘exp(
$k)=Ai,\
Var[ r] -A:,
x inL, T, R,
(5.2)
J
and E,[T] -A,[
I-exp(
?)I, x inB,
Var,[T] -A:[
1 - exp( ?)I.
(5.3)
1 Model 2
The diffusion approximation is characterized by m(x)=rx(l-x/k), (5.4)
u(x)=px+cXx*,
i.e., e(x)=cx, a constant, in (3.2). Models with environmental randomness as reflected in the term ox* were discussed in [23], [7], and [26] without demographic randomness and can be interpreted as the result of a Gaussian white noise perturbation of r with r/k fixed using Ito calculus. This model including demographic randomness was studied in [20] and [15]. Using (4.6)-(4.9), we find that for x in L, T, or R
E[rl
-
Var[ ~1 -A;, where
TP Nkr3(l +ak//3)
1 I/*
exp[2N@2(0)]=A2,
(5.5)
PERSISTENCE IN STOCHASTICPOPULATIONS
107
while for x in B E,[T]
Var,[7] Model 3
-A2[1-exp(-2A+x/P)19
(5.6)
-AQl-exp(-4NrxlP)l.
The diffusion approximation is characterized by m(x)=rx(l-x/k),
(5.7)
o(x)=Px+(yx2(1-x/k)‘.
Here e(x) = o( 1 - x/k)2 in (3.2) can be interpreted as the result of introducing Gaussian white noise in r with k fixed. This model was discussed by Keiding [ 151. Using (4.6)-(4.9) we find that for x in L, T, or R E[ 71 -[ ?rp/Nkr3]“2exp[2N@3(0)]
=A,,
Var[ 7 ] -A:,
(5.8)
where
For x in B,
-A3[l-exp(-2~rx/P)1,
EB[7]
Var,[r] Model 4
-A’,[l-q$--4Nrx/kQl.
(5.10)
The diffusion approximation is characterized by
m(x)=rx(l
-x/k),
v(x)=px+ax4/k2.
(5.11)
This model was mentioned by Keiding [15]. Without demographic stochasticity (/3= 0), this diffusion was studied by Feldman and Roughgarden [7] and Turelli [26], and can be interpreted as the result of introducing Gaussian white noise in r/k with r fixed. Using (4.6)-(4.9), we find that for x in L, T, or R,
E[rl-
4 Nr3k(l
Var[ 71 -A’4
+ak/P)
1 l/2
exp[2NQ4(0)]=A4,
(5.12)
CHARLES
108
TIER AND FLOYD
B. HANSON
where (5.13) For x in B, E8[7] Var,[r]
-A4[l-exp(-2A+x/P)19
(5.14)
-A24[1-exp(-4Nrx/P)le
We can now compare the persistence of the four models described above as measured by the mean extinction time E [ T]. First we note that if cx= 0, models 2,3, and 4 reduce to model 1 since the environmental randomness enters the models only through the term containing CLThe parameter cx is the same in models 2,3, and 4 by our scaling and measures the variance of the environmental noise as in (2.12). In Fig. 2, we sketch E[T] versus x for N>>l. The graphs were sketched qualitatively using EB[7] in (5.3), (5.6), (5. lo), and (5.14), since these results are uniformly valid for x in all regions. The rapid change in the mean as x varies through region B is illustrated by enlarging region B as in Fig. 2. In addition, we obtain an ordering of the models in terms of persistence. This ordering is the same as the ordering of
i =l
@I
I I I i=4
p----_ A---l ’
0
i=2
X
k
FIG. 2. The mean extinction time I?[ T] for models 1 to 4 sketched qualitatively from the asymptotic results (5.3), (5.6), (5.10) and (5.14) for N>l. The regions B, L, T, and R are given.
PERSISTENCE
IN STOCHASTIC
109
POPULATIONS
5O.Or
0.04 0.0
5.0
10.0
Q Fro. 3. log Ai is calculated from the asymptotic a function of (I with r=k= 1, N= 100, and fi= 1.
results (5.2),(5.5),(5.8),
and (5.12) aa
the terms Ai, i= 1, . . . ,4. For the case in Fig. 2, we have A,>A,>A,>A,.
(5.15)
This ordering can be shown to be true if N> 1, or more precisely, N//3>>1, since @i(O) = rk/2/3 > ad(O) > aa, > Qz(0). In Figs. 3-5, we graph Ai versus a for different values of p. Here hi, i= 1, . . . ,4, is used to distinguish the mean extinction times of the four models using (5.2), (5.5), (5.9, and (5.12). For these figures, we have chosen r= k= 1 and N= 100. We see from Figs. 3-5 that Ai decreases, and hence the mean extinction time decreases as a increases. As a increases, the size of the random fluctuations increases, which in turn increases the chance that the sample paths of the population reach extinction. As B increases, but with the asymptotic results still valid, the ordering A,>A,>A,>A,
(5.16)
is possible for certain ranges of a. Figure 3 with fi= 1 has the ordering (5.15). In Fig. 4 with /3=5, the ordering changes from (5.15) to (5.16) at about a= 6.7. When p= 10 as in Fig. 5, the ordering is given by (5.16) for all a. We note the scale of the ordinate is different in Figs. 3-5.
110
CHARLES
TIER AND FLQYD
B. HANSON
i=l
6.0-
FIG. 4. log Ai is calculated from the asymptotic results (5.2), (5.5), (5.8), and (5.12) as afunctionofawithr=k=l,N=lOO,and~=5.
0.01 0.0
5.0
10.0
Q FIG. 5. log Ai is calculated from the asymptotic results (5.2), (5.5), (5.8), and (5.12) as afunctionofawith/Pk=l,N=100,andB=10.
PERSISTENCE
IN STOCHASTIC
POPULATIONS
111
50.0 r i =l
Fro. 6. log Ai is calculated from the asymptotic results (5.2), (5.5), (5.8), and (5.12) as afunctionof/3withr=k=l,N=lCQ,andn=5.
For N sufficiently large,- our asymptotic approximations involving the Ai given by (5.2), (5.5) (5.8), and (5.12) imply the ordering (5.15) should be valid. However, for fixed large N with specific values of a,j3, r, and k, the numerical evaluation of our asymptotic approximation can also yield the ordering (5.16) when fl is sufficiently large. One reason is that as /3becomes large the mean extinction time decreases significantly and the difference between the four models becomes small. This is shown in Fig. 6, where Ai is plotted versus 8. Further, N always appears in the ratio N//I, as is clearly seen in (5.2), (5.5), (5.9, and (5.12). Thus an increase in p is equivalent to a decrease in N, which reduces the importance of the exponent 2N@(O) and increases the importance of the coefficient multiplying the exponential function in Ai. This can cause a change in the ordering of the models. The ordering of the models in (5.15) for N//3>>1 coincides precisely with the ordering of the respective values of u(x) in the four models in the neighborhood of x=0. Hence for N/P>> 1, the size of the fluctuations near extinction (X =0 ‘) is crucial for distinguishing the four models. As B increases, the asymptotic inequality N/P>> 1 is weakened and the size of the random fluctuations near x = k becomes more important. Figures 3- 5 show
CHARLES
112
TIER AND FLOYD
B. HANSON
that the values of A, and A4 are close and an increase in /I can cause the ordering between them to change. The ordering (5.15) can easily be seen from the tangent line approximation (4.12) of the four models, giving
2iWi(0)+[
l_$Qi].
i=l,...,
4.
(5.17)
Here the factors Qi are defined as Qi = 0 for model 1, Q2 = 1 + l/(1 +ak/P) for model 2, Qs = 1 for model 3, and Q4 = l/(1 +cuk/P) for model 4. Clearly,
which gives the ordering (5.15) when Ai, for N>> 1, is dominated by the exponent 2 N@(O). In general, using our results, we can conclude that under similar demographic circumstances, populations in a constant environment have the greatest persistence, as in model 1. The persistence of populations subject to environmental randomness depends upon the manner in which this noise is introduced. Model 2 was the least persistent, which is reasonable, since the noise as measured by u(x) is largest of all the models for each x. The difference in persistence of models 3 and 4 is less clear, since near x = 0 the noise in model 3 is greater than in model 4, but the ordering changes near x= k. The asymptotic results help to distinguish these models in order of increasing persistence where intuition fails.
APPENDIX
1.
THE EXTINCTION
TIME PROBLEM
The extinction problem for a stochastic population is the first passage problem in probability theory (Darling and Siegert [5]; Ricciardi [25]). Let r be the first passage time to the zero population size or the extinction of the population defined on (0, co); then 7(X)=SUp{t~0
(Al.l)
of r is given by
F(x,t)=Pr{7
where p(x, t, y) is the probability
density
function
(Al ,2) for X(t)=y
given that
PERSISTENCE
X(0)=x
IN STOCHASTIC
113
POPULATIONS
and X(t*)
From diffusion theory,p(x,t,y)
P,=;~wPx,+wx)Px~
o
jy=qx-Y),
satisfies (Al .3)
with the proper boundary conditions. The probability density function of r, f(x, t), is then given by (Al .4)
f(x,I)=F;(x,t)=-lga~,(x,t,y)dy with a moment generating function t((x,s)=
s0
mf(x,t)e-“dr.
(Al .5)
Using (A1.4) with (A1.3) we see that u(x, s) satisfies
tfqx)%,
+M(x)u,
=su,
O
(Al .6)
The proper boundary conditions depend upon whether or not the process can exit from the boundary points x= 0 and x = cc. For M(x) and V(x) of the form (1.5), the boundary points x= 0 and x= cc are singular points. These points must be classified using Feller’s theory of singular boundary points [9]. For d(x) > 0 and e(x) bounded, x = 0 is an exit boundary and hence the proper boundary condition is u(O,s)=l.
(Al .7)
We find the boundary point x = cc to be an entrance boundary, and the proper boundary condition is lim u(x,s)
APPENDIX
2.
ASYMPTOTIC
(Al .8)
MATCHING
We now give a sketch of the matching of the asymptotic approximations UB, UL, UT9 and uR. Complete details for the more general problem in a finite domain are in Hanson and Tier [12]. In order for u to be correctly represented by the asymptotic approximations us, uL, uT, and uR, there must be some overlap in the asymptotic validity of each neighboring pair of approximations. The matching consists of equating neighboring approximations in the overlap regions. This not only insures that we have a valid
114
CHARLES
TIER AND FLOYD
B. HANSON
approximation to the solution for all X, but allows us to write all constants in terms of a,, since b, =O. THE MATCHING
OF uR AND uT
In order to match uR and uT we require that uR as x-+k + asymptotically match ur as {-+ + co. Hence as x-+k + , (3.10) with (3.11) and (3.13) gives uR(x)-aR[NT(x-k)]-S’r.
W.1)
For ur. we use the asymptotic form for the standardized parabolic cylinder functions U and I/ from [ 11, for I++ co, U(a,S)-_S-“-0.5exp(-~2/4), V(a, S)-
Ir-
fl
(42.2)
+“-0.5exp( +12/4),
(A2.3)
and we obtain u*({)-I&
“/‘exp({‘/2).
--s’r
The comparison of (A2.4) with (A2.1) and I= A$.( x -k) and aT =uR as given in (3.16). THE MATCHING
(A2.4)
requires that b, = 0
OF uT AND uL
We require that ur as {+ - co match uL as x-+k -. As x+k -, we must retain the most significant parts of the expansion of @ and @’ which appear in the function A defined in (3.11), namely Q(x)-@“(k)(x-k)‘/2=r(xk)‘/2u(k) and a’(x)-@“(k)(x-k)= +r(k-x)/o(k). Therefore (3.10) becomes uL(x)-aJNr(k-x)1-“”
+~~[%(‘=)I””
[
&
1
r/Zexp[ rN(x-k)2/u(k)] r(k-x)/u(k)
* (A2.5)
Since the expansions for U and V given in (A2.2) and (A2.3) are not valid as 14 - co, we must use the connection formulas (see [l]),
U(a,S)=
nV(a,-l) l?(a+0.5)
-sin(aa)U(a,-1)
w4
PERSISTENCE
IN STOCHASTIC
115
POPULATIONS
and V(a,~)=~(a+0.5)cos2(sa)U(a,-S)/a+sin(au)V(u,-S), (fQ.7) so that when I+ - co, (A2.2) and (A2.3) may be used with 1 replaced by -l for -I-++ co. The result, using the fact that b, -0, is
m
uT(.+u,
ITi+;“~:““”
+cos(““) r
s r
The comparison NT = d2rN/o( k) in (3.17).
,S, +].
(AM)
of (A2.8) and (A2.5) with { = Nr(x - k) and yields al.-u~cos(as/r) and b,-u,r/T(s/r) as given
THE MATCHING OF uL AND us
Finally, we require that uL as x-+0 + match uB as .$+ + co. As x+0 + , it follows from (3.10) with (3.11) that
udx)-uL[
q]“‘&exp[
&]“‘“,,[
-s
],
(A2.9)
where we have used a( x)-@(O) - rx/u’(O), W(x)- r/u’(O) as x+0 + and A, = lim,, h(x; N). For the leading term asymptotic expansion of uB, we need the well-known formula for the Whittaker function, as
Wr(k, 5) = W,, 0.5(5)-5ke-f/2
.++cc,
(A2.10)
and the “well-defined” complete asymptotic expansion, w,( k, O-5
-ke +‘I2
as
.++oo,
(A2.11)
which was selected for only having exponentially large terms in its complete expansion. Hence from (3.21) and (3.22), we have (A2.12)
~,(~)--a,E-“/‘e-6+b,~+~/~.
The comparison of (A2.12) and (A2.9) along with (3.17), E=Nex, and (3.19) for Nsr yields u,rA,[
k2NBNT] ‘/I
a,-
b,-a,cos
r(s/r) (
y
which are identical to those in (3.23).
>
’ [ k2NBNT]-s’r
116
CHARLES
TIER AND
FLOYD
B. HANSON
REFERENCES 1 2
M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Nat. Bur. Stand. Math. Series No. 55, Washington, 1964. K. B. Athreya and S. Karlin, On branching processes with random environments I, Ann. Math. Stat. 42:1499-1520 (1971).
3
A. D. Barbour, Quasi-stationary distributions ounces in AppI. Probability 8:296-314 (1976).
4
A. D. Barbour, Second order limit theorem for the Markov branching process in random environments, Stochastic Processes Appl. 4~33-40 (1976). D. A. Darling and A. J. E. Siegert, The fit passage problem for a continuous Markov process, Ann. Math. Stat. 24~624-639 (1953).
5 6 7
in Markov
population
processes,
Ad-
S. N. Ethier and M. F. Norman, An error estimate for the diffusion approximation to the Wright-Fisher model, Proc. Nat. Acad. Sri. U.S.A. 74:50965098 (1977). M. W. Feldman and J. Roughgarden, A population’s stationary distribution and chance of extinction in a stochastic environment with remarks on the theory of species packing, Theoret. Population BioIogv 7: 197-207 (1975).
8
W. Feller, Diffusion processes in genetics, in Proceedings of the Second Berkeley Symposium, Univ. of Calif. Press, Berkeley, 1951, pp. 227-246.
9
W. Feller, The parabolic differential equations and transformations, Ann. Math. 55:468-519 (1952).
10 11 12 13 14 15
the associated
semi-groups
of
J. H. Gillespie, The effects of stochastic environments on allele frequencies in natural populations, Theoret. Population BioIogv 3:241-248 (1972). N. S. Goel and N. Richter-Dyn, Stochastic Models in BioIogv, Academic, New York, 1974. F. B. Hanson and C. Tier, An asymptotic solution of the first passage time problem for singular diffusion in population biology, SIAM J. AppI. Math., 40 (1981). S. Karlin and B. Levikson, Temporal fluctuations in selection intensities: case of small population size, Theoret. Population BioIogv 6:383-412 (1974). N. Keiding, Extinction and exponential growth in random environments, Theoret. Population BioIogv 8:49-63 (1975). N. Keiding, Population growth and branching processes in a random environment, in Proceedings of the 9th Biometric IntemationaI Conference, Invited Pqoers, Vol. II, 1976, pp. 149- 165.
16
T. G. Kurtz, Strong approximation theorems Stochastic Processes AppI. 6:223-240 (1978).
17
C. Lipow, Limiting diffusions for population size dependent branching processes, J. AppI. Probabiliry 14: 14-24 (1977). D. Ludwig, Stochastic PopuIation Theories, Springer, New York, 1974. D. Ludwig, Persistence of dynamical systems under random perturbations, SIAM Reu. 17:605-640 (1975). D. Ludwig, A singular perturbation problem in the theory of population extinction, SIAM-AMS Proc. 10:87- 104 (1976). D. R. McNeil and S. Schach, Central limit analogies for Markov population processes, Proc. Roy. Sot. Stat. Ser. B 35: l-23 (1973). B. J. Matkowsky and 2. Schuss, The exit problem for randomly perturbed dynamical systems, SIAM J. AppI. Math. 33:365-382 (1977). R. M. May, Stabi& and Cornplexiry in Model Ecosystems, Princeton U. P., Princeton, 1973.
18 19 20 21 22 23
for density
dependent
Markov
chains,
PERSISTENCE
IN STOCHASTIC
POPULATIONS
117
24
R. E. O’Malley, Jr., Introduction 1974.
25
L. Ricciardi, Di@sion Processes and Related Topics in Biology, Springer, New York, 1977. M. Turelli, A re-examination of stability in random varying versus deterministic
26
to Singular Perturbations,
Academic,
New York,
environments with comments on the stochastic theory of limiting similarity, Theoret. Population Biology 13:244-267 (1978).