TIIEORETICAL

POPULATION

Random

12, 140-l 78 (1977)

HIOLOGY

Environments MICHAEL

and Stochastic

Calculus

TLIRELLI*

Center for Quuntitative Science in Forestry, Fisheries and Wildlifr, University of Washington, Seattle, Washington 98195 Received October

22, 1976

This paper investigates the use of heuristically derived stochastic differential equations (SDEs) as models in population biology. It is stressed that these equations are best viewed as approximations for more realistic, but often analytically intractable, models. A criterion is presented for determining which interpretation (e.g., Ito or Stratonovich) is likely to serve as the most useful approximation. Several limit theorems are presented which illustrate the use and implications of this criterion. In particular, it is shown that the solutions to sequences of models which approach a given SDE may converge to a diffusion process which corresponds to no solution of the SDE. However, arguing that in population biology the SDEs are generally serving as approximations to stochastic difference equations with autocorrelated noise, it is shown for a variety of models that the Ito calculus may provide a more useful approximation than the Stratonovich. Important limitations of this result and on the use of SDEs are indicated. These findings and observations are compared with those of several papers in the recent literature.

1.

INTRODUCTION

There IS currently a great deal of interest in incorporating randomly varying parameters into the standard models of theoretical ecology and population genetics (see, for instance, Feldman and Roughgarden, 1975; Kciding, 1975; Karlin and Leiberman, 1974; Gillespie, 1975; and their bibliographies). The biological questions which motivate this interest involve the effects of environmental unpredictability on population sizes, species interactions, and gene frequencies. This paper concentrates on a particular class of models, heuristically derived stochastic differential equations (SDEs), which have become fairly popular as tools for investigating these effects. WC use the term SDE in a narrow sense, referring only to models of the form dS/dt

= f(X)

- g(X) A(t),

where A(.) is a so-called Gaussian white noise process. *Present 95616.

address:

Department

of Genetics,

Copyright Q 1977 by Academic Press, Inc. All rights of reproduction in any form reserved.

University

of California,

Davis,

Calif.

140 ISSS

OIMO-5809

STOCHASTIC

141

CALCULUS

The standard “derivation” of such models is as follows. deterministic differential equation

We begin with a

d-wt = flu7 + rfz.(W in which the parameter y is assumed to reflect the state of the environment. To incorporate environmental variation, y is replaced by a mean value 7, which represents the average state of the environment, plus a scaled white noise variation term ~/l(t) in which the scale parameter a reflects the level of fluctuation. The resulting model has the form of (1.1) with f = fi + Tf2and g = ofi . The difficulty with and interest in such models can be traced to the nature of the white noise term fl. In the derivation above, it is treated as an ordinary stochastic process modeling environmental fluctuation. However, this “process” is not a stochastic process in the usual mathematical sense and it has no direct physical interpretation. Instead it is a “generalized stochastic process” (cf. Arnold, 1974, pp. 50-56) and its precise definition is closely linked to that of the Dirac a-function (which is really not a function but a “generalized function”). In technical terms, .4 is a “generalized derivative” of the standard Brownian motion process, W, and is sometimes symbolically denoted dW/dt. In crude terms, /l can be described as a one-parameter collection of independent Gaussian random variables with mean zero and infinite variance, i.e., a stationary, zero mean Gaussian process whose covariance kernel is the Dirac delta function. To illustrate the connection between the precise and crude descriptions recall that W has independent Gaussian increments (i.e.. for fixed times 0 < t, < “’ < t, the random variables W(tJ,..., W(tn) - W(t& are independent and Gaussian) and that for s < t, W(t) - W(s) has mean zero and variance t - S. Hence for s < t and h < t - s the difference quotients w(t + h) h

qg

and

W(s + h) h

W(s)

(1.2)

are independent Gaussian random variables with mean zero and variance l/h. The heuristic description of /l is motivated by letting h --j 0 in (1.2). The complications which surround fl and hence (1.1) arise from the fact that this limit cannot be usefully defined within the context of ordinary stochastic processes. The difficulty of assigning a meaning to dW/dt is implicit in the more conventional mathematical notation for (1.1) dX =:

f(X)

dt f g(X) dW.

(1.3)

Equations of this form were introduced by Ito (see references in McKean, 1969). His purpose was not to consider stochastic analogs of deterministic equations

142

MICHAEL

TURELLI

but rather to represent general diffusion processes as transformations of a Brownian motion, W. To accomplish this he proposed a particular interpretation for (1.3). However as discussed in Section 2, (1.3) has many other possible interpretations and which diffusion process it represents depends on which interpretation is employed. Whereas the work on (1.3) by Ito and other mathematicians revolves around representing and manipulating a known diffusion process, the heuristics leading to (1.1) are aimed at finding stochastic counterparts for deterministic processes. The purpose of adding Gaussian white noise to the deterministic model rather than adding an ordinary stochastic process is so that the stochastic output will be a diffusion and hence amenable to some standard analytical techniques. (This issue is discussed in Wong, 1971, Chap. 4.) The price that is paid for this convenience is that having arrived at (1.3) d irectly by inserting a noise process which cannot be physically realized, it is not obvious which if any interpretation of it will provide useful information concerning the effects of environmental stochasticity. Despite the facts that the physical assumptions implicit in the heuristic derivation of (1.1) are not at all clear and that it is mathematically ambiguous, models of this form enjoy a good deal of attention in recent theoretical literature. Certainly the most ambitious (and controversial) application is May’s construction of a theory to predict the limiting degree of similarity, consistent with coexistence, between competitors whose limiting resource fluctuates randomly in abundance (May and MacArthur, 1972; May, 1973; May, 1974). Mathematical deficiencies in the foundation of this theory have been discussed by Feldman and Roughgarden (1975), Abrams (1976), and most recently elaborated by Turelli (1977). SDE growth models have also been studied by Capocelli and Ricciardi (1974) Tuckwell (1974), Feldman and Roughgarden (1975), Ladde and Siljak (1975) Ludwig (1975) and Levikson (1976). Also, Tuckwell (1976) has used SDEs to reanalyze the gene frequency dynamics of a single diallelic locus for which fitnesses vary stochastically. In most of these papers, the authors have pointed out the difficulties associated with the ambiguity of these models. The most frequently used interpretation of (1.3) other than Ito’s is due to Stratonovich (1966). Both interpretations have been employed in the biological literature. The crux of the difficulty surrounding the use of heuristically derived SDEs is that the biological conclusions reached often seem to depend critically on which mathematical convention is employed. This dilemma has been clearly illustrated by Feldman and Roughgarden (1975) who point out that the results of May’s theory of limiting similarity can be completely overturned by using the Stratonovich interpretation of his equations rather than the Ito interpretation which he chose. Despite the widespread awareness of this difficulty, there has been no systematic discussion of how one might decide which, if any, interpretation of the SDE (1.3) is likely to provide useful information as an approximate description

STCCHASTIC

CALCULUS

143

of a particular biological phenomenon. Certainly if these heuristically derived models are to serve any productive purpose, some precise information regarding their possible utility and limitations is necessary. It is our intention to provide some of the necessary information. We begin in the next section by discussing briefly the origin of the difficulty in interpreting (1.3) and proposing a criterion to decide which interpretation is likely to be most useful in a particular context. Some fairly general convergence theorems are presented which are pertinent to the application of this criterion. Unfortunately, these theorems do not cover the caseswe feel are most likely to be relevant in population biology. In Section 3, the simple but interesting equation dX = rXdt

+ aXdW

modeling exponential growth in a random environment is discussed in detail. Limit results are presented which indicate that even in this very simple setting no interpretation of the SDE may be wholly adequate. However, particular properties of interest, namely conditions for population extinction and explosion, are correctly prescribed (according to the criterion of Sect. 2) by the Ito solution. Generalizations of this result are given in Section 4 which imply that of the two standard interpretations of (1.3), the Ito is likely to be more useful in population biology. In Section 5, this conclusion is compared to those of other investigators and the reasons for discrepancies are examined. The last section attempts to put these results into perspective. Implicit in posing the problem of deciding on an appropriate interpretation of (1 .l) is the assumption that this model is meaningful. Given the nonphysical nature of white noise, this is not self-evident. Some severe limitations of (1.1) as a heuristically derived diffusion approximation are indicated. We also discuss some obvious limitations of the criterion proposed in Section 2 and of our applications of this criterion. In all but the last section, the SDEs considered are models of population growth. In Section 6, we comment briefly on Tuckwell’s (1976) genetic application. It should be noted that heuristically derived SDEs had been used and discussed extensively by physicists and electrical engineers well before their appearance in population biology. In particular, much has been written concerning the problem of interpreting these ambiguous models (see, for instance, Gray and Caughey, 1965; Saaty, 1967, Chap. 8; Mortensen, 1969; Wong, 1971, Chap. 4; and McShane, 1974). However, because of basic differences in the phenomena being modeled (see Sect. 2), these discussions have inherent biases which detract from their usefulness to population biologists. Thus although the criterion presented in Section 2 has been discussed previously, the earlier treatments have omitted one or more of the limit results discussed below which are pertinent to its application. Although all of the mathematics employed is straightforward, some is rather technical. The basic ideas and results are presented heuristically in the text.

144

MICHAEL

TURELLI

In particular, we do not specify the modes of convergence utilized in the various convergence results discussed. However, all of the theorems we refer to which concern convergence of stochastic processes are or can be stated in terms of “convergence in distribution.” This concept is discussed in the Appendix.

2. A CRITERION Stochastic population

FOR CHOOSING

A

STOCHASTIC

CALCULUS

growth models of the form dX =:: f(X)

dt + g(X) dW

(2.1)

which are heuristically derived, as described in the Introduction, from deterministic differential equations are the focus of this paper. An essential difficulty with such models is that they are mathematically ambiguous. The SDE (2.1) is actually shorthand for the random integral equation

X(t) - X(t,) = [lf(X(s)) - ‘0

ds -+ if ,q(X(s)) dW(s). . to

The ambiguity arises from the fact that Brownian rapidly to define

motion,

W, oscillates

too

(2.2) in the usual Riemann-Stieltjes (or Lebesgue-Stieltjes) sense. The result is an infinity of definitions for (2.2) each of which gives rise to a different solution for (2.1). Only two of these definitions, those attributed to Ito and Stratonovich, yield solutions which are likely to be useful as models. Associated with each of these standard integration procedures for (2.2) is a “differentiation rule” for the stochastic processes obtained as solutions of (2.1). Following the terminology of Feldman and Roughgarden (1975), the combined rules of integration and differentiation will be referred to as the Ito and Stratonovich calculi. The properties which distinguish these calculi as modeling tools are discussed below. One approach to sorting out the ambiguity inherent in (2.1) is to keep in mind that it is a highly idealized model. Its principal attraction is that once a suitable definition for (2.2) is chosen, the solution will be a well-defined diffusion process whose properties can be investigated by standard analytical methods. Recall that apart from boundary behavior a diffusion process, X(.), is characterized by two functions its infinitesimal mean, M(x), and infinitesimal variance, V(x). In most cases, they can be calculated according to

STOCHASTIC

145

CALCULUS

and

V(r) = FT (IPZ)E[(X(t + h) - x(t))” I X(t) = x]. Intuitively, M(x) describes the systematic (deterministic) pressure the process experiences when its value is x while I/(x) reflects the amount of stochastic jostling it undergoes. If the Ito interpretation of (2.2) is employed, the diffusion which satisfies (2.1) has infinitesimal mean

JfI(4 = f(4 and infinitesimal

(2.3)

variance V(x) = g”(x).

The Stratonovich mal mean is

solution has the same infinitesimal

variance, but its infinitesi-

Despite the mathematical convenience of working with diffusion processes, the phenomenon being modeled by (2.1) may in fact be more accurately described by a stochastic difference equation rather than a differential equation, and the random perturbations to which it is subject may be autocorrelated rather than being white noise. The difficulty, of course, is that the more realistic model is liable to be analytically intractable. Nonetheless, it is important to recognize that if the stochastic differential equation (2.1) is replaced by a model involving differences rather than differentials, or autocorrelated rather than white noise, the new model will be mathematically unambiguous since the troublesome “dW” of (2.1) will have vanished. This observation suggests a criterion for choosing a solution procedure for (2.1). Let P denote a definition for (2.2) and XP the corresponding solution of (2.1). For XP to be considered a reasonable approximation to the behavior of a particular phenomenon, it seems natural to require that if (2.1) is replaced by a more realistic (and hence unambiguous) model which is “close” to it, the resulting solution should be “close” to X, . An obvious way to apply this criterion is to consider a sequence of models, indexed by N, with solutions denoted X,v . Each model in the sequence should involve differences and/or correlated noise depending on the nature of the phenomenon being described. As N --) CO the models should approach (2.1), i.e., differences should converge to differentials and autocorrelation in the noise should vanish. If a solution to the idealized model (2.1) is to be a useful approximation, its behavior should agree with the behavior of XV for large iv. Some limit results which characterize

146

MICHAEL

TURELLI

conditions under which the Ito and Stratonovich solutions of (2.1) provide acceptable approximations according to this criterion are presented below. These examples also clarify our intuitive notion of convergence of a sequence of models to (2.1). Suppose the process of interest is more accurately described by a difference equation than a differential equation but that the noise it experiences can be reasonably considered uncorrelated. According to the proposed criterion, we can determine which, if any, solution of (2.1) is likely to be useful by examining the solutions to a sequence of stochastic difference equations, involving uncorrelated noise, which converge to (2.1). The most obvious such sequence (indexed by N = 1, 2,...) is (N) Xt+(l/N)

+ g(X!N))(Wt+(l,N) - XcN) I = ft~~“‘)(W) t = 0, l/N, 2/N ,..., for

- wt)

which is obtained by simply replacing the differentials of (2.1) by differences (cf. Arnold, 1974, Sect. 10.3). However, in order to obtain a sequence of models that behave sensibly, a bit more care is generally required. In particular, in order to ensure that the supposedly more realistic difference equation models do not predict negative population sizes, it will often be necessary to modify the forms of f(x) and g(x) for large values of x and to restrict the range of the noise process. These considerations, which are elaborated in Section 4, lead to sequences of the form x(N)

t+(llN)

-

XcN) = fN(XiN))(l/N) t

