- Email: [email protected]
The regulation of inhomogeneous populations
J. theor. Biol. (1975) 52, 441-457
The Regulation of Inhomogeneous Populations W. S. C. GURNEY AND R. M. NISBET Department of Applied Physics, University of Strathclyde, Glasgow G4 ONG, Scotland (Received 29 May 1974) A number of authors have argued that dispersal may play a key role in the regulation of populations of some species. In this paper we examine, on an ecological time scale, mechanisms whereby a population existing
in an inhomogeneous environment may use dispersal to regulate its sire below the carrying capacity set by the supply of available nutrients. We show that dispersal produced by wholly random motion is incapable of exerting any stabilizing influence, but that the introduction of a suitable non-linearity into the dispersal behaviour of a species whose characteristics are otherwise wholly linear can lead to stabilization under a wide range of conditions.
1. Introduction A number of author (e.g. Lidiker, 1962; Howard, 1960; Murray, 1967) have argued that emigration or dispersal may play a key role in the regulation of populations of some species. A persuasive example of this suggestion occurs in a paper by Carl (1971), whose observations on a population of arctic ground squirrels showed that this species achieves population control by the dispersal of animals from densely populated areas of favourable habitat into unfavourable areas where burrow sites are not available, and where they are subject to intensive predation. It seems plausible that adaptative advantage might accrue to a species with a stable strategy for limiting its numbers below the maximum carrying capacity of its environment, since it is thereby provided with considerable protection against fluctuations in the supply of available nutrients. Lidiker (1962) has argued that such a strategy might well operate through dispersal, and would then bring with it the additional advantages of increased heterozygosity and more probable genetic recombination. That dispersal can indeed operate to control a population at a level below the nominal carrying capacity of its locality has recently been most elegantly demonstrated for several vole species by Krebs, Keller & Tamarin (1969). 441
442
W.
S. C. GURNEY
AND
R.
M.
NISBET
In this paper we examine, on an ecological time scale the mechanisms whereby a population in a spatially variable environment might regulate its numbers through dispersal. The problem we shall consider is that of a single mobile species existing in an area within which large variations in the local supply of some critical resource? produce large differences in the local birth (B(r)) and death (4(r)) rates of the species. Since we seek mechanisms by which the population can regulate its numbers well below the “nutrient defined carrying capacity” we shall assume that its activities cannot affect the supply of the critical resource at any locality sufficiently to produce a significant change in either the local birth or death rates. For simplicity, we shall assume that both p(r) and 4(r) are wholly independent of population density p(r). The effect on our population of the totality of relevant abiotic factors is thus described by a single growth function G(r) = B(r) - 444 which depends only on position. We now ask how p(r), the population density at some point r, changes with time. Animals can enter the population at r only by birth or immigration and can leave it only by death or emigration. Keeping this in mind, a simple piece of “book-keeping” is sufficient to show that if the dispersal processes which are the subject of this paper produce a local population current density j(r, t), then: -$ (r, t) = G(r)p(r, t) - V .j(r, I).
2. Models
To solve equation (1) for the time development of our population we need a detailed description of the transport behaviour of the species concerned, so as to provide an explicit form for the current density j(r, r). In this paper we examine three quite distinct microscopic models. (A) THE RANDOM MOTION (RM) MODEL
The simplest assumption one can make is that the individual members of the population move entirely at random, with a mean free path which is small compared to the natural distance scale of the problem. The transport behaviour of the population as a whole will then be the simple linear diffusion normally observed in physical systems, and we will have: j = -DVp, (2) t We include refugefrom predation within our definition of a critical resource.
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
which, together with equation (I), tells us immediately
that:
ap
at = G(r)p +DV’p. Although the assumption of a randomly moving, non-interacting population leading to linear diffusional behaviour is very common in work on inhomogeneous populations,t we believe that it is only plausible in a small number of special cases (such as an algal population in turbulent water where random motion is impressed upon the population). In general we believe that interactions between different individuals, and between single animals and external factors will render the behaviour more complex than that represented by equation (2). (B)
THE
BIASED
RANDOM
MOTION
(BRM)
MODEL
If we imagine the movements of individuals within the population to be largely random but modified to a small extent by the gross distribution of their fellows, then we have a model which seems appropriate to a nonterritorial species, whose members travel continuously in search of food and interact only weakly with each other. In Appendix A we consider in detail the problem of a population whose members walk pseudo-randomly through a rectangular grid with the directional probability distribution for each step distorted a little by the local population density gradient.$ We show that provided the step length is small, the population current density is: j = -DVp-ApVp, (4) where rZis a positive constant depending upon the proportionality between bias and density gradient. It follows that in the case of biased random motion our population will develop according to: ;$ = G(r)p+DV’p+IzV.(pVp). (C)
DIRECTED
MOTION
If a species is strongly territorial, then altRough in the field we observe each animal in apparently “random” motion, such motion plays no part in population transport since each animal regularly returns to its starting point. The movements which are of interest to us in the present context are movements of territorial centre. Such movements will in general be made either t See, for example, Fisher (1937), Skellam (1991), Bradford & Philip (1970a,b), Yoshizawa (1970), Segal &Jackson (1972), Gurtin (1973). $ It is possible to construct alternative microscopic models for BRM and DM in which the frequency of the random moves depends on the local population density-see Gurney & Nisbet (1975).
444
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
by mature animals driven out by invaders or by young animals just reaching maturity moving out of their parental territory to establish breeding territory of their own. In both cases, rather than the movements being regarded as in any sense random, it is much more plausible to suppose that they will be directed towards nearby vacant territory. In this model, therefore, movement will take place almost exclusively “down” the population density gradient, and will be much more rapid at high population densities than at low ones. In an attempt to model this situation we consider in Appendix A a walk through a rectangular grid, in which at each “step” an animal may either stay at its present location or may move in the direction of lowest population density, the probability distribution among these two possibilities being determined by the magnitude of the population density gradient at the grid site concerned. We show that in this case: j
= -KpVp,
(6)
and hence that
ap at = G(r)p + KV .(pVp). The microscopic model which we have used to arrive at equation (7) is clearly not a good representation of the details of the “spacing” or territorial behaviour described by Watson & Moss (1970) and Wynne-Edwards (1962) in the context of population regulation. We have been unable to evolve a more realistic dynamic model which retains the elegant simplicity and tractability of equation (7). However, despite its lack of detailed realism we believe equation (7) reproduces many of the consequences of territoriality which are important in the context of population regulation. (D)
NORMALIZATION
As we are concerned with a one-species system, we are flee to regard the “natural” scales of our problem as properties of that species. By doing so we immediately turn many of the apparent differences between species with the same type of transport behaviour into simple differences in scale. An elegant demonstration of the power of this technique for systematizing apparently disparate experimental data is given by Gold (1973). Before proceeding with the analysis of our three models it is desirable to eliminate from them all “scaling” differences and leave only genuine qualitative differences in behaviour. We do this by making explicit our choice of scale (po, re, to) for each variable and writing: P’ = PIP09
r’
= r/ro,
t’ =
t/t,.
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
In terms of our new scaled variables, equation becomes :
445
(3) for the RM model
ag=(t,G(r))p’ +rJ$] v’2p’ where V’ is the operator
Provided that we restrict ourselves to growth functions with a single, well-defined positive maximum value, we can allow the growth function to determine the time scale by demanding that: toWmax = 1, and writing t,G(r) = G'(r'j G 1. We now let the diffusion coefficient set our distance scale by putting rf = Dt 0’
This choice of scales turns equation (3) into: -w = G’(r’)p’ + V2p’.
at’
Analysis of the properties of this equation as a function of the various possible shapes of the normalized growth function G’ will now answer all our questions about the behaviour of the RM model. Performing exactly the same process on equation (5) causes us to choose: toG(r),,,
= 1,
r-0”= Dt,,
p. = A,-‘&
thus reducing our BRM model to
w=
z
G’(r’)p’ + V”p’ + V’ . (p’v’p’).
In the case of our third (DM) model we choose: toWmax = 1, to reduce equation (7) to:
w = G’(r’)p’
z
6 = Kto po,
+V’ .(p’v’p’).
(11)
446
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
This normalization process has left the scale of population density unrestrained in the case of the RM and DM models, whereas in the BRM model the presence of a second adjustable parameter (A) has forced us also to choose the population scale if we are to achieve the maximally simple representation. (E)
THE
GROWTH
FUNCTION
We shall restrict our attention to a largely hostile (G’(r’) < 0) “universe” containing a single region of viable habitat (G’(r’) > 0). If the total population of our “universe” is to remain finite we must insist that p’(r’) + 0 as Ir’l becomes very large. To this end we demand that the normalized growth function satisfies the condition: G’(r’) + - 00 as lr’l + co. The qualitative parts of our analysis have been carried out for growth functions of the general form shown in Fig. 1 and satisfying the above
-4L
FIG.
’ -6
I -4
I / 1 I -2 0 2 4 Normalized distance (x/x0)
! 6-
1. A typical normalized growth function.
condition. In order to carry out any quantitative analysis we must provide an explicit form for G’(r’). Our computer studies indicate that if the growth function has the form of Fig. 1, the main changes in behaviour occur as the region of net growth (NGR) changes width, and that the detailed shape of the curve is of only secondary importance. We therefore limit discussion to growth functions of the form: G’(r’) = 1 --Ad”.
(12)
REGULATION
OF
INHOMOGENEOUS
3. Quditative
POPULATIONS
447
Analysis
The development with time of a population obeying equation (9), (10) or (11) is uniquely determined for a given growth function if we specify the normalized population p’(r’, tb) at some initial time tb. However, only for one of our models [the RM model defined by equation (9) which is linear in p’] were we able to obtain analytic solutions. For the other non-linear models, explicit solutions of the initial value problem had to be obtained numerically, but it was possible in addition to analyse rigorously certain qualitatioe properties of these solutions (i.e. properties such as existence, positivity and stability). This qualitative analysis is now given; the reader uninterested in the mathematical aspects may jump to the final paragraph of this section where we have summarized our principal qualitative results. (A)
THE
RANDOM
MOTION
MODEL
As we work with scaled variables throughout the rest of section 3, in the interests of clarity we shall drop the primes used to denote scaled quantities. The scaled RM model is now defined by: m
$ = G(r)p +I+,
(9)
which, apart from a factor i on the left-hand side, has the form of the Schriidinger equation of quantum mechanics (see e.g. Schiff, 1955). Guided by this similarity, we define an operator A(r) by A(r) = G(r) + V2. (13) Provided G(r) + - 00 as lrl + co, this operator has a discrete set of eigenfunctions {Ui(r)}, with a corresponding set of eigenvalues {A,}. The solution of equation (9) is then p(r, t) = C Ci e”” ul(r) (14) where the constants Ci are simply the coefficients in the expansion of the initial distribution in terms of the eigenfunctions ul(r). Provided G(r) is of the form shown in Fig. 1, the maximum eigenvalue IzO is non-degenerate and the corresponding eigenfunction uO(r) may be chosen to be everywhere positive. It is then easy to show that the coefficient CO is non-zero, and hence that after any transient behaviour has died out both the local population density and the total population ultimately grow or decay exponential&f t There is, of course, the pathological case I = 0 which leads ultimately to a stationary solution. We regard this solution as of little biological interest since it is structurally unstable in the sense that an infinitesimal change in any coefficient causes a qualitative change in the solution.
448
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
(depending on whether I, is positive or negative). Moreover, for many simple G(r), the eigenvalues of A(r) are readily accessible in the quantum mechanics literature. It is thus frequently possible to write down explicitly the rate of this exponential growth or decay. It is the form of the normalized growth function G(r) which determines the sign of &, and hence the fate of the population. For future reference we call a growth function which leads to an exponentially increasing population in the RM model a “random motion increaser” (RMI) and analogously define a “random motion decrease? (RMD). (B)
THE
BIASED
RANDOM
MOTION
MODEL
After scaling, this model is defined by equation (10) g = G(r)p + V’p + V . (pVp).
(10)
If G(r) is a RMD growth function, it is not difficult to demonstrate that the only stationary solution of this equation is the trivial solution p(r) = 0 for all r. In Appendix B we show that if G(r) is a RMI growth function, there exists in addition a non-trivial, non-negative, stationary solution p,(r). We can prove that this solution p,(r) is stable by examining the behaviour of small fluctuations in the population density. If we set (15) p(r) = d9 + 44 with /e(r)] 6 [P,(r)1 for all r, a sufficient condition for the stability of p,(r) is that the mean square fluctuation decreases with time, a result which implies that?
zj(s(r))2 dr<0
(16)
with equality only if s(r) = 0 for all r. In practice, when a system is perturbed, it is often found that the condition (16) is only satisfied after a finite period of transient behaviour. In this cast, the ultimate decay of the fluctuations is guaranteed provided there exists a strictly positive weighting function w(r) such that i j w(r>[&(r>]’ dr < 0 with equality only when a(r) = 0 for all r. f The integral in equation (16) is over the entire “universe” we are considering. As our subsequent analysis is very general it is unnecessary to distinguish between Z- and 3-dimensional systems.
REGULATION
OF INHOMOGENEOUS
POPULATIONS
449
To prove the stability of the stationary solution p,(r) of equation (IO), we first note that for the fluctuations c(r) defined in equation (15), at = G(r)& + pJ% + 2Vp,V& + &V2p, + V’E, to first order in e. We then note that as long as p,(r) is strictly positive, G(r) may be calculated in terms of p,(r) from the steady state condition Wh
+ V’P, + V .(PJP~ = 0.
(19)
Substitution from (19) into (18) now gives as/& in terms of ps, e, and their derivatives. To complete the proof of the stability of p,(r) with respect to small fluctuations we must choose the weighting function. While there is no well-defined recipe for obtaining a weighting function suitable for a particular problem, Glansdoe & Prigogine (1971) give a number of examples of suitable weighting functions for thermodynamic systemls with non-linear constitutive equations. Guided by these examples we choose w(r) = 1 +p,(r).
(20)
It is now a matter of tedious and rather uninstructive algebra to show that provided e(r) is everywhere sufficiently small for us to use (18), and provided c(r), p,(r), and their gradients vanish on the boundaries of the “universe”, then?
- j E’(V~~)’ dr < 0 (21) with equality only when c(r) = 0 for all r. This completes the proof of the stability of the steady state solution p,(r). (C) THE DIRECTED MOTION MODEL
The analysis of this model is very simple once the case of BRM is understood. In Appendix B we prove that for any un-normalized growth function which is positive in some finite region, there exists a non-trivial, non-negative, stationary solution of equation (11). The stability of the solution is again proved by considering small fluctuations c(r); the ensuing algebra is very similar to that for BRM except that we choose w(r) = p,(r).
(22)
t The key to deriving (21) is the result that if a(r), b(r) are scalar fields which vanish on the boundary of a volume V then 6 nVZbdr = - j Va. Vbdr.
450
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
The final result is that f $ j p,e2 dr = - j [E’(V~~)* +p,z(V~)~] dr < 0, with equality only when s(r) = 0 for all r. This proves that the stationary state is stable. (D)
SUMMARY
OF THE
QUALITATIVE
ANALYSIS
We have proved that in the RM model, the local population density and the total population eventually grow or decay exponentially. We call a normalized growth function which leads to an exponentially increasing population a “random motion increaser” (RMI) and analogously define a “random motion decrease? (RMD). For the BRM model we show that with a RMD growth function the population proceeds exponentially to extinction. With a RMI growth function there is however a non-trivial, non-negative steady-state solution which is stable with respect to small fluctuations in population, For the DM model, there is a non-trivial, non-negative, steady-state solution provided the growth function is positive in some finite region and negative elsewhere. This steady state is also stable under small fluctuations in the population. 4. Numerical
Analysis
The qualitative analysis described in section 3 is able to tell us the conditions under which the equations describing our three models have stable, non-trivial, stationary solutions, but it can tell us nothing about the transient behaviour of our system. Such transient behaviour is of clear interest since stable solutions to equations (9)-(11) are of no biological interest unless they are reached within a reasonably short space of time. In our normalized equations the natural unit of time (to) is directly related to the reproductive rate of the species concerned, and we thus simply ask if the stable state is reached within a reasonable number of “natural time units”. To gain at least a partial answer to this question we investigated in detail the behaviour of a number of special cases of our three models, by using a simple predictor-corrector technique (Kelly, 1967) to solve a one-dimensional version of the relevant equation C(9), (10) or (1 1)] numerically. A typical set of results is reproduced in Fig. 2, which shows the time-development of the total number (n) of individuals? in populations of similar initial size t For computational convenience we have normalized the total population n by a factor n, = pox0 and it is the normalized total population (n/n,) which is displayed in Figs 2 and 3.
REGULATION
OF INHOMOGENEOUS
451
POPULATIONS
operating in two different “universes”; one described by a RMI growth function, and the other by a RMD growth function. In the RMD “universe” [Fig. 2(a)] we see that after an initial transient, populations described by both RM and BRM models decrease exponentially with a time constant of about 2*5t,,t while that described by the DM model takes a total of about Sr, to relax to its stable level from a population initially five times as large. In the RMI “universe” [Fig. 2(b)] the randomly moving population increases exponentially with a time constant -2t0, while both those which show some
Directed
motion Directed
Normalized
motion
time (t/t,)
RIO. 2. Time development of model populations “universe”.
in (a) a RMD
“universe”,
(b) a RMI
kind of density-dependent motion achieve stable states after about lo&,. In these special cases the detailed solution of equations (9)-(11) has yielded results whose qualitative features agree well with the general results in section 3(~). It is clear that in those cases where a stable state exists, it is reached, even for a grossly perturbed initial state, within - lot,. More detail of the behaviour of populations exhibiting BRM and DM in RMI “universes” is shown on a linear scale in Fig. 3, where we display the build-up of a population from a very low initial level, as a function of the width of the region (NGR) in which G’(C) is positive. As one would expect, t We have checked the rate of population decay shown by our numerical solutions of equation (9) against the maximum eigenvalue of the operator defined by equation (13) with the appropriate growth function. We find the two in good agreement. T.B. 29
452
..,
n
w.
2 EL 3
#.
b.
L.
^..-.l-.l
1
.I-..
-
ANIJ
CIUKNEY
.~
.-
K.
.
.
M.
NGR
width
5.7
NGR
width
4.0
-.
.._^---
NlStlhl
~~__..___
250
g 200 z
. 1
1.50 INGR
width
5.7
NGR
width
4.0
100 50
0
5
IO
15 Normalized
(b)
FIG. 3. Stabilization directed motion.
of
increasing
20
25
--30
35
- 1
time (t/r,,)
populations
exhibiting
(a)
biased
random
motion
the stable population level increases with the NGR width, but the period of transient behaviour is both relatively short and only weakly dependent on the growth function. 5. Discussion The analyses described above show that in a constant inhomogeneous environment a single species population cannot be stabilized by dispersal achieved simply by random motion. This is not surprising since postulating both density-independent birth and death rates trrld a linear dispersal process leads to a model containing no non-linearity suitable for providing control. This becomes clearer if we write the population current density of equation (1) as the product of a population density and a drift velocity j = pv. (23)
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
453
In the case of random motion we then see that the drift velocity, which determines the transport rate for any individual, depends only upon the reZative population density gradient thus : v=D-.
VP P
This makes it intuitively clear that although for any given growth function the population distribution will have a stable 4rhrzpe,its absolute size is not subject to any constraint. It is thus plausible that after an initial transient, during which this stable shape is attained, the amplitude of the population distribution will either grow or decay according to the nature of the growth function. In fact, for a particular functional form of normalized growth function the change from increasing to decreasing behaviour is simply a matter of the width of the region of net growth. In the case of the DM model it is clear from a cursory examination of equation (11) that the size of the dispersal contribution to dp/St goes to zero along with p. It is very easy to see how such a system will develop a suitable stable distribution for any growth function which has a region of net growth, no matter how small. The key to understanding the biased random motion model is to realize that the essential effect of introducing a density-dependent bias into the system is to increase the diffusion rate at high population densities. At low population densities the motion is indistinguishable from RM. Thus if the growth function is such that a randomly moving population would tend to increase, then the bias will increase the diffusion rate in densely populated regions where births exceed deaths, thus stabilizing the population. However, if the growth function is a random motion decreaser, the bias simply makes the decrease more rapid at high population levels, and as the population decays its rate of decrease relaxes back towards its “random motion” value.
6. Conclusion It has long been believed that the kind of sigmoid growth curves shown in Fig. 3 are associated with the presence of a density-dependent factor controlling the growth of population. The classic assumption of theoretical population dynamics is that this density dependence occurs in the birth or death rates of the species concerned. We have shown that the introduction of a suitable non-linearity into the dispersal behaviour of a species which behaves in an otherwise linear way can, in an inhomogeneous environment, lead to an exactly similar regulatory effect. Just such a mechanism is required
454
W.
S.
C‘.
GURNEY
AND
K.
M.
NISBkr
to explain the observation that populations of certain species appear to be able to regulate their numbers below the carrying capacity set by the supply of nutrients available to them. REFERENCES BRADFORD, E. W. & PHILIP, J. R. (19704. BRADFORD. E. W. & PHILIP. J. R. (197Ob). CARL, E. A. (1971). Eculog; 52, 39j. ’ FISHER, R. A. (1937). Ann. Eugen. 7, 355. GLANSWRFF, P. & PRIGOGINE, I. (1971).
J. theor. Biul. ZY, 13. J. theor. Biol. 29. 28.
Thermodynumic T/wry of Strncture, Stability and Fluctuations. London: Wiley-Interscience. GOLD, A, (1973). Science, N. Y. 181, 275. GURNEY, W. S. C. & NISBET, R. M. (in press). J. theor. Biul. GURTIN, M. E. (1973). J. theor. Biol. 40, 389. HOWARD, W. E. (1960). Am. Midl. Nat. 63, 152. KELLER, H. B. (1970). J. Diff. Eqs 7,417. KELLY, L. G. (1967). Handbook of Numerical Methods and Applicufiuns. Reading, Massachusetts: Addison Wesley. KREBS, C. J., KELLER, B. L. & TAMARIN, R. H. (1969). Ecology 50, 587. LIDIKER. W. Z. (1962). Am. Nat. 96.29. MONTR~LL, E. w. &WEST, B. J. (19?3). In Synergetics (H. Haken, cd.). Stuttgart: Teubner. MURRAY. B. G. (1967). Ecology 48,975. SCHIFF, L. I. (1955). Quantum Mechanics. New York: McGraw-Hill. SEGAL, L. A. & JACKSON, J. L. (1972). J. theor. Biol. 37, 545. SKELLAM, J. G. (1951). Biometrika 38, 196. STACKGOLD, I. (1971). SIAM Rev. 13, 289. WATON, A. & Moss, R. (1970). Animal Populatiun~ ill Relation ID their Food Resuarces. Oxford and Edinburgh: Blackwell Sci. Pub. WYNNE-EDWARDS, V. C. (1962). Animal Dispersiun in Relutiun to Suciul Behaviorrr. Edinburgh and London : Oiiver & Boyd. YOSHIZAWA, S. (1970). Math. Biosci. 7, 291.
APPENDIX
A
Microscope Models for Biased Random Motion and Directed Motion It is well known that in a population of random walkers the macroscopic current density is given by j = - DVp. (Al) Montroll & West (1973) recently modelled populations of “biased” random walkers who prefer to move up or down a population gradient; our present interest is in the case of walkers who dislike congestion and prefer to move down the population gradient. We imagine a “walk” in a d-dimensional rectangular grid (where d = I, 2 or 3). We let a be the distance between adjacent points in the grid and assume that during each time interval z every walker moves to one of its nearest neighbour points. If the motion were random, the probability of jumping
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
455
from a point r to a neighbouring point r +aB would be 2-‘; we bias this result by setting the probability of such a jump to P(r, ati) where ~(r, ae) = 2-d--+q.e, 642) and p > 0. With this definition the component of microscopic current density in the direction B is simply j,(r) = UT- ‘(p(r)P(r, d>-p(r + di)P(r+ ae, - ae)). (A31 If we assume that (I and r are small, and set D = 2-du2z-1 and I = 2&r-’ then we find: j,(r) = ( - DVp - ApVp) .A. (A41 It follows that the macroscopic current density is simply j = - DVp - IpVp, 645) a result in disagreement with Montroll & West (1973) who appear to have dropped some vital second-order terms in the1 derivation of their continuity equation. The concept of “directed motion” was introduced in the text to cover the important case of a territorial animal whose random motion makes no contribution to the macroscopic current density. The only choices offered to the walker are to stay put or to move down the population gradient. We choose the probability that an individual will jump to a nearest-neighbour point to be P(r, ae) = -ptavp . e H( - vp . e), WI where H(x) is the step function defined by H(x) = 1 if x > 0, 647) = 0 otherwise. The microscopic current density in direction A is then j, = 4pvp.e (fw where a is small and K = pa2z-‘, while the corresponding macroscopic current density is j = -KpVp. (A9
APPENDIX
B
The Existence of Positive, Stationary Solutions to the\ BRM and DM Equations Let 1 be a real Banach space with norm I( 11,and let A: g + W be a non-linear mapping which satisfies A(0) = 0. Consider the non-linear eigenvalue equation A(P) = PP9
PEW,
PE%
(AlO)
456
W.
S.
C.
GURNEY
AND
R.
M.
NISBE’I
which obviously has the trivial solution p = 0. A blench point of (AlO) is a number no with the property that in every neighbourhood of (PO, 0) in g x a there exists a solution of (AIO) with llyj # 0. The theory of branch points of (AlO) involves “linearizing” the operator A. If, when /i.ll is small, we can write A(p+1)-A(p)
= L,;1+om
(Al 1) then we say that A is linearizable at p and call L,, the FrPchet dericatice of A at p. We write Lo for the Frtchet derivative of .4 at p = 0. Stackgold (1971) in a lucid expository article quotes two theorems which enable us to locate the branch points of (AlO) provided A is linearizable atp = 0: (i) /lo can only be a branch point of (AIO) if it is an eigenvalue of L,, (ii) if A is completely continuous and p” ( # 0) is an eigenvaluc of Lo with odd multiplicity, then /I’ is a branch point of (AIO). To apply these theorems to our models it is convenient to start with the mnormafized BRM equation (5), a steady-state solution of which must satisfy the non-linear eigenvalue equation 0)
= PP?
(A121
where ,4(p) = - 1V. (pVp) - DV’p + &r)p.
(Al3)
The Frcchet derivative of ! at p = 0 is Lo = - DV’ +4(r),
(Al4)
which is essentially the operator we used to analyse the RM model. When the growth function is of the form shown in Fig. I, all the eigenvalues of I., are positive and there is a simple minimum eigenvalue ~1’. When j3 < /lo the only non-negative solution of (A12) is the trivial solution p = 0. To prove this we choose as norm in J IIPII = S p2 dr(Al? and use the original differential equation (5) to establish that for all strictly positive functions p
As the condition /j < no is equivalent to the condition that the normalized growth function be a RMD, we have proved that there is no positive stationary solution with a RMD growth function. However, the two theorems quoted above guarantee that the point /I = p” is a branch point of (A12). A simple perturbation analysis about the branch point reveals that the “new”
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
457
solution of (A12) which emanates from /3 = p” is initially (i.e. when /I-p0 is small and positive) proportional to the eigenfunction corresponding to PO. For growth functions of the form of Fig. 1 we know that this eigenfunction is of fixed sign and in fact we may take it to be positive (see Keller, 1970), for analysis of this point in a similar example). It remains to be shown that this solution persists for all /3 > PO. This is done by considering functions of very small and very large norm. It is easy to show that if p = au with llull = 1, then for all positive functions U, there exists an a0 (dependent on u) such that $ ilp/ < 0
for a > cto.
(Al7)
It is also easy to show that for all positive u, there exists an a, such that after some initial transient T(U)
; llpll ’ 0
acal,
t> T(u).
GW
As the inequalities (A17) and (A19) hold for dip, the existence of a stationary solution is guaranteed, although of course its stability and the possibility of periodic behaviour remains to be investigated. This was done in the text. We can prove the existence of a positive stationary solution to the DM equation by exactly the same strategy as was used in the previous paragraph, i.e. by considering functions of large and small norm.
The Regulation of Inhomogeneous Populations W. S. C. GURNEY AND R. M. NISBET Department of Applied Physics, University of Strathclyde, Glasgow G4 ONG, Scotland (Received 29 May 1974) A number of authors have argued that dispersal may play a key role in the regulation of populations of some species. In this paper we examine, on an ecological time scale, mechanisms whereby a population existing
in an inhomogeneous environment may use dispersal to regulate its sire below the carrying capacity set by the supply of available nutrients. We show that dispersal produced by wholly random motion is incapable of exerting any stabilizing influence, but that the introduction of a suitable non-linearity into the dispersal behaviour of a species whose characteristics are otherwise wholly linear can lead to stabilization under a wide range of conditions.
1. Introduction A number of author (e.g. Lidiker, 1962; Howard, 1960; Murray, 1967) have argued that emigration or dispersal may play a key role in the regulation of populations of some species. A persuasive example of this suggestion occurs in a paper by Carl (1971), whose observations on a population of arctic ground squirrels showed that this species achieves population control by the dispersal of animals from densely populated areas of favourable habitat into unfavourable areas where burrow sites are not available, and where they are subject to intensive predation. It seems plausible that adaptative advantage might accrue to a species with a stable strategy for limiting its numbers below the maximum carrying capacity of its environment, since it is thereby provided with considerable protection against fluctuations in the supply of available nutrients. Lidiker (1962) has argued that such a strategy might well operate through dispersal, and would then bring with it the additional advantages of increased heterozygosity and more probable genetic recombination. That dispersal can indeed operate to control a population at a level below the nominal carrying capacity of its locality has recently been most elegantly demonstrated for several vole species by Krebs, Keller & Tamarin (1969). 441
442
W.
S. C. GURNEY
AND
R.
M.
NISBET
In this paper we examine, on an ecological time scale the mechanisms whereby a population in a spatially variable environment might regulate its numbers through dispersal. The problem we shall consider is that of a single mobile species existing in an area within which large variations in the local supply of some critical resource? produce large differences in the local birth (B(r)) and death (4(r)) rates of the species. Since we seek mechanisms by which the population can regulate its numbers well below the “nutrient defined carrying capacity” we shall assume that its activities cannot affect the supply of the critical resource at any locality sufficiently to produce a significant change in either the local birth or death rates. For simplicity, we shall assume that both p(r) and 4(r) are wholly independent of population density p(r). The effect on our population of the totality of relevant abiotic factors is thus described by a single growth function G(r) = B(r) - 444 which depends only on position. We now ask how p(r), the population density at some point r, changes with time. Animals can enter the population at r only by birth or immigration and can leave it only by death or emigration. Keeping this in mind, a simple piece of “book-keeping” is sufficient to show that if the dispersal processes which are the subject of this paper produce a local population current density j(r, t), then: -$ (r, t) = G(r)p(r, t) - V .j(r, I).
2. Models
To solve equation (1) for the time development of our population we need a detailed description of the transport behaviour of the species concerned, so as to provide an explicit form for the current density j(r, r). In this paper we examine three quite distinct microscopic models. (A) THE RANDOM MOTION (RM) MODEL
The simplest assumption one can make is that the individual members of the population move entirely at random, with a mean free path which is small compared to the natural distance scale of the problem. The transport behaviour of the population as a whole will then be the simple linear diffusion normally observed in physical systems, and we will have: j = -DVp, (2) t We include refugefrom predation within our definition of a critical resource.
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
which, together with equation (I), tells us immediately
that:
ap
at = G(r)p +DV’p. Although the assumption of a randomly moving, non-interacting population leading to linear diffusional behaviour is very common in work on inhomogeneous populations,t we believe that it is only plausible in a small number of special cases (such as an algal population in turbulent water where random motion is impressed upon the population). In general we believe that interactions between different individuals, and between single animals and external factors will render the behaviour more complex than that represented by equation (2). (B)
THE
BIASED
RANDOM
MOTION
(BRM)
MODEL
If we imagine the movements of individuals within the population to be largely random but modified to a small extent by the gross distribution of their fellows, then we have a model which seems appropriate to a nonterritorial species, whose members travel continuously in search of food and interact only weakly with each other. In Appendix A we consider in detail the problem of a population whose members walk pseudo-randomly through a rectangular grid with the directional probability distribution for each step distorted a little by the local population density gradient.$ We show that provided the step length is small, the population current density is: j = -DVp-ApVp, (4) where rZis a positive constant depending upon the proportionality between bias and density gradient. It follows that in the case of biased random motion our population will develop according to: ;$ = G(r)p+DV’p+IzV.(pVp). (C)
DIRECTED
MOTION
If a species is strongly territorial, then altRough in the field we observe each animal in apparently “random” motion, such motion plays no part in population transport since each animal regularly returns to its starting point. The movements which are of interest to us in the present context are movements of territorial centre. Such movements will in general be made either t See, for example, Fisher (1937), Skellam (1991), Bradford & Philip (1970a,b), Yoshizawa (1970), Segal &Jackson (1972), Gurtin (1973). $ It is possible to construct alternative microscopic models for BRM and DM in which the frequency of the random moves depends on the local population density-see Gurney & Nisbet (1975).
444
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
by mature animals driven out by invaders or by young animals just reaching maturity moving out of their parental territory to establish breeding territory of their own. In both cases, rather than the movements being regarded as in any sense random, it is much more plausible to suppose that they will be directed towards nearby vacant territory. In this model, therefore, movement will take place almost exclusively “down” the population density gradient, and will be much more rapid at high population densities than at low ones. In an attempt to model this situation we consider in Appendix A a walk through a rectangular grid, in which at each “step” an animal may either stay at its present location or may move in the direction of lowest population density, the probability distribution among these two possibilities being determined by the magnitude of the population density gradient at the grid site concerned. We show that in this case: j
= -KpVp,
(6)
and hence that
ap at = G(r)p + KV .(pVp). The microscopic model which we have used to arrive at equation (7) is clearly not a good representation of the details of the “spacing” or territorial behaviour described by Watson & Moss (1970) and Wynne-Edwards (1962) in the context of population regulation. We have been unable to evolve a more realistic dynamic model which retains the elegant simplicity and tractability of equation (7). However, despite its lack of detailed realism we believe equation (7) reproduces many of the consequences of territoriality which are important in the context of population regulation. (D)
NORMALIZATION
As we are concerned with a one-species system, we are flee to regard the “natural” scales of our problem as properties of that species. By doing so we immediately turn many of the apparent differences between species with the same type of transport behaviour into simple differences in scale. An elegant demonstration of the power of this technique for systematizing apparently disparate experimental data is given by Gold (1973). Before proceeding with the analysis of our three models it is desirable to eliminate from them all “scaling” differences and leave only genuine qualitative differences in behaviour. We do this by making explicit our choice of scale (po, re, to) for each variable and writing: P’ = PIP09
r’
= r/ro,
t’ =
t/t,.
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
In terms of our new scaled variables, equation becomes :
445
(3) for the RM model
ag=(t,G(r))p’ +rJ$] v’2p’ where V’ is the operator
Provided that we restrict ourselves to growth functions with a single, well-defined positive maximum value, we can allow the growth function to determine the time scale by demanding that: toWmax = 1, and writing t,G(r) = G'(r'j G 1. We now let the diffusion coefficient set our distance scale by putting rf = Dt 0’
This choice of scales turns equation (3) into: -w = G’(r’)p’ + V2p’.
at’
Analysis of the properties of this equation as a function of the various possible shapes of the normalized growth function G’ will now answer all our questions about the behaviour of the RM model. Performing exactly the same process on equation (5) causes us to choose: toG(r),,,
= 1,
r-0”= Dt,,
p. = A,-‘&
thus reducing our BRM model to
w=
z
G’(r’)p’ + V”p’ + V’ . (p’v’p’).
In the case of our third (DM) model we choose: toWmax = 1, to reduce equation (7) to:
w = G’(r’)p’
z
6 = Kto po,
+V’ .(p’v’p’).
(11)
446
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
This normalization process has left the scale of population density unrestrained in the case of the RM and DM models, whereas in the BRM model the presence of a second adjustable parameter (A) has forced us also to choose the population scale if we are to achieve the maximally simple representation. (E)
THE
GROWTH
FUNCTION
We shall restrict our attention to a largely hostile (G’(r’) < 0) “universe” containing a single region of viable habitat (G’(r’) > 0). If the total population of our “universe” is to remain finite we must insist that p’(r’) + 0 as Ir’l becomes very large. To this end we demand that the normalized growth function satisfies the condition: G’(r’) + - 00 as lr’l + co. The qualitative parts of our analysis have been carried out for growth functions of the general form shown in Fig. 1 and satisfying the above
-4L
FIG.
’ -6
I -4
I / 1 I -2 0 2 4 Normalized distance (x/x0)
! 6-
1. A typical normalized growth function.
condition. In order to carry out any quantitative analysis we must provide an explicit form for G’(r’). Our computer studies indicate that if the growth function has the form of Fig. 1, the main changes in behaviour occur as the region of net growth (NGR) changes width, and that the detailed shape of the curve is of only secondary importance. We therefore limit discussion to growth functions of the form: G’(r’) = 1 --Ad”.
(12)
REGULATION
OF
INHOMOGENEOUS
3. Quditative
POPULATIONS
447
Analysis
The development with time of a population obeying equation (9), (10) or (11) is uniquely determined for a given growth function if we specify the normalized population p’(r’, tb) at some initial time tb. However, only for one of our models [the RM model defined by equation (9) which is linear in p’] were we able to obtain analytic solutions. For the other non-linear models, explicit solutions of the initial value problem had to be obtained numerically, but it was possible in addition to analyse rigorously certain qualitatioe properties of these solutions (i.e. properties such as existence, positivity and stability). This qualitative analysis is now given; the reader uninterested in the mathematical aspects may jump to the final paragraph of this section where we have summarized our principal qualitative results. (A)
THE
RANDOM
MOTION
MODEL
As we work with scaled variables throughout the rest of section 3, in the interests of clarity we shall drop the primes used to denote scaled quantities. The scaled RM model is now defined by: m
$ = G(r)p +I+,
(9)
which, apart from a factor i on the left-hand side, has the form of the Schriidinger equation of quantum mechanics (see e.g. Schiff, 1955). Guided by this similarity, we define an operator A(r) by A(r) = G(r) + V2. (13) Provided G(r) + - 00 as lrl + co, this operator has a discrete set of eigenfunctions {Ui(r)}, with a corresponding set of eigenvalues {A,}. The solution of equation (9) is then p(r, t) = C Ci e”” ul(r) (14) where the constants Ci are simply the coefficients in the expansion of the initial distribution in terms of the eigenfunctions ul(r). Provided G(r) is of the form shown in Fig. 1, the maximum eigenvalue IzO is non-degenerate and the corresponding eigenfunction uO(r) may be chosen to be everywhere positive. It is then easy to show that the coefficient CO is non-zero, and hence that after any transient behaviour has died out both the local population density and the total population ultimately grow or decay exponential&f t There is, of course, the pathological case I = 0 which leads ultimately to a stationary solution. We regard this solution as of little biological interest since it is structurally unstable in the sense that an infinitesimal change in any coefficient causes a qualitative change in the solution.
448
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
(depending on whether I, is positive or negative). Moreover, for many simple G(r), the eigenvalues of A(r) are readily accessible in the quantum mechanics literature. It is thus frequently possible to write down explicitly the rate of this exponential growth or decay. It is the form of the normalized growth function G(r) which determines the sign of &, and hence the fate of the population. For future reference we call a growth function which leads to an exponentially increasing population in the RM model a “random motion increaser” (RMI) and analogously define a “random motion decrease? (RMD). (B)
THE
BIASED
RANDOM
MOTION
MODEL
After scaling, this model is defined by equation (10) g = G(r)p + V’p + V . (pVp).
(10)
If G(r) is a RMD growth function, it is not difficult to demonstrate that the only stationary solution of this equation is the trivial solution p(r) = 0 for all r. In Appendix B we show that if G(r) is a RMI growth function, there exists in addition a non-trivial, non-negative, stationary solution p,(r). We can prove that this solution p,(r) is stable by examining the behaviour of small fluctuations in the population density. If we set (15) p(r) = d9 + 44 with /e(r)] 6 [P,(r)1 for all r, a sufficient condition for the stability of p,(r) is that the mean square fluctuation decreases with time, a result which implies that?
zj(s(r))2 dr<0
(16)
with equality only if s(r) = 0 for all r. In practice, when a system is perturbed, it is often found that the condition (16) is only satisfied after a finite period of transient behaviour. In this cast, the ultimate decay of the fluctuations is guaranteed provided there exists a strictly positive weighting function w(r) such that i j w(r>[&(r>]’ dr < 0 with equality only when a(r) = 0 for all r. f The integral in equation (16) is over the entire “universe” we are considering. As our subsequent analysis is very general it is unnecessary to distinguish between Z- and 3-dimensional systems.
REGULATION
OF INHOMOGENEOUS
POPULATIONS
449
To prove the stability of the stationary solution p,(r) of equation (IO), we first note that for the fluctuations c(r) defined in equation (15), at = G(r)& + pJ% + 2Vp,V& + &V2p, + V’E, to first order in e. We then note that as long as p,(r) is strictly positive, G(r) may be calculated in terms of p,(r) from the steady state condition Wh
+ V’P, + V .(PJP~ = 0.
(19)
Substitution from (19) into (18) now gives as/& in terms of ps, e, and their derivatives. To complete the proof of the stability of p,(r) with respect to small fluctuations we must choose the weighting function. While there is no well-defined recipe for obtaining a weighting function suitable for a particular problem, Glansdoe & Prigogine (1971) give a number of examples of suitable weighting functions for thermodynamic systemls with non-linear constitutive equations. Guided by these examples we choose w(r) = 1 +p,(r).
(20)
It is now a matter of tedious and rather uninstructive algebra to show that provided e(r) is everywhere sufficiently small for us to use (18), and provided c(r), p,(r), and their gradients vanish on the boundaries of the “universe”, then?
- j E’(V~~)’ dr < 0 (21) with equality only when c(r) = 0 for all r. This completes the proof of the stability of the steady state solution p,(r). (C) THE DIRECTED MOTION MODEL
The analysis of this model is very simple once the case of BRM is understood. In Appendix B we prove that for any un-normalized growth function which is positive in some finite region, there exists a non-trivial, non-negative, stationary solution of equation (11). The stability of the solution is again proved by considering small fluctuations c(r); the ensuing algebra is very similar to that for BRM except that we choose w(r) = p,(r).
(22)
t The key to deriving (21) is the result that if a(r), b(r) are scalar fields which vanish on the boundary of a volume V then 6 nVZbdr = - j Va. Vbdr.
450
W.
S.
C.
GURNEY
AND
R.
M.
NISBET
The final result is that f $ j p,e2 dr = - j [E’(V~~)* +p,z(V~)~] dr < 0, with equality only when s(r) = 0 for all r. This proves that the stationary state is stable. (D)
SUMMARY
OF THE
QUALITATIVE
ANALYSIS
We have proved that in the RM model, the local population density and the total population eventually grow or decay exponentially. We call a normalized growth function which leads to an exponentially increasing population a “random motion increaser” (RMI) and analogously define a “random motion decrease? (RMD). For the BRM model we show that with a RMD growth function the population proceeds exponentially to extinction. With a RMI growth function there is however a non-trivial, non-negative steady-state solution which is stable with respect to small fluctuations in population, For the DM model, there is a non-trivial, non-negative, steady-state solution provided the growth function is positive in some finite region and negative elsewhere. This steady state is also stable under small fluctuations in the population. 4. Numerical
Analysis
The qualitative analysis described in section 3 is able to tell us the conditions under which the equations describing our three models have stable, non-trivial, stationary solutions, but it can tell us nothing about the transient behaviour of our system. Such transient behaviour is of clear interest since stable solutions to equations (9)-(11) are of no biological interest unless they are reached within a reasonably short space of time. In our normalized equations the natural unit of time (to) is directly related to the reproductive rate of the species concerned, and we thus simply ask if the stable state is reached within a reasonable number of “natural time units”. To gain at least a partial answer to this question we investigated in detail the behaviour of a number of special cases of our three models, by using a simple predictor-corrector technique (Kelly, 1967) to solve a one-dimensional version of the relevant equation C(9), (10) or (1 1)] numerically. A typical set of results is reproduced in Fig. 2, which shows the time-development of the total number (n) of individuals? in populations of similar initial size t For computational convenience we have normalized the total population n by a factor n, = pox0 and it is the normalized total population (n/n,) which is displayed in Figs 2 and 3.
REGULATION
OF INHOMOGENEOUS
451
POPULATIONS
operating in two different “universes”; one described by a RMI growth function, and the other by a RMD growth function. In the RMD “universe” [Fig. 2(a)] we see that after an initial transient, populations described by both RM and BRM models decrease exponentially with a time constant of about 2*5t,,t while that described by the DM model takes a total of about Sr, to relax to its stable level from a population initially five times as large. In the RMI “universe” [Fig. 2(b)] the randomly moving population increases exponentially with a time constant -2t0, while both those which show some
Directed
motion Directed
Normalized
motion
time (t/t,)
RIO. 2. Time development of model populations “universe”.
in (a) a RMD
“universe”,
(b) a RMI
kind of density-dependent motion achieve stable states after about lo&,. In these special cases the detailed solution of equations (9)-(11) has yielded results whose qualitative features agree well with the general results in section 3(~). It is clear that in those cases where a stable state exists, it is reached, even for a grossly perturbed initial state, within - lot,. More detail of the behaviour of populations exhibiting BRM and DM in RMI “universes” is shown on a linear scale in Fig. 3, where we display the build-up of a population from a very low initial level, as a function of the width of the region (NGR) in which G’(C) is positive. As one would expect, t We have checked the rate of population decay shown by our numerical solutions of equation (9) against the maximum eigenvalue of the operator defined by equation (13) with the appropriate growth function. We find the two in good agreement. T.B. 29
452
..,
n
w.
2 EL 3
#.
b.
L.
^..-.l-.l
1
.I-..
-
ANIJ
CIUKNEY
.~
.-
K.
.
.
M.
NGR
width
5.7
NGR
width
4.0
-.
.._^---
NlStlhl
~~__..___
250
g 200 z
. 1
1.50 INGR
width
5.7
NGR
width
4.0
100 50
0
5
IO
15 Normalized
(b)
FIG. 3. Stabilization directed motion.
of
increasing
20
25
--30
35
- 1
time (t/r,,)
populations
exhibiting
(a)
biased
random
motion
the stable population level increases with the NGR width, but the period of transient behaviour is both relatively short and only weakly dependent on the growth function. 5. Discussion The analyses described above show that in a constant inhomogeneous environment a single species population cannot be stabilized by dispersal achieved simply by random motion. This is not surprising since postulating both density-independent birth and death rates trrld a linear dispersal process leads to a model containing no non-linearity suitable for providing control. This becomes clearer if we write the population current density of equation (1) as the product of a population density and a drift velocity j = pv. (23)
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
453
In the case of random motion we then see that the drift velocity, which determines the transport rate for any individual, depends only upon the reZative population density gradient thus : v=D-.
VP P
This makes it intuitively clear that although for any given growth function the population distribution will have a stable 4rhrzpe,its absolute size is not subject to any constraint. It is thus plausible that after an initial transient, during which this stable shape is attained, the amplitude of the population distribution will either grow or decay according to the nature of the growth function. In fact, for a particular functional form of normalized growth function the change from increasing to decreasing behaviour is simply a matter of the width of the region of net growth. In the case of the DM model it is clear from a cursory examination of equation (11) that the size of the dispersal contribution to dp/St goes to zero along with p. It is very easy to see how such a system will develop a suitable stable distribution for any growth function which has a region of net growth, no matter how small. The key to understanding the biased random motion model is to realize that the essential effect of introducing a density-dependent bias into the system is to increase the diffusion rate at high population densities. At low population densities the motion is indistinguishable from RM. Thus if the growth function is such that a randomly moving population would tend to increase, then the bias will increase the diffusion rate in densely populated regions where births exceed deaths, thus stabilizing the population. However, if the growth function is a random motion decreaser, the bias simply makes the decrease more rapid at high population levels, and as the population decays its rate of decrease relaxes back towards its “random motion” value.
6. Conclusion It has long been believed that the kind of sigmoid growth curves shown in Fig. 3 are associated with the presence of a density-dependent factor controlling the growth of population. The classic assumption of theoretical population dynamics is that this density dependence occurs in the birth or death rates of the species concerned. We have shown that the introduction of a suitable non-linearity into the dispersal behaviour of a species which behaves in an otherwise linear way can, in an inhomogeneous environment, lead to an exactly similar regulatory effect. Just such a mechanism is required
454
W.
S.
C‘.
GURNEY
AND
K.
M.
NISBkr
to explain the observation that populations of certain species appear to be able to regulate their numbers below the carrying capacity set by the supply of nutrients available to them. REFERENCES BRADFORD, E. W. & PHILIP, J. R. (19704. BRADFORD. E. W. & PHILIP. J. R. (197Ob). CARL, E. A. (1971). Eculog; 52, 39j. ’ FISHER, R. A. (1937). Ann. Eugen. 7, 355. GLANSWRFF, P. & PRIGOGINE, I. (1971).
J. theor. Biul. ZY, 13. J. theor. Biol. 29. 28.
Thermodynumic T/wry of Strncture, Stability and Fluctuations. London: Wiley-Interscience. GOLD, A, (1973). Science, N. Y. 181, 275. GURNEY, W. S. C. & NISBET, R. M. (in press). J. theor. Biul. GURTIN, M. E. (1973). J. theor. Biol. 40, 389. HOWARD, W. E. (1960). Am. Midl. Nat. 63, 152. KELLER, H. B. (1970). J. Diff. Eqs 7,417. KELLY, L. G. (1967). Handbook of Numerical Methods and Applicufiuns. Reading, Massachusetts: Addison Wesley. KREBS, C. J., KELLER, B. L. & TAMARIN, R. H. (1969). Ecology 50, 587. LIDIKER. W. Z. (1962). Am. Nat. 96.29. MONTR~LL, E. w. &WEST, B. J. (19?3). In Synergetics (H. Haken, cd.). Stuttgart: Teubner. MURRAY. B. G. (1967). Ecology 48,975. SCHIFF, L. I. (1955). Quantum Mechanics. New York: McGraw-Hill. SEGAL, L. A. & JACKSON, J. L. (1972). J. theor. Biol. 37, 545. SKELLAM, J. G. (1951). Biometrika 38, 196. STACKGOLD, I. (1971). SIAM Rev. 13, 289. WATON, A. & Moss, R. (1970). Animal Populatiun~ ill Relation ID their Food Resuarces. Oxford and Edinburgh: Blackwell Sci. Pub. WYNNE-EDWARDS, V. C. (1962). Animal Dispersiun in Relutiun to Suciul Behaviorrr. Edinburgh and London : Oiiver & Boyd. YOSHIZAWA, S. (1970). Math. Biosci. 7, 291.
APPENDIX
A
Microscope Models for Biased Random Motion and Directed Motion It is well known that in a population of random walkers the macroscopic current density is given by j = - DVp. (Al) Montroll & West (1973) recently modelled populations of “biased” random walkers who prefer to move up or down a population gradient; our present interest is in the case of walkers who dislike congestion and prefer to move down the population gradient. We imagine a “walk” in a d-dimensional rectangular grid (where d = I, 2 or 3). We let a be the distance between adjacent points in the grid and assume that during each time interval z every walker moves to one of its nearest neighbour points. If the motion were random, the probability of jumping
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
455
from a point r to a neighbouring point r +aB would be 2-‘; we bias this result by setting the probability of such a jump to P(r, ati) where ~(r, ae) = 2-d--+q.e, 642) and p > 0. With this definition the component of microscopic current density in the direction B is simply j,(r) = UT- ‘(p(r)P(r, d>-p(r + di)P(r+ ae, - ae)). (A31 If we assume that (I and r are small, and set D = 2-du2z-1 and I = 2&r-’ then we find: j,(r) = ( - DVp - ApVp) .A. (A41 It follows that the macroscopic current density is simply j = - DVp - IpVp, 645) a result in disagreement with Montroll & West (1973) who appear to have dropped some vital second-order terms in the1 derivation of their continuity equation. The concept of “directed motion” was introduced in the text to cover the important case of a territorial animal whose random motion makes no contribution to the macroscopic current density. The only choices offered to the walker are to stay put or to move down the population gradient. We choose the probability that an individual will jump to a nearest-neighbour point to be P(r, ae) = -ptavp . e H( - vp . e), WI where H(x) is the step function defined by H(x) = 1 if x > 0, 647) = 0 otherwise. The microscopic current density in direction A is then j, = 4pvp.e (fw where a is small and K = pa2z-‘, while the corresponding macroscopic current density is j = -KpVp. (A9
APPENDIX
B
The Existence of Positive, Stationary Solutions to the\ BRM and DM Equations Let 1 be a real Banach space with norm I( 11,and let A: g + W be a non-linear mapping which satisfies A(0) = 0. Consider the non-linear eigenvalue equation A(P) = PP9
PEW,
PE%
(AlO)
456
W.
S.
C.
GURNEY
AND
R.
M.
NISBE’I
which obviously has the trivial solution p = 0. A blench point of (AlO) is a number no with the property that in every neighbourhood of (PO, 0) in g x a there exists a solution of (AIO) with llyj # 0. The theory of branch points of (AlO) involves “linearizing” the operator A. If, when /i.ll is small, we can write A(p+1)-A(p)
= L,;1+om
(Al 1) then we say that A is linearizable at p and call L,, the FrPchet dericatice of A at p. We write Lo for the Frtchet derivative of .4 at p = 0. Stackgold (1971) in a lucid expository article quotes two theorems which enable us to locate the branch points of (AlO) provided A is linearizable atp = 0: (i) /lo can only be a branch point of (AIO) if it is an eigenvalue of L,, (ii) if A is completely continuous and p” ( # 0) is an eigenvaluc of Lo with odd multiplicity, then /I’ is a branch point of (AIO). To apply these theorems to our models it is convenient to start with the mnormafized BRM equation (5), a steady-state solution of which must satisfy the non-linear eigenvalue equation 0)
= PP?
(A121
where ,4(p) = - 1V. (pVp) - DV’p + &r)p.
(Al3)
The Frcchet derivative of ! at p = 0 is Lo = - DV’ +4(r),
(Al4)
which is essentially the operator we used to analyse the RM model. When the growth function is of the form shown in Fig. I, all the eigenvalues of I., are positive and there is a simple minimum eigenvalue ~1’. When j3 < /lo the only non-negative solution of (A12) is the trivial solution p = 0. To prove this we choose as norm in J IIPII = S p2 dr(Al? and use the original differential equation (5) to establish that for all strictly positive functions p
As the condition /j < no is equivalent to the condition that the normalized growth function be a RMD, we have proved that there is no positive stationary solution with a RMD growth function. However, the two theorems quoted above guarantee that the point /I = p” is a branch point of (A12). A simple perturbation analysis about the branch point reveals that the “new”
REGULATION
OF
INHOMOGENEOUS
POPULATIONS
457
solution of (A12) which emanates from /3 = p” is initially (i.e. when /I-p0 is small and positive) proportional to the eigenfunction corresponding to PO. For growth functions of the form of Fig. 1 we know that this eigenfunction is of fixed sign and in fact we may take it to be positive (see Keller, 1970), for analysis of this point in a similar example). It remains to be shown that this solution persists for all /3 > PO. This is done by considering functions of very small and very large norm. It is easy to show that if p = au with llull = 1, then for all positive functions U, there exists an a0 (dependent on u) such that $ ilp/ < 0
for a > cto.
(Al7)
It is also easy to show that for all positive u, there exists an a, such that after some initial transient T(U)
; llpll ’ 0
acal,
t> T(u).
GW
As the inequalities (A17) and (A19) hold for dip, the existence of a stationary solution is guaranteed, although of course its stability and the possibility of periodic behaviour remains to be investigated. This was done in the text. We can prove the existence of a positive stationary solution to the DM equation by exactly the same strategy as was used in the previous paragraph, i.e. by considering functions of large and small norm.