Stochastic models for toxicant-stressed populations

Bulletino/ MathematicalBiologyVol. 54, No. 5, pp. 827 837, 1992. Printed in Great Britain.

0092-8240/9255.00 + 0.00 Pergamon Press Ltd © 1992 Society for Mathematical Biology

STOCHASTIC MODELS FOR TOXICANT-STRESSED POPULATIONS THOMAS

C.

GARD

Department of Mathematics, The University of Georgia, Athens, GA 30602, U.S.A. (E.mail: [email protected].)

We obtain conditions for the existence of an invariant distribution on (0, oe) for stochastic growth models of Ito type. We interpret the results in the case where the intrinsic growth rate is adjusted to account for the impact of a toxicant on the population. Comparisons with related results for ODE models by Hallam et al. are given, and consequences of taking the Stratonovich interpretation for the stochastic models are mentioned.

1. Introduction. Dynamical models of population growth with time-dependent parameters are considered more realistic than constant parameter versions. Intrinsic growth rates and carrying capacities, for example, naturally vary with time-evolving environments. However, such nonautonomous models are generally less tractable mathematically, especially when complete time histories of the parameters are required to decide important qualitative issues. So a pertinent question is, to what extent does partial information about the parameters, namely parameter statistics, suffice to determine qualitative behavior? Do average growth rates, for example, characterize persistence? Recently, Vance (1990) and Vance and Coddington (1989) have addressed the general question for deterministic ordinary differential equation population models: d x / d t = xf(t, x).

(1)

They specified conditions under which the dynamics given by the nonautonomous equation (1) could be mimicked by that of the associated autonomous model: d x / d t = xg(x)

(la)

where g is the long-time average off. Hallam and Ma Zhien (1986) and Ma Zhien et al. (1989) have considered a more specific question; models of stressed population growth of the form: d x / d t = x/g(x ) Jr(t)-fix)x]

(2)

where the intrinsic growth rate r is given by: 827

828

T.C. GARD

r(t)=r o --rlCT(t) with r o representing the constant unstressed intrinsic growth rate, ca-, some time-varying measure of the toxicant stress on the population (e.g. the concentration of the toxicant in a typical individual in the population), and r I , a units conversion parameter. They established threshold values for persistence of the population in (2) in terms of long-time averages of cT(t). In this paper we consider stochastic analogues of (2) of the form:

dX=X/g(X) [ ( a - f ( X ) X ) dt +b dW]

(3)

which incorporate explicit r a n d o m variation in the measure of toxicant stress cT and hence in the intrinsic rate r(t). Equation (3) is a stochastic differential equation which arises from (2) when r(t) is taken as:

a+bN(t), where N denotes standard white noise, the derivative, in a generalized sense, of Brownian motion W; in terms of the toxicant stress, we are taking here:

cT(t ) =

CO -- ClN(t

)

(4)

and so:

a=ro-rxCo,

and

b = r l c 1.

In (4), c o represents the average value of the toxicant stress, and c 1 denotes the intensity of r a n d o m fluctuations of the stress about the average; since for each t, N(t) is a standard normal r a n d o m variable, Cr(t ) is also normal but with mean c o and variance cf. By the law of large numbers, c o also represents the long time average of cv. Equation (3) constitutes a generic model for the effects of a toxicant on a population in a r a n d o m environment; it allows all basic growth models treated in the literature as deterministic counterparts (2) (see Hallam, 1986, for example), and it postulates the effect of the toxicant simply as biased time-varying r a n d o m adjustment to the population's intrinsic growth rate. In contrast, for example, is the specific stochastic toxicant-population interaction model treated in Gard (1990), where the main source of the r a n d o m toxicant stress fluctuation is assumed to be variation in the toxicant washout rate from the environment and the interaction is taken as an uptake mechanism of the toxicant in environment. Our goal in this paper is to characterize persistence of the population under stress of the toxicant in terms of the statistics c o and c 1 of the stress measure. The results obtained here are stochastic analogues of those given for the deterministic counterpart models (2) in Hallam and M a Zhien (1986) and M a Zhien et al. (1989).

TOXICANT-STRESSED