+ gN(XiN)) [I”’

for t = 0, l/N,...,

(2.6)

where for each N, 61”’ is a sequence of independent, identically distributed (i.i.d.) random variables with mean zero and variance l/N (which is the variance of [W(t + (l/N)) - W(t)]). In (2.6), the coefficients fN and g, agree with f and g of (2.1) for x less than some large value C, , and these transition points, C, , approach infinity as N---f co. Noting that the passage of a unit of time in (2.1) corresponds to N iterations of (2.6), (2.6) can be replaced by the notationally more convenient form x(N)

n+1 -

xzNLN) = fN(x;Nv))(

1 /N) + &.,(x:N,N’)(q;:;/N”2),

(2.7)

where for each N, 71I”’ is a sequence of i.i.d. random variables with mean 0 and variance 1. Our interest is in the limiting behavior of the sequence of processes{X,(~)}~=, defined in terms of the solutions to (2.7) by XN(t) = X[cj, for each t > 0. (Here [x] denotes the greatest integer 0. Roughly speaking, this implies that extinction is more likely in a more uncertain environment. By contrast, for X, extinction is certain only if v < 0 and us plays no role in determining whether the population persists. Assuming the dynamics of the populations of interest can be adequately approximated by a stochastic difference equation with autocorrelated noise, the criterion of the previous section suggests that we investigate the limiting behavior of a sequence of processes {XN}E=l obtained from models approaching (3.2) which involve differences and correlation. We set XN(t) == XiEil , where XI”‘) denotes the solution of x~~‘, - xLN) _ ()cJj,r) X(N) n $. ax(N)(p) n n+lpvly

(3.5)

x(N) y= ‘y ,) . 0

For each N, we assume that ~1.~’ is a stationary process with mean zero, variance one, and COV(~~~~,qL$) =- K&s), nonzero for at least some values of s. Note that in these models, dX/dt is approximated by N(XkTi - XAy)) and dW/dt is approximated bv N1’z~~;y). Also as indicated by the definition of XN(t), the passage of a unit of time in (3.2) corresponds to N interations of (3.5). Hence, in order that the autocorrelation in (3.5) vanish as N-+ “o, we require that for each fixed positive t

i.e., K.V(N) = 0(1/N) as N-t crj. This condition parallels the one discussed above for continuous time. In order to prove convergence results for the processes X, , it suffices to assume that the sequences q!‘“’ satisfy mixing conditions such as the discrete analog of (2.9) together with a condition on the rate at which 4(s) 4 0. If (2.9) is assumed, then we also have that for each fixed 0 ,< s < t the random variables which are the noises the Nth model experiences at times s ?;“,l] and $$j, and t, respectively, are asymptotically independent in the sense that for fixed Bore1 sets A and B P(T#$J E -4 ! 7&

E B) - P(T&

E A)\ + 0

as

N + cx3.

So, at least intuitively, if these conditions are satisfied, the models (3.5) approach (3.2) as N + a. According to (3.5) we have that for each N and n

STOCHASTIC

153

CALCULUS

thus Lvtl

Xnr(t) = x0 n

(1 + (f/N)

t (u/N”“)

TjJI’“‘).

i=l

