J. theor. Biol. (1977) 68, 143-160

An Evaluation of Linear Models of Population Fluctuations R. M. NISBET, W. S. C. GURNEY AND M. A. PETTIPHER Department of Applied Physics, University of Strathclyde, Glasgow G4 ONG, Scotland (Received 15 December 1976, and in revisedform 11 March 1977) The analysis of non-linear stochastic population models is notoriously difficult and there is little or no data against which they may he tested. It is therefore important to investigatethe range of validity of modelsof fluctuating populationswhich are linear in the deviation of the population from its mean value. In this paper we demonstrate the robustnessof somewell-known, linear methods of analysing fluctuations, the power of the methods being evaluated by comparison with numerical results for somenon-linear models involving demographic or environmental stochasticity.

1. Introduction The central problem of population dynamics is to extract information on the mechanisms of population regulation from a time series of population values. The implicit assumption in much of this work is that the time series may be represented as a superposition of deterministic and random components, the aim of statistical work being to separate these and to produce a functional representation of the deterministic part. This problem of

separating “signal” from “noise” is notoriously intractable and we know of no natural populations and of few laboratory populations where it has been possible to meaningfully fit a non-linear? model to existing data.

On the other hand, non-linear representations of population fluctuations have received considerable attention from those who study model populations. In such models fluctuations may arise from both “demographic stochasticity” (caused by the fact that such populations contain a discrete number of individuals) and “environmental

environmental

fluctuations

which ultimately

stochasticity”

(caused by

lead to birth and death rates

1 In this paper we define a linear model to be a model with equations which are linear in the deviation of a populationfrom its meanvalue. 143

144

R.

M.

NISBET,

W.

S.

C.

GURNEY

AND

M.

A.

PETTIPHER

which are functions of time). May (1973a,b) gives a lucid introduction to the analysis of environmental stochasticity while Bartlett (1957, 1960) and Bartlett, Gower & Leslie (1960) give an authoritative analysis of demographic stochasticity. Considerable mathematical subtlety is required to analyse non-linear stochastic models (see, for example, Mortensen, 1969), and if the effects of non-linearity are experimentally unobservable, it seems reasonable to ask if these mathematical difficulties may normally be circumvented? In this paper we answer that question in the affirmative by demonstrating the robustness of some simple, well-known, linear methods of analysing fluctuations, the utility of the methods being assessed by comparison with numerical resuIts for several non-linear models. We find that the approximate results obtained may be good even for rather large fluctuations. This work is an extension and generalization of the work of May (19736) who emphasized the utility of simple approximate treatments of population fluctuations, and used one such method to analyse otherwise intractable competition models with environmental stochasticity. A similar approach has also been used by Roughgarden (1975) on a simple discrete-time model.

2. Theory (A)

SPECTRAL

ANALYSIS

Suppose we know a sequence of values, N(t), of a fluctuating population over a long time interval, and that this time series shows no obvious deterministic trend. It is natural to attempt to model the population time series as a realization of some stationary stochastic process; indeed many statistical forecasting techniques rely on this procedure (see, for example, Box & Jenkins, 1970; Poole, 1976). The population variance 0; defined by? CT;f (N(t)2) - (N(f))2

(1)