POPULATIONS

829

2. Preliminaries and Deterministic Results. Persistence is the model analogue for survival of the population species represented. Although there are, in general, a n u m b e r of different mathematical definitions for persistence (see Butler et al., 1986; M a Zhien et al., 1989), for a u t o n o m o u s population growth models such as:

dx/dt=x/9(x)[a-f(x)x]

(5)

persistence a m o u n t s to global (relative to the positive axis, 0 < x < oo) stability of a positive equilibrium, or equilibria configuration. In (5),fand 9 are assumed at least to be positive continuous functions defined on the positive axis, so that (5) represents a generalization of the well-known logistic equation. If: g(O) > 0 and a-f(x)x < O, for all sufficiently large x,

(6)

then (5) has at least one positive equilibrium, and persistence holds. F o r stochastic versions of (5):

dX=X/g(JO [(a-fiX)X) d r + b dW]

(7)

the analogue of persistence is positive recurrence. This implies that for every e > 0 , there is an M > 0 , such that: lim inf p X { M - 1 <~X(t) <. M } >~ 1 -- E

(8)

t--* oo

where px{x(t)eB} denotes the probability that the solution X of (7) with X(0) = x is in the set B at time t. Positive recurrence is equivalent to the existence of an invariant distribution p with nowhere zero density on the interval (0, oo) which is stable, i.e.: limit P x { X ( t ) ~ B} = p(B) for any Boreal set B; the density p of the distribution # is a non-negative solution to the F o k k e r - P l a n c k equation: 1 d2 d 2 dx 2 [(bx/g(x))2P] +dx [(x(a-f(x)x)/g(x))p] = 0

(9)

on the interval (0, oo) that satisfies:

f oop(x) dx = 1. 0

The condition that such a solution exists is that the expression:

(10)

830

T . C . GARD

S(x) = (g(x)/bx) 2 exp (b2i I f ~ -ag(y) d y - f x f(y)g(y)dy 1) Y

(11)

is integrable on (0, oo). In this case, if we take:

S=

~0°°

S(x) dx

we have:

p(x) = S(x)/g. On the other hand, if:

R(x)=exp(-

If Xag(Y)dY- f Xf(Y)g(y)

is integrable near zero, but not integrable near oo, then any solution X(t) of (7) with X(0)=xe(0, oo) tends to zero with probability one, i.e. the population goes extinct almost surely. This follows from the fact that the probability Pe that X(t)-,O, as t ~ oo is given by: &

_

T(O)--

where T(x) = fx R(y) dy.

(Gard, 1988, for example, contains details and basic results on stochastic differential equations.) In the next two sections we establish persistence thresholds for the stochastic model (7), by giving conditions o n f a n d g which guarantee each of the above alternatives. Then we interpret the results in the toxicant-stress situation. 3. The g ( x ) - I Case. Most of the deterministic continuous population models (5) proposed in the literature exhibit approximate linear growth at small population levels. But general observations, such as that population carrying capacity and steady state may be distinct due to food consumption varying not only with population density but also with population growth rate (see Hallam, 1986), have suggested models with nonconstantf(x) and 9(x). The function f(x) in (5) in particular can represent the proportion of the food consumption rate at density x to the rate consumed at equilibrium. With fix) = a/K, K denoting the carrying capacity, the equation:

dx/dt = x(a -fix)x)

(12)

TOXICANT-STRESSED P O P U L A T I O N S

831

corresponds to the well-known Verhulst or logistic equation. But (12) has also been proposed with: f i x ) = (log x)/x (Gompertz, 1925), f ( x ) = x ~, - 1 < 7 < 0

(Goel et at., 1971),

f i x ) = (e-~/X)/x (Gallopin, 1971), and f ( x ) -- yfl++a~x (DeAngelis et al., 1975).

In the latter two cases, it is assumed that a < 1 and a < 6, respectively, so for each of these models, (6), and hence persistence holds. Since xf(x) is monotone in each of the first four cases cited above, there is a unique positive equilibrium in these cases, but the last example permits multiple positive equilibria for some choices of parameters. We also note that for each of the above examples: f i x ) > a/x, for sufficiently large x;

