Equilibrium and stability in populations whose interactions are age-specific

J. theor. BioL (1975) 54, 207-224

Equilibrium and Stability in Populations whose Interactions are Age-specific MANUEL ROTENBERG

University of California San Diego La Jolla, California 92037, U.S.A. (Received 19 July 1974, and in revisedform 21 January 1975) The continuous equations for a closed set of interacting populations are examined with regard to equilibrium and stability. The novel feature that is introduced is the age-specificity of the interactions between populations as well as the age-specificity of the effects of limited resources. First, single self-limiting populations are studied. Next, a set of interacting populations is examined for equilibrium points. Finally, the formal theory of stability in an age-specific description is detailed. Numerous simple examples are given.

1. Introduction

The notion of the Leslie matrix (Leslie, 1945, 1948) has given population biologists a facile conceptual and computational tool to deal simultaneously with time and age variables in ecological problems (Usher, 1972). The great appeal of the Leslie method, however, seems to have diverted attention from the fundamental differential equations from which it derives. In consequence, the age-dependent problem has escaped the analysis that the single-variable Volterra & Verhulst systems have been subjected to (God, Maitra & Montroll, 1971; May, 1973; Rescigno & Richardson, 1973). In this paper we examine the equilibrium and stability of interacting populations whose growth is described by differential equations rather than matrices. The analytic niceties of the continuous description are hard to ignore, as are the generalities it affords; the matrix description, for all its usefulness, is mathematically speaking a finite difference approximation to the differential equations. We shall first develop the basic ideas of equilibrium for a single species. These are then generalized to interacting systems, and finally the problem of local stability is examined. 207

208

M. ROTENBERG

2. Single S~eies (A) GENERALEQUATIONS

Imagine a closed single-species female population whose age distribution is described by the function n(a, t), where n(a, t)da is the number of individuals between the ages of a and a+da at time t. The integral of n(a, t) over all ages gives N(t), the total number of individuals at time t. The equations that govern the development of n(a, t) are (yon Foerster, 1959; Rotenberg, 1972)

On(a, t) Ot

On(a, t) - = -Do{a; In(a, t)]}n(a, t); a > 0 Oa B(t) = ~ re{a; In(a, t)]}n(a, t) da

(la) (lb)

where Do is the death rate (number of deaths per individual per unit time) and m is the fertility rate, analogously defined. In a limited envffonment both death and fertility rates are functionals of the standing population; this is indicated by the square-bracket notation and will be discussed further in the next section. In general, an explicit time dependence also appears in the vital rates, but since our interest centers about equilibrium populations, we consider only those vital rates that have no explicit dependence on time. The time appears in equations (1) only implicitly through the functionals of the population. B(t) is the number of births at time t; it is a convenient notation for n(0, t) and is therefore a density function. The infinite limit of age integration in equation (lb) is also a convenience since the oldest age of reproduction is controlled by the fertility function m. The interpretation of equation (la) becomes clear when it is recognized that the left side is nothing but the total time derivative; the total rate of charge of the population is due to deaths so long as births (a = 0) are avoided. The a = 0 case is covered by equation (lb) which sums the contributions from all elements of the population to the number of births at any instant of time t. A low order finite difference approximation to equations (1) yields the Leslie model. Equation (la) can be integrated. The result is

n ( a - t , O) exp -

a--[

Do(~; In(t-a)]) de ; a > t

(2a)

n ( a , t) =

