Nonlinear Differential Equations, Electron. J. Diff. Eqns., Conf. 05, 2000, pp. 81–90 http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)

Transient effects of stochastic multi-population models ^{∗}

Thomas C. Gard

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.

1 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)

dX = D(X)[f(X) dt + G(X) dW].  (1.1)

Here

X = (X_{1}, X_{2}, ..., X_{n})

with X_{i} representing the ith-species population density, and

D = D(X) = diag{X_{1}, X_{2}, ..., X_{n}};

∗Mathematics Subject Classifications: 92D25, 60H10.

Key words: multi-population models, stochastic differential equations, transient behavior.

\begin{matrix}c \\ ◯ \end{matrix}2000 Southwest Texas State University.

Published October 25, 2000.

more precisely,

X = X(t, ω) = {X_{1}(t, ω), X_{2}(t, ω), ..., X_{n}(t, ω)}

and

W = W (t, ω) = {W_{1}(t, ω), W_{2}(t, ω), ..., W_{m}(t, ω)}

with the W_{j} denoting independent Brownian motion processes, t > 0, and ω ∈ Ω, 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),

F_{i}(x) = \frac{1}{x_{i}} \frac{dx_{i}}{dt} ,  (1.2)

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, ω) = {N_{1}(t, ω), N_{2}(t, ω), ..., N_{m}(t, ω)}

with the N_{j}, t > 0, and ω ∈ Ω, specifying independent white noises.

Specifically, we are interested in the solutions X = X(t, ω; x) of (1.1) which satisfy the initial conditions

X(0, ω; x) = x ∈ B,

for almost all ω, 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

ξ^{T} D(x)G(x)G^{T} (x)D(x)ξ ≥ m∥ξ∥^{2}  (H)

for all ξ 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 (H) in B, then any such solution X must exit B in finite (random) time

τ = τ_{x} = inf{t ≥ 0 : X(t) ∉ 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(τ_{x})

is known ([8]) to solve the following boundary value problem:

Lv(x) = v̇(x) + \frac{1}{2} trace(D(x)G(x)G^{T} (x)D(x)v_{xx}(x)) = −1, x ∈ int B (1.3)
v(x) = 0, x ∈ ∂B  (1.4)

where above

v̇(x) = ∑_{i=1}^{n} x_{i}f_{i}(x)\frac{∂v(x)}{∂x_{i}} and v_{xx}(x) = \{\frac{∂^{2}v(x)}{∂x_{i}∂x_{j}} \}_{i,j=1}^{n}.

(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 (H) and that at least an estimate of E(τ) 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) = ν_{j}},

where the ν_{j} are positive constants. Such functions, sometimes referred to as average Liapunov functions, have the property that, if ν_{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) = ∏_{i=1}^{n} x_{i}^{α_{ki}},

where the α_{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.

Theorem 1 Suppose V is a non-negative C^{2} function defined on the bounded set B with

η = inf{V (x) : x ∈ B} and μ = sup{V (x) : x ∈ B},

and assume that the level surface

S = {x : V (x) = η}

forms part of the boundary of B. Suppose also that condition (H) holds in B. For any x ∈ B, let X(t) = X(t, ω; x) be the solution of (1.1) with X(0) = x, and let

τ = τ_{x} = inf{t ≥ 0 : X(t) ∉ B},

the first exit time of X starting from x in B. If there is a positive constant c, such that

LV = V̇ + \frac{1}{2} tr(DGG^{T} DV_{xx}) ≥ c  (1.5)

then

Prob\{X(τ) ∈ S\} ≤ \frac{μ − V (x) − cE(τ)}{μ − η}.  (1.6)

If the level surface

T = {x : V (x) = μ}

forms part of the boundary of B, if (H) holds in B, and if there is a positive constant c such that

LV = V̇ + \frac{1}{2} tr(DGG^{T} DV_{xx}) ≤ −c  (1.7)

then

Prob\{X(τ) ∈ T \} ≤ \frac{V (x) − η − cE(τ)}{μ − η}.  (1.8)

Proof: For x ∈ B, η ≤ V (x) ≤ μ. Let

S = {x ∈ ∂B : V (x) = η}, T = {x ∈ ∂B : V (x) = μ}, R = ∂B − S − T,

and further, let

p = Prob{X(τ) ∈ S}, q = Prob{X(τ) ∈ T }, r = Prob{X(τ) ∈ R}.

Since (H) holds in B, p + q + r = 1. Therefore,

E(V (X(τ))) = pη + ∫_{R} V (x)dP + qμ  (1.9)
≤ pη + rμ + qμ = pη + (1 − p)μ

On the other hand, by Dynkin’s Formula and (1.5), we get

E(V (X(τ))) − V (x) = E ∫_{0}^{τ} LV (X(s)))ds ≥ cE(τ)  (1.10)