so, in particular for any number b, there is an e > 0 such that: f i x ) >i (a - b 2/2 + E)/x, for sufficiently large x.

(13)

For corresponding stochastic models: d X = X[(a -f(JOX) dt + b d IT],

(14)

now we have the following result. THEOREM 1. Let f(x) be a positive continuous function defined on (0, oo) which satisfies (13), where a - b 2 / 2 > 0 , and assume f is integrable near O. Then the solution X(t) of the stochastic differential equation (14) with X(0)= xe (0, oo) is positive recurrent on (0, oo ) and (14) has an invariant distribution. If, on the other hand, a - b 2 / 2 < O, and the other conditions on f hold, then: limit X(t)=0, w.p. 1. t~0

(15)

Before giving the proof, we interpret Theorem 1 in the toxicant-stress context. Recall from the introduction that: ,

a = r o - - r l C o , and b = r l c 1 ,

where r o, Co, c I and r~ represent the population's (unstressed) intrinsic growth rate, the average value of the toxicant stress measure, the intensity of random stress fluctuations, and a units conversion factor, respectively. The theorem

832

T . C . GARD

predicts survival in the sense of (8) if the unstressed intrinsic growth rate exceeds a weighted combination of the m e a n and variance of the measure of toxicant stress: r 0 :> r 1(cO +

r~c2/2),

and almost certain extinction if this inequality is reversed. We emphasize that the persistence/extinction criterion depends on both the mean and variance of the toxic stress. We will c o m m e n t further on this feature in the discussion section at the end of the paper. Proof. The function S given by (11) here has the form:

for some positive constants C and x 0. Since f is integrable near 0, and a-b2/2>O, there are positive constants C o and 6 such that:

S(x) <~Cox - 1+~, for x near 0

(17)

and so S is integrable near 0. O n the other hand, for large x we have, by (13):

f

x' f(y ) dy >1(a - b 2/2 + e) log(x/Xo), o

so that

S(x)<~C1 x - l - ~ , for large x

(18)

with C 1 and y positive constants. So S is integrable near oo. Combining (17) and (18) with the fact that S(x) is continuous on (0, oo), we see that S is integrable on (0, oo) and the result is established. In fact, letting:

= f o S(x) dx we have p(x)= S(x)/S. For the proof of (15), as noted in the previous section:

(

px _ lim X(t) = 0 \t~oo

where

-

r(x)-r(oo) T(O)

T(oo)'

(19)

TOXICANT-STRESSED POPULATIONS

T(x):fxR(y)dY:fxC2y-2a/b2 exp (~--~fY f(z) dz) dy.

833

(20)

Since f satisfies (13), (20) implies that: T(oo) = lim

R(y) = - or.

x---~ oo

Because f is integrable near zero, and a - b 2/2 < 0:

fx

R(Y) dy~x

2./b2+1 for x n e a r 0

and so T(0) -= x-~o+limfx

R(y) dy < + ~.

Therefore, from (19), we have:

px (lim X(t)=0) = 1. t--~ oO

This completes the proof of the theorem. 4. The Nonconstantg(x) Case. Nonconstant choices for the function g(x) can arise from the consideration that per capita intrinsic growth rates are density dependent. Specific candidates g(x)=g+hx, and g(x)=x 1/2 have been suggested by Smith (1963) and Rosenweig (1971), respectively. In this section we treat two classes of functions g(x) which include these examples as special cases. For the first case the result is a direct generalization of Theorem 1 which requires the following assumptions on the functions fix) and g(x):

f(x) is a non-negative continuous function on the interval [0, ~ ]

(G1)

g(x) is a positive continuous function on the interval [0, ~ ) , with g(0)=g (G2) for some ~>0, and 0 < M < 1/(1-2~/b2), for sufficiently large x

f(x) >~(a-- b 2/2 + E)/x

(G3)

and

ff(x)<~M.

(G4)

834

T.C.

GARD

THEOREM 2. If the functions f(x) and g(x) satisfy assumptions (G1)-(G4), and if ag - b2/2 > 0, then equation (7) has an invariant distribution. On the other hand, if (G1)-(G4) hold, and if a g - b Z / 2 < 0 , then every solution X(t) of (7) with X(0)=xe (0, oo) tends to zero at t--.oo with probability 1. Proof. The proof is similar to that for Theorem 1. The condition ag b2/2 > 0 gives integrability of S(x) near zero, just as in Theorem 1 using the trivial estimate 9 for the function g(x) near zero. And (G3) and (G4) guarantee integrability near aD. Indeed, suppose these conditions hold for all x >/x 1 : then there is a positive constant K such that for all such x: -

S(x) <~Kx ~I(1- 2~/b2)-2

-

(21)

since

S(x)=Kl(g(x)/x)2 exp(~--y fx I [a/y-f(y)]g(y)dy)

<;;

<~KI(M/x) 2 exp ~ -

[a/y--fly)] dy

