\documentclass[twoside]{article}
\usepackage{epsf} % To input PostScript figures
\pagestyle{myheadings}
\markboth{ Transient effects of stochastic multi-population models }
{ Thomas C. Gard }
\begin{document}
\setcounter{page}{81}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
Nonlinear Differential Equations, \newline
Electron. J. Diff. Eqns., Conf. 05, 2000, pp. 81--90\newline
http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)}
\vspace{\bigskipamount} \\
%
Transient effects of stochastic \\ multi-population models
%
\thanks{ {\em Mathematics Subject Classifications:} 92D25, 60H10.
\hfil\break\indent
{\em Key words:} multi-population models,
stochastic differential equations, transient behavior.
\hfil\break\indent
\copyright 2000 Southwest Texas State University.
\hfil\break\indent Published October 25, 2000. } }
\date{}
\author{ Thomas C. Gard }
\maketitle
\begin{abstract}
We give some estimates for exit probabilities through specific
portions of the boundary of bounded subsets of the feasible region for
solutions of stochastic population interaction models. These exit
probability estimates can indicate initial tendencies for survival or
extinction for the modeled populations. When the subset boundaries are
given by level curves of multiple Liapunov type functions, the estimates
are more tractable.
\end{abstract}
\newtheorem{theorem}{Theorem}
\catcode`@=11
\@addtoreset{equation}{section}
\def\theequation{\thesection.\arabic{equation}}
\catcode`@=12
\section{Introduction and main result}
Transient effects are often
overlooked in the mathematical analysis of population dynamics models.
Generally the focus of qualitative investigations is on long term
properties
such as asymptotic stability, while for the short term picture,
simulations
are carried out. Even after fairly exhaustive efforts of the latter type,
it may be still difficult to conclude anything of very general nature
concerning short-to-intermediate time horizon events. Such information
could have crucial practical significance. If the location of a stable
equilibrium, periodic solution or strange attractor is close to the
boundary
of the feasible region, or if trajectories typically sweep close to the
boundary before settling down near some special trajectory, the transient
behavior of the model may be more important than some of these other
qualitative properties. The situation may be even more critical if random
features ([1]) are incorporated in the model.
We address this situation with a result for a generic such model - an Ito
stochastic $n$-dimensional multi-population interaction model
(in Kolmogorov form)
\begin{equation}
dX=D(X) [f(X)\,dt+G(X) \,dW].
\end{equation}
Here
$$
X=(X_1,X_2,...,X_n)
$$
with $X_i$ representing the ith-species population density, and
$$
D=D(X)=\mathop{\rm diag}\{X_1,X_2,...,X_n\};
$$
more precisely,
$$
X=X(t,\omega)=\{X_1(t,\omega),X_2(t,\omega),...,X_n(t,\omega)\}
$$
and
$$
W=W(t,\omega)=\{W_1(t,\omega),W_2(t,\omega),...,W_m(t,\omega)\}
$$
with the $W_j$ denoting independent Brownian motion processes, $t>0$, and
$\omega\in\Omega$, a probability sample space. Formally, (1.1) can arise
if
one models the impact of environmental stochasticity by characterizing the
net per capita growth rates of the populations $F(x)$,
\begin{equation}
F_i(x) = \frac{1}{x_i}\frac{dx_i}{dt},
\end{equation}
as random noise fluctuations about some average values
$f(x)$ with intensities $G(x)$ dependent on the population size $x$
$$
F(x)=f(x)+G(x)N
$$
where
$$
N=N(t,\omega)=\{N_1(t,\omega),N_2(t,\omega),...,N_m(t,\omega)\}
$$
with the $N_j$, $t>0$, and $\omega\in\Omega$, specifying independent white
noises.
Specifically, we are interested in the solutions
$X=X(t,\omega;x)$ of $(1.1)$
which satisfy the initial conditions
$$
X(0,\omega;x)=x\in B,
$$
for almost all $\omega$, where $B$ is a bounded subset of the the positive
cone in $R^n$
$$
R^n_{+}=\{x=(x_1,...,x_n):x_i>0,i=1,...,n\},
$$
the interior of the feasible region for (1.1). The result in this paper
gives estimates for exit probabilities of $X$ through certain portions of
the boundary of $B$. The result generalizes slightly an earlier result
([9],[10]) of the author, and indicates the use of
multiple Liapunov functions in its application in the study of practical
persistence for stochastic population interaction models([11]).
The set-up thus far is nearly identical to the setting of the
Wentzell-Freidlin exit problem ([6]). Our goal here also is
similar to that in the W-F problem - to obtain features of the exit
probability of $X$ from $B$. However, unlike the W-F problem, we do not
consider the noise necessarily as a small perturbation, and the set $B$ is
not required to be an attractor or a basin of attraction of an equilibrium
of
a corresponding deterministic model. Also, instead of trying to obtain
the
complete exit distribution - which may not be needed to deduce important
biological implications, we will be satisfied, for now, with just
determining some relevant characteristics of the exit probability.
We will assume that (1.1) is non-degenerate in $B$, i. e., that
there is a
positive constant $m$ such that
$$
\xi^TD(x)G(x)G^T(x)D(x)\xi\geq m\|\xi\|^2\eqno(H)
$$
for all $\xi$ in $R^n $and $x$ in $B$ (superscript $T$ denoting
transpose).
It is well-known (see [8], for example) that if the noise
intensity matrix G satisfies {\it (H)} in $B$, then any such solution $X$
must
exit $B$ in finite (random) time
$$
\tau=\tau_x=\inf\{t\geq 0:X(t)\notin B, X(0)=x\}
$$
almost surely (no matter what the form of the drift term $f(x)$). In fact,
the expected exit time
$$
v(x)=E(\tau_x)
$$
is known ([8]) to solve the following boundary value problem:
\begin{eqnarray}
&\mathcal{L} v(x)=\dot{v}(x)+
\frac{1}{2}\mathop{\rm trace}(D(x)G(x)G^T(x)D(x)v_{xx}(x))=-1, \;
x\in \mathop{\rm int}B&\\
&v(x)=0, \quad x\in\partial B&
\end{eqnarray}
where above
$$
\dot{v}(x)=\sum_{i=1}^nx_if_i(x)\frac{\partial
v(x)}{\partial x_i} \quad \mbox{and} \quad
v_{xx}(x)=\Bigl\{\frac{\partial^2v(x)}{\partial x_i\partial
x_j}\Bigr\}^n_{i,j=1}.
$$
(Equation (1.3) is known as Dynkin's equation ([5]).) The
expected exit time has been suggested as characterizing relative
persistence
in the case of models of the form (1.1) and computation techniques have
been proposed ([14],[15]). For simplicity and
clarity here we will assume {\it (H)} and that at least an estimate of
$E(\tau)$
has been determined. The question then is where (through which portion of
the boundary of $B$) does the exit take place. The answer may indicate if
modeled populations are at risk in the near future.
The main point in this paper is that it may be useful to consider
sets $B$ defined by certain multiple Liapunov type functions $V_k$.
In particular we assume the boundary of $B$ is given by pieces of
certain level surfaces of the functions $V_k$
$$
S_{kj}=\{x:V_k(x)=\nu_j\},
$$
where the $\nu_j$ are positive constants. Such functions, sometimes
referred to as
average Liapunov functions, have the property that, if $\nu_j$ is small,
some component of $x$ is small for every $x$ in $S_{kj}$. Construction of
Liapunov
functions has been a principal technique in the investigation of
permanence or
uniform persistence ([13]), in particular, practical persistence, in
deterministic
dynamical system models of population interactions ([2],[3],[4],[7],[12]).
The typical form for such functions is
$$
V_k(x)=\prod_{i = 1}^nx^{\alpha_{ki}}_i ,
$$
where the $\alpha_{ki}$ are constants, at least one of which being
positive
for each $k$. We are ready to state the main result. For simplicity, we
state
the result for a single Liapunov type function $V$.
\begin{theorem} Suppose $V$ is $a$ non-negative $C^2$ function
defined on the bounded set $B$ with
$$
\eta =\inf\{V(x):x\in B\} \quad\mbox{and}\quad
\mu =\sup\{V(x):x\in B\},
$$
and assume that the level surface
$$
S =\{x:V(x)=\eta\}
$$
forms part of the boundary of $B$. Suppose also that condition (H) holds
in $B$.
For any $x\in B$, let $X(t)=X(t,\omega;x)$ be the solution of (1.1) with
$X(0)=x$,
and let
$$
\tau=\tau_x=\inf\{t\geq 0:X(t)\notin B\},
$$
the first exit time of $X$ starting from $x$ in $B$. If there is a
positive constant $c$, such that
\begin{equation}
\mathcal{L} V=\dot{V}+\frac{1}{2}\mathop{\rm tr}(DGG^TDV_{xx}) \geq c
\end{equation}
then
\begin{equation}
\mathop{\rm Prob}\{X(\tau)\in S\}\leq\frac{\mu-V(x)-cE(\tau)}{\mu-\eta}.
\end{equation}
If the level surface
$$
T =\{x:V(x)=\mu\}
$$
forms part of the boundary of $B$, if (H)
holds in $B$, and if there is a positive constant $c$ such that
\begin{equation}
\mathcal{L} V=\dot{V}+\frac{1}{2}\mathop{\rm tr}(DGG^TDV_{xx}) \leq-c
\end{equation}
then
\begin{equation}
\mathop{\rm Prob}\{X(\tau)\in T\}\leq\frac{V(x)-\eta-cE(\tau)}{\mu-\eta}.
\end{equation}
\end{theorem}
\paragraph{Proof:} For $x\in B$, $\eta\leq V(x)\leq\mu$. Let
$$
S =\{x\in\partial B:V(x)=\eta\},\quad T =\{x\in\partial B:V(x)= \mu\},
\quad R =\partial B-S-T,
$$
and further, let
$$
p=\mathop{\rm Prob}\{X(\tau)\in S\}, \quad
q=\mathop{\rm Prob}\{X(\tau)\in T\}, \quad
r=\mathop{\rm Prob}\{X(\tau)\in R\}.
$$
Since {\it (H)} holds in $B$, $p+q+r=1$. Therefore,
\begin{eqnarray}
E(V(X(\tau))) &=& p \eta +\int_RV(x)dP+q\mu \\
&\leq& p \eta+r\mu+q\mu=p \eta+(1-p)\mu\nonumber
\end{eqnarray}
On the other hand, by Dynkin's Formula and (1.5), we get
\begin{equation}
E(V(X(\tau)))-V(x)=E\int_0^\tau {\cal L} V(X(s)))ds\geq cE(\tau)
\end{equation}
Now from (1.9) and (1.10), we have
\begin{equation}
p \eta+(1-p)\mu\geq V(x)+cE(\tau)
\end{equation}
Solving (1.11) for $p$, we obtain
$$
p\leq\frac{\mu-V(x)-cE(\tau)}{\mu-\eta}
$$
which is the first conclusion of the Theorem.
Similarly
\begin{eqnarray}
E(V(X(\tau))) &=& p \eta + \int_RV(x)dP+q \mu \\
&\geq& p \eta+r \eta+q \mu=(1-q) \eta+q\mu \nonumber
\end{eqnarray}
which, similarly to (1.10) and (1.11), together (1.7) leads to
\begin{equation}
(1-q)\eta+q\mu\leq V(x) - cE(\tau).
\end{equation}
Solving for $q$ in (1.13) gives
$$
q\leq\frac{V(x)-\eta-cE(\tau)}{\mu-\eta}
$$
which is the second conclusion, and the proof is
complete.\hfill$\diamondsuit$\medskip
Before giving a population example, we will consider the
following simple
example - the simplest of all SDEs - to get a feeling about what this
result gives us.
In this case we know both the expected exit time $E(\tau)$ and the
probability of exit $q$ explicitly, and so we can compare with the
estimate
(1.8)
obtained in the result above.
\paragraph{Example 1.} We consider $dX=dW, X(0)= x\in (0,1]$. For
some $r$ in $(0, 1)$, take
$$V(x)=x^r.
$$
In this example, $W$ is a scalar Brownian motion (with
$W(0)=0$ a.s.), and thus $X$ is Brownian motion starting at $x$. The set
$B$ is the unit interval $[0,1]$. The boundary of $B$ is made up of the
level sets for $V$ corresponding to $\eta=0$ and $\mu=1$. That $V$ is not
differentiable at $0$ causes no problems. Indeed, for $x\in(0,1],$
\begin{equation}
{\cal L}V(x)=\frac{1}{2}r(r-1)x^{r-2}\leq-\frac{1}{2}r(1-r).
\end{equation}
The probability of exit for $X$ through the top boundary
$$
u(x)=\mathop{\rm Prob}\{X(\tau)=1\}(=q)
$$
solves the boundary-value problem
$$
0={\cal L}u(x)=\frac{1}{2}u''(x),\quad u(0)=0,u(1)=1.
$$
(See ([8]), for example.) The solution is
\begin{equation}
u(x)=x.
\end{equation}
The expected exit time $v(x)=E_x(\tau)$
on the other hand solves the problem:
$$
-1={\cal L}v(x)=\frac{1}{2}v^{ \prime\prime}(x), v(0)=0,v(1)=0.
$$
Its solution is
\begin{equation}
v(x)=x-x^2.
\end{equation}
Now, the conclusion (1.8) of the theorem is equivalent to
\begin{equation}
u(x)\leq\frac{V(x)-0-cv(x)}{1-0}
\end{equation}
Using (1.15) and (1.16) and taking the value of $c$ from (1.14) in
(1.17) results in the estimate of $x$ by a concave function:
\begin{equation}
x\leq x^r-\frac{1}{2}r(1-r)(x-x^2),\ \ 0\leq x\leq 1.
\end{equation}
Inequality (1.18) could be called a ``bow saw inequality''.
\begin{figure}[ht]
\begin{center}
\epsfxsize=9cm
\epsffile{fig1.eps}
\end{center}
\caption{}
\end{figure}
\section{A predator-prey example}
We consider the following predator-prey example
to illustrate that the key conditions (1.2) and (1.3) are obtainable in
the population interaction situation.
\paragraph{Example 2.}
\begin{eqnarray}
dX&=&X[ a-bX-Yh(X,Y)] dt+g_1X dW_1\\
dY&=&Y [-k+mXh(X,Y)] dt+g_2Y dW_2\nonumber
\end{eqnarray}
In (2.1), $a, b, k, m, g_1$ and $g_2$ are positive constants and the
function $h$ satisfies
$$
h(x,y) >0 \quad \mbox{in } R^2_{+}=\{(x,y):x>0,y>0\}.
$$
We consider the multiple Liapunov functions
$$
V_0=mx+y, \quad V_1=xy^{-\alpha},\quad V_2=x^\beta y
$$
where $\alpha$ and $\beta$ are positive constants to be determined
possibly
from the analysis of the corresponding deterministic predator-prey model
\begin{eqnarray}
\frac{dx}{dt}&=&x [a-bx-yh(x,y)]\\
\frac{dy}{dt}&=&y [-k+mxh(x,y)] \nonumber
\end{eqnarray}
For (2.2) the derivative of $V_0$ along solution trajectories satifies
\begin{equation}
\dot{V}_0=mx(a-bx)-ky\leq c-kmx-ky=c-kV_0,
\end{equation}
for a sufficiently large constant $c$. If we take
$$
K=\frac{c}{k},
$$
from (2.3) it follows that all trajectories of (2.2) must eventually
reside in the set
$$
C=\{(x,y)\in\overline{R}^2_{+}:V_0(x,y)=mx+y\leq K\}.
$$
It is of interest, therefore, to choose positive constants
$\eta_1,\eta_2,\mu_1$,and $\mu_2$ so that
$$
B=\{(x,y)\in\overline{R}^2_{+}:\eta_1\leq V_1\leq\mu_1,\eta_2\leq
V_2\leq\mu_2\}\subseteq \mathop{\rm int} C.
$$
The boundary $\partial B$ of $B$ is contained in the set $C$ and is
composed
of the four segments $S_{\eta_1}, S_{\mu_1}, S_{\eta_2},$ and $S_{\mu_2},$
where
$$
S_{\eta_1}=\{(x,y)\in\overline{R}^2_{+}:V_1=\eta_1, \eta_2\leq
V_2\leq\mu_2\}
$$
with the other segments given analogously.
\begin{figure}
\begin{center}
\epsfxsize=9cm
\epsffile{fig2.eps}
\end{center}
\caption{}
\end{figure}
Now, for the choices of $V_1$ and $V_2$ above, we have
\begin{eqnarray}
\dot{V}_1&=&V_1 [a-bx+\alpha k-h(x,y)(y+\alpha mx)]\\
\dot{V}_2&=&V_2 [\beta(a-bx)-k-h(x,y)(
\beta y-mx)].\nonumber
\end{eqnarray}
The noise intensity matrix here is
$DGG^TD$ = $\mathop{\rm diag}\{g_1^2x^2,g_2^2y^2\}$
and thus we obtain
\begin{eqnarray}
\mathop{\rm tr}(DGG^TDV_{1xx})&=& V_1 [\alpha(\alpha+1)g_2^2] \\
\mathop{\rm tr}(DGG^TDV_{2xx}) &=& V_2 [\beta(\beta-1)g_1^2]\nonumber
\end{eqnarray}
Putting together (2.4) and (2.5) we get
\begin{eqnarray}
{\cal L}V_1 &=& V_1 [a-bx+\alpha k-h(x,y)(y+\alpha
mx)+\frac{1}{2}(\alpha(\alpha+1)g_2^2)] \\% (2.6)
{\cal L}V_2 &=& V_2 [\beta(a-bx)-k-h(x,y)(\beta
y-mx)+\frac{1}{2}(\beta(\beta-1)g_1^2)].
\end{eqnarray}
For the special case of Michaelis-Menten dynamics,
$$
h(x,y) = \frac{r}{s + x}
$$
where $r$ and $s$ are positive constants, we have
\begin{eqnarray*}
\lefteqn{ a-bx + \alpha k-h(x,y)(y+\alpha mx) }\\
&=& a-bx+\alpha k-r(y+\alpha mx)/(s+x)\\
&=&a + \alpha k - [b+r\alpha m/(s+x)]x -ry/(s+x)
\end{eqnarray*}
and
\begin{eqnarray*}
\lefteqn{ \beta(a - bx)-k-h(x,y)(\beta y-mx) }\\
&=&\beta( a-bx)-k-r(\beta y-mx)/(s+x)\\
&=&\beta a-k- [\beta b-rm/(s+x)]x - r\beta y/(s+x)
\end{eqnarray*}
both of which are bounded in $B$. Therefore, if
$g_2^2$ sufficiently large, there is a $c_1>0$, such that in $B$
$$
{\cal L}V_1\geq c_1.
$$
Thus, from the theorem
$$
p_1=\mathop{\rm Prob}\{(X(\tau),Y(\tau))\in
S_{\eta_1}\}\leq\frac{\mu_1-V_1(x)-c_1E(\tau)}{\mu_1-\eta_1}.\eqno{(2.7)}
$$
Similarly, if $\beta<1 $ and $ g_1^2$ sufficiently large, there
is a $c_2>0$, such that in $B$,
$$
{\cal L}V_2\leq-c_2,
$$
and so
$$
q_2=\mathop{\rm Prob}\{(X(\tau),Y(\tau))\in
S_{\mu_2}\}\leq\frac{V_2(x)-\eta_2-c_2E(\tau)}{\mu_2-\eta_2}.\eqno{(2.8)}
$$
Inequalities (2.7) and (2.8) give estimates for first tendencies of
predator extinction and prey explosion respectively.
\begin{thebibliography}{00} {\frenchspacing
\bibitem{a1} L. Arnold, Random Dynamical Systems. Springer-Verlag, New
York, 1998.
\bibitem{c1} R. S. Cantrell and C. Cosner, Practical persistence in
ecological models
via comparison methods. Proc. Royal Soc. Edinburgh 126A(1996), 247-272.
\bibitem{c2} R. S. Cantrell and C. Cosner, Practical persistence in
diffusive food
chain models. Natural Resource Modeling 11(1998), 21-34.
\bibitem{c3} Y. Cao and T. C. Gard, Practical persistence for
differential delay
models
of population interactions. Electronic J. Differential Equations (1998)
Conference 01, 1997, 41-53.
\bibitem{d1} E. B. Dynkin, Markov Processes I, II. Springer-Verlag, New
York, 1965.
\bibitem{f1} M. I. Freidlin and A. D. Wentzell, Random Perturbations of
Dynamical
Systems. Springer-Verlag, New York, 1984.
\bibitem{g1} T. C. Gard, Uniform persistence in multispecies population
models. Math.
Biosci. 85(1987), 93-104.
\bibitem{g2} T. C. Gard, Introduction to Stochastic Differential
Equations. Marcel
Dekker, New York, 1988.
\bibitem{g3} T. C. Gard, Stochastic population models on a bounded
domain: exit
probabilities and persistence. Stochastics and Stochastics Reports
33(1990),
63-73.
\bibitem{g4} T. C. Gard, Exit probabilities for stochastic population
models: initial
tendencies for extinction, explosion, or permanence. Rocky Mountain J.
Math.
20(1990), 917-932.
\bibitem{g5} T. C. Gard, Practical persistence in stochastic population
models.
Fields
Institute Communications 21(1999), 205-215.
\bibitem{h1} V. Hutson and K. Mischaikow, An approach to practical
persistence. J.
Math. Biol. 37 (1998), no. 5, 447-466.
\bibitem{h2} V. Hutson and K. Schmitt, Permanence and the dynamics of
biological
systems. Math. Biosci. 111(1992), 1-71.
\bibitem{l1} D. Ludwig, Persistence in dynamical systems under random
perturbations.
SIAM Rev. 17(1975), 605-640.
\bibitem{r1} H. Roozen, An asymptotic solution to a two-dimensional exit
problem
arising in population dynamics. SIAM J. Appl. Math. 49(1989), 1793-1810.
}\end{thebibliography}\medskip
\noindent{\sc Thomas C. Gard}\\
Department of Mathematics\\
The University of Georgia \\
Athens GA 30602 USA \\
e-mail: tgard@uga.edu
\end{document}