\documentclass[12pt,twoside]{article}
\usepackage{expdist}
\usepackage{epsf}
%TwoSide:
% Modified JMB Web copy Aug 1997 to correspond to galleys.
% Modified JMB Submit May 1997 to incororate figs.
% Revised JMB Submit January 1997
\pagestyle{myheadings}
\markboth {\rm F. B. Hanson, H. C. Tuckwell}{\rm Population growth with 
randomly distributed jumps}
%
\begin{document}
\thispagestyle{empty}
\NODENT
{\Large\bf Population growth with randomly distributed jumps}%%
\footnote{This paper is to appear in {\it J. Math. Biol.}, vol. 36 (2),
pp. 169-187, December 1997.}
~\\*[1em]
{\large\bf Floyd B. Hanson$^1$,  Henry C. Tuckwell$^2$}\\*[0.1in]
{\small $^1$Department of Mathematics, Statistics, and Computer Science
(M/C 249), University of Illinois at Chicago, 851 S. Morgan, Chicago, IL
60607-7045, USA\\e-mail:~ hanson@uic.edu}\\
{\small $^2$Epidemiologie et Sciences de l'Information,
Universite Paris 6 (INSERM U444),
Faculte de Mede\-cine St Antoine,
27 rue Chaligny, 75571 Paris Cedex 12,
FRANCE\\e-mail:~ tuckwell@b3e.jussieu.fr}\\*[1.2em]
%
Received 21 June 1996; received in revised form 17 February 1997\\*[2.5em]
%
{\bf Abstract.}
The growth of populations with continuous deterministic and random jump
components is treated.
Three special models in which random jumps occur at the time of
events of a Poisson process and admit formal explicit solutions are considered:
A)  Logistic growth with random disasters having exponentially
distributed amplitudes;
B)  Logistic growth with random disasters causing
the removal of a uniformly distributed fraction of the population
size;  and  C)  Exponential decay with sudden increases (bonanzas)
in the population and with each increase being an exponentially
distributed fraction of the current population.  Asymptotic and
numerical methods are employed to determine the mean extinction
time for the population, qualitatively and quantitatively.
For Model~A, this time becomes exponentially large
as the carrying capacity becomes much larger than the mean
disaster size.  Implications for colonizing species for Model~A are discussed.  
For Model~B and C, the practical notion of
a small, but positive, effective extinction level is chosen, and in these
cases the expected extinction time rises rapidly with population size, yet
at less than an exponentially large order.
\\*[1em]
{\bf Key words:}  Population extinction -- Randomly distributed jumps -- 
Bonanzas~~-- Disasters \\*[0.2in]
\section{Introduction}
%
We are interested in determining the time to extinction of a
population which is subjected to randomly occurring disasters
when these disasters have a continuous distribution of sizes.
In the time intervals between disasters the population evolves
deterministically.
\REDENT

Previously we have analyzed population growth which is
logistic and suffers disasters at the times of occurrence of
events in a simple Poisson process.  In the first case \cite{HT78},
the disasters were of constant size.  
In a subsequent investigation \cite{HT81},
the disasters were density independent in the sense that they were
a constant multiple of the current population size.  In both cases
the expected time to extinction satisfied a singular
differential-difference equation.
These equations were solved by the method of steps modified by a singular
decomposition for the equations associated with biological extinction.  
This singular decomposition was rigorously analyzed in
\cite{Hanson83} within the framework of a general discrete jump model.
It was shown that the singular decomposition
makes the singularities removable by transforming the
original singular boundary value problem to two nonsingular problems, so that
numerical procedures were well conditioned.
The extinction time grew exponentially with the carrying capacity of the
environment.

Hanson and co-workers \cite{RHMB85,RHJMB86,CJFAS87,NRM88} have applied this
Poisson disaster model to the optimal harvesting in an uncertain environment
with a focus on fisheries applications.
In \cite{RHMB85}, Ryan and Hanson treated the problem of exponential growth in
the presence of random disasters and bonanzas.
In \cite{RHJMB86}, they considered a logistically growing population with
density independent (i.e., linear in population size) random disasters.
Optimal harvesting in the presence of density dependent (i.e., nonlinear in
amplitude) disasters is studied in \cite{NRM88}.
A multispecies model of a Lake Michigan predator-prey system was presented and
computationally analyzed in \cite{CJFAS87} using fishery catch data 
for alewives undergoing occasional massive mortalities to estimate disaster 
parameters.

In most previously studied stochastic population models it was
assumed that populations change size continuously, as for example in
diffusion processes.  Tuckwell and Koziol \cite{TuckKoz87} have reviewed the
logistic model with diffusion type noise, and have also modeled a stochastic
logistic process with dispersal according to either Bose-Einstein or
Maxwell-Boltzmann statistics.

However, all real populations experience 
catastrophic disasters or bountiful bonanzas of various magnitudes, ranging in
size from very small to large.  Adverse phenomena can have many causes:
floods, tsunami, sea level changes, earthquakes, volcanoes, impacts of
extraterrestrial objects, plate tectonics, climate changes, famines,
predation, competition, war, fires, accidents, infectious diseases and
storms or other adverse weather effects.  Examples of favorable effects are
bumper harvests, destruction of predators, destruction of competitors,
sudden abundance of prey or other food, sudden availability of colonizable
patches, favorable ocean currents and beneficial weather.

Flessa \cite{FLES90} lists several contrasting
interpretations of the fossil record for the cause or causes of
mass extinctions that occurred at the Cretaceous-Tertiary boundary.
Jablonski \cite{JABL86} discusses the plausibility of many potential causes
of mass extinctions.

Ludwig et al.\cite{LHW93} strongly emphasize the importance of
uncertainty and how it can mask overexploitation effects, leading
to unwise exploitation policies, especially those based upon sustainable
yield theory.

Previously, Tsurui and Osaki \cite{TO76} have treated the first
passage time to exceed a certain threshold in the case of exponentially
distributed, but positive Poisson amplitudes
(i.e., bonanzas exponentially distributed in both time and amplitude) in a
background of deterministic exponential decay.  Mangel and Tier
\cite{MT93a,MT93b}
derive a birth and death model with catastrophes, presenting solution by
both asymptotic and numerical methods.  
In \cite{MT93a}, they emphasize a simple direct method to compute
persistence times for birth and death models with catastrophes.
Pakes et al. \cite{PTB79} modeled the problem of mass
emigration of an exponentially growing population with randomly distributed
emigration occurring according to a population size dependent jump process.
In a related paper, Trajstman \cite{Traj81} formulated the problem of growth
with emigration more succinctly as a nonstationary compound 
Poisson process superimposed upon a linear bounded deterministic component.
Brockwell el al. \cite{Brock82} treat the births, immigration and
catastrophes as discrete state-space Markov models using distributions that
result in explicit solutions, with comparisons to simple {\it semi-stochastic}
continuous state-space linear models.
In a subsequent paper \cite{Brock83}, they treat catastrophes in
continuous state-space in the present of upward growth according to a compound
Poisson process with an arbitrary distribution of the upward jumps using the
methods of their discrete state-space paper \cite{Brock82}.
Brockwell \cite{Brock85} presents another
variation on the emigration (catastrophe or disaster) model using a birth,
death and catastrophe process, while also treating a corresponding jump
diffusion (stochastic diffusion with Poisson distributed disasters) model.
Bartoszynski et al. \cite{BBCP89} analyze disaster and bonanza models based
upon Markov branching processes that include the birth, death and
population size independent Poisson disaster or bonanza processes.