1

)

and

fx

~ [a/y-f(y)] dy ~<(b212-E)ln(xlxl). 1

Thus (G4) implies that the exponent in (21) exceeds - 1, and integrability at oo follows. The extinction part of the theorem also follows as in Theorem 1. • We note that the continuity at zero condition on 9(x) in the above theorem is not essential; the proof makes it clear that the result stands for g(x) continuous on (0, oo) with 9(0) replaced by: lim inf g(x) and lim sup g(x) x~O +

x~O +

in the first and the last parts of the theorem, respectively. By taking into account this observation, examples such as g(x)= x - 1/2 can be treated. In fact in this situation we have ag-b2/2>O and (G4) satisfied for any choices of positive a and b, and Theorem 2 reduces to: COROLLARY 3. If g(x)= x -° with p > 0 , and/ff(x) satisfies (G1) and (G3), then (7) has an invariant distribution on (0, oo). T h e o r e m 2 does not apply to either the Smith choice g(x)=g+hx or Rosenzweig's model g(x) = x 1/2 which violate the boundedness condition (G4).

TOXICANT-STRESSED POPULATIONS

835

However if we restrict f(x) further, we can relax (G4) to include some unbounded choices for 9(x). The following is a sample of possible results in this direction. THEOREM 4. The conclusion of Theorem 2 holds if f(x) and g(x) satisfy the conditions of Theorem 2 except that (G3) and (G4) are replaced by

for some constants C, e, 6 > O, 0 <%p < 1, and all sufficiently large x f(x) >/(a + e)/x

