- Email: [email protected]

Available at

MATHEMATICAL

www.EIsevierMathematics.com l ONIISD .I SCICNCC @ DINCDT-

COMPUTER MODELLING

Mathematical

and

Computer

Modelling

38 (2003)

1437-1442 www.elsevier.com/locate/mcm

Population Growth as a Nonlinear Stochastic Process K. SOBOLEVA

AND A. B. PLEASANTS AgResearch, Ruakura Research Centre Private Bag 3123, Hamilton, New Zealand

T.

Abstract-The evolution of the probability density of a biological population is described using nonlinear stochastic differential equations for the growth process and the related Fokker-Planck equations for the time-dependent probability densities. It is shown that the effect of the initial conditions disappears rapidly from the evolution of the mean of the process. But the behaviour of the variance depends on the initial condition. It may monotonically increase, reaching its maximum in the steady state, or have a rather complicated evolution reaching the maximum near the point where growth rates (not population size) is maximal. The variance then decreases to its steady-state value. This observation has implications for risk assessments associated with growing populations, such as microbial populations, which cause food poisoning if the population size reaches a critical level. @ 2003 Elsevier Ltd. All rights reserved. Keywords-Nonlinear evolution,

Risk

stochastic

processes,

Population

growth,

Fokker-Planck

equations,

Mean

assessments.

1. INTRODUCTION The study of growth is a central problem in biology. To describe growth, biologists have used a differential equation in the form dx = (a - ,x,-l)

xdt,

m 1 2.

(1)

Here z is the size of an organism or a biological population and LY is the Malthusian growth parameter. The second term in equation (1) describes the restriction in growth, which occurs due to the effects of crowding and competition for resources, typical in biological growth. An example of this formulation is the logistic equation (m = 2). Another equation, which incorporates a limit on size, is the Gompertz equation dx = (a - plnx)xdt.

(2)

These equations are deterministic, and ignore the effect of random fluctuations on the growth of the population. It is generally assumed that such random fluctuations are independently and identically distributed about the growth path described by equation (1) or (2). In each case it is assumed that the fluctuations do not interact with the growth law described by the deterministic equation. The

help

of K. Vetharaniam

0895-7177/03/L - see front doi: lO.l016/SOS95-7177(03)00358-3

is gratefully matter

@ 2003

acknowledged. Elsevier

Ltd.

All rights

reserved.

Typeset

by .A.-‘QX

1438

T. K.

SOBOLEVA

AND A. B. PLEASANTS

A more general approach is to recognise that the fluctuations are an intrinsic part of the growth process and include them into the formulation of the equations. This approach does not involve any unreasonable assumptions of independence of the growth path and previous fluctuations. These are naturally included in dynamics. This paper investigates the situation when the deterministic growth equations are replaced by their stochastic forms, noting particularly the relationship between the mean and the variance of the process as growth proceeds.

2. THE

MODEL