Our task is to describe the limiting behavior of these processes as N+ co. Not surprisingly, the result depends on the structure of the processes ylN’. In particular, it depends on whether or not as Nco the covariance functions K,v approach the Kronecker delta function s lJ,m := 0

=1

if

m#O

if

m = 0,

which is the covariance function for an i.i.d. sequence with variance 1. This distinction is precisely analogous to the one made in discussing the limit theorems of Khas’minskii (1966) and Papanicolaou and Kohler (1974) for sequences of differential equations with correlated noise. In that context, if the noise processes [,v in (2.8) were generated from a single stationary process v with covariance function K, the limit, Xy, of the solutions XcN) to (2.8) corresponded to neither the Ito nor the Stratonovich solution of the SDE being investigated unless s K(x) dx = 1. However, if the fN’s were constructed from a sequence of stationary processes whose covariance functions converged to the Dirac S-function, then the Stratonovich solution was obtained in the limit. The two cases below provide difference equation counterparts to these results. Case I. Suppose that the noise processes 77I”’ in (3.5) do not change with N. Let 7. be a stationary process with mean 0, variance 1, and covariance function (K(s)}. Setting phi) = yn for each N and n, (3.6) becomes X,(t)

= X,, n

(1 + (y/N) f (u!N~/~) TV).

(3.7)

i=l