is normally used as a measure of the intensity of the population while the autocorrelation function (ACF) C(z), defined by C(z) = ((NO - = IwH(f>l

2.

on the popu(20

The quantity T(f) is known as the transfer function and its properties are well documented in the signal analysis and control theory literature (e.g. Lynn, 1973). The simplest interpretation of the transfer function is reached by considering the special case of a small sinusoidal input r(t). The population will eventually fluctuate sinusoidally at the same frequency, and the transfer function is then a complex function whose modulus is the ratio of the amplitude of the population oscillations to that of the input oscillations, and whose argument is the phase difference between these oscillations. The modulus of the transfer function has characteristic forms for overdamped, critically damped, and underdamped systems which we have illustrated in Fig. 1. It is important to note that at high frequencies the transfer function is always very small, the system showing no significant response to perturbations at frequencies much greater than the reciprocal of the smallest ecological time constant in the system. We may thus safely represent random environment fluctuations as white noise if (and only if) the relevant spectral density is flattish over the range of frequencies for which the magnitude of the transfer function is significant. Where the input noise is not white its spectrum must be measured or guessed at before we can calculate the population spectral density from (21). However, no new principles are involved in the spectral analysis of non-white noise, a point which has been recently emphasized by Roughgarden (1975). When the environmental noise is white the transfer function(s) of a system alone determine the nature of the population fluctuations. For example,

MODELS

OF

POPULATION

v

FLUCTUATIONS

149

. I (

Frequency

(orbitrory

units)

FIG. 1. Typical forms for the modulus of the transfer function for (a) overdamped, (b) critically damped, and (c) underdamped systems.

with a single randomly fluctuating parameter in a population demographic stochasticity ai = (x2) = 2a2 f\7’(j)i2

df

with negligible (22)

and C(r) = 2a2 rcos 27rfrlT(f)12 df,

(23

where a2 is the spectral density of the white noise. Where there are several sources of environmental noise, we can (still only within the linear approximation) obtain the population spectral density as a sum of terms of the type CW. (D)

CYCLING

POPULATIONS

Considerable interest has been generated over the years by the small number of animal populations which exhibit “quasi-cyclic” fluctuations. In a recent paper (Nisbet & Gurney, 19766) we argued that these oscillations might arise from the “filtering” of environmental or demographic fluctuations by a system whose transfer function had a sharp peak. We argued that provided the spectral density of the environmental fluctuations is flattish over the range of frequencies for which the magnitude of the transfer function is appreciable, then quasi-cyclic fluctuations will result. We now use the theory from the previous sections to develop a more detailed description of these quasi-cyclic fluctuations. We shall define

150

R.

M.

NISBET,

W.

S.

C.

GURNEY

AND

M.

A.

PETTIPHER

quasi-cyclic fluctuations to be those for which the autocorrelation function (ACF) consists of a series of damped oscillations. We define the coherence number n, as the number of cycles required for the amplitude of the ACF to be reduced by a factor e. The long-term pattern of such fluctuations consists of irregular bursts of “coherent” oscillation each containing two to three times n, cycles (since this is the time scale over which “memory” of phase information is lost), interspersed with short periods of incoherent noise. We demonstrate an example of this pattern of behaviour in section 3(B). Thus if we are to perceive a truly cyclic character in the fluctuations we require n, Iv > 1. For example, in many systems the ACF has the form C(T) = A e-“’ cos WT,

(24)

n, = w,/27d.

(25)

for which In such a system the fluctuations

will appear cyclic if i. 5 w,/(2n).

(26)

The appearance of the fluctuations at very high values of n, depends somewhat on the details of the underlying population dynamics. We have not explored this point in detail, but it is clear that if n, 2 100 the fluctuations observed in most normal situations would appear almost perfectly regular. 3. Tests of the Linear Theory (A)

DEMOGRAPHIC

STOCHASTICITY

We have tested the previous approximate treatment of demographic stochasticity on a simple birth and death model with j?(N) = aN and 6(N) = uN2/K. Then No = K and from equation (17) the population variance is also K. By solving (1 l), it can be shown that the steady-state probability distribution is, in the linear approximation, Gaussian and hence T,, the mean time to extinction is of order of magnitude eKi2. Each of the above quantities? (mean, variance and extinction time) may be computed numerically from the exact differential equations describing the birth and death process. The results are shown in Table 1 for a number of values of K, The following points are evident: (i) the mean population is approximately K- 1 for all but very small K; (ii) the predicted variance is valid for K > 10; (iii) the prediction of extinction time is poor. t Strictly speaking, we evaluated the mean and not lzecoming extinct.

variance

conditional on the population

MODELS

OF

POPULATION TABLE

151

FLUCTUATIONS

1

Exact values of the population mean, variance, and logarithm of mean extinction time for a “logistic” birth and death model ___-

K

09

2 3 5 7 10 15 20 40 60 80 100 I20

1,871 2473 4.004 5.853 8.556 13.91 18.94 38.97 58.98 78.99 98.99 119-o

According

modifications In

(arc)

to the linear in Appendix

theory B change

UN2 -

1.128 2.055 4.416 6.895 IO.16 15.11 20.07 40.02 60.01 80.00 100~0

120-o = crNa = Kand the predictions

--

ln (ard

-.

144 2.23 3.84 5.42 7.95 12.4 17.1 36.4 56.0 75.6 95.4 115.2 In (ma) = fK. The non-linear to (N) = K - 1, aNz =

K,

== K.

Thus even the very crude approximate equation (11) is satisfactory in predicting the population mean and variance for all but the smallest populations. It is not surprising that the estimate of the order of magnitude of the extinction time is poor since this will obviously be affected by the detailed shape of the distribution near N = 0. However, a small increase in sophistication of treatment improves the accuracy of the estimate of extinction time and explains the pull-down of the mean population to K- 1. As this is only of mathematical interest we leave the demonstration of these results to Appendix B. Having demonstrated the robustness of the linear estimates of population mean and variance, we tested the predictions of the dynamic characteristics of the population fluctuations (i.e. the spectral density and the autocorrelation function). Using our previous expressions for j?(N) and 6(N) we obtained numerical realizations of the non-linear stochastic differential equation (9) over a range of values of K. Because of the well-known difficulties in estimating spectral densities from realizations of a stochastic process [see section 2(A)], we chose to calculate the autocorrelation function for each run and to compare the result with the linear estimate given by equation (18). Typical examples are shown in Fig. 2 for a large value of K (K = 100) and a small value (K = 10). The robustness of the linear predictions is obvious.

152

R.

M.

NISBET,

W.

S.

C.

GURNEY

I

AND

M.

I

I

I

2

3

4

A.

PETTIPHER

5

or

FIG. 2. Predicted and observed autocorrelation demographic stochasticity. (B)

ENVIRONMENTAL

STOCHASTICITY LOGISTIC

functions for a logistic model with

IN

THE

TIME-DELAYED

MODEL

We tested the linear expression (22) for the intensity of fluctuations due to environmental stochasticity on the time-delayed logistic model, chosen as the simplest model of a single population in which the return to a stable equilibrium may be either overdamped or underdamped. For stochastic analysis it is convenient to write this equation in the form

WO __ = rN(t){l-flN(t--r)}, dt

where r is the intrinsic growth of the species, fi-’ is the carrying capacity, and t a time delay. We model a fluctuating carrying capacity by setting PO> = BOO + r(O) (28) and define our units of population such that j?,, = 1. The popuIation fluetuations may be derived in the linear approximation from equation (21) with the transfer function

T(f) = - l/[cos 27fz - i(sin 27cfT- 27cf/r)].

(29

We discussed the behaviour of this transfer function in a previous paper (Nisbet & Gurney, 1976a) but for completeness we recall that the system

MODELS

is underdamped

OF POPULATION

FLUCTUATIONS

153

if e-l < r5 < 42

(30)

and overdamped if rz < e-‘. We examined the robustness of the linear theory as an approximate description of the population fluctuations resulting from large fluctuations in the carrying capacity. With a white noise input the fluctuations should be given by (22) which for this model implies (31) an integral which is easily evaluated numerically. We compared the population mean and standard deviation observed in realizations of the model with that predicted by the linear theory. In all the realizations we kept Y = 1, and varied Q and r. In Fig. 3 we show the variation in observed standard deviations as c?, the variance of the noise in carrying capacity, is increased, for an overdamped system (T = 0.25) and an underdamped system (r = l-2). The very encouraging feature of these results is that the linear approximation is reasonable right up to the point where the probability

FIG. 3. Predicted and observed standard deviations in a population described by the time-delayed logistic equation with (a) rr = 0.25, an overdamped case, and (b) rz = l-2, an underdamped case. The standard deviations are averages from four realizations, each of which ran to a time 5OOr-‘. c2 is the variance of the white noise in 8, the reciprocal of the carrying capacity. 7.B. II

154

R.

M.

NISBET,

W.

S.

C.

GURNEY

AND

M.

A.

PETTIPHER

of extinction becomes appreciable on an ecological time scale. The approximation is particularly strong for the overdamped case, where even with o‘ = 0.50 there is only a 4% discrepancy between the size of predicted and observed fluctuations. The approximation is clearly poorer with an underdamped system, the predicted and observed results diverging steadily for B > 0.20 and disagreeing by 12% at rr = 040. However, even this accuracy compares well with the precision attainable in field or laboratory measurements of the intensity of population fluctuations. For this model we also tested the robustness of the prediction of the autocorrelation function, the linear estimate of which was derived by using (29) to obtain the spectral density and then numerically integrating the expression in equation (23). The effect of non-linearities was assessed by comparison of this prediction with autocorrelation functions derived from realizations with small values of c (for which the linear theory should be very accurate) and with large values of cr (where non-linearities might be important). We performed our tests on an underdamped example (r7 = 1.4) for which the autocorrelation function consists of a set of damped oscillations. We eliminated the effect of sampling error which (see, for example, Bartlett, 1966, p. 302) may cause an observed autocorrelation function to exhibit less damping than a theoretical one, by comparing runs in which the realization of the environmental noise was generated by a single sequence of random numbers which were then multiplied by the appropriate value of 6. Every test yielded results like those in Fig. 4, from which we conclude that very large non-linearities will increase the damping in the autocorrelation function while the effect of smaller non-linearities will normally be negligible in comparison with sampling error. For this example (r-c = 1.4) we were also able to demonstrate the utility of the coherence number [defined in section 2(~)] in characterizing quasicyclic population fluctuations. For this case n, - 2.8 and examination of realizations like that in Fig. 5 reveals precisely the “bursts’ of oscillations we should expect with this coherence number. (C)

ENVIRONMENTAL

STOCHASTICITY

IN

A PREDATOR-PREY

MODEL

As a second example of environmental stochasticity, we consider a twospecies Lotka-Volterra predator-prey model with self-limitation of the prey (see, for example, Maynard Smith, 1974, p. 19). With a suitable choice of units for both populations and time we may write the equations as P=P(I-/P-C),

(32)

c = -C(b--P).

(33)

MODELS

OF

POPULATION

I 5

FLUCTUATIONS

I IO

7

155

I

I

15

20

25

FIG. 4. Predicted and observed autocorreiation functions for the time-delayed logistic model with white noise in the reciprocal of the carrying capacity. For each graph rr = 1.4, an undamped case. -, Linear theory; -,withv=O.Ol;-----,withg-0.1; *....a., with cr = 0.25.

I

I

I

20

I

I

40

I 60

I

I 80

I

Time

FIG. 5. An example of “quasi-cycles” an undamped case.

in the time-delayed logistic model with rT = 1.4,

156

R.

M.

NISBET,

W.

S.

C.

GURNEY

AND

M.

A.

PETTIPHER

We now assume that there are random fluctuations in S, the predator death rate, and set 6 = 6,+y(t) (34) and t = P-PO,

q=c-c

0’

where PO, Co are respectively the steady-state populations predator. Then F(f) = TIUMf)

and

(35) of prey and

V(f) = T,(fW)~

(36)

with (37)

WOO -Do> + W1 - PSo)f 4n f - i27$6,f+ S,(j?S, - 1)’

7-2(f) = -2-y--

(38)

When r(t) is white noise of spectral density 02, the mean square fluctuations in each population may be calculated by substituting for T,(f) and T,(f) in equation (22) to get -mob2

(39)

-pso>(p2s,-ps,+1)~2.

(40)

e2> = W’U and (a2> = (2p&J-‘(I

We tested the effect of non-linearities on these approximations by running realizations for an overdamped and an underdamped case, with various values of 0. The results are shown in Fig. 6. As with the previous model, our conclusion is that the linear theory yields pretty robust estimates of fluctuation intensity. For this model the linear estimates of the ACF for prey (C,(z)) and for predator (C,(7)) can be obtained analytically by substituting for Tr(f) and T2(f) in equation (23) to get in the underdamped case (which is the case of interest as regards the occurrence of quasi-cyclic behaviour) C&7) = (t2> exp (-j/3S07)(cos o1 7 + (/36,/2co,) sin o1 z), C,(T) = (#) exp (-3~So7>(cos 0, s-(j?So/ol) sin 0, z),

(414 (41b)

where 01 = {60(1- BJO) - tB2G)+, (42) while (c2) and (1’) are given in equations (39) and (40). It is evident that significant oscillatory behaviour only occurs if PS, $ cur in which case the

MODELS

OF

POPULATION

157

FLUCTUATIONS

_--.

r

o’4G--

n

/

O-I

02

o-3

FIG. 6. Predicted and observed standard deviations in prey and predator populations described by a damped L.otka-Volterra model. The parameter values are (a) /I = So = 0.5, an underdamped case, and (b) j = 2, A, = 0.35, an overdamped case.

coherence number n, is easily shown to be n, = w,/(7cg8,) ” ?r-‘P-‘s,+(l --/Is,>+. (43) We have previously (Nisbet & Gurney, 1976b) shown an example of large amplitude quasi-cycles for the case fl = O-1, 6 = 0.91 (i.e. n, = 3.2). 4. Conclusions The aim of this work was to test the accuracy of some simple linear methods of analysing large population fluctuations. We have demonstrated how to estimate the population variance, spectral density and autocorrelation function (ACF) for systems with demographic and/or environmental stochasticity. We tested the approximate treatment of demographic stochasticity on a simple birth and death model chosen to resemble the logistic model. Except with very small populations the linear theory gave good predictions of mean, variance and ACF, but a non-linear treatment was necessary to get accurate predictions of extinction times. We tested the approximate treatment of environmental stochasticity on two models-the time-delayed logistic with a fluctuating carrying capacity and a damped Lotka-Volterra predator-prey with a fluctuating predator death rate. In both cases there was acceptable agreement between the intensities of predicted and observed population fluctuations right up to the point where the fluctuations were large enough for there to be an appreciable probability of extinction on an ecological time scale. For one of the above models (the time-delayed logistic) we demonstrated the reliability of the linear estimate of the ACF and showed that if the system is underdamped, the effect of

158

R.

M.

NISBET,

W.

S.

C.

GURNEY

AND

M.

A.

PETTIPHER

non-linearities is to slightly increase the damping in the ACF. For the special case of an underdamped system we introduced a “coherence number” which provides a measure of the extent to which the population fluctuations will be “quasi-cyclic” in appearance. The message to emerge from these numerical investigations is thus one of reassurance. From linear approximations we have obtained representations of population fluctuations which prove very robust when tested on data generated from familiar non-linear models. We hope these results give encouragement to those workers who adopt simple, linear, statistical representations of field and laboratory data. We thank J. McKinstry for detailed comments on an early draft of this paper. One of us (M.A.P.) held an S.R.C. studentship while carrying out this work. REFERENCES BARTLETT, BARTLETT,

M. M.

S. (1957).Biometrika 44, 27. S. (1960). Stochastic Population

Models. London and Colchester: Spottiswoode, Bailantyne& Co. Ltd. BARTUIT, M. S. (1966).An Introduction to Stochastic Processes, 2nd ed. Cambridge: CambridgeUniversityPress. BARTLETT, M. S., GOWER, J. C. & LESLIE, R. H. (1960).Biometrika 47, 1. Box, G. P. & JENKINS, G. M. (1970). Time Series Analysis-Forecasting and Control. San Francisco:Holden-Day. GURNEY, W. S. C. & NISBET, R. M. (1976). Am. Nut. (in press) JENNISON, R. C. (1961).Fourier Transforms and Convolutions for the Experimentalist. Oxford and New York: PergamonPress. LYNN, P. A. (1973). An Introduction to the Analysis and Processing of Signals. Londonand Basingstoke:Macmillan. MAY, R. M. (1973a).Stability and Complexity in Model Ecosystems. Princeton:Princeton UniversityPress. MAY, R. M. (19736). Am. Nat. 107,621. MAYNARD SMITH,J. (1974).Models in Ecology. Cambridge:CambridgeUniversityPress. MORTENSEN, R. E. (1969). J. statist. Phys. 1,271. NISBET, R. M. & GURNEY, W. S. C. (1976a). .Z. theor. Biol. 56,459. NISBET, R. M. & GURNEY, W. S. C. (19766). Nature, Land. 263, 319. PLATT, T. & DENMAN, K. L. (1975). A. Rev. Ecol. System. 6, 189. POOLE, R. W. (1976). Theor. Pop&. ROUGHGARDEN, J. (1975). Am. Nat.

Biol. 9 25.

109,713. WAX, N. (1954).Selected Papers on Noise and Stochastic Processes. New York: Dover.

APPENDIX A

Combined Environmental

and Demographic Stochasticity

Consider a single species birth and death model, in which both the birth t) and the death rate 6’(N(t), t) are explicit functions of time.

rate B’(N(t),

MODELS

OF

POPULATION

FLUCTUATIONS

159

Let (Al)

P’W(t)v 0 = P(N) + e1tt>, ww), 0 = W)+e,tt),

&‘I

and e(t) = q(t)-+(t). (43) Suppose cl(t), e,(t), and any population fluctuations are small. Then the fluctuation spectrum of the population, in the linear approximation, is easily shown to be Z(f)

= (r+2b~f)-~{Z(f)+(2Q)~

e’“),

where 4 is a phase factor which will vary from realization to realization. Provided the environmental noise e(t) is totally uncorrelated with the demographic population fluctuations, the resulting population spectral density will be (in the notation used in the text) (A4) which is of the form of a sum of terms representing environmental and demographic stochasticity. This type of analysis is easily extended to the case of many species and also the case of discrete or distributed time-delays in birth or death rates (by using the convolution theorem). No new principles are involved and the final expression for the spectral density may always be written as a sum of terms representing environmental and demographic stochasticity.

APPENDIX

B

Further Comment on Demographic Stochasticity in a Logistic Model In the main text we considered a linear treatment of a single species birth and death mode1 with p(N) = c&’ and 6(N) = c@/K. Strictly the evolution of this system is controlled by the non-linear stochastic differential equation dIV(t) = ctN(1 -N/K) dt+ {crN(l +N/K)}” dZ(t), (A5) where dZ(t) is a Poisson variable with mean zero and variance dt. When we linearized this equation we predicted (N)=K, o~=(N)-(N)~=K, 2’ 7E cc efK . C-46) Bartlett (1960, p. 32) has shown how to compute the next order of approximation for these quantities. It turns out for our model that the mean popu-

160

R.

M.

NISBET,

W.

S.

C.

GURNEY

AND

M.

A.

PETTIPHER

l&ion should be pulled down by an amount &(N) = 1, i.e. the second-order prediction of mean population is (N)=K-1. (A7) This is in very close agreement with the exact results in Table 1. Our poor estimate of the extinction time was due to the replacement of the Poisson increment dZ(t) with the Wiener increment dw(t). The assumption of a Wiener process led us to a normal distribution with mean and variance given by (A6). When we realize that the final distribution should be much closer to a Poisson distribution with mean K, extrapolation to N = 0 suggests that the extinction time should be proportional to eK (and not eK”). This is consistent with the exact results in Table 1.