(G3')

6<~g2(x)<~CxP e x p ( ~ fX g(Y) dY)

(G4')

and

Proof. (GY) and the second part of (G4') imply, similarly to the proofs of Theorems 1 and 2, that, for some positive constant K: S ( x ) <...Kx p - 2

for sufficiently large x. Thus S is integrable on (0, oo). Also (G4') guarantees that R(x) is not integrable at oo, a condition that is needed to get almost sure extinction when ag-b2/2O in Theorems 2 and 4 is just the persistence condition in Theorem 1 with a and b replaced by a/g and b/g, respectively. The interpretation of these results for the toxicant-stress situation then is basically the same as indicated in Section 3. Under fairly mild conditions on the functionsf(x) and g(x), survival [-in the sense of the existence of a stable invariant distribution on (0, oo)] occurs if:

ro>q(Co+rtc~/2g)

(22)

and almost sure extinction results if (22) is reversed. We see in (22) that the value of the density dependent adjusted intrinsic growth rate at zero can offset or enhance the effect of variation in the toxicant stress; larger values of g allow larger values of the stress variance c~. In any case, the persistence threshold is given by:

836

T.C. GARD

ro/q = co + rlc2/2g; it depends on both the m e a n and variance of the stress. The results here are based on the Ito interpretation of the stochastic equation (7). It is well k n o w n that this is not the only interpretation possible for (7). More precise model formulation is necessary to distinguish whether Ito or some other interpretation of (7), such as Stratonovich, is appropriate. In some cases this argument amounts to establishing confidenceqn a related m o r e basic difference or differential equation model, from which (7) usually follows via some limiting procedure. For example, if one has a good deal of confidence in the differential equation model (5) in the unstressed situation, rather than some even more fundamental model, the result of W o n g and Zakai (1965) suggests that the Stratonovich interpretation of (7) m a y be appropriate. Fortunately our results can be easily reinterpreted in this situation. Indeed, one can transform the Stratonovich interpreted equation (7) into the corresponding Ito equation:

dX=(X/g(X))I(a-Xf(X)+ (b2/2) (g(X~--\ -g2-~Xg'(X)"~))dt+b dW 1 which can be rewritten in the form: dX=

(J(/g(X)) [(a +

b 2/(2g) -

XF(X) ) dt + b d 142J

(7s)

where

g(X)=g+ Xh(X) and

( h )2 (h(X)+(X/2)(h'(X)+h2(X)/g)]"

F(X)=f(X)+ g+Xh(X)

The results for the Stratonovich version of (7) are obtained by applying our results to the corresponding I t o equation (7S), i.e. by replacing a with a+b2/(2g) and f(x) with F(x). In particular we notice that the persistence condition becomes:

0 < g(a + b 2/(2g))

--

b 2/2

=

ga

or just a>0 if g = g(O) is assumed positive. Returning to the toxicant-stress model we recall that

TOXICANT-STRESSED P O P U L A T I O N S

a =

r o -- rlc

837

o

so t h a t positivity of a is e q u i v a l e n t to c o < ro/r I ;

the t h r e s h o l d for p e r s i s t e n c e for the S t r a t o n o v i c h m o d e l is c h a r a c t e r i z e d b y the m e a n of the stress alone. T h i s m a k e s sense to s o m e e x t e n t intuitively b e c a u s e a c c e p t i n g the S t r a t o n o v i c h m o d e l a m o u n t s to a s s u m i n g n e a r l y valid c o n t i n u ous d a t a . T h e results are n o t s u r p r i s i n g l y in a c c o r d with the d e t e r m i n i s t i c results of H a l l a m a n d of V a n c e , b e c a u s e these results are for differential equation models.

LITERATURE Butler, G., H. I. Freedman and P. Waltman. 1986. Uniformly persistent systems. Proc. Am. math. Soc. 96, 425-430. DeAngelis, D. L., R. A. Goldstein and R. V. O'Neill. 1975. A model for trophic interaction. Ecology 56, 881-982. Gallopin, G. C. 1971. A generalized model of a resource-population system: I. General properties. II. Stability analysis. Oecologia 7, 382-413; 7, 414-432. Gard, T. C. 1988. Introductions to Stochastic Differential Equations. New York: Marcel Dekker. Gard, T. C. 1990. A stochastic model for the effects of toxicants on population. Ecol. Modellin 9 51,273-280. Goel, N. S., S. C. Maitra and E. W. Montroll. 1971. On the Volterra and other nonlinear models of interacting populations. Rev. Mod. Phys. 43, 231-276. Gompertz, B. 1925. On the nature of the function expressive of the law of human mortability. Phil. Trans. 115, 513 585. Hallam, T. G. 1986. Population dynamics in a homogeneous environment. In Mathematical Ecology, T. G. Hallam and S. A. Levin (Eds). Berlin: Springer-Verlag. Hallam, T. G. and Ma Zhien. 1986. Persistence in population models with demographic fluctuations. J. math. Biol. 24, 327 339. Ma Zhien, Song Baojun and T. G. Hallam. 1989. The threshold of survival for systems in a fluctuating environment. Bull. math. Biol. 51, 311-323. Rosenzweig, M. 1971. The paradox of enrichment: destabilization of exploitation ecosystems in ecological time. Science 171,385-387. Smith, F. E. 1963. Population dynamics in Daphnia magna and a new model for population growth. Ecology 44, 651-663. Vance, R. R. 1990. Population growth in a time-varying environment. Theor. Pop. Biol. 37, 438-454. Vance, R. R. and E. A. Coddington. 1989. A nonautonomous model of population growth. J. math.Biol. 27, 491-506. Wong, E. and M. Zakai. 1965. On the convergence of ordinary integrals to stochastic integrals. Ann. math. Stat. 36, 1560-1564. R e c e i v e d 27 D e c e m b e r 1990