\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}