B ( t - a ) exp ( - oi D°('; [n(t-a)])d~X);

a < t.

(2b)

EQUILIBRIUM

209

AND STABILITY

Neither equation (2a) nor (2b) can be regarded as a solution since n appears on the right sides and the birth rate is unknown. In equation (2a), n(a-t, O) is to be regarded as the given initial age distribution. It is equation (2b) that accounts for the regeneration of the population. Equations (2) are substituted in equation (lb) to obtain an equation for births for all time 8 ( 0 = F(0 + ~ B ( t - a)m(a; [,,(t)]) exp o

( - IoDo(C; [ n ( t - . ) ] ) de }

(3)

where

Equation (4) represents the births due to the initial population; its effect vanishes after all the members of the initial population pass the age of reproduction. Since we shall be concerned with asymptotic behavior, the basic equations do not involve equations (2a); they are equations (2b) and (3) at large times:

n(a, t) = B(t-a) exp ( - i Do(~; [n(t-a)] d~} B(t)=

~oB(t-a)m(a; In(t)])exp ( - i D o ( ~ ; In(t-a)])d~}

(5) da.

(6)

(13) THE VITAL RATES

Equations (5) and (6) are rigorously correct for a large, spatially homogeneous population. Ironically, "the behavior of the population with respect to time depends most sensitively on quantities which are the most difficult to observe: Do(a)and m(a). A population is "modeled" once these functions are chosen. The simplest non-trivial example of a death rate that depends on the population itself is one in which the population enters as a linear functional,

Do(~; [n(a, t)]) = d(~) + ~ u(~, a)n(a, t) da.

(7)

o

The function d(O is the usual age specific death rate. The effects of the standing population is embodied in the interaction function u(~, a). The quantity u(¢, a)n(a, t) da is the increment to the death rate on the population of age ~ due to the population between the ages of a and a+da. Suppose that the increased death rates were due to limited resources, then u(~, a) serves to transform the cause (the population around age a) into the effect T.B.

14

210

M.

ROTENBERG

(increased deaths in the population of age O. Since different age groups tap the resources at different rates, and since the lack of resources effect different age groups in the population differently, u in general is a function of both and a as shown in equation (7). Much analytical progress can be made when it is realized that age-specific correlation between cause and effect in u is, for a wide class of problems, either weak or non-existent. The rate of consumption of resources may be a-dependent, and the effect on deaths of the lack of resources may be ~dependent, but the ages of the deaths in the population may not be correlated with the ages of the consumers in the population. A counter-example would be a situation in which mature hunters succeed in creating a food shortage by feeding themselves, leaving the very young and very old to starve. Barring such correlations, the interaction function can be written in the product form u(~, a) = q(Op(a) and the interaction function is called separable. It is assumed henceforth that the interaction is separable. To see how use is made of the separability assumption, suppose a population is limited by resources so that it asymptotically reaches a time independent age distribution ~(a) = n(a, co). (8) (Asymptotic quantities will be denoted by a circumflex.) We now introduce the age-specific consumption rate p(a) (calories per individual per age interval per time interval) such that the total consumption rate is asymptotically

P1 = ~ p(a)ft(a) da.

(9)

O

Consumption of limited resources limits the population by inducing an increase in the death rate. Call the induced death rate q(~), such that at equilibrium Do(l; If~(a)]) =--1~o(0 = d(O + q(OPl.

(10a)

It will appear that whenever a separable form is used, then at equilibrium the energy consumption P'~ will. be the solution of an eigenvalue problem, The linear form indicated above, while simple, may not be adequate for some purposes. Once the separable assumption is made almost any function of P1 can be used. An example worked out in the next section is the "catastrophe model"

~o(~) = d(~) + q(~)P1/(P-P~)

(lOb)

which behaves like the linear model for P~ ~ P; the lack of resources becomes emphasized as the consumption approaches the availability P.

EQUILIBRIUM

AND

STABILITY

211

The functional dependence of fertility rates on the ambient population is usually ignored in mathematical population analysis, but it need not be. Depending on the species, fertility rates can either be enhanced or inhibited by crowding. An example exhibiting inhibitory response is m(a; [~t(a)] -~ ~(a) = (M - #/~)a a e- ra (11 a) where/~ is the total population at equilibrium, = ~ ft(a) da

(lib)

o

and # and M are constants. When the stability of interacting populations is discussed, it will be seen that the precise forms taken for Do and m are important. (C) THE EQUILIBRIUM POPULATION

If the population reaches an asymptotic distribution ~i(a) which is independent of time, it does so as the number of births reaches a constant -

B(oo).

Likewise, death rates and fertility rates become functions independent of time, but dependent upon functionals of fi(a). At equilibrium, therefore, equation (6) becomes 1=

r~(a) exp

- f/5o(~)d~

o

da

(12)

0

which will be an eigenvalue equation when specific forms for 130 and ~ are substituted. Equation (12) does not involve /~. The asymptotic births are obtained from equation (5) at ~quilibrium after multiplying it through by p(a) and integrating over all ages,

p(a)~(a) da B =

o

(13)

I p(a) exp o

/3o(~)d~ da o

where the value of the integral in the numerator of equation (13) is obtained as an eigenvalue from equation (12). Finally, the equilibrium age distribution is obtained from the asymptotic form of equation (5) ~(a) = / ~ exp

I - I/)o(~) ° d~] ,

(14)

o

and the carrying capacity/~ is found by integrating equation (14) over all ages.

212

1~f. R O T E N B E R G

(D) EXAMPLES Some simple algebraic examples follow. In order to avoid computers we do not take full advantage of the method and use constants in the vital rates wherever possible. We go through the equations as if such simplifying assumptions did not provide shortcuts; the examples are intended only to serve as mathematical travel guides. Take the death rate to be

(15)

bo(~) = d + q P l

where d and q are constants and PI is defined by equation (9). We take p to be constant also so that P1 = P-~, where ~ is the total asymptotic population. This is not unlike the Verhulst form. For the fertility rate take an inhibitory form m ( a ) = ( M - p l ~ ) 2 a e -ra (16) where M,/z and ~ are constants characteristic of the species,Equation (12) can now be written as 1

(M-~/~)2a exp [ - ( ~ + d + q P 1 ) a ]

da

(17)

o

where PI is to be regarded as an eigenvalue. The real eigenvalue is sought; in general there is only one (Coale, 1972, p. 64). The solution is P~ =

M - ~_,~ - ~ - d

q

(18)

where we pretend that/~ is as yet unknown. The births are given by equation (13) B =--P1 ( d + q P 1 )

(19)

P which becomes = ( M - 1a1~- y - tO(M - la2~ - ~,)/(pq)

(20)

when equation (18) is substituted. With P1 and ~ thus determined, the age distribution is given by equation (14) exp [ - ( M - / z ~ - d ) a ] .

ft(a) = ( M - / z R - r - a ) ( M - l a ~ - d ) Pq

(21)

The total population _~, or carrying capacity, is obtained from equation (21) by integrating it over all ages: /~

= M-- y - d

It + Pq

(22)

EQUILIBRIUM AND STABILITY

213

and can be substituted back into equations (20) and (21). T h e equilibrium birth equation (20) is particularly important because its going to zero marks extinction of the species:

B = ( \# + pqPq/~ 2 ( M - , - d ) ( M - ,

# d) . - -~

(23)

Another example is provided by the catastrophe model given in equation (10b). The same assumptions used in the previous example, except that (10b) is taken as the death rate, provides, through equation (12)

( M-#l~-~,-d~ P1 = P \ M _ l z t ~ _ ~ _ d + q ] .

(24)

The energy consumed is always a fraction of the energy provided, as expected. Equation (13) yields the equilibrium births

P ( M - #~ - ~,- d)( M - Izl~ - y) = -p (M-IM~-y-d+q)

(25)

and the carrying capacity ~ is given by integrating equation (14). The result is a quadratic which is solved for N:

P ( M - , - d ) ] 112}

.

(26)

The ambiguity in sign is resolved by requiring that/~ vanish as the available energy P vanishes. As a final example we relieve the restriction that p(a) is a constant and take p(a) = pa so that

el =- P S a~(a) da. That is, consumption increases linearly with age. We continue to assume that the rate of starvation is independent of age so that q remains constant. The form equation (16) is retained for the fertility rate. Equations (15), (17) and (18) remain the same, but equation (19) is replaced by equation (13) which now reads

=

~

Pt =~(d+qPt)2. P P S a exp [-P+qPl)a] da o

214

M: ROTENBERG

This equation, isreturned to~equation.(!4), for.:the..-age:distribution~.

~(a) = ~ (d + qPl) 2 exp [ - ( d

+ qPi)a].

(27)

P1 is given by equation(18) but is still unknown because ~ is unknown. To determine ~ , equation (27) is integrated over all ages and equation (18) is substituted for PI: q This quadratic is solved

_

2.

2./

?

7

and now Px is known through equation (18), ~(a) is known through equation (27) and/~ through equation (19). We see that the form ofp(a) only complicates the solution for ~ . It may seem that the tractability of the foregoing analysis rests on two assumptions. First, that the interaction function u(¢; a) is separable, and second, that the energy consumption is always a linear functional of ~(a). In fact, the simplicity rests only on the first assumption. The second one may be physically true, but mathematically it is little more than a convenience. Suppose, for example, through some curious eating habits, consumption depends on the square of the population so that the death rate appears as

SP(a)[A(a)] 2 da

O0

~0(¢) = d(¢)+ q(¢)

(28)

o

instead of equation (10a). But since the integral is only. a number, it is regarded as an eigenvalue to be obtained from equation (12), and the analysis proceeds essentially as before. Instead of equation (I3) the equation for equilibrium births will turn out to be

SOP(a)l'~(a)]2 da p(a)

112

¢] d

where the numerator, as before, is the eigenvalue obtained from equations (12) and (28).

EQUILIBRIUM AND STABILITY

215

3. Interacting Species (A) 6ENFa~L EQtrAr,ONS

We consider the differential equations for v interacting populations without distinguishing between competing or prey-predator relationships. The generalization of equations (la) and (lb) for the ith population is (indices run from 1 to v unless explicitly restricted)

Oni(a, t) + Oni(a, t_~)= _Dos(a; [n(t)])n~(a, t) dt Oa

(29a)

and

Bi(t) = ~ mi(a;

In(a,

t)])ni(a, t)da.

(29b)

0

The bold face n in the death and fertility rates signifies a functional dependence upon all populations. It is assumed that Do~and mi are such that stable o r unstable equifibrium can result even though such equilibrium may involve the extinction of one or more populations. The procedure that led to equations (5) and (6) can be used to integrate equation (29a) at equilibrium. Equilibrium quantities are indicated by a circumflex, as before.

~(a) = B, exp [ - i 3o,(~) d~]

(30)

oo

Bi = ~ m~a)~i(a ) da.

(31)

0

Substitution of equation (30) into equation (31) results in the equation

l = i mi(a) exp [ - i Do,(~) d~] da.

(32)

(B) SPECIFICFORMS The strategy used in solving equations (31) and (32) depends on the form adopted for Dos and rn~. As a useful and simple example, we choose an additive model: co

Doi(~) = di(~) + ~ S ui~(~, a)ttj(a) da, j o

03)

where di is the death rate for the ith population with ample resources in the absence of other ,population. Define oo

wo(~) -- ~ u~j(~, a)~j(a) da. 0

(34)

216

M. ROTENBERG

The equations which must be solved self-consistently are then equation (32) 1 = o~n~(a)exp [-D~(a)] exp [ - ~ i wu(~)d~] da,

(35)

and the equations obtained by multiplying equation (30) through by uu and integrating over all ages: w,j(~) = Bj i uj(~, a)exp

[-Di(a)] exp [ - ~ io wjt,(~)de] da.

(36,

We have used the notation i1

Dj(a) ~ $ dj(~) de. o

(37)

Equations (35) and (36) are to be solved for wu(~) and/~j. These quantities are then used in equation (30) to obtain the equilibrium age distribution. We have not attempted a general study of the system of equations (35) and (36) for uniqueness, conditions for existence, algorithms for solution, and so on; it would appear to be a non-trivial task. But once again, if uu(~, a) is assumed to be separable, the equations reduce to eigenvalue problems. Take

uu(~, a) = qu(~)pu(a)

(38)

where pij(a) is the age-specific consumption rate by species j of species i (cannibalism is permitted). The age-specific death rate induced by species j in species i of age ~ is qtj(~). For convenience, denote a

Qu(a) --- S qu(~) d~

(39)

Pti = ~o Po(a)ftJ(a)da

(40)

i Dol(~)d~ = Dt(a) + ~ Q~j(a)Ptj.

(41)

o

and

so that o

j

Since ~j is unknown, the Pij remain to be determined; there are v2 of them. Equation (35) appears as

1= ~o n~,(a)exp [-Di(a)- ~j Qu(a)P,Jl da. There are v such equations.

(42,

EQUILIBRIUM AND STABILITY

217

Equation (36) becomes

P'J= ~J i P'J(a)exp [--Dj(a)-- ~'QJ~(a)Pjk]

(43)

The unknown/~j are eliminated by taking ratios of successive pairs of Plj in equation (43)

i = 1. . . . . v - l ; j

= I .....

v.

(44)

The range of indices shows that there are v(v-1) such equations. These, plus the v equations from (42) complete the set of v2 equations required for the Pw The appropriate Pij are returned to equation (43) to obtain the ~i, and the equilibrium age distributions are obtained from equation (30). (C) EXAMPLE For reasons of algebraic simplicity we consider a system of two populations. We further resort to mi(a ) -- M2a exp (-),ia)

Dt(a) = Dia

(45)

Qij(a) = Qija pij(a) = ni~ where quantities without the functional dependence explicitly indicated are stipulated to be constants. It will be seen that these constants are not completely at our disposal. Equation (42) integrates to

1 = M2(~+D~+Q~IP~I+Qi2 P 12)- ,"

(46)

and equation (44) is simply

Pljw=~.~lJ P2j

(47)

~2j

The ensuing algebra results in Pxl =

7~11~22Q22 × 7~12n21Q12Q21 - 7~117~22Q1:tQ22 7~12Q12 ×[×(M2"7)~2-D2)-(Ml-Yl-D1)]

(48a)

218

M. R O T E N B E R G

P22 =

7~117~22Q11

x

7~127~21QlzQ21 - 7~117~22Q11Q22

x[7~21Q~21(MI--~l-DI)-(M2--~2--D2)] LrcIIQI i

(48b) (48c)

PI2 ---"P22 7~12 ~22

P21 ----Pil

7~21 7Zl I

(48d)

-

The number of births is

Ptj J~j ----- - (Dj + Qj tPjl -i- Qj2Pj2) = 7tij

(My --

~j).

(49)

The first equality comes from equation (43), the second from equation (46). Note that the right side of equation (49) must be independent of i, which it is because of equation (47). We now discuss equations (48) and (49) in the context of a prey-predator system. Imagine that species 1 is an herbivor and that species 2 preys solely upon species 1; no cannibalism takes place. The following significance is attached to the various quantities: nl 1 is the rate at which species 1 consumes energy per individual in the form of herbage; n12 is the rate at which the predator consumes energy per individual predator in the form of prey; n21 is the rate at which the prey yields energy to the predator; and n22 = 0 since there is no cannibalism. In Qll is embodied information regarding the availability of herbage. A small Q11 indicates that large consumption is possible. Q12 reflects the food value of the prey. If Qt2 is small, then each individual prey represents a large amount of energy to the predator; that is, the smaller Q~2 is, the smaller the kill rate is for a given predator population. Q2~ represents the sensitivity of the predator to the prey population. Since the predator by hypothesis is dependent solely upon the prey, it must be that Q21 < 0. That is, increasing numbers of species I induces a decrease of the death rate of species 2. The birth and death rate of species 2 must be such that extinction results if species 1 is absent. This will become apparent shortly. Equations (48a)-(48d) become, since rc22 is zero, 7~11

Pxl -- 7~2~21 ( M 2 -- r2 -- 9 2 )

(50a)

EQUILIBRIUM AND STABILITY P22

(50b)

= 0

PI2 =

219

(MI_y,_Di)_(M2_r2_D2) ]

~llQil [n21Q2____~l n21Q12Q21 LffI1QI1 1

P21 = ~ (M2 -~2 -D2). ~g2l

(50c)

(50d)

All energy rates and populations are positive, so that all P~s > 0. Since Q21 < 0, it follows from equations (50a) and (50d) that M 2 - Y 2 - D 2 is negative. Furthermore, the square bracket in equation (50c) must be negative so the constants must be chosen to satisfy ~21[Q211 (M 1- Y l - D I ) > I(Ma-y2-Dz)I. ~IIQll

(51)

A species is defined as being in the process of extinction when its equilibrium births B~ becomes zero. Referring to equation (49) we find that I

Bl = (M1-Yl)(M2-~'2-D2) 7c21Q21

(52a)

and ¸ B2 =

(M,_yl_D1)_(M2_~2_D2)].(52b)

zqlQ____2£ [~21Q21

Thus we see that by increasing the death rate D1 or decreasing the birth rate Mt it is possible to reduce B2 to zero, but because M 1 - 71- D1 will become zero before M x - h does, Bx will remain positive when B2 = 0. Similarly, by manipulating ntj and Qlj B2 will be reduced to zero without extinguishing species 1. Further, if it is attempted to reduce B I to zero by increasing 71, for example, it will be seen that Be becomes zero before/~ does. Conservation of energy imposes its own restriction through the calorie budget. The source of energy is through species 1, and species 2 receives the energy only after it has been degraded; therefore Pit > Pt2.

(53)

In addition, the energy given up by species 1 to species 2 must approximate closely that received by species 2 from species 1, the difference being wasted:

PI2 ~ P21

(54)

220

M. zoz~Nn~go

or, from equations(50c)and (50d) P,2-P,t =

1

[I

~2iQ,2xq'O" {Mz-,z-Dz[-(M,-y,-D,)]

~

O.

(55)

Equations (53) and (54)imply that P1, > P2,

(56)

which, from equations (50a) and (50d), requires that

~,

>

~2~.

(57)

4. Stability (A) GENERAL

We shall now investigate the behavior of an age specific trophic system that has somehow reached equilibrium and is subsequently perturbed by a small but arbitrary fluctuation in population starting at t =- 0. It is necessary to return to the time dependent equations,

n,(a, t) = B~(t-a) exp { - i Do,(~; [n(t-a)]) d~}

(58)

B,(t) = i mi(a; [n(t)])n~(a, t).

(59)

o

The time dependence in the vital rates is limited to that of the population. For example (use will be made of these later) o0

Doi(~; ['n(t-a)])

= di(~) + ~ qu(O S pij(x)nj(x, j

En(t)J)=

o

t-a + x) dx

exp

(60)

(61)

Substitution of equation (58) into equation (59) gives the equation governing the births as a function of time

8~(t)-- B~O-a)m,(a; [n(O]) exp -

Do~(~; [nO-a)]) d~ da.

(62)

o

It is the time delayed regeneration that make for interesting stability problems, and at least mathematically, they are quite distinct from the Volterra problem. The births and the populations are perturbed from their equilibrium values by small amounts beginning at t = 0: B~(t) = B~+~B~(t) t > 0 (63a) nj(a, t) = ~l(a)+ 6nl(a, t) t > 0. (63b)

EQUILIBRIUM

AND

STABILITY

221

Substitution into equation (62) and expansion to lowest order gives the equilibrium equations studied in the last section. It will be assumed that the equilibrium equations are solved and that the equilibrium quantities can be used in the first-order perturbation equations. It is difficult to write the first-order perturbation equations in general without bringing to bear a rather involved notational system for the necessary functional derivatives. It will be seen, however, that one is led to complete generality by examining the case in which the vital rates are given by equations (60) and (6I). The substitution of equations (63) into (60) and (61) gives co

Do,(~; l-n(t- a)]) =/)o,(0 + ~ qij(~) ~ P~j(x)rni(x, t - a + x) dx 2

(64)

0

and

m,(a; [n(t)]) = ~ ( a ) - 2M i ~. piia #' exp ( - yta) 6Nj(t)

(65)

d

where O3

cSNj(t) - S 6nj(a, t) da.

(66)

o

Equations (63)-(65) are substituted in equation (62). The first order result is 6Bi(t) = 1 6B,(t-a)rhl(a ) exp - I/)ol(O d~ da + 6Ai(t), o

(67)

o

where

~At(t) = Bt So r~l(a) exp -

D~t(~)d~ x o

x j

I

I p,j(xl n/x, t-a+x) dx+

o

o

oxp <-,,o>] × x exp -

~ot(~) d~ da 6Nj(t).

(68)

The upper limit of integration in equation (67) is t rather than infinity because it has been specified that the perturbation is zero for t < 0, that is, 6Bi is zero when its argument is negative. Two points are worth noting regarding equation (67). First, the first term on the right is entirely independent of the form of the perturbation ~n~(t).

222

M. ROTENBERG

Second, .t h e second term, 6At(t), involves only known quantities; further, it depends not only on the form of the perturbation, but on the functional derivatives of the vital rates~with respect t o n. Put another way, the first term on the right is characteristic of the system, the second term is the driving term characteristic of fluctuations from without the system. We now proceed to find the natural modes of the system. Taking the Laplace transform of equation (67), it is found that

~SA*(s)
(69)

where starred quantities are Laplace transforms

y*(s) -- ~ y(x) e -'x dx

(70)

0

and, in particular, q~*(s) =

~,(x) exp

-

0

Dol(~) d~ e -'x dx.

j"

(71)

0

The inverse Laplace transform of equation (69) gives the time behavior of the birth due to the perturbation:

~B~(t) - -,~o+,S~SB*(s) e~ ds = -,~o+, ~ 1 -q~*(s) e" ds

(72)

where e is the appropriate positive quantity (Morse & Feshbach, 1953). The integral is performed by contour integration in the usual way, with the result ~Bj(t)

=

-- ~

~SA*(sYk)

~

"

~.~k,

(73)

where the s~k are the complex roots of q~*(sj~) = 1 = exp (2xki);

(74)

there are an infinite number of them. The time behavior of the perturbation is a function of where the s~k fall on the complex plane. If the real part of any sik is positive, the perturbation grows; if it is negative, the perturbation damps. If sjk is pure imaginary, the perturbation oscillates without damping. It is now seen why it is unnecessary to delve into the details of a functional expansion as we did in equation (68). 6A;(sik) merely gives the amplitude of each complex eigenmode sii. These amplitudes are generally uninteresting quantities, for i t is the eigenfrequency with the largest positive real part that controls how the system goes unstable, not the amplitude with which

EQUILIBRIUM AND STABILITY

223

that mode is excited. The term "unstable" must, as usual, be understood within the context of small perturbation and linear expansions. There is no guarantee that in general an instability will continue to grow exponentially. (B) EXAMPLE

We treat a system of two species, and we again use constants wherever possible for ease of manipulation. Species 1 is oviparous. Its numbers are controlled by species 2 which is ovivorous. Members of species I, once born, enjoy ample resources. Thus we take for the equilibrium vital rates #/l(a) = (M 1 --/z/~2)a p' exp (--yla)

/5ol(a ) = Dla th2(a ) = M2ap" exp (-r2a)

(75)

/3o2(a) = (D2-0~/~l)a. The constant ~ plays the role of q21P21 in previous examples. The total equilibrium populations Rx and/~2 are treated as eigenvalues of equation (32). The results are N1 = lo~{D2 +Y2

- [ M 2 F ( f l 2 ) ] 11(#2+ 1))

(76a)

and /~2=~[ MI

(Dl+vx)a~+x]F-~x) J

(76b)

where F(x) is the gamma function. We are now prepared to find q~7(s) from (71) and solve equation (74) for the eigenmodes. By putting equations (75) into (71) we have (MI -/aR2)F(flt ) q~*(s) = (01 +?t +s) a' + 1

(77a)

and M2F(fl2)

~ ( s ) = '(D2+72'--a!~1 +s) p~+I

(77b)

and equation (74) gives the roots, after using equation (76),

k/h+:/

S2k = [exp ( - 2 ~ k i ) - 1] [M2F(fl2)]'/(''+').

(78a) (78b)

224

M. ROTENBERG

Assuming that p~ are irrational numbers, all roots are distinct. In case all roots are not distinct, a special discussion is necessary (Feller, 1941). We notice that the real parts of all the s~k are negative with the exception of sto, both of which are zero. With these exceptions then, the system is unconditionally stable. T h e Sto are neutral modes in that once excited, they do not oscillate, nor does the system return to the equilibrium point. The source of this mode is from the fact that the vital rates we have chosen [equation (75)] are sensitive only to aggregate populations, and the aggregates are indifferent to the kind o f perturbation that involves the simple interchange of the numbers of individuals in different age groups. REFERENCES COALE,A. J. (1972). The Growth and Structure o f Human Populations. Princeton: Princeton University Press. FELLER,W. (1941). Ann. math. Statist. 12, 243. GOEL,N. S., MAtTaA,S. C. & MONTROLL,E. W. (1971). Rev. mod. Phys. 43, 231. L~r~E, P. H. (1945). Biometrika 33, 183. L~t.m, P. H. (1948). Biometrika 35, 213. MAY, R. M. (1973). Model Ecosystems. Princeton: Princeton University Press. MORSE,P. M. & F~HS^CH, H. (1953). Methods in TheoreticalPhysics. New York: McGrawHill. Rr:scaoNo, A. & RaCHARDSON,I. W. (1973). In Foundations o f Mathematical Biology, Vol. 3 (R. Rosen, ed.). New York: Academic Press. ROT~B~G, M. (1972). J. theor. Biol. 37, 291. USHER, M. B. (1972). In Mathematical Models in Ecology. 12th symposium of the British Ecological Society O. N. R. Jeffers, ed.). London: Blackwells. 'CoN FoERsr~rt, H. (1959). In Kinetics of Cellular Proliferation (F. Stohlman, Jr, ed.) New York: Grune and Stratton.