The precise assumptions on the process 7. and the details of a convergence proof for the processes X, will not be given. As discussed below, a proof can be found in Guess and Gillespie (1976). W e p resent here a heuristic development which contains the essence of a proof. Applying the logarithm to (3.7) we obtain log XN(t) - log x0 = [f’

log (1 + $

+ &$

Vi)

i=l

j&?'

=E(y&+

as

N+co,

-

;$li2

)

+o

k&4

154

MICHAEL

TURELLI

where the last term results from expanding

the logarithm.

Note that as N +

co,

(3.8a) [Ntl

(l/N)

1 7i2+t

(3.8b)

i=l

by the Ergodic

Theorem,

and [Ntl

(l/N’/2)

c

7ji +

1+ 2 f L

i=l

s-1

K(s) II2 W(t), I

(3.8~)

according to a generalized version of the Central Limit Theorem, where W is a Brownian motion process. The expression for the variance parameter in (3.8~) is crucial so we motivate it by noting that Var

(1/N1/2) %:I TV,) = 1 + (2/N) Ni’ 2 K(j i=l j=i+l

- i)

N-l =

1 +

(2/N)

=

1 +

(2/N)

+

1 + 2 i K(s) s=l

C *=l

(N - s) K(s)

N-l

Assembling

1

s=1

where

Z(s)

as

N -

Z(s) = i K(i) i=l 00

by Cesaro summation.

the pieces from Eqs. (3.8) we see that

log Xiv(t) - log x, + A - &Jat + (T 1 + 2 5 K(s) 1’2 W(t). I 8x1 Thus, K(s)

as

N-co.

(3.9) As in the case of the “y-limits” for solutions to the sequences of differential equations described in the previous section, this limit process is in general neither the Ito nor the Stratonovich solution of the SDE, (3.2), being considered. In fact because of the appearance of the covariance function K, this limit cannot in general be obtained from any of the one parameter family of stochastic integration procedures described at the end of the previous section. However, the conditions for extinction, r < &s2, and exponential growth, F > $a”, of this limiting process agree with those of Xi , the Ito solution of (3.2).

STOCHASTIC CALCULUS

155

After constructing a direct proof of (3.9) (see Turelli, 1977, Appendix II), we became aware of independent work by Gillespie and Guess concerning limits for solutions to sequences of stochastic difference equations involving autocorrelated noise (Guess and Gillespie, 1976; Gillespie and Guess, 1976). Their motivation was the investigation of particular genetic models rather than the problem of interpreting heuristically derived SDEs. Nevertheless, in the first of their papers, they have rigorously proved limit theorems for the solutions to a class of discrete time stochastic models which essentially includes (3.5). In their second paper, they heuristically extend their result to a wider class of models which includes diploid random selection models. Their procedure and its relationship to the application of SDEs is discussed in Section 6. Case II. Here we assume that 7.W) has covariance function KN and that as N---f co, KN approaches the Kronecker delta function. Proceeding heuristically as in Case I, the limits corresponding to (3.8a) and (3.8b) are unchanged. However, retracing the calculation following (3.8) we expect that [Ntl

(l/N”‘)

c

~i(~)--t W(t).

i=l

(An appropriate assumption is that each 71”” satisfies a uniform mixing condition and that the same +-function applies for every N; see Turelli, 1977, Appendix III.) Substituting this limit for (3.8), we see that in this case the sequence of solutions converges to precisely Xr . From these results we see that although the processes X, may not converge to either X1 or Xs , the limit processes share with Xr the same conditions for extinction, i.e., [X(t) -+ 0 as t -+ co], and explosion, i.e., [X(t) + 0~)as t --f co]. This parallels the situation for the sequences of differential equation models discussed in Section 2. There, although the limit processes did not always agree completely with the Stratonovich solution, it was asserted (and is demonstrated in the next section) that they had the same extinction and explosion behavior. Since in many contexts, most notably the theory of niche overlap, the biological question of interest turns on whether a species persists, it may be useful to know that the extinction behavior of the solutions to a broad class of stochastic difference equations of the form (2.7) involving correlated noise converges to that of the Ito solution for (2.1). This issue together with an examination of conditions under which various diffusion process growth models predict that population explosion might occur is taken up in the next section.

4. ASYMPTOTIC BEHAVIOR OF GROWTH MODELS The asymptotic behavior of the Ito and Stratonovich solutions (denoted X, and Xs) of the SDE (2.1), as well as the diffusion limits, Xy, obtained in Sec-

156

MICHAEL

TURELLI

tion 2, is investigated and compared to the asymptotic behavior of solutions (denoted X,v) to models of the form (2.7) which involve differences and correlated noise. The criterion of Section 2 suggests comparing the asymptotic behavior of X1 and X, to that of X,,, , for large N, to determine which, if either, is likely to provide a useful approximation. In the previous section, it was shown that for the exponential growth model the asymptotic behavior of the limits of such sequences agreed with that of X, . Results are given below which generalize this observation to a large variety of SDEs which might arise by randomizing the parameters of standard ecological models. These results center on probabilities of population extinction and explosion. Let {X(t)}tao denote a stochastic process modeling population size. As in the discussion of exponential growth, we associate the event [X(t) ---f 0 as t ---f so] with extinction. Similarly, the event [X(t) ---f CO as t - 001 is referred to as a population explosion. (This last definition differs from the mathematical convention (cf. McKean, 1969, p. 50) w h’ic h associates explosion with [X(t) == CC for some t < co]; but as we see shortly, in many cases of interest both events have the same probability.) To discuss the probabilities assigned to these events by the diffusions Xi , X, and Xy, it is useful to introduce some terminology. Let X be a diffusion process with state space E, an interval of the real line with endpoints (possibly infinite) ri < Y$ . The boundary point yi is said to be attracting if P[X(t) + ri as t + CO] > 0; otherwise it is called repelling. An attracting boundary ri is attainable if P[X(t) = ri for some t < co] > 0; otherwise it is unuttuinubZe. Simple integral criteria involving the infinitesimal mean and variance of X are available to determine the nature of a boundary (see Prohorov and Rozanov, 1969, pp. 265-270). For growth models of the form (2.1) d erived from deterministic equations, the state space will typically be (0, r) where Y is possibly $ co. The boundary classification machinery will be employed to determine conditions on the parameters of the original deterministic model and of the noise process which ensure that population extinction and/or explosion (as defined above) are certain or are impossible. We begin with an example which illustrates a disturbing feature of the Stratonovich solutions to SDEs obtainable from density-dependent growth equations. Following Feldman and Roughgarden (1975) we consider the SDE dX = YX(I - (X/K,))

dt + arX2 dW,

(4.1)

which is formally derived from the logistic equation by replacing l/k’ with (l/K,,) + a(dW/dt). As Feldman and Roughgarden point out, both the Ito and Stratonovich solutions have the property that extinction is impossible (i.e., 0 is repelling) as long as Y > 0. Let M, and M, denote the infinitesimal means of the Ito and Stratonovich solutions of (4. I). A ccording to formulas (2.3) (2.4) and (2.5),

STOCHASTIC

CALCULUS

157

MI(X) = 41 - (xiql)), M,(x)

= rx(1 -

(X/K,))

+ O”raXa,

and for both solutions the infinitesimal variance is V(x) = 02y2x4. Applying the classification procedure to the boundary point +co, it is easily verified to be a repelling boundary under the Ito calculus. Hence according to this model, the population size never converges to co, as we would expect from an equation involving density dependence. However, when the behavior of the Stratonovich solution is checked, + co is not only an attracting boundary, it is also attainable! So if Stratonovich calculus is employed, model (4.1) predicts that, with probability I, the population size becomes infinite in a finite time as long as Y > 0. The reason is straightforward. Whereas, M, represents logistic growth pushing the population size toward the deterministic equilibrium K, , Ms consists of this quadratic density-regulation term plus the positive cubic term u2r*x3. For large X, this cubic dominates and provides a strong deterministic push toward + 00 which makes it attainable. In spite of this, Feldman and Roughgarden display a function which they call the stationary probability density for the Stratonovich solution of (4.1) (see their Eq. (24)). Indeed the function they display, call it p, is a stationary solution of the Kolmogorov forward equation associated with Ms and V and it satisfies .$‘p(x)dx = 1. H owever, since 00 is an attainable boundary (in fact, it is a regular boundary according to Feller’s 1952 classification scheme), the stationary distribution given implies a particular sort of behavior when GO is reached. The implicit assumption is that it is a reflecting boundary so that once the process reaches co it “bounces off” in a continuous fashion. If instead, 03 was assumed to be an absorbing boundary, the “stationary distribution” would consist of unit probability mass concentrated at the point +co. The biological interpretation of the reflecting process is tenuous at best, and the bizarreness of its behavior is reflected in the fact that its stationary distribution has an infinite mean. We now consider a class of SDEs, obtainable from density-dependent growth models, whose Ito and Stratonovich solutions display the same explosion behavior as seen in the example above. Since boundary behavior is a local phenomenon, classifying the boundary at co requires specifying the infinitesimal mean and variance only for large values of x. Hence, it suffices to describe the coefficients f and g of (1 .l) only for large values of X. For the sake of simplicity, we consider models of form (1 .I) satisfying the conditions

g”(x) > 0

for

E > 0,

f(x) N cx* and g(z) N bxA as x - co, where c, b and h are constants satisfying c < 0, A > 1,

(4.2a)

(4.2b)

158

MICHAEL

TURELLI

and f and g are linear combinations

of (possibly

noninteger)

powers of x.

(4.2~)

These conditions restrict the form of the original deterministic model as well as the way randomness is incorporated into it. Condition (4.2a) ensures that irrespective of which calculus is used the solution process can attain arbitrarily large values so that the asymptotic expressions in (4.2b) take effect. Conditions (4.2b) and (4.2~) allow us to perform the boundary classification simply. The examples below illustrate that these conditions are satisfied by a large variety of SDEs that can be obtained from density-dependent growth equations by incorporating randomness only into the “regulation” term. Ayala et al. (1973) have argued for the utility of competition equations based on the single-species growth models dX/dt

= YX(I -

(X/K)“),

dX/dt

= rX[l

(X/K)(l

where

0 > 0,

(4.3)

and -

+ /IX)].

(4.4)

Both of these are generalizations of the standard logistic equation. If we assume that environmental unpredictability is manifested only in variation of the carrying capacity K, and model this by substituting (l/K$ + u(dW/dt) for l/K0 in (4.3) and (l/KJ + a(dW/dt) for 1/K in (4.4), the resulting SDEs dX = rX(1 -

(X/KJe)

dt + c~rXl+~ dW

(4.3a)

dX = rX[l

(X/K,,)(l

+ /3X)] dt - arX2(1 + flX) dW

(4.4a)

and -

satisfy (4.2) as long as we assume /3 > 0 in (4.4). The result of interest is that cc is a repelling boundary for the Ito solution to any SDE satisfying conditions (4.2), whereas co is both attracting and attainable for the Stratonovich solution as well as the diffusion limit Xy specified by (2.10) and (2.11). Recall that the infinitesimal mean, My, of the process Xy is f(x) + hg(x) g’(x) and that M, corresponds to the case y = 1. As in the logistic example, the mathematical origin of the “explosions” of Xs and Xy is clear. From (4.2b) and (4.2~) we see that M,(x) N cxA,

where

c < 0,

as

.r -

cc,

whereas My(x)

N cxA + ~yA62x2A-1 = &yhb2~“A-1 + o(x~~-~)

as

x-

r*).

STOCHASTIC

159

CALCULUS

Hence when x is large, the Ito diffusion process experiences a systematic drift toward the origin, whereas Xs and Xy are subject to even stronger drift out toward co. This coupled with the fact that the infinitesimal variance of X does not grow too much more quickly than My as x -+ co makes co attainable. In the present context, it hardly seems necessary to employ the criterion of Section 2 to decide which interpretation of (2.1) is more reasonable. Nevertheless for the sake of completeness, we discuss the behavior of the solutions, denoted X, , to sequences of stochastic difference equations of the form (2.7) with correlated noise which are derived from an SDE whose coefficients f and g satisfy (4.2). First observe that it is trivially true that co is unattainable for all of these processes since in a finite time each undergoes only a finite number of finite jumps. In order to evaluate P[X,(t) + co as t - co], it is necessary to digress and discuss how the coefficientsfN and g,v of (2.7) are determined. Since this specification is also essential to the discussion of extinction which follows, the coefficients f and g of the original model are not restricted b.y (4.2). We assume that f(x) -

ax

as

x-+0;

f(x)

L-xAl

as

x-+

g(x) -

bx’

as

g(x) -

dxA”

as

-

co,

where

c < 0,

A, > 1

x--,0,

where

y > 1;

x--too,

where

1 < A, < Ai .

(4.5)

These conditions are met by all the density-dependent growth models that we consider. To see the necessity of modifying f and g when forming stochastic difference equations from the SDE (2.1), let us consider how models (2.7) would behave if fN and g, were equal to f and g for all N. The basic difficulty is that since the solutions move by jumps rather than continuously, the fact that the coefficients f and g are zero at the origin does not prevent these processes from becoming negative. This point was overlooked in Levikson’s (1976) discussion of discrete versions of (4.3). The difficulty arises from two sources: (1) the behavior of the deterministic portion of the models x(N)

n+1 - XLN’ = (l/N)f(XLN’),

and (2) the action of the noise g(X~~))(~~~)l/N1/z). The first point is illustrated by the following sequence of discrete analogs of the logistic X t+(l/,v) - Xt = (GO

X4

-

&G/W).

Note that if the initial population size is very large, in particular X0 > K(1 + (N/r)), this model predicts that the population size will be negative

160

MICHAEL

TURELLI

in the next generation, then will become increasingly negative with a quadratic growth (shrinkage ?) rate (assuming Y > 0). Although this model is nonsensical for such large starting values, it will be well behaved if begun at smaller values, i.e., X0 < K(l + (N/Y)), if N is so large that (r/N) < 2 (see May and Oster, 1976). Thus, these deterministic models will behave nicely as long as the initial values are suitably restricted; and the permissible upper bound goes to infinity as AJ+ co. However, once stochasticity is introduced, say by replacing l/K by (l/K,) + ON*$(‘V)t 1 there is nothing to prevent the process from wandering out to the large values where the deterministic model is no longer sensible or from jumping directly to a negative value due to the action of the noise. Since the purpose of introducing stochastic difference equations is that they should provide a more realistic model than that ambiguous SDE, it is pointless to consider such models whose behavior is biologically nonsensical. A particularly straightforward way to modify the models so that they will behave reasonably (i.e., remain positive) is to set

define g, in the same way, and restrict ~1~~’ so that / ?tN’ 1 < B, . A little algebra shows that under conditions (4.3, these truncation constants, BN and C, , can be chosen so that they go off to infinity as N -+ co, allowing the discrete models to approach the continuous one being investigated. (In particular, B,v and CN can be chosen to be O(N1i4).) it becomes very simple to verify that If this procedure is followed, P[X,(t) -+ co as t --f co] = 0 for the discrete processes obtained from an SDE satisfying (4.2). It follows from observing that once X, exceeds the level C, , its governing equation becomes x(N) nt1

_ X(N) ?I

= (l/N)f(C

N

) +- g(C N )(7+“‘/P), n+1

where f(CN) is a negative constant. So if X,v were to stay above C, , it would be a sum of identically distributed random variables with negative expectation. Thus, if the qtN” s satisfy any sort of mixing condition, such as ergodicity, which ensures that they are not completely dependent, X, could not tend to + co. The fact that F’[X,(t) ---f m as t + 031 = 0 is not tied to the particular scheme we have chosen to specify fN and gN . Indeed, if these coefficients were not leveled off entirely but instead allowed to continue growing (with f becoming more negative) as x --f co, the processes XN would experience an even stronger drift back toward the origin once they reached large values. After this unfortunately long winded but necessary digression, the immediate message is simply that the “explosion” behavior of these discrete time processes agrees with that of Xr rather than Xs .

STOCHASTIC

161

CALCULUS

Now we turn our attention to probabilities of extinction and consider a class of SDEs which is complementary to the one specified by conditions (4.2). As discussed previously, conditions (4.2) are likely to apply to SDEs obtained by introducing stochasticity only into the “regulation” term of a density-dependent growth equation. In contrast, the results to be presented on extinction depend on conditions (see (4.6) below) which are likely to be satisfied by SDEs obtained by introducing stochasticity in the “intrinsic growth rate” parameter irrespective of whether the noise affects regulation parameters. The conditions we impose are: f(x)

-

ax and g(x) -

bx as x + 0, where a and b are constants,

(4.6a)

and f and g are linear combinations

of (possibly noninteger)

powers of x.

(4.6b)

Applying the criterion for 0 to be a repelling boundary, it is easily verified that according to the Ito calculus extinction never occurs (i.e., P[X(t) -+ 0 as t -+ co] = 0) if a - $b” > 0; whereas the Stratonovich calculus as well as the “Y-limits” of Section 2 give the condition a > 0. These conditions are precisely analogous to those given for the exponential growth model (3.2) with r replaced by a and u replaced by b. The similarity is hardly surprising since condition (4.6a) says that near the origin the model we are examining looks like the exponential growth model. In the present context, unlike the discussion of population explosion, it is not obvious which of the calculi gives a more reasonable prediction. However, if it is agreed that a stochastic difference equation with correlated noise would provide a suitable model for the growth process of interest, the criterion of Section 2 can be applied. Consider a sequence of models of form (2.7) which have been derived from an SDE satisfying (4.6). We assume that fN, gN and the range of the noise processes 7.ok) have been chosen, as specified above, so that the solution processes X, remain positive. For each A’, 7 1”’ is assumed to be a stationary, ergodic process with mean 0 and variance 1 and to satisfy 1 + (u/N) + (b/W/?) T(N) > 6, where 0 < 6 < 1. With these assumptions, it is easy to mimic the proof of Theorem 1 in Karlin and Lieberman (1975) to show that P[X,(t) ---f 0 as t-+co] =Oif E[N ln(l + (u/iv) + (b/Nl”)

#‘)I

> 0.

(See Turelli, 1977, Appendix IV, for details.) Letting N-t co so that the sequence of difference equations approaches the SDE, this condition becomes a -

&b2 > 0

162

MICHAEL

TUFtELLI

which agrees with the Ito condition for the origin to be repelling. Conversely, by following the proof of Karlin and Lieberman’s Theorem 2, it can be shown that for every E > 0 there exists a value x,v such that if Xh”’ < x, , P[X,(t)-+Oast+co]>l--if E[N ln(1 + (a/N)

+ (b/P”)

vjN’)] < 0.

We further conjecture that since the process X, cannot escape to infinity, it cannot stay out of any neighborhood of origin indefinitely. This would imply that if the condition above is met, P[X,(t) --, 0 as t - co] = 1 irrespective of the starting value Xi’“‘. Letting N ---f CO, this condition agrees with the Ito condition for 0 to be attracting. Hence, as in the discussion of probabilities of explosion, the behavior of the Ito solution of the SDE model is consistent with the behavior of solutions to biologically more realistic equations whereas the behavior of the Stratonovich solution is not. As stated earlier, the class of SDEs specified by (4.6) is complementary to of that specified by (4.2). I n many contexts the biological interdependence growth and regulation parameters makes it unrealistic to assume that environmental variation affects only the regulation parameter and leaves the growth parameter fixed. As an example, consider a population undergoing logistic growth in an unpredictable environment. However, instead of the standard form of the logistic in which r and K appear as seemingly independent parameters, consider Schoener’s (1973) parametrization

dX/dt = RX[e(T - AX) - c - yX]. Here, R is a conversion factor of energy into offspring, T is the proportion of time spent in feeding and intraspecific interaction, E is the net energy obtainable per unit foraging time, and fl, c, and y measure costs of metabolic maintenance and intraspecific interaction. (Here fl does not denote white noise.) In terms of these parameters, the standard r and K are given by

r = R(eT - c)

and

K = (ET - c)/(y + c/l).

(4.7)

Following the standard heuristic procedure, we might model the effects of fluctuating resource abundance by replacing t with E,,+ o(dW/dt). Doing this, we get the SDE

dX = r,X(l

- (X/K,))

dt -t aRX(T - AX) dW,

where Y,, and K, are given by (4.7) with E replaced by E,, . Note that at X = T/A, the noise coefficient vanishes and the deterministic motion is downward (since T/A > K,,). Hence if the process begins inside the interval [0, T/A], it will stay

163

STOCHASTIC CALCULUS

there. Thus the results on explosion do not apply. Clearly, however, the model does satisfy (4.6). More generally, if we begin with a growth model of the form

dXjdt = X a, + f aiXY” ,

where

yi > 0

for each

i = 1, 2 ,..., n,

i=l

and introduce stochasticity into a rate parameter so that a, varies, the resulting SDE will satisfy (4.6). A s indicated in Section 6 below, such models may be more biologically meaningful than those derived by introducing white noise directly into the carrying capacity. In any case, for a very large variety of SDEs likely to arise by introducing stochasticity into a deterministic growth model one or the other of the results in this section will apply. Both imply that the Ito solution is likely to be more useful than the Stratonovich one contrary to various statements in the recent literature.

5. WHAT

HAPPENED

TO STRATONOVICH?

In recent papers, Capocelli and Ricciardi (1974) and Tuckwell (1974, 1976) have applied Stratonovich calculus to SDEs heuristically derived from standard ecological and genetic models. Their justifications for this choice have implicitly (e.g., Tuckwell’s references to Gray and Caughey (1965) and Wong and Zakai (1965)) or explicitly (Capocelli and Ricciardi, 1974, discussion on p. 32) involved the assumption that the “true model,” for which the SDE is serving as a convenient approximation, is a differential equation subject to autocorrelated noise. Under this assumption, the appropriateness of the Stratonovich calculus would follow from the various limit theorems for differential equation models cited in Section 2. However, we would assert that most biological populations fall, more or less, into one of the three categories discussed in Section 2 for which stochastic difference equations are more realistic models. If this is indeed the case, the results of the two preceding sections suggest that the use of the Ito calculus may be appropriate. Contrary to this, Feldman and Roughgarden (1975) have argued, using the SDE exponential growth model (3.2) as an example, that the appropriateness of the Stratonovich calculus is supported by considering “... the continuous form [i.e., (3.2)] . . . as a limiting representation of [a] discrete process.” In the hope of clarifying at least a small part of what has become an exceedingly confusing area, we discuss this discrepancy in detail. In their discussion, Feldman and Roughgarden consider the model

dN=pNdt+uNdW, They point out that according 653/12/2-4

N(0) = No.

to the Stratonovich

calculus In(N(t)/N,,)

(5.1) has a

164

MICHAEL

TURELLI

normal distribution with mean pt and variance a’t while the Ito calculus leads to ln(N(t)/Na) having mean (p - $9) t. These assertions follow directly from (3.3) and (3.4). They then noted that I‘... the exact discrete time analysis of the process K,+, = Z,N, with E[ln(Z,)] = X and Var[ln(Z,)] = y” gives ln[NJN,,] normal with mean At and variance yPf. In this simple case . . . the Stratonovich calculus seems to be physically appropriate.” To see what is going on here, let us forget stochasticity for a moment and compare the differential equation dN/dt = pN with the recurrence relation N,,~, = IN, In order for their solutions to coincide, we require 1 en or equivalently In(Z) = p, Now introduce white noise variation into the differential equation. Feldman and Roughgarden imply that to get a “corresponding” stochastic difference equation we should require E[ln(Z,)] However,

= p.

the assumed analogy is arbitrary.

(5.2)

\Vhp not instead require

E(Z,) == ep ?

(5.3)

Indeed, if we assume that the sequence (2,) is composed of independent lognormally distributed random variables with mean eo and variance ua, we find that In(n;,/lv,) is normal with mean (p - 47”) t, where 72 == ln(1 + u2e-*~). Note that if u2 and p are small, 72 N u2 so that ln(N,/N,,) is approximately normal with mean (p - &r”) t and variance u2t. From this we would conclude that the Ito “... calculus seems to be physically appropriate.” This dilemma arises at least in part from the ambiguity of the phrase “mean growth rate.” If we agree to call the ratio N,+,/N, the growth rate of the population, Feldman and Roughgarden’s specification (5.2) says that the population has a geometric mean growth rate of ep. (The geometric mean of a positive random variable X is defined as eEdnX).) On the other hand, our condition (5.3) says that the arithmetic mean growth rate of the population is en. Since the geometric mean ln(N,/N,) = cf=, ln(ZJ, it is E[ln(Z,)] or equivalently growth rate that determines the fate of the population. By the law of large numbers, the population will grow beyond any bound if E[ln(Z,)] > 0, and will go extinct if E[ln(Z,)] < 0. Under assumption (5.2), extinction occurs if p < 0. But when (5.3) is assumed, the approximation E[ln(Z,)] = p -

$T~,

(5.4)

which is valid if p and ~9 are small, implies that extinction is sure if p < $0’. In this context the origin of the difference between the conclusions about the effects of environmental variation on the likelihood of extinction is clear. It is not a subtle mathematical point involving white noise but simply a matter of

STOCHASTIC

CALCULUS

165

keeping straight the interpretation of the parameters. If p is defined as the log conclusion that of the geometric mean growth rate, then the “Stratonovich” extinction requires p < 0 is appropriate. However, if p is defined as the log of the arithmetic mean growth rate, the “Ito” conclusion that extinction occurs if p < &a is appropriate. It is important to recognize that these results are not contradictory; in fact as shown by (5.4) they are (at least approximately) the same. The approximation which enters into this discussion would vanish if instead of comparing the behavior of the continuous time model (5.1) to that of a single difference equation, we considered a sequence of stochastic difference equations converging to (5.1). Although these manipulations with arithmetic and geometric means indicate that both the Ito and Stratonovich extinction criteria can be “justified” by considering discrete time models, the paradox presented is transparent. If we follow the standard convention and interpret “mean” as mathematical expectation, i.e., arithmetic mean, then the Ito criterion follows, in concordance with the discussions in Sections 2-4 of convergence theorems for discrete time models. This conclusion also agrees with the results of Keiding (1975), who explicitly calculated a diffusion approximation for a discrete time model of unrestricted population growth.

6.

DISWSSIOX

To avoid leaving the impression that the basic conclusion of this paper is that the use of heuristically derived SDE models by population biologists is proper as long as the Ito calculus is employed, we point out some obvious limitations on the application of these models as well as weaknesses in the criterion we have proposed for determining an appropriate solution. An intuitive restriction on the heuristic derivation of SD& is that the parameter which is randomized be capable of assuming any value between +30 and -co. The rationale behind this restriction is that after replacing the parameter y of (1 .I) by y(t) = 7 + d(t) th e new stochastic parameter will take on all of these values since (1 does. This can be glibly justified by noting that (1 is a generalized Gaussian process hence its range is the entire real line. An alternative argument which is constructive and more in line with the theme of this paper involves regarding the SDE (1.3) as a limit of stochastic difference equations. To keep track of the behavior of the randomized parameter, rewrite (1.3) as

dX = (fi + j&2) dt + uJ2dW. As in Section 2, we introduce x(N)

fl(l’N)

ON) -

xt

the sequence of difference equation approximations

166

MICHAEL

TURELLI

where the random variables Q(“) all have mean 0 and variance 1. This equation can be rewritten as x(N)

t+(l/N)

-

y(N)

1 t

=

&

[fl(XjN’)

+

(?/ -

ON1’27)r(N))f2(Xt(N))],

which displays more obviously the way in which the randomized parameter is assumed to behave. Since the qr’ ‘s have mean zero and variance one, it follows that as N---f CC the range of the random variable 7 $ ,N1+~N) must become unbounded above and below. This observation raises questions about the meaning of SDE models which are derived by randomizing a parameter whose range is intrinsically restricted. The most obvious ecological examples are the various randomized logistic equations in which the carrying capacity, K, or some function of it is replaced by a mean value plus white noise. The resulting stochastic equations have negative carrying capacities at times. This turns the model into one with undamped quadratic growth rather than linear growth with quadratic damping. The effect of this can be seen in the Ito solution of Feldman and Roughgarden’s (1975) model dx

= rX(1 -

(X/K,))

dt - a(r/K,)

X2 dW.

As noted in Section 3, because of the form of the infinitesimal mean associated with Stratonovich integration, the Stratonovich solution of this equation reaches, infinity in a finite time. This behavior is clearly nonsensical. Although the Ito solution does not explode, its stationary distribution is p,(x) = F

exp [$;

(1 - +$)I,

which possesses only two finite moments (i.e., E(Xk) = cc for k > 3). This reflects the presence of a lot of probability associated with large population sizes and can be attributed to the strong push out to infinity the model experiences when K(t) becomes negative. Because of this property, the usefulness of this model for discussing the effects of environmental variation in resource abundance is rather uncertain. Hence, it seems that one reasonable restriction to impose on the derivation of SDEs from deterministic models is that only parameters of potentially unbounded range should be randomized. For instance, there is no logical difficulty with the SDE (3.2) for exponential growth since negative growth rates are quite acceptable. Similarly, the randomized logistic (4.8), derived by varying one of the rate parameters introduced by Schoener (1973), does not impose unrealistic values on its parameters. Moreover, since the solution of this SDE lives on a finite interval all of its moments will be finite. Although

STOCHASTIC

CALCULUS

167

Schoener’s parameterization is based on some rather strong assumptions, the corresponding SDE may provide a more reasonable model for evaluating the effects of fluctuating resource abundance than the models discussed by May (1973), and Feldman and Roughgarden (1975). The properties of this model as well as a difference equation version of it will be discussed in a later publication on the problem of limiting similarity in a stochastic environment. As the “stochastic K” example illustrates, a danger associated with the casual derivation of SDEs is that it conceals the assumptions inherent in the final model. As we have stressed, it is helpful to view these idealized models as limits of more realistic models. Very often the more realistic models will be stochastic difference equations. Yet, the randomization procedure discussed in the Introduction begins with a differential equation and the parameter of interest is assumed to enter linearly. If the true model is a stochastic difference equation and the varying parameter enters nonlinearly, application of the SDE formalism can lead to a good deal of confusion. Tuckwell’s (1976) study of the effects of stochastically varying selection coefficients provides an example. Tuckwell begins with continuous time models to avoid “the necessity of carrying out limiting operations to generate them.” Yet he uses the analytical results he obtains to analyze data on an organism, Punaxiu dominula, with nonoverlapping (annual) generations (cf. Ford, 1965, p. 108). Also, he contrasts his findings with those of Kimura (1954), w h o considered a diffusion approximation of a discrete generation model. The model which is central to Tuckwell’s investigation involves a single diallelic locus with additive fitnesses. The differential equation he employs is dX/dt = $X(1 - X).

U-5.1)

The parameter s is described as “the measure of the difference in the fitnesses of the two homozygous classes.” To explore the situation in which the selection differential fluctuates randomly but no genotype is favored on the average, Tuckwell replaces the constant parameter s of (6.1) by a(dW’/dt) and interprets the resulting SDE in the Stratonovich sense. His principal qualitative conclusion is that irrespective of the initial frequency, each allele has probability 0.5 of being quasifixed (i.e., arbitrarily near 1 in frequency). For the case where one allele has a mean selective advantage, he concludes that it will become quasifixed irrespective of its initial frequency. A good part of the difficulty involved in evaluating the accuracy of Tuckwell’s results arises from his ambiguous use of the term “fitness.” In particular, it is not at all clear whether he intends the parameter s to be a measure of the difference between the Malthusian parameters of the homozygotes or the difference between their “Wrightian” fitnesses (cf. Crow and Kimura, 1970, Sects. 1.2 and 5.3). The fact that he uses a differential equation formulation of the selection process suggests that s be considered as a selection differential on the Malthusian

168

MICHAEL

TURELLI

scale. However, in applying his result he estimates the mean selection coefficient for the medionigra allele in Panaxia by equating the expected change per generation of the allele to the observed change. Also, he compares his findings to those of Kimura (1954), who explicitly states that his model assumes discrete generations, no dominance, and that s is the difference in the Wrightian fitnesses of the homozygotes. These applications suggest that Tuckwell intends his SDE model to serve as a diffusion approximation of a discrete time model, with additive (Wrightian) fitnesses, and fitness differential s between the homozygotes. Considering s initially as a constant, we see that there is a one-parameter family of such models. Denoting the genotypes AA, Aa, and aa, note that the fitness triples (1 + (c + 4) s, 1 + cs, 1 + (c - +) s) satisfy the specified conditions for each value of c. Assuming random mating and infinite (i.e., very large) population size, the frequency of allele A evolves according to AX = gql

-- X)/(1 $ s(c - $ A--X)).

(6.2)

If the fixed parameter s is replaced by a sequence of stationary, stochastic fitness values s, , the resulting nonlinear stochastic difference equation is of a type that has been much studied recently. In particular, it is well known that in the absence of dominance the outcome of selection (i.e., the probability of quasifixation for each allele) is determined by the geometric mean fitnesses of the genotypes (cf. Gillespie, 1973; Karlin and Lieberman, 1975). The situation discussed by Kimura (1954) and referred to by Tuckwell requires E(s,) == 0. For convenience we assume that the approximation E[ln(l + Ks~)] m &!?(s,) - $KZE(s,;3)

(6.3)

is valid. From these assumptions, it follows that the outcome of the stochastic selection process varies according to the value of the parameter c as follows: if-$ 3

as

N-+

cc.

170

MICHAEL

TURELLI

Under appropriate assumptions on the function f, it can be shown that the processes X!“) converge to a diffusion with infinitesimal mean, M(x) = i(d/?s)f(x, 0) + &T*(a2/ai~)f(x, O), and infinitesimal variance, V(X) = u”[(Z/&)f(x, O)]“. Karlin and Levikson (1975), and Levikson and Karlin (1975) discuss such approximations for stochastic selection models. The diffusions obtained in this way from (6.2) and (6.4) have precisely the same qualitative behavior (with respect to quasifixation) as the discrete processes being approximated. It is to be expected that they would model the processes more accurately than Tuckwell’s procedure since they incorporate the form of the nonlinear dependence on the selection parameter. Gillespie and Guess (1976) have heuristically derived a diffusion limit for the sequence of stochastic difference equations (6.5) which permits the ,iN” s to be autocorrelated. Their procedure, like the one above, assumesf(x, 0) == 0. On the other hand, the SDE procedure allows .f(x, 0) #: 0 but requires that f be a linear function of s. When both techniques are applicable, they lead to the same diffusion and hence are complementary. In his 1974 monograph, Ludwig proposed a heuristic procedure for obtaining diffusion approximations for differential equations of the form

dX/dt = F(X, R) (see Ludwig, 1974, Sect. IV.l) in which the parameter R is subject to “small” random perturbations. The point we are trying to make by mentioning these various approximations is that the SDE procedure is just one of several ways of deriving diffusion approximations for processes generated by varying the parameters of deterministic models. In many contexts, it will be necessary to develop the relevant approximation de novo since the standard procedures are limited in scope. In particular, to study the effects of autocorrelation on the transient behavior or stationary distribution of stochastic difference equations in which the randomized parameter enters linearly, it would be useful to establish convergence results for sequences of processes such as (2.7) with correlated noise. As illustrated by the exponential growth example in Section 3, the resulting diffusions would explicitly incorporate autocorrelation parameters. Hence in general they would not correspond to any solution of an SDE heuristically derived via white noise randomization. These approximations would provide general discrete time counterparts to the y-limits of Section 2. The criterion proposed in Section 2 for choosing an interpretation of an SDE emphasizes that the SDE is to be viewed as providing a diffusion approximation of a more realistic model. To apply the criterion, one has merely to specify the form of the more realistic model (i.e., whether it involves differences or differentials, white or correlated noise) then appeal to the pertinent general limit theorem of Section 2 or to one of the results in Section 3 or 4 to decide which interpretation is likely to be most useful. The difficulty, however, is

STOCHASTIC

CALCULUS

171

that the “true model” is never really specified and so its relationship to the SDE is rather ambiguous. In this regard, we should make clear that the criterion of Section 2 is a very lax one. For instance, it would “permit” approximating a difference equation model by a differential equation since it is known that the sequence of solutions to the difference equations x(N) tt(1lN

- xp

= (1/N)f(X:“‘)

converges to the solution of the differential equation dX/dt =f(X) (cf. Birkoff and Rota, 1969, Sect. 7.2). Yet as May has recently emphasized (May, 1975; May and Oster, 1976) difference analogs of standard differential equation models can exhibit surprisingly complex behavior which is absent from their differential equation counterparts. Hence, we certainly cannot expect the diffusion approximations arising from an SDE to give a completely accurate description of the behavior of the “true model” even when it is clear which stochastic calculus is applicable. Yet it is reasonable to suppose (until proved otherwise!) that when carefully applied, formally derived SDEs can lead to qualitatively correct conclusions concerning the effects of randomly varying parameters. Monte Carlo tests of this conjecture are currently under way. Last, we point out an obvious flaw in the use of the term “extinction” in this paper. It has been used as a synonym for the event [X(t) -+ 0 as t -+ co]. Yet it is clearly unrealistic to assert that a population whose dynamics are approximated by a diffusion process X will not go extinct if P[X(t)-+O

as t+

co] = 0

(f-5.6)

(cf. Capocelli and Ricciardi, 1974; Ludwig, 1975). The problem is that population sizes take on only integer values whereas the diffusion models are only restricted to being nonnegative. This indicates that an approximation of the form P[population size equals rz.at time t] = P[n - 4 < X(t) < n + &] for n = 1, 2,..., and P[population size is zero at t] = P[X(t) < &] is implicit in these continuous state space models. In light of this, it would seem reasonable to say that the model predicts extinction for the population when X(t) falls below 4. Suppose the event [X(t) < $1 or, more generally, any event of the form [X(t) < c], where c > 0 is fixed, is equated with extinction. It is easy to check that the Ito solution of any growth model incorporating density dependence predicts that extinction is certain in a finite time irrespective of the values of its parameters. Similarly, if the Stratonovich solution does not explode, it too would yield extinction in a finite time. In spite of this observation, we believe that the “extinction” result of Section 4 is valuable. Although it is stated in terms of condition (6.6), this result indicates, at least heuristically, how the solutions of models involving differences and correlation will behave near the origin. In this regard, it is important to recognize

172

MICHAEL

TURELLI

that when both the Ito and Stratonovich solutions of an SDE predict that absorption at a particular boundary will occur in a finite time, their mean absorption times may be quite different. The reason is that the associated diffusion processes move at different rates. In particular, the diffusion corresponding to the Stratonovich solution has a larger infinitesimal mean near the origin, so the expected absorption time for this process will be greater than that of the Ito solution. Ludwig (1975) h as suggested that the mean time until absorption at a boundary associated with “extinction” is a useful measure of the effect of random perturbations on a dynamical system. He presents a general technique for approximating this value when the underlying model is a vector valued stochastic differential equation with “small” noise terms. His figure 6.1.2 shows that this mean persistence time can vary appreciably for a given model depending on which integration procedure is used. To the extent that such quantitative differences separate the behavior of the Ito and Stratonovich solutions of a particular equation, results such as ours which bear on the appropriateness of these solutions are useful.

APPENDIX

This appendix is devoted to stating and proving a convergence result for a sequence of stochastic processes. The mathematical setting, notation, and techniques come primarily from Billingsley (1968). The continuous time processes whose limiting behavior we want to investigate are obtained from discrete time processes by holding the value of the discrete process fixed between sampling periods. Because of this, the realizations of the processes all live in the space D[O, 00) of functions f: [O, co) --t R which are right continuous and have left-hand limits. For each fixed N, the process X,(.) is considered as a random element of this function space. Our interest is centered on the probability measure, PN , that it induces on a Bore1 field, denoted 9, generated by a suitable metric on D[O, co) (cf. L o&e, 1973, Chap. 1, Sect. 2). Since the limit processes we obtain all have continuous sample paths and the jumps of the processes X,(.) occur at fixed times, it suffices to take for @ the Bore1 sets of the topology generated by a metric, say d, for uniform convergence on finite time intervals (cf. Billingsley, 1968, Sect. 18; Whitt, 1970). Hence if (xn} and x are elements of D[O, co), d(x, , x) - 0 as rz --f co if and only if for each finite --f 0 as n -+ co. (If the discontinuity points of T > 0, SUP, 0 be fixed. Iff andg were bounded, from Kushner (1974, Theorem 1). we consider the processes X,,,, , X,,, respectively, with f replaced by

x 0 is fixed)

:=- smooth (Cl) transition

for

M 0 be given. We show that for all N sufficiently

large

P,(F) - P,(F) < E. From (i) and (ii), we see that 3Mo 3

From (iii), 3Nr 3 VN > Nr ,

hwOW

- &,(F)

< c/3.

From (iv), 3Na 3 N >, N, implies

I p(7~.~o > T) - p(,,

> T)! < c/6.

This together with the bound on P(T~, > T) implies that

(v)

P(%,M, > q 2 1 - 43.

Let N,, = max(iVr , NJ. We assert that P,(F) - P,(F) < E VN 2 No . Take @ 3 No and define B = (~20,~~ > T). From (v),

P(Bc) < 43.

176

XIICHAEL TURELLI

Thus we have Pfi(F) ~ P,q,Mo(F) ~==P(Xzc E F, B) f P(Xfl E F, B’) - P(Xfi,MO t F, B) -- P(X& = P(Xfi since {X, Hence,

M E F, Bc)

EF, BC) - P(Xfi,A,O EF, B”),

E F, B) = (X, ,.,,n E F, B).

PAP)

- P,(F) = P#)

So, PR(F) -

- Pip,,,@)

-ALP,“(F)

Pfi,MB(F) < P(B”)

+ f’fi.~,JF)

- P,(F)

< c/3.

- %JF)

< E.

Y.E.U.

To take account of the fact, discussed in Section 2, that if f and g are not modified for large values of X, the solutions to (A2) map not behave sensibly, we state the following. COROLLARY. Suppose f and g of (A2) are replaced successively by iV and g, where fN = f and g,” = g for x < C,V and C, ---f cx, as iV + CO. Then if me denote the corresponding processes obtained from the modi$ed (A2) by X,Vx, S,” 5 X in (D[O, OS), a’).

Proof. Identical to that of the Proposition, constants M should be replaced by C,M - 1.

except

that

the truncation

ACKNOWLEDGMENTS The author gratefully acknowledges many productive and exhausting discussions with his doctoral advisor Joseph Felsenstein. Conversations with and encouragement from T. W. Schoener, K. B. Erickson, J. Gillespie, J. Roughgarden, R. Blumenthal, M. Slatkin, M. Feldman, and B. Levikson gave rise to many of the ideas presented. Comments on an earlier draft by S. Karlin and an anonymous reviewer were also very helpful. This research was supported by fellowship funds provided by the Ford Foundation under Grant 680-0183 A and by the National Institutes of Health under Grant 5 TO1 G&101269.

REFERENCES ABRAMS, P. 1976. Niche overlap and environmental variability, Moth. Biosc. 28, 357-372. ARNOLD, L. 1974. “Stochastic Differential Equations: Theory and Applications,” WileyInterscience, New York. AYALA, F. J., GILPIN, M. E., .~ND EHRENFELD, J. G. 1973. Competition between species: Theoretical models and experimental tests, Theor. Pop. Biol. 4, 331-356. of Probability Measures,” Wiley, Sew York. BILLINGSLEI', P. 1968. “Convergence

STOCHASTIC

CALCULUS

177

P. 1971. “Weak Convergence of Measures: Applications in Probability,” Philadelphia. BIRKHOFF, G., AND ROTA, G. C. 1969. “Ordinary Differential Equations,” 2nd ed., Blaisdell, Wltham, Mass. CAPOCELLI, R. M., AND RICCIARDI, L. M. 1974. A diffusion model for population growth in a random environment, Theor. Pop. Bid. 5, 28-41. CROW, J. F., AND KIMURA, M. 1970. “An Introduction to Population Genetics Theory,” Harper & Row, New York. FELDMAN, M. W., AND ROUGHGARDEN, J. 1975. A population’s stationary distribution and chance of extinction in a stochastic environment with remarks on the theory of species packing, Theor. Pop. Biol. 7, 197-207. FELLER, W. 1952. The parabolic differential equations and the associated semigroups of transformations, Ann. Math. 55, 468-519. FORD, E. B. 1965. “Ecological Genetics,” Methuen, London. GRAY, A. H., AND CAUGHEY, T. K. 1965. A controversy in problems involving random parametric excitation, J. Math. Phys. 44, 288-296. GIHMAN, I. I., AND SKOROKHOD, A. V. 1969. “Introduction to the Theory of Random Processes,” Saunders, Philadelphia. GILLESPIE, J. H. 1973. Polymorphism in random environments, Theov. Pop. Biol. 4, 193-19s. GILLESPIE, J. H. 1975. The role of migration in the genetic structure of populations in temporally and spatially varying environments. I. Conditions for polymorphism, Amer. Nutur. 109, 127-135. GILLESPIE, J. H., AND GUESS, H. 1976. The effects of environmental autocorrelation on the progress of selection in a random environment, Amer. Natur., to appear. Guess, H. A. 1973. On the weak convergence of Wright-Fisher models, Stochastic PYOC. Appl. 1, 287-306. GUESS, H. A., AND GILLESPIE, J. H. 1976. Diffusion approximations to linear stochastic difference equations with stationary coefficients, J. Appl. Probability, to appear. IBRAGIMOV, I. A., AND LINNIK, Yu. V. 1970. “Independent and Stationary Sequences of Random Variables,” Wolters-Noordhoff, Groningen. KARLIN, S., AND LEIBERMAN, U. 1974. Random temporal variation in selection intensities: Case of large population size, Theor. Pop. Biol. 6, 355-382. KARLIN, S., AND LEIBERMAN, U. 1975. Random temporal variation in selection intensities: One-locus two-allele model, 1. Math. Biol. 2, I-17. KARLIN, S., AND LEVIKSON, B. 1974. Temporal fluctuations in selection intensities: Case of small population size, Theor. Pop. Biol. 6, 383-412. KEIDING, N. 1975. Extinction and exponential growth in random environments, Theor. Pop. Biol. 8, 49-63. KHAS’MINSKII, R. Z. 1966. A limit theorem for the solutions of differential equations with random right-hand sides, Theor. Probability Appl. 3, 390-406. KIMURA, M. 1954. Process leading to quasifixation of genes in natural populations due to random fluctuation of selection intensities, Genetics 39, 280-295. KUSHNER, H. J. 1974. On the weak convergence of interpolated Markov chains to a diffusion, Ann. Probability 2, 40-50. LADDE, G. S., AND SILJAK, D. D. 1975. Stability of multispecies communities in randomly varying environment, J. Math. Biol. 2, 165-178. LEVIKSON, B. 1976. Regulated growth in random environments, J. Math. Biol. 3, 19-26. LE~IKSON, B., AND KARLIN, S. 1975. Random temporal variation in selection intensities acting on infinite diploid populations: Diffusion method analysis, Theor. Pop. Biol. 8, 292-300. BILLINGSLEY,

SIAM,

178

MICHAEL

TURELLI

LINDVALL, T. 1973. Weak convergence of probability measures and random functions in the function space D[O, 01, J. Appl. Probability 10, 109-121. LOEVE, M. 1963. “Probability Theory,” 3rd ed., Van Nostrand, New York. LUDWIG, D. 1974. “Stochastic Population Theories,” Lecture Notes in Biomathematics, Vol. 3, Springer-Verlag, Berlin. LUDU’IG, D. 1975. Persistence of dynamic systems under random perturbations, SIAM Rev. 11, 605-640. MAY, R. M. 1973. “Stability and Complexity in Model Ecosystems,” Princeton Univ. Press, Princeton, N. J. MAY, R. M. 1974. On the theory of niche overlap, Theor. Pop. Biol. 5, 297-332. MAU, R. M. 1975. Biological populations obeying difference equations: Stable points, stable cycles, and chaos, J. Theor. Biol. 49, 511-524. MAY, R. M., AE.‘DMACARTHUR, R. H. 1972. Niche overlap as a function of environmental variability, Proc. iVat. Acad. Sci. U.S.A. 69, 1 109-I 1 13. MAY, R. M., AND OSTER, G. F. 1976. Bifurcations and dynamic complexity in simple ecological models, Amer. Latur. 110, 573-599. MCKEAN, H. P. 1969. “Stochastic Integrals,” Academic Press, New York. MCSHANE, E. J. 1974. “Stochastic Calculus and Stochastic Models,” Academic Press, New York. MORTENSEN, R. E. 1969. Mathematical problems of modeling stochastic nonlinear dynamic systems, J. Statist. Phys. 1, 271-296. MORTON, J. B., AND CORRSIN, S. 1969. Experimental confirmation of the applicability of the Fokker-Planck equation to a non-linear oscillator, J. Math. Phys. 10, 361-368. Processes and Learning Theory,” Academic Press, NORMAN, M. F. 1972. “Markov New York. PAPANICOLAOU, G. C., AND KOHLER, W. 1974. Asymptotic theory of mixing stochastic ordinary differential equations, Comm. Pure Appl. Math. 27, 641-668. PROHOROV, Yu. V., AND ROZANOV, Yr.. A. 1969. “Probability Theory,” Springer-Verlag, New York. SAATY, T. L. 1967. “Modern Nonlinear Equations,” McGraw-Hill, Ne&- York. SCHOENER, T. W. 1973. Population growth regulated by intraspecific competition for energy or time: Some simple representations, Theor. Pop. Biol. 4, 56-84. SKOROKHOD, A. V. 1958. Limit theorems for Markov processes, Theov. Prob. Appl. 3, 202-246. STRATONOVICH, R. L. 1966. A new representation for stochastic integrals and equations, SIAM J. Control 4, 362-371. TUCKWELL, H. C. 1974. A study of some diffusion models of population growth, Theor. Pop. Biol. 5, 345-357. TIICKWELL, H. C. 1976. The effects of random selection on gene frequency, Math. Biosci. 30, 129-l 37. TLJRELLI, M. 1977. “Random Environments, Stochastic Calculus, and Limiting Similarity,” Dissertation, University of Washington, Seattle. WHITT, W. 1970. Weak convergence of probability measures on the function space C[O, m], Ann. Math. Statist. 41, 939-944. WONG, E. 1971. “StochasticProcesses in Information and Dynamical Systems,” McGrawHill, New York. WONG, E., AND ZAKAI, M. 1965. On the convergence of ordinary integrals to stochastic integrals, Ann. Math. Statist. 36, 156@-1564.