There are different ways to incorporate stochastic components into deterministic models (1) or (2). In the process under consideration the fluctuations are assumed to be faster than rmacro = l/o, which characterises the time scale of the evolution of the macroscopic variables in these equations. Random fluctuations, such as rapid environmental changes, affect the system through external parameters. Suppose that the parameter Q in equation (1) or (2) is regarded as a random variable in the form at = a+&, (3) where CYis the mean, & is Gaussian white noise, and (T is the intensity of the noise. Then equations (1) and (2) are replaced by stochastic differential equations (SDE) for the random process Xt . dX, = (a - PXF-‘) X, dt + aXt dW, (4 and dt + ax, dW,,

(5)

respectively. W, in equations (4) and (5) is a Wiener process. For convenience, ,Ll = (Ylnx* has been substituted in equation (5) and it is assumed that the variable x is dimensionless, i.e., the equations describe changes in relative population size. Aspects of these models have been considered by different authors (see, for example, [l] and references therein, [2,3]). The choice of an appropriate diffusion process, i.e., the choice between the Ito and Stratonovich interpretations of the above SDEs depends on the nature of the process being modelled. In the case considered here there are no qualitative differences in any of the results between either approach. However, for definiteness, the Ito stochastic calculus is used here.

3. THE

RESULTS

Each solution of an SDE describes one realisation, or growth path. The ensemble of realisations generated by an SDE is described by a transitional probability density P(x, t), which satisfies the corresponding Fokker-Planck equation

-ap(x, t) = -;

at

f(x)P(x, t) + ; g

x2P(x, t),

where f(x) = (CY- Pxmml )x or ax In x*/x, respectively. The steady-state probability Ps(x) d ensity can be evaluated from this. For equation P,(x) = N8x2(alua-1)

exp

- 0~~~~:,> (

And for the Gompertz

(4) this is (7)

model, p&)

= Nsx+ln~‘l~w

exp(-+),

(8)

Population

1439

Growth

where N, is a normalising constant. Note that the probability density (8) is normalised and finite at any positive value of the parameters. The probability density (7) is normalised (i.e., the steady-state probability exists) for cr > u2/2. The condition cr > a2/2 coincides in the Ito interpretation with the condition that z = 0 is a natural boundary for the stochastic process. At a2/2 < cx < g2 the probability density (7) is divergent at z = 0. This means that the class of models (4) has a noise induced phase transition at the point a = g2. Horsthemke and Lefever [l] give details on this phenomena. Henceforth, we suppose that cr > u2 and that the process evolves toward the steady-state density (7) or (8). The solution of the Fokker-Planck equation (6) describes the evolution of the probability density through time. The class of models given by equation (1) all demonstrate the same qualitative behaviour. All allow the substitution of y = zVrn+l, which linearize the SDEs. In the case when m = 2, the substitution y = l/z changes the logistic SDE into the linear SDE dY, = [(a2 - a) Yt + p] dt + aY, dW,

(9)

with solution Y,=exp((G-a)ttoW,){Y~+fi~exp[(-(G-a)s-oWt)]

ds}.

(10)

In this solution the Wiener process is in the exponential function, so that the Gaussian fluctuations undergo a nonlinear transformation. Thus, the solution of this linear SDE is not a Gaussian process. The linearization procedure cannot assist in the determination of the time dependent solution of the Fokker-Planck equation (6). However, it can be found numerically. The results of a numerical solution of equation (6) will be presented below. For the Gompertz SDE the substitution y = lnz transforms the nonlinear process into the Ornstein-Uhlenbeck process. This means that the substitution can be used to find the time dependent solution of the Fokker-Planck equation for the Gompertz stochastic process. When the initial population size is known exactly, i.e., P(z, 0) = 6(z - Q), the solution is e,t

I zo) =

f-2(lnz - ln2* + a2/2a - eVat ln20)2 u2 (1 - e-2at) az~/7r (1 T e-2at) /a exp ( ).

(11)

Both the mean and variance from this time-dependent probability density increase monotonically towards the steady-state mean and variance. These statistics depend on the parameters of the SDE only. If only the probability of a range of initial population sizes is known the analytical solution of the Fokker-Planck equations has a complicated form of an infinite series over orthogonal polynomials. It could be useful to analyse the behaviour of the first two moments from the time-dependent probability density P(z, t). Ordinary differential equations for the mean mt and variance st can be derived from the Fokker-Planck equation (6) or directly from the SDE. When the instantaneous fluctuation is proportional to the state Xt, then different approximation procedures lead to the same equations dst -dt = 2ast - 4Pmtst + a2 (772: + st) , (12) dmt - dt =amt -P(m:+a) for the logistic stochastic process, and %=2ost(ln($)-2)+o’(mf+s,),

(13)

1440

T. K. SOBOLEVA AND A. B. PLEASANTS

.18

Figure 1. Numerical solution of equation (12) for the mean (m) and the variance (s) of the logistic stochastic process at (Y = 0.04; p = 0.005; CJ = 0.01. Bold curves correspond to the initial conditions m(O) = 0.5, s(0) = 0.0001; thin curves correspond to the initial conditions m(O) = 0.5, s(0) = 0.005.

1.8 8-

1.6

c

'

_

76s

Figure 2. Numerical solution of equation (13) for the mean (m) and the variance (s) of the Gompertz stochastic process at a = 0.03; I’ = 7.45; Q = 0.03. Bold curves correspond to the initial conditions m(O) = 0.07, s(0) = 0.0007; thin curves corre spond to the initial conditions m(O) = 0.07, s(O) = 0.007.

for the Gompertz case. The solutions of the set of equations (12) and (13) are shown in Figures 1 and 2, respectively. In each case, the mean and variance evolve to a steady state, which is unique for each system. Note that the values for the steady-state mean and variance could also be obtained from probability densities (7) or (8). But the path by which the variance reaches the steady state differs depending on the initial conditions. The main difference occurs in the behaviour of the variance over time. While the mean always monotonically evolves toward the steady state, the variance could quickly increase above its steady-state value, reach some maximum at a time when the rate of growth is near its maximum, and then decrease to the steady state. This situation for the logistic model is illustrated in Figures 1 and 2 for the

Population

Growth

Figure 3. Example numerical solutions of the Fokker-Planck equation for the logistic stochastic process, using the set of parameters a = 0.045; fl = 0.003; cr = 0.005. Curves l-3 show the evolving probability density at times tr < t2 < ts, respectively; curve 4 corresponds to the steady state.

Figure 4. Example numerical solutions of the Fokker-Planck equation for the Gomperte stochastic process, using the set of parameters (I = 0.03; z* = 15; u = 0.01. Curves l-3 show the evolving probability density at times tr < tz < t3, respectively; curve 4 corresponds to the steady state.

Gompertz process. When initially the ratio of the standard deviation to the mean of the initial population size is small, the variance increases monotonically to the steady-state variance. But when this ratio is initially close to or greater than this ratio at the steady state, then the variance increases to a point greater than the steady-state variance before shrinking back to the steadystate variance. It could be that this behaviour is due to the approximations that are made in deriving equations (12) and (13). To resolve this, direct numerical solutions of the Fokker-Planck equation for both the Gompertz and logistic processes were obtained. These solutions are plotted in Figures 3 and 4. The Gompertz stochastic process shows the same behaviour, suggesting that the behaviour of the variance does not depend on the approximations. The evolution of the probability density is shown for increasing time (from left to right). The last curve corresponds to the steady-state

1442

T. K. SOBOLEVA

AND A. B. PLEASANTS

distribution. The nonmonotonic changes of the variance of the stochastic growth processes are clearly expressed by nonmonotonic changes of the amplitudes of the probability density curves. For example, it is clear that the steady-state variance is greater than the process variance at t = ti, but less than the variance at t = t2 (tl < t2). A possible explanation for this behaviour is that the scale of the response of the mean and the variance differ as growth proceeds. If the ratio of the initial standard deviation to the initial mean is greater than the ratio of the steady-state standard deviation to the steady-state mean, then this difference in scale produces the behaviour observed. While the mean can only increase, the variance can increase above the steady-state variance at some time before the process gets to steady state, and then shrink back as growth rates shrink to zero. 4.

CONCLUDING

REMARKS

Growth equations are assumed to provide an adequate description of biological growth. Dependencies between the mean and the variance have been observed, particularly in cases of animal growth [4]. These results show that the relationship between the mean and variance may be more complex when nonlinear growth processes are involved. In particular it should be expected that circumstances will exist where the variance of the process will increase to a level above the steady-state variance before falling back to the steady-state variance. Intuitively, one might expect that the process would reach maximum variance at the steady state, and in this case it would be expected that a measurement of the variance made at the steady state would establish an upper bound. This analysis makes clear the danger of such an approach if, for example, the goal was to measure the risk to food safety or to water quality through micro-organism growth. Note that the behaviour described is a specific feature of nonlinear growth. Both the mean and variance of growth processes that can be described by linear differential equations (for example, limited or decreasing exponential growth) with a random growth rate, always monotonically evolve to their steady-state values. This work also suggests that the statistical fitting of growth equations to data would be enhanced by using a model which allowed for the variance to change with time.

REFERENCES 1. W. Horsthemke and 2. P. Holgate, Varieties Biology 139, 369-378, 3. C. Wang, Stochastic Modelling, Lyngby, 4. A.B. Plessants, W.H. strategies for control,

Ft. Lefever, Noise-Znduced ‘i%xasitions, Springer-Verlag, Berlin, (1984). of stochastic model: A comparative study of the Gompertz effect, Journal of Theoretical (1989). differential equations and a biological system, Ph.D. Thesis, Institute of Mathematical Denmark, (1994). McMillan and R.A. Barton, The evolution of liveweight variance in Angus steers and In Proceedings of the New Zealand Society of Animal Production (to appear).