An element of realism is introduced if multiple disasters or bonanzas are
incorporated in deterministic models of population growth.  
Models will be examined for a population which grows logistically
except when it undergoes sudden changes at times of events of a
Poisson process and whose sizes are random with continuous distributions.
Three special cases are presented in Sect. 3 in
which formal explicit solutions for the extinction time are
obtained.  In the first case disasters occur whose magnitude is
exponentially distributed, but is independent of the population size.
Also, the asymptotic behavior is given for large values of the limiting
population size. 
In the second case the disaster size is density independent
(i.e. proportional to the population size) and uniformly
distributed.  Asymptotic solutions are given for large initial population size.
In the third case the exponentially distributed,
density-independent bonanza is treated.
%
\section{The general models}
%
A population that starts out at size  $N(0)=x$  and grows
logistically to $N(t)$ at time $t$ satisfies the deterministic process,
$dN(t)/dt=F(N(t))$,
in the absence of stochastic effects.

If the population, instead, grows deterministically
only between jumps in population size which arrive at
the times of events of a Poisson process and with a magnitude
that is continuously distributed, then the growth of the
population satisfies the following stochastic differential equation,
\eq\label{SDEGen}
dN(t) = F(N(t)) dt + \int_{-\infty}^{\infty} J(N(t),u){\cal P}[du,dt],
~~ N(0) = x ~~ \in (a,b].
\nq
Here ${\cal P}$ is a Poisson random measure (Gihman and Skorohod \cite{GS72}).
The symbol ${\cal P}[du,dt]$
counts the number of jumps of a compound Poisson process with
amplitudes $J(x,u)$ generated in the mark interval $(u,u+du)$ during
the time interval $(t,t+dt)$, while
$\int_{-\infty}^{\infty}{\cal P}[du,dt]$
is the number of jumps in $(t,t+dt)$.
The mean rate of arrival of jumps with amplitudes in $(u,u+du)$ is
$\mbox{E}[{\cal P}[du,dt]]=\phi(u)du~dt$.
Thus, the total jump rate is
$\lambda=\int_{-\infty}^{\infty}\phi(u)du > 0$,
$\phi(u)/\lambda$ is the probability density function for the marks $u$.
The expected jump size conditioned on the population size at $t$ is
$\mbox{E}[J(N(t),u)|N(t)=x]=\int_{-\infty}^{\infty}\phi(u)J(x,u)du/\lambda$.

In our first paper \cite{HT78}, the special case of a 
downward jump with fixed magnitude $\eps$,
$\phi(u)=\lambda\delta(u-\eps), ~J(N,u)=-u, ~F(N)=r N (1-N/K)$
was treated, and
\eq\label{HT78N}
dN(t) = r N(t)\left(1-N(t)/K\right) dt-\eps d\Pi(t;\lambda),
\nq
where $\Pi(t;\lambda)$ is a simple Poisson process with rate
parameter $\lambda$,
$r$ is the intrinsic per capita rate of increase, 
$K$ is the carrying capacity of the environment,
$N(0)\in(0,K]$ and $\eps>0$.  
The Poisson random integral in (\ref{SDEGen}) is replaced by 
$-\eps d\Pi(t;\lambda).$

In \cite{HT81} we considered the particular choices
$\phi(u)=\lambda\delta(u-\eps), ~J(N,u)=-N u$ 
and $F(N)=r N (1-N/K)$ so that  
\eq\label{HT81N}
dN(t) = r N(t)\left(1-N(t)/K\right) dt - \eps N(t) d\Pi(t;\lambda),
\nq
where $N(0)\in(\Delta,K]$ for some $0<\Delta<K$ and $0<\eps<1$.
The case $\eps\geq 1$ is trivial, since the first jump leads to certain
extinction.  Note that in this case the single mark $u=\eps$ is the
proportion lost in this density independent disaster, whereas in
the previous case the mark represents the absolute amount lost in a population
size independent disaster.
The choice of a positive 
{\it effective extinction level} $\Delta$ permits the
population to go extinct in finite time in analogy with the discrete population
case where extinction occurs when the population gets below reproducing 
pair, rather then zero or below.

In the above models, the quantity of immediate interest is the
extinction or persistence time of the population size on some
interval $(a,b)$.
Let
\eq\label{T_ext}
T_{\rm ext}(x) = \inf_t\{N(t) \notin (a,b) ~|~ N(0) = x \in (a,b)\}
\nq
be the random variable denoting this persistence time and its expectation be
$T(x)=\mbox{E}[T_{\rm ext}(x)]$.

For the model of equation (\ref{SDEGen}) with $F(N)$ appropriate for 
logistic growth, this quantity satisfies the integro-differential equation 
\cite{HT76}
\eq\label{DE-T-Log}
F(x) T\PRIME(x) +
\int_{-\infty}^{\infty} T(x+J(x,u))\phi(u)du - \lambda T(x) = -1,
\nq
when $x$ is in $(a,b)$.  The choice of the minimum population $a$
depends on the choice of the drift $F(x)$ and the jump size $J(x,u)$.
In those cases where the value zero is attainable, $a=0$ is chosen 
to identify extinction with the passage of the population size to
zero or below.  If zero is not attainable in a finite time, then
for practical considerations an effective extinction should be defined as the
first passage to some small positive population size $a=\Delta$.
%
\subsection{Boundary conditions}
%
The general boundary condition for the
integro-differential equation (\ref{DE-T-Log})
is that
$T(x)=0 ~{\rm for}~ x\leq a ~{\rm or}~ x\geq b$.
This is a simple consequence of the fact that the first exit time from
$(a,b)$ is zero if the process starts outside that interval.