Now from (1.9) and (1.10), we have

pη + (1 − p)μ ≥ V (x) + cE(τ)  (1.11)

Solving (1.11) for p, we obtain

p ≤ \frac{μ − V (x) − cE(τ)}{μ − η}

which is the first conclusion of the Theorem.

Similarly

E(V (X(τ))) = pη + ∫_{R} V (x)dP + qμ  (1.12)
≥ pη + rη + qμ = (1 − q)η + qμ

which, similarly to (1.10) and (1.11), together (1.7) leads to

(1 − q)η + qμ ≤ V (x) − cE(τ).  (1.13)

Solving for q in (1.13) gives

q ≤ \frac{V (x) − η − cE(τ)}{μ − η}

which is the second conclusion, and the proof is complete. ♢

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(τ) and the probability of exit q explicitly, and so we can compare with the estimate (1.8) obtained in the result above.

Example 1. We consider dX = dW, X(0) = x ∈ (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 η = 0 and μ = 1. That V is not differentiable at 0 causes no problems. Indeed, for x ∈ (0, 1],

LV (x) = \frac{1}{2}r(r − 1)x^{r−2} ≤ −\frac{1}{2}r(1 − r).  (1.14)

The probability of exit for X through the top boundary

u(x) = Prob{X(τ) = 1}(= q)

solves the boundary-value problem

0 = Lu(x) = \frac{1}{2}u′′(x), u(0) = 0, u(1) = 1.

(See ([8]), for example.) The solution is

u(x) = x.  (1.15)

The expected exit time v(x) = E_{x}(τ) on the other hand solves the problem:

−1 = Lv(x) = \frac{1}{2}v′′(x), v(0) = 0, v(1) = 0.

Its solution is

v(x) = x − x^{2}.  (1.16)

Now, the conclusion (1.8) of the theorem is equivalent to

u(x) ≤ \frac{V (x) − 0 − cv(x)}{1 − 0}  (1.17)

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:

x ≤ x^{r} − \frac{1}{2}r(1 − r)(x − x^{2}), 0 ≤ x ≤ 1.  (1.18)

Inequality (1.18) could be called a “bow saw inequality”.

Figure 1:

2 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.

Example 2.

dX = X[a − bX − Y h(X, Y)]dt + g_{1}XdW_{1}  (2.1)
dY = Y [−k + mXh(X, Y)]dt + g_{2}Y dW_{2}

In (2.1), a, b, k, m, g_{1} and g_{2} are positive constants and the function h satisfies

h(x, y) > 0 in R_{+}^{2} = {(x, y) : x > 0, y > 0}.

We consider the multiple Liapunov functions

V_{0} = mx + y, V_{1} = xy^{−α}, V_{2} = x^{β}y

where α and β are positive constants to be determined possibly from the analysis of the corresponding deterministic predator-prey model

\frac{dx}{dt} = x[a − bx − yh(x, y)]  (2.2)
\frac{dy}{dt} = y[−k + mxh(x, y)]

For (2.2) the derivative of V_{0} along solution trajectories satifies

V̇_{0} = mx(a − bx) − ky ≤ c − kmx − ky = c − kV_{0},  (2.3)

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) ∈ \overline{R}_{+}^{2} : V_{0}(x, y) = mx + y ≤ K\}.

It is of interest, therefore, to choose positive constants η_{1}, η_{2}, μ_{1},and μ_{2} so that

B = \{(x, y) ∈ \overline{R}_{+}^{2} : η_{1} ≤ V_{1} ≤ μ_{1}, η_{2} ≤ V_{2} ≤ μ_{2}\} ⊆ int C.

The boundary ∂B of B is contained in the set C and is composed of the four segments S_{η_{1}}, S_{μ_{1}}, S_{η_{2}}, and S_{μ_{2}}, where

S_{η_{1}} = \{(x, y) ∈ \overline{R}_{+}^{2} : V_{1} = η_{1}, η_{2} ≤ V_{2} ≤ μ_{2}\}

with the other segments given analogously.

Now, for the choices of V_{1} and V_{2} above, we have