In addition to the above boundary condition, the
integro-differential equation is solved subject to the constraints that
$T(x)$ be bounded and continuous in its domain, deriving natural boundary
conditions in each special case where appropriate.
%
\section{Special examples}
%
The integro-differential equation (\ref{DE-T-Log})
can be reduced to quadratures in special cases.
%
\subsection{Model A:~ logistic growth with exponentially distributed,
N-independent disasters}
%
For a population that has a limited environment, a logistic model with fixed
carrying capacity serves as a simple deterministic, continuous model.  However,
it is unrealistic to expect a population 
to increase continuously without being interrupted
by occasional disasters of various magnitudes. We will
assume, in accordance with intuition, that large disasters
occur less frequently than small ones. Hence we choose
an exponential distribution for the magnitude of the
downward jumps. Thus we choose
\eq
F(x)&=&r x\left(1-x/K\right),~~
J(x,u) = -u,~~ {\rm and}~~
\phi(u) = \left\{\begin{array}{ll}
\lambda/m e^{-u/m}, &  u \geq 0\\*[0.5ex]
0, & \mbox{otherwise}
\end{array} \right. ,
\nonumber\\&~&
\nq
where $r>0$ and $K>0$, so that the drift is logistic and the amplitude of
the disasters is independent of population size.
The conditional mean jump amplitude is $\mbox{E}[J(N(t),u)|N(t)=x]=m$
and the conditional variance is $\mbox{Var}[J(N(t),u)|N(t)=x]=m^2$.
Choosing $a=0$ and $b=+\infty$, (\ref{DE-T-Log}) becomes
\eq\label{DE-T-A}
r x\left(1 - x/K\right) T\PRIME(x)
+ \frac{\lambda}{m}\int_0^{\infty}  e^{-u/m} T(x-u) du
- \lambda T(x) = -1.
\nq

Applying the general boundary condition and changing variables
in the integrand, 
\eq\label{DE-T-A2}
\left[r x\left(1 - x/K\right)T\PRIME(x) - \lambda T(x) + 1
\right] e^{x/m}
+ \frac{\lambda}{m}\int_0^x e^{v/m}T(v) dv = 0.
\nq
A {\it natural boundary condition} at extinction follows from taking the limit
$x\ARW 0^+$ in (\ref{DE-T-A2}) which yields
$T(0^+)=1/{\lambda}$.
This condition means that sufficiently small populations
would go extinct at the first disaster.
Note that $T(0^+)$ bounded implies that $T\PRIME(x)$ is less
than order $1/x$ as $x\ARW 0^+$.

A single differentiation in $x$ converts the integro-differential
equation (\ref{DE-T-A2}) to a first order ordinary differential
equation in $T\PRIME(x)$,
\eq\label{DE2ModA}
m r x \left(1-x/K\right) T''(x)
+ \left[r x\left(1-x/K\right)
+ m r\left(1-2 x/K\right)-m\lambda\right]T\PRIME(x) = -1.
\nq
The first integral of (\ref{DE2ModA}) yields
\eq\label{DE-T-A3}
\lambda T\PRIME(x) = x^{q-1}(K-x)^{-q-1} e^{-x/m}
\left[C_1 + q \frac{K}{m}\int_x^K e^{z/m}
z^{-q} (K-z)^q dz\right],
\nq
where
$q\equiv\lambda/r$
is the {\it bio-disaster growth ratio} defined in our previous papers
\cite{HT78,HT81}.  The parameter $q$ is the ratio of the
deterministic return (e-folding) time, $1/r$, to the 
expected waiting time, $1/\lambda$, between disasters.

    The constant, $C_1$, of integration in (\ref{DE-T-A3}) must be
zero when $q >0$ ($\lambda>0$ and $r>0$),
because $T\PRIME(x)$ must be less than order
$1/(K-x)$ for $T(x)$ to be bounded as $x\ARW K^-$.  Thus,
$\lambda T\PRIME(K^-)=q/((q+1)m)$
is a second natural condition.
Another integration with this
natural condition converts (\ref{DE-T-A3}) into closed form:
\eq\label{TModA}
\lambda T(x) = 1+\frac{q K}{m}\int_0^x e^{-y/m} y^{q-1}(K-y)^{-q-1}
\int_y^K e^{z/m}z^{-q}(K-z)^q dz dy .
\nq
It is important to note that this formal solution has been obtained
from the singular functional differential equation (\ref{DE-T-A2}) using
natural boundary conditions that are deduced from (\ref{DE-T-A2}) and 
boundedness.  No externally imposed Dirichlet conditions
have been or can be imposed, assuming $q>0$.

The second order differential equation (\ref{DE2ModA}) 
for mean first exit time $T(x)$
has an interesting interpretation as arising from a stochastic diffusion.
This equation would govern
the first passage time for a diffusion process governed by the
It\^{o} stochastic differential equation
\eq\label{SDE-ADiff}
d\widetilde{N}(t)&=&\left[r\widetilde{N}(t)
\left(1-\widetilde{N}(t)/K\right)
+ m r\left(1-2\widetilde{N}(t)/K\right)-m\lambda\right] dt
\nonumber\\
&~& +\sqrt{2 m r\widetilde{N}(t)\left(1-\widetilde{N}(t)/K\right)}dW(t),
\nq
where $W(t)$ is the Wiener or Gaussian process.
The diffusion coefficient, $2m r\widetilde{N}(1-\widetilde{N}/K)$, is
proportional to the logistic growth function.  The drift in the square bracket
of (\ref{SDE-ADiff}) is basically the logistic growth
perturbed by the mean jump rate, $-m \lambda$, and by the
derivative of the diffusion coefficient.  The boundary conditions
as a set are natural for the original expected extinction time of
the exponentially distributed disaster model.

There appears to be no simpler form for (\ref{TModA}) in terms of
elementary functions.  When $K$ is very large compared to the mean jump size
$m$, the double integral in (\ref{TModA}) can be
simplified asymptotically by a 2-dimensional extension of Laplace's method
\cite{HT82}.  The method relies on the fact that the dominant contribution
of the exponentials in the double integral occurs at the maximum
of the exponent of the integrand at $(y,z)=(0,K)$.
Expanding the integrals about this maximum and extending the limits of
integration of the inner integral to infinity and
retaining the outer integral limits to preserve the uniform validity
at $x=0$, 
\eq\label{TModAAs}
\lambda T(x) &\sim& 1 + q\Gamma(q+1)R^{-2q} e^R
\left[1+\frac{q(q+1)}{R}(1+\frac{\DS(q+1)(q+2)}{\DS 2 R})
\right.
\nonumber \\*[0ex]
&~&
+ \left.\oh\left(\frac{\DS 1}{\DS R^3}\right)\right]
\left[\left(1+\frac{q(q+1)}{R}
+\frac{\DS q(q+1)^2 (q+2)}{\DS 2 R^2}\right)
\gamma\left(q,\frac{x}{m}\right)\right.
\\*[0ex] \nonumber
&-& \left.\frac{(q+1)}{R}\left(1
+\frac{\DS (q+2)(q+1+x/m)}{\DS 2 R}\right)
\left(\frac{x}{m}\right)^q e^{-x/m}
+ \oh\left(\frac{1}{R^3}\right)\right]
\nq
$\mbox{as}~ R\equiv K/m\ARW\infty$,
where $\Gamma(q)$ is the gamma function
and $\gamma(q,x)$ is the incomplete gamma function.
The approximation (\ref{TModAAs}) is
uniformly valid when $x=0^+$
because $\gamma(q,0)=0$.  In (\ref{TModAAs}), only
the first three exponentially large terms with exponent $R$ have been
retained.  When
$x=K^-$ and $\gamma(q,R)\sim\Gamma(q)$, the approximation
(\ref{TModAAs}) simplifies to
\eq\label{TKModAAs}
\lambda T(K^-) \sim 1 + \Gamma^2(q+1) R^{-2q} e^R
\left[1 + \frac{2q(q+1)}{R}
\left(1+\frac{(q+1)^2}{R}\right)
+\oh\left(\frac{1}{R^3}\right)\right]
\nq
as $R\equiv K/m\ARW\infty$.  The simplified asymptotic approximation
(\ref{TKModAAs}) indicates that $\lambda T(K^-)$ will be exponentially
large when $R$ is large, i.e., its behavior is exponentially
dominated by the factor, $\exp(R)$, when the carrying capacity $K$ is a large
multiple of the mean disaster amplitude $m$.  There is also a weaker, algebraic
dependence on the bio-disaster growth parameter $q$ due to the factor
$R^{-2q}$.  The approximation,
(\ref{TModAAs}) or (\ref{TKModAAs}), should be useful for very large $R$.

Extinction according to this model is certain, but the population will
persist for extremely long times, provided $x/m$ and $R$ are
sufficiently large.  It is well known that almost all species that
have ever existed have become extinct.  Raup \cite{Raup81b} gives a summary
of mean species duration times for many fossil groups, the range being
between one and eleven million years.

In contrast, small values of $x$ are of interest in studying the
relative dangers in going extinct from a small number of colonizers.
Very close to the extinction level, as $x/m\ARW 0^+$ and using
(\ref{TModA}), the composite dependence is found
\eq\label{TModAExtAs}
\lambda T(x) \sim 1 + \left\{ \begin{array}{ll}
\frac{\DS q}{\DS (q-1)} \frac{\DS x}{\DS m}, & q > 1 \\*[0.4ex]
-\frac{\DS x}{\DS m} \ln\left(\frac{\DS x}{\DS m e}\right), & q = 1\\*[0.4ex]
R\left(\frac{\DS x}{\DS K}\right)^q
\int_0^1 z^{-q}(1-z)^q e^{R z} dz, & 0 < q < 1
\end{array} \right. ,
\nq
due to the singularities in the $z$-integral.
Approximation (\ref{TModAExtAs}) implies that when the bio-disaster growth
ratio $q$ is less than one, $\lambda T(x)$ will increase most
rapidly; in fact, with a vertical asymptote as $x/m\ARW 0^+$ and also
with exponential dependence on $R$ when $R$ is large.  When $q$ is
greater than one the increase in $x$ is only linear in $x$.  Hence
$q<1$ or $r>\lambda$ is advantageous to the colonizing population,
because very long persistence time will be most probable for a
relatively small initial population size.

The same singular features of the solution (\ref{TModA}) that limit the
usefulness of the asymptotic approximation (\ref{TModAAs}) also hinder the
numerical approximation of (\ref{TModA}) for general values of $q$ and
moderate values of $R$. 
We here first isolate 
the singular powers containing $q$ and approximate only the
nonsingular terms.  Let $u(x)=x/(K-x)$ when $x<K/2$ and
$v(x)=(K-x)/x=1/u$ when $x>K/2$, so that (\ref{TModA}) is
transformed to a more manageable form:
\eq\label{ModANum}
\nonumber
\lefteqn{\lambda T(x)}\\
&=&\left\{\begin{array}{ll}
\begin{array}{l}
1 + q R e^R \int_0^{u(x)} U^{q-1} A(U)\\*[0.2ex]
~~~~ \times \left[\int_0^1 V^q B(V) dV
+ e^{-R/2}\int_U^1 W^{-q} C(W) dW \right] dU,
\end{array}
&x \leq \frac{K}{2}
\\*[0.2ex]
 \lambda T\left(\frac{K}{2}\right)
+ q R e^{R/2} \int_{v(x)}^1 V^{-q-1} D(V)
~ \int_0^V W^q B(W) dW dV ,
&x \geq \frac{K}{2}
\end{array}\right. ,
\nq
where
$A(u)=\exp(-R u/(1+u))$, $B(u)=A(u)/(1+u)^2$,
$D(u)=\exp(-R/2)/A(u)$, and $C(u)=D(u)/(1+u)^2$.
The relatively nice functions $A$, $B$, $C$, and $D$ are approximated
by linear interpolation with sufficient accuracy using parabolically
spaced nodes and the remaining integrals can be done exactly, because
their integrals contain only powers, some singular.

Some numerical results using (\ref{ModANum}) for
{\it Model~A} with $\lambda T(x)$ are given
in the semilog plot of Fig.~\ref{FIG1}.  
Here, $K$ has the value 10, while $q = 0.5$,
1.0 and 2.0.  These values of $q$ are just representative of the three
distinct asymptotic cases, which are $0<q<1$, $q=1$ and $q>1$.
Taking into account the logarithmic scale, it can be
seen that the differences 
between the values of $\lambda T(x)$
are large 
at small and moderate values of $x$, for various values of $q$.
The curves for $q=0.5$ (the smallest bio-disaster growth ratio
considered) rise very sharply away from $x=0^+$ in comparison with
the curve for $q = 2.0$.  Note the tremendous advantage that 
populations with large per capita growth rates $r$ have because their
mean persistence times are the longest and are achieved with a small
colonization effort.  Another qualitative biological implication of the results
is that the population can reach nearly maximal persistence outside of
a sufficiently small neighborhood of extinction, but strong colonization
properties are critical within this neighborhood.
%
\begin{figure}[htb]
\begin{center}
~\\*[-20pt]
\hspace*{10pt}
\epsfysize350pt
\epsffile{expdistfig1w.eps}
\caption{\label{FIG1}
Logarithm of the scaled, expected extinction time,
$\log\{\lambda T(x)\}$, for {\it Model~A} with exponentially distributed,
size independent disasters as a function of the initial population size, $x$,
for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$}
\end{center}
\end{figure}
%

In Fig.~\ref{FIG2},
the asymptotic (\ref{TModAExtAs}) and numeric (\ref{ModANum}) approximations
to the exact solution for the expected extinction time are compared in terms
of the percent relative difference.  The difference is quite small for values
of the bio-disaster ratio of $q=0.5$ and $1.0$, but quite large for $q=2.0$,
reaching about $30\%$ for small, colonization level, population sizes.  Thus,
for small population sizes and higher bio-disaster growth ratios, the
asymptotic approximation significantly underestimates the extinction time, and
so overestimates the danger of extinction.
%
\begin{figure}[htb]
\begin{center}
~\\*[-20pt]
\hspace*{10pt}
\epsfysize350pt
\epsffile{expdistfig2w.eps}
\caption{\label{FIG2}
The percent relative difference between the asymptotic
and numerical approximations to expected extinction time,
$(T_{\rm as}(x)/T(x) -1)*100$, for {\it Model~A} with exponentially distributed,
$N$-independent disasters as a function of the initial population size, $x$,
for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$}
\end{center}
\end{figure}
%
%
\subsection{Model B:~ density independent disasters with uniformly
distributed jump size}
%
When the effect of a disaster is proportional to population size
the disaster is said to be {\it density independent},
because the relative disaster size is independent of $N$.
The density independent disaster is usually due to abiotic causes which
affect the population as a whole, whereas biological causes
due to species interactions, such as epizootics, may lead to nonlinear or other
density independent disasters \cite{CJFAS87}.
For convenience, the size of the disaster per unit population size is
assumed to be uniformly distributed on $(0,1)$.
Note that any disaster greater than the population size has the same effect
as a disaster wiping out the whole population.
An example of possible density independent disasters with uniform distribution
is the occasional massive die-offs of alewives in Lake Michigan, but in
\cite{CJFAS87} only a simple Poisson disaster in a multi-species
bioeconomic model is used due to limited field data.

Hence,
\eq\label{JModB}
J(x,u) = - x u,
~~F(x) =  r x \left(1-x/K\right)
~~\mbox{and}
~~\phi(u) = 
\left\{\begin{array}{ll}
\lambda, & 0 < u < 1 \\
0, & \mbox{otherwise}
\end{array} \right. .
\nq
The conditional mean jump size
is $\mbox{E}[J(N(t),u)|N(t)=x]=-x/2$
with conditional variance
$\mbox{Var}[J(N(t),u)|N(t)=x]=x^2/12.$
Since the disaster gets small with the population size, ideal extinction
at $x=0^+$ will take an infinite amount of time.
Therefore, let $a=\Delta>0$ be the {\it effective extinction
size}, the lowest viable population size, in analogy with a discrete population
being extinct when there is less than one mating pair.
The expected extinction time
equation (\ref{DE-T-Log}) with general boundary condition reduces to
\eq\label{DEModB}
r x \left(1-x/K\right)T\PRIME(x)
+ \lambda \int_0^{1-\Delta/x} T(x - u x) du
- \lambda T(x) = -1
\nq
on $(\Delta,K)$.
The effective extinction boundary condition for (\ref{DEModB}) does
not have the simple form as the extinction natural boundary condition of 
{\it Model~A}, but instead has a mixed condition,
$r\Delta(1-\Delta/K)T\PRIME(\Delta^+)=\lambda T(\Delta^+)-1$,
which arises because (\ref{DEModB}) is not singular at $x=\Delta^+>0$.

As in the prior case, the 
integro-differential equation can be converted to a first order
differential equation in $T\PRIME(x)$ similar to the stochastic diffusion-like
expected extinction time equation (\ref{DE2ModA}) for {\it Model~A}.
Transforming $x$ out of the integrand, multiplying by $x$ and differentiating,
\eq\label{DE2ModB}
r x^2\left(1-x/K\right)T''(x)
+ \left[\left(r x^2\left(1-x/K\right)\right)'
- \lambda x\right] T\PRIME(x) = -1.
\nq

Otherwise, the two integrations with $T(K^-)$ bounded work the same as
before and
\eq\label{TModB}
\lambda T(x) = 1
   &+& \Delta^{q-1}(K-\Delta)^{-q}\int_{\Delta}^K z^{-q}(K-z)^q dz
\nonumber\\*[0.2ex]
   &+& q K \int_{\Delta}^x y^{q-2}(K-y)^{-q-1}
\int_y^K z^{-q}(K-z)^q dz dy
\nq
is the formal result.  For numerical purposes, (\ref{TModB}) can be transformed
to isolate potential troublesome arbitrary powers related to the
disaster growth ratio $q$ as in {\it Model~A} to the more manageable form:
\eq\label{TModBNum}
\lambda T(x)
=\left\{\begin{array}{l}
1+u^{q-1}(\Delta)~A(\Delta)\left[\int_{u(\Delta)}^1
~U^{-q} B(U) dU+\int_0^1 V^q  B(V)dV
\right],\\*[0.2ex]
~~~~~x=\Delta^+\\*[0.2ex]
\lambda T(\Delta^+)+q
\int_{u(\Delta)}^{u(x)} U^{q-2} A(U)
\left[\int_U^1 W^{-q}
+\int_0^1 W^q \right] B(W) dW dU,
\\*[0.2ex]
~~~~~\Delta<x\leq K/2\\*[0.2ex]
\lambda T(K/2)+q
\int_{v(x)}^{1} V^{-q-1} A(V) dV
\int_0^V W^q B(W) dW dV,
\\*[0.2ex]
~~~~~K/2<x\leq K
\end{array}\right.
\nq
but with new definitions,
$A(u)=1+u$ and $B(u)=1/(1+u)^2$.

For large carrying capacity $K$ relative to the effective extinction
level $\Delta$, the formal solution (\ref{TModB}) has the asymptotic
simplification as follows,
\eq\label{TModBAs}
\lambda T(x) \sim \left\{ \begin{array}{l}
\frac{\DS 1}{\DS 1-q}\left[\left(\frac{K}{\Delta}\right)^{1-q}
- q\left(\frac{K}{\DS x}\right)^{1-q}-q\ln\left(\frac{\DS ex}{\Delta}\right)
\right] + \oh(1), 
\\*[0.2ex]
~~~~~0 < q < 1,~ x\geq \Delta \\*[0.2ex]
\left[\half\ln\left(\frac{\DS K^2}{{\DS \Delta x}}\right) -1\right]
\ln\left(\frac{\DS x}{\Delta}\right) + \ln\left(\frac{K}{\Delta}\right),
\\*[0.2ex]
~~~~~q = 1,~  x\geq \Delta \\*[0.2ex]
\frac{\DS q}{\DS q-1}\ln\left(\frac{\DS ex}{\Delta}\right) +\oh(1),
\\*[0.2ex]
~~~~~q > 1,~ x\geq \Delta \end{array} \right. ,
\nq
as $K/\Delta\ARW\infty$.  This result has been obtained by a
somewhat different method 
through simply factoring out the most dominant
parts of the integrals in (\ref{TModB}).  
The non\-exponential integrands lack the sharp maxima needed to apply 
Laplace's method.  
The approximation (\ref{TModBAs}) indicates that expected persistence time 
grows algebraically or logarithmically as $K/\Delta\ARW\infty$, while
in {\it Model~A} the asymptotic growth is exponential in $R\equiv K/m$.

The numerical results for {\it Model~B}, using similar techniques to those
employed in {\it Model~A}, are exhibited in Fig.~\ref{FIG3} for $q=0.5,
1.0$ and $2.0$.
Again, the advantage for a population described by {\it Model~B} occurs 
when $q$ is less than unity, but
it is less pronounced than in {\it Model~A}.  
%
\begin{figure}[htb]
\begin{center}
~\\*[-20pt]
\hspace*{10pt}
\epsfysize350pt
\epsffile{expdistfig3w.eps}
\caption{\label{FIG3}
Logarithm of the scaled, expected extinction time,
$\log\{\lambda T(x)\}$, for {\it Model~B} of uniformly distributed, density
independent disasters as a function of the initial population size, $x$,
for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$}
\end{center}
\end{figure}
%
Figure~\ref{FIG4} shows the relative
difference between the asymptotic approximation and the numerical computation
for {\it Model~B} for the same set of bio-disaster values.  The asymptotic
approximation is reasonable for the case $q=1.0$,
as long as the initial population size is not near extinction,
somewhat less reasonable for $q=2.0$, and very bad for $q=0.5$.
%
\begin{figure}[htb]
\begin{center}
~\\*[-20pt]
\hspace*{10pt}
\epsfysize350pt
\epsffile{expdistfig4w.eps}
\caption{\label{FIG4}
The percent relative difference between the asymptotic
and numerical approximations to expected extinction time,
$(T_{\rm as}(x)/T(x) -1)*100$, for {\it Model~B} with exponentially distributed,
$N$-independent disasters as a function of the initial population size, $x$,
for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$}
\end{center}
\end{figure}
%
The qualitative biological implications of the results are similar to those
of {\it Model~A}, because sufficiently far away from extinction, nearly maximal
persistence is achieved.  Due to a combination of the uniformity
of the disaster amplitude distribution (the larger disasters are just as
likely as the smaller ones) and the density independent property of the
disaster amplitudes (the larger disasters can be comparable in size to the
population size $N$), the nearly maximal persistence times are not
exponentially large, but are only algebraic-logarithmically large.
%
\subsection{Model C:~ exponentially distributed, N-independent bonanza}
%
Some populations such as scrub oaks, sea palms, according to Paine
\cite{Paine79}, benefit from disasters to other species because such disasters
provide open patches which are bonanzas for them as colonizing species.
In the absence of the availability of these patches these populations decline
in numbers.
Hence, let the decline be a linear decaying function, $F(x)=-r x$.

{\it Model~C} is also motivated by the simple model of certain fish
stocks that remain at low levels most of the time, but have occasional
bonanzas in recruitment.  Hennemuth et al. \cite{Henne80} offer examples
from the North Atlantic fisheries.
The term, $-r x$,
corresponds to the usual assumption that adult fish populations
decay exponentially between recruitments (Beverton and Holt \cite{BH57}).
A model of optimal harvesting in the presence of discrete bonanzas is
given by Ryan and Hanson \cite{RHMB85}.

In this model, a {\it bonanza}
causes a rapid upsurge in the population size.  However, it is assumed
that the benefits of the upsurge are limited by an amount $K-x$ so that the
upsurged population size will not exceed some limit $K$.
Hence, the jump size is
\eq\label{JModC}
J(x,u) = \left\{ \begin{array}{ll}
u, & u \in (0, K-x) \\*[0.2ex]
0, & u \notin (0, K-x)
\end{array} \right. ,
\nq
and the distribution is exponential as in {\it Model~A}.
The conditional mean jump size is
$\mbox{E}[J(N(t),u)|N(t)=x]=m [1-(1+\xi)\exp(-\xi)]\in[0,m)$,
where $\xi\equiv (K-x)/m$, and the conditional variance is
$\mbox{Var}[J(N(t),u)|N(t)=x] =
m^2[1-\xi^2\exp(-\xi)-(1+\xi)^2\exp(-2\xi)]\in[0,m^2)$.

In this case, zero cannot be reached in a finite
time, since zero can only approached like $O(\exp(-r t))$,
so a minimum population size $a=\Delta$ is selected.

The equation governing the mean persistence time reduces to
\eq\label{DEModC}
&~&-r x T\PRIME(x) + \frac{\lambda}{m} \int_0^{K-x}
e^{-u/m} T(x+u) du
\nonumber\\&~~~~~& + 
\frac{\lambda}{m} T(x) \int_{K-x}^{\infty}
e^{-u/m} du -\lambda T(x) = -1,
\nq
provided $x$ is in $(\Delta,K]$.

As with {\it Models A} and {\it B},
Equation (\ref{DEModC}) may be converted by differentiation to a first
order ODE for $T\PRIME(x)$ in the form like that for a stochastic diffusion,
\eq\label{DE2ModC}
m r x T''(x) + \left[-r (x-m) + \lambda m\left(1
- \exp\left(-\frac{(K-x)}{m}\right)\right)\right] T\PRIME(x) = -1.
\nq
The general boundary condition and the fact that at
$x=\Delta^+$ the population will drift below the effective extinction
level $\Delta$ with rate $-r$
before the first upsurge with probability one imply the initial
condition, $T(\Delta^+)=0$.
In addition, allowing $x\ARW K^-$ in (\ref{DEModC}) yields the
non\-singular, final condition,
$T\PRIME(K^-)=1/(r K)$.
The solution of (\ref{DE2ModC}) from quadratures is
\eq\label{TModC}
\lambda T(x) = q R^q e^{-R}\int_{\Delta/m}^{x/m}
\frac{e^y y^{-1-q}}{h(y)} \left[1 + R^{-q}e^R \int_y^R e^{-z}z^q h(z) dz
dy \right] ,
\nq
where the variables of integration have been scaled by the mean jump amplitude
$m$, $R\equiv K/m$, $q\equiv\lambda/r$ is now the {\it bio-bonanza ratio},
and
$h(x)=\exp[q e^{-R}\int_x^R(e^t/t)dt]$
is the integrating factor for (\ref{DE2ModC}).

For the numerical evaluation of multiple integral form of the solution
(\ref{TModC}), a multiple integral modification of the standard Simpson's
rule is used to achieve sufficient accuracy.
The principal multiple integration modification arises in the double integral,
\eq\label{Simp2}
\int_{x_{j-1}}^{x_j} G(y)\int_y^{x_n}H(z) dz dy &\simeq& \frac{x_j-x_{j-1}}{6}
~\left(\left(G_{j-1}+4G_{j-0.5}+G_j\right)\int_{x_j}^{x_n}H(z) dz \right.
\nonumber\\
&~~~~~& + \left.G_{j-1}\int_{x_{j-1}}^{x_j} H(z) dz
+4G_{j-0.5}\int_{x_{j-1}}^{x_{j-0.5}} H(z) dz\right),
\nonumber\\&~&
\nq
where $G_j\equiv G(x_j)$ for $j=0$ to $n$.
For the last inner integral,
modification of the standard Simpson's rule by using the same parabolic
approximation over the half-subinterval
$(x_{j-1},x_{j-0.5})$ with compatible
accuracy of standard Simpson's, i.e.,
\eq\label{SimpHalf}
\int_{x_{j-1}}^{x_{j-0.5}} H(z) dz \simeq \frac{x_j-x_{j-1}}{24}~
\left(5H_{j-1}+8H_{j-.05}-H_j\right),
\nq
although the extra accuracy of Simpson's rule for cubics must be sacrificed.
Hence, using the combined full and half-interval Simpson's rule
yields a numerical method for the double integrals with quadratic accuracy.
Similarly, the full solution (\ref{TModC})
can be approximated as triple and double integrals similar to the double 
integral case.

As $R\ARW\infty$, a simple asymptotic approximation does not
appear to work, in spite of the appearance of the exponentially
small factor $e^{-R}$ and exponentially large factor $e^R$ in (\ref{TModC}),
which are cancelled out nonuniformly by exponentials in the integrands.

Some examples of the numerical results obtained for {\it Model~C} are
shown in Fig.~\ref{FIG5}, just scaled by the disaster rate $\lambda$ and not 
in logarithm form, in contrast to figures for the previous models.
%
\begin{figure}[htb]
\begin{center}
~\\*[-20pt]
\hspace*{10pt}
\epsfysize350pt
\epsffile{expdistfig5w.eps}
\caption{\label{FIG5}
Scaled, expected extinction time,
$\lambda T(x)$, for {\it Model~C} with exponentially distributed,
$N$-independent bonanzas as a function of the initial population size, $x$,
for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$}
\end{center}
\end{figure}
%
Here, $K=10$, $m=1$ and $\Delta=1$, so that $R=10$ and
$\Delta/m=1$.  The {\it bio-bonanza growth ratio}, $q \equiv\lambda/r$, has
the values $q=0.5, 1.0$ and $2.0$ for the curves starting from the bottom to
top, respectively.  Apparently, $q$ should be greater than one for a population
to attain reasonable extinction time levels, i.e., so that the bonanzas
more than compensate for the deterministic exponential decay.  The qualitative
biological implications of the results are somewhat similar to those of
{\it Models~A} and~{\it B}, considering that the results of {\it Model~C} are 
not presented in logarithmic form, in that, away from effective extinction, 
extinction times
take on larger values, though less dramatically than the former models.
{\it Model~C} is quite different from the other models in terms of the roles 
played by deterministic and stochastic processes, because extinction is
caused by exponential deterministic decay to an effective extinction
(true extinction to zero population is not attainable), while persistence is
strengthened by random bonanzas.

Note that the expected extinction time is an increasing
function of the bio-bonanza ratio in Fig.~\ref{FIG5} for {\it Model~C}, 
but it is a decreasing
function of the bio-disaster ratio in the disaster examples of 
Fig.~\ref{FIG1} for
{\it Model~A} and Fig.~\ref{FIG3} for {\it Model~B}, since the bonanza is 
beneficial while the disaster is adverse to the population.

Even though the relative carrying capacity $R=10$, the curves do not exhibit
near horizontal asymptotes as the relative population size $x/m$ approaches
$R$, confirming the difficulty in finding an asymptotic behavior for large $R$.
%
\section{Conclusions}
%
Two models of growth with randomly occurring disasters
and one model of decay with randomly occurring bonanzas have 
been presented for which the equations reduce
to first order  ordinary differential equations 
as in Eqs. (\ref{DE2ModA}), (\ref{DE2ModB}),
and (\ref{DE2ModC}).  These equations would correspond to the expected
extinction time for a stochastic diffusion with drift, except that the
boundary conditions, natural or otherwise,  as a set do not correspond to the
usual diffusion.

Due to the nature of the reduced equation, the expected
extinction time can be further reduced to quadratures, i. e., to a formal
solution in terms of explicit integrals.
The reduction to quadratures still leaves nontrivial
explicit integrals, so special numerical techniques have been used in each
case, adapted to some critical form of the integrals for better conditioning.

For {\it Models~A} and {\it B}, which included disasters, 
the expected extinction time
$T(x)$ quickly levels off outside of a neighborhood of extinction or
effective extinction,  but the fast growth near extinction is exponential
for the exponential distributed disaster of {\it Model~A}, but 
algebraic-logarithmic for the uniformly distributed disaster of {\it Model~B}.
For the exponentially distributed bonanza case of
{\it Model~C}, the growth with respect to scaled, initial population size $x$ is
much more varied in dependence.  Also, the expected extinction time is an
decreasing function of the bio-disaster ratio $q=\lambda/r$ in the disaster
models, but an increasing function of the bio-bonanza ratio $q=\lambda/r$
in the bonanza model.

In all the models, the common qualitative biological implication of
the results is that sufficiently far away from the colonizing important
neighborhood of extinction or effective extinction,
the larger persistence times are achieved.
However, the quantitative levels of the persistence times varies considerably
for the three models.   The size and distribution of the random disasters
undoubtedly provides the difference between {\it Models~B} and~{\it A}, 
with more dangerous and larger disasters being more likely for {\it Model~B}.  
{\it Model~C} is quite different from the other two models, in that the 
deterministic process makes the population tend toward extinction, while the 
stochastic process
make it tend away from extinction, the opposite tendencies in contrast
to the other two models.

In the case of the models with disasters, {\it A} and {\it B}, 
asymptotic simplifications
have been found useful for particular ranges of the bio-disaster ratio, away
from extinction.  The asymptotic results are useful in determining the
qualitative behavior of the expected extinction time.
Asymptotic simplification was not fruitful for the third, bonanza model,
{\it Model~C}.
%
\\*[0.1in]
%
{\small\it Acknowledgements.}~~ 
{\small This research supported in part by National Science
Foundation  Grants  DMS 91-02343, 93-01107 and 96-26692
at the University of Illinois at Chicago.  
We performed part of this research in
the Department of Biomathematics at the University of California at Los
Angeles and in the Department of Mathematics at Monash University.
We thank two reviewers for their constructive criticisms which led to a
more concise paper.}
%
\begin{thebibliography}{99}
%
\bibitem{BBCP89} Bartoszynski, R., B\"{u}hler, Chan, W., Pearl, D. K.:
Population processes under the influence of disasters occurring independently
of population size.  J. Math. Biol. {\bf 27}, 167-178 (1989)
%
\bibitem{BH57} Beverton, R. J. H., Holt, S. J.:
On the dynamics of exploited fish populations.
Fisheries investigations, {series II, vol. 19}.  London:  Ministry of
Agriculture, 1957
%
\bibitem{Brock85} Brockwell, P. J.:
The extinction time of a birth, death and catastrophe process and of a
related diffusion model.  Adv. Appl. Probab. {\bf 17}, 42-52 (1985)
%
\bibitem{Brock82} Brockwell, P. J., Gani, J., Resnick, S. I.:
Birth, immigration and catastrophe processes.  Adv. Appl. Probab.
{\bf 14}, 709-731 (1982)
%
\bibitem{Brock83} Brockwell, P. J., Gani, J. M., Resnick, S. I.:
Catastrophe processes with continuous state-space.  Austral. J. Statist.
{\bf 25}, 208-226 (1983)
%
\bibitem{FLES90}  Flessa, K. W.:
The ``facts'' of mass extinctions.  In: V. L. Sharpton, P. D. Ward:
Global catastrophes in earth history (Special Paper 247, pp. 1-7)
Boulder: Geological Society of America 1990
%
\bibitem{GS72}  Gihman, I. I., Skorohod, A. V.:
Stochastic differential equations.  New York: Springer 1972
%
\bibitem{Hanson83} Hanson, F. B.:
Analysis of a singular functional differential equation
that arises from stochastic modeling of population growth.
J. Math. Anal. Applic. {\bf 96}, 405-419 (1983)
%
\bibitem{CJFAS87}  Hanson, F. B.:
Bioeconomic model of the Lake Michigan alewife fishery.
Can. J. Fish. Aquat. Sci. {\bf 44}, suppl. 2, pp. 298-305 (1987)
%
\bibitem{NRM88}  Hanson, F. B., Ryan, D.:
Optimal harvesting with density dependent random effects,
Natural Resource Modeling {\bf  2}, 439-455 (1988)
%
\bibitem{HT82}  Hanson, F. B., Tier, C.:
A stochastic model of tumor growth.  Math. Biosci. {\bf 61}, 73-100 (1982)
%
\bibitem{HT78}  Hanson, F. B., Tuckwell, H. C.:
Persistence times of populations with large random fluctuations.
Theor. Popul. Biol. {\bf 14}, 46-61 (1978)
%
\bibitem{HT81}  Hanson, F. B., Tuckwell, H. C.:
Logistic growth with random density independent disasters.
Theor. Popul. Biol. {\bf 19}, 1-18 (1981)
%
\bibitem{Henne80}  Hennemuth, R. C., Palmer, J. E., Brown, B. E.:
A statistical description of recruitment in eighteen selected fish
stocks.  J. Northw. Atl. Fish. Sci. {\bf 1}, 101-111 (1980)
%
\bibitem{JABL86}  Jablonski, D.:
Causes and consequences of mass extinctions:  A comparative approach.
In: D. K. Elliot: Dynamics of extinction (pp. 183-229)  New York: Wiley 1986
%
\bibitem{LHW93}  Ludwig, D., Hilborn, R., Walters, C.:
Uncertainty, resource exploitation, and conservation: Lessons from
history.  Science {\bf 260}, 17 \& 36 (1993)
%
\bibitem{MT93a}  Mangel, M., Tier, C.:
A simple, direct method for finding persistence times of populations
and application to conservation problems.  Proc. Nat. Acad. Sci. USA {\bf 90},
1-28 (1993)
%
\bibitem{MT93b}  Mangel, M., Tier, C.:
Dynamics of metapopulations with demographic stochasticity and
environmental catastrophes.  Theor. Popul. Biol. {\bf 44}, 1-31 (1993)
%
\bibitem{Paine79}  Paine, R. T.:
Disaster, catastrophe, and local persistence of the sea palm
{\it Postelsia palmaeformis}.  Science {\bf 205}, 685-686 (1979)
%
\bibitem{PTB79} Pakes, A. G., Trajstman, A. C., Brockwell, P. J.:
A stochastic model for a replicating population subjected to mass
emigration due to population pressure.  Math. Biosci. {\bf 45},
137-157 (1979)
%
\bibitem{Raup81b}  Raup, D. M.:
Extinction:  Bad genes or bad luck?.
Acta Geologica Hispanica, Concept and Method in Paleontology {\bf 16},
25-33 (1981)
%
\bibitem{RHMB85}  Ryan, D., Hanson, F. B.:
Optimal harvesting with exponential growth in an environment with
random disasters and bonanzas.
Math. Biosci. {\bf 74}, 37-57 (1985)
%
\bibitem{RHJMB86} Ryan, D., Hanson, F. B.:
Optimal harvesting of a logistic population in an environment
with stochastic jumps. J. Math. Biol. {\bf 24}, 259-277 (1986)
%
\bibitem{Traj81}  Trajstman, A. C.:
A bounded growth population subjected to emigrations due to population
pressure.  J. Appl. Probab. {\bf 18}, 571-581 (1981)
%
\bibitem{TO76}  Tsurui, A., Osaki, S.:
On a first-passage problem for a cumulative process with exponential
decay.  Stoch. Processes Applic. {\bf 4}, 79-88 (1976)
%
\bibitem{HT76} Tuckwell, H. C.:
On the first exit time problem for temporally homogeneous Markov
processes.  J. Appl. Prob. {\bf 13}, 39-48 (1976)
%
\bibitem{TuckKoz87} Tuckwell, H. C., Koziol, J. A.:
Logistic population growth under random dispersal.
Bulletin Math. Biol. {\bf 49}, 495-506, 1987.
%
\end{thebibliography}
%
%\newpage
%\begin{center}
%{\bf Figure Captions}
%\end{center}
%\begin{description}
%\item [Figure 1.]  Logarithm of the scaled, expected extinction time,
%$\log\{\lambda T(x)\}$, for {\it Model~A} with exponentially distributed,
%size independent disasters as a function of the initial population size, $x$,
%for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$.
%\item [Figure 2.]  The percent relative difference between the asymptotic
%and numerical approximations to expected extinction time,
%$(T_{\rm as}(x)/T(x) -1)*100$, for {\it Model~A} with exponentially distributed,
%$N$-independent disasters as a function of the initial population size, $x$,
%for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$.
%\item [Figure 3.]  Logarithm of the scaled, expected extinction time,
%$\log\{\lambda T(x)\}$, for {\it Model~B} of uniformly distributed, density
%independent disasters as a function of the initial population size, $x$,
%for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$.
%\item [Figure 4.]  The percent relative difference between the asymptotic
%and numerical approximations to expected extinction time,
%$(T_{\rm as}(x)/T(x) -1)*100$, for {\it Model~B} with exponentially distributed,
%$N$-independent disasters as a function of the initial population size, $x$,
%for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$.
%\item [Figure 5.]  Scaled, expected extinction time,
%$\lambda T(x)$, for {\it Model~C} with exponentially distributed,
%$N$-independent bonanzas as a function of the initial population size, $x$,
%for several values of the bio-disaster growth ratio, $q = 0.5, 1.0, 2.0$.
%\end{description}
\end{document}