V̇_{1} = V_{1}[a − bx + αk − h(x, y)(y + αmx)]  (2.4)
V̇_{2} = V_{2}[β(a − bx) − k − h(x, y)(βy − mx)].

88 Transient effects of stochastic multi-population models

Figure 2:

The noise intensity matrix here is DGG^{T} D = diag{g_{1}^{2}x^{2}, g_{2}^{2}y^{2}} and thus we obtain

tr(DGG^{T} DV_{1xx}) = V_{1}[α(α + 1)g_{2}^{2}]  (2.5)
tr(DGG^{T} DV_{2xx}) = V_{2}[β(β − 1)g_{1}^{2}]

Putting together (2.4) and (2.5) we get

LV_{1} = V_{1}[a − bx + αk − h(x, y)(y + αmx) + \frac{1}{2}(α(α + 1)g_{2}^{2})] (2.6)
LV_{2} = V_{2}[β(a − bx) − k − h(x, y)(βy − mx) + \frac{1}{2}(β(β − 1)g_{1}^{2})]. (2.7)

For the special case of Michaelis-Menten dynamics,

h(x, y) = \frac{r}{s + x}

where r and s are positive constants, we have

a − bx + αk − h(x, y)(y + αmx)

= a − bx + αk − r(y + αmx)/(s + x)
= a + αk − [b + rαm/(s + x)]x − ry/(s + x)

and

β(a − bx) − k − h(x, y)(βy − mx)

= β(a − bx) − k − r(βy − mx)/(s + x)
= βa − k − [βb − rm/(s + x)]x − rβy/(s + x)

both of which are bounded in B. Therefore, if g_{2}^{2} sufficiently large, there is a c_{1} > 0, such that in B

LV_{1} ≥ c_{1}.

Thus, from the theorem

p_{1} = Prob\{(X(τ), Y (τ)) ∈ S_{η_{1}}\} ≤ \frac{μ_{1} − V_{1}(x) − c_{1}E(τ)}{μ_{1} − η_{1}}.  (2.7)

Similarly, if β < 1 and g_{1}^{2} sufficiently large, there is a c_{2} > 0, such that in B,

LV_{2} ≤ −c_{2},

and so

q_{2} = Prob\{(X(τ), Y (τ)) ∈ S_{μ_{2}}\} ≤ \frac{V_{2}(x) − η_{2} − c_{2}E(τ)}{μ_{2} − η_{2}}.  (2.8)

Inequalities (2.7) and (2.8) give estimates for first tendencies of predator extinction and prey explosion respectively.

References

[1] L. Arnold, Random Dynamical Systems. Springer-Verlag, New York, 1998.

[2] R. S. Cantrell and C. Cosner, Practical persistence in ecological models via comparison methods. Proc. Royal Soc. Edinburgh 126A(1996), 247-272.

[3] R. S. Cantrell and C. Cosner, Practical persistence in diffusive food chain models. Natural Resource Modeling 11(1998), 21-34.

[4] 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.

[5] E. B. Dynkin, Markov Processes I, II. Springer-Verlag, New York, 1965.

[6] M. I. Freidlin and A. D. Wentzell, Random Perturbations of Dynamical Systems. Springer-Verlag, New York, 1984.

[7] T. C. Gard, Uniform persistence in multispecies population models. Math. Biosci. 85(1987), 93-104.

[8] T. C. Gard, Introduction to Stochastic Differential Equations. Marcel Dekker, New York, 1988.

[9] T. C. Gard, Stochastic population models on a bounded domain: exit probabilities and persistence. Stochastics and Stochastics Reports 33(1990), 63- [10] T. C. Gard, Exit probabilities for stochastic population models: initial ten-

dencies for extinction, explosion, or permanence. Rocky Mountain J. Math. 20(1990), 917-932.

[11] T. C. Gard, Practical persistence in stochastic population models. Fields Institute Communications 21(1999), 205-215.

[12] V. Hutson and K. Mischaikow, An approach to practical persistence. J. Math. Biol. 37 (1998), no. 5, 447-466.

[13] V. Hutson and K. Schmitt, Permanence and the dynamics of biological systems. Math. Biosci. 111(1992), 1-71.

[14] D. Ludwig, Persistence in dynamical systems under random perturbations. SIAM Rev. 17(1975), 605-640.

[15] H. Roozen, An asymptotic solution to a two-dimensional exit problem arising in population dynamics. SIAM J. Appl. Math. 49(1989), 1793-1810.

Thomas C. Gard Department of Mathematics The University of Georgia Athens GA 30602 USA e-mail: tgard@uga.edu