Murray J.D. Mathematical Biology: I. An Introduction

Подождите немного. Документ загружается.

36 1. Continuous Population Models for Single Species

1.7 Population Model with Age Distribution

For some populations, one of the deﬁciencies of the above ordinary differential equation

models is that they do not take into account any age structure which, in many situations,

can inﬂuence population size and growth in a major way. Although most natural popu-

lations have some structure, such as age, stage or whatever, it is not always of primary

importance but, when it is, we must know how to incorporate it in a model. So, we

consider here a ﬁrst extension to include age dependence in the birth and death rates.

One way of incorporating structure in a population is by Leslie matrices after

Leslie

(1945). For example, these can incorporate different age classes, such as juvenile

and adult, and quantify movement from one to the other. The Leslie matrix model then

relates the adult and juvenile population at a later time in terms of those at the earlier

time: the terms in the matrix incorporate, for example, the data on births and survival.

The books by Charlesworth (1980), Metz and Diekmann (1986) and Kot (2001) give a

good survey as well as the wide spectrum of applicability of age-structured models.

Let n(t, a) be the population density at time t in the age range a to a +da.Letb(a)

and µ(a) be the birth and death rates which are functions of age a: for man, for example,

they qualitatively look like the curves in Figure 1.21. For example, in a small increment

of time dt the number of the population of age a that dies is µ(a)n(t, a) dt. The birth

rate only contributes to n(t, 0); there can be no births of age a > 0. The conservation

law for the population now says that

dn(t, a) =

∂n

∂t

dt +

∂n

∂a

da =−µ(a)n(t, a) dt.

The (∂n/∂a) da term is the contribution to the change in n(t, a) from individuals getting

older. Dividing this equation by dt and noting that da/dt = 1sincea is chronological

(a) (b)

b(a)

1010 50 70

µ(a)

a (years)a (years)

Figure 1.21. Qualitative birth (a) and death (b) rates for man as functions of age in years.

Patrick Leslie was a mathematician who worked with the founding father of animal ecology, Charles

Elton, in his Bureau of Animal Population in Oxford in the mid-1930’s. In spite of the increasing use of

mathematical modelling, Elton never really embraced such interdisciplinary research and was on occasion

outspokenly critical. I met Elton once at a Fellows’ garden party in Corpus Christi College, Oxford, in the

1970’s. When I expressed an interest in animal ecology and mentioned his seminal work on the data from

the Hudson Bay Company (see Chapter 2) he started to talk enthusiastically about oscillatory behaviour

in populations. When I said I was a mathematician working in biology there was a notable cooling of his

enthusiasm and he said, ‘Oh, you’re one of them’ and added, ‘I thought you were somebody else.’

1.7 Population Model with Age Distribution 37

age, n(t, a) satisﬁes the linear partial differential equation

∂n

∂t

∂n

∂a

=−µ(a)n, (1.54)

which holds for t > 0anda > 0. For example, if µ = 0, it reduces to a conservation

equation which says that the time rate of change of the population at time t and age a,

∂n/∂t, simply changes by the rate at which the population gets older, namely, ∂n/∂a.

Equation (1.54) is a ﬁrst order partial differential equation which requires a condi-

tion on n(t, a) in t and in a. The initial condition

n(0, a) = f (a), (1.55)

says that the population at time t = 0 has a given age distribution f (a). The other

boundary condition on a comes from the birth rate and is

n(t, 0) =



∞

b(a)n(t, a) da, (1.56)

where, for mathematical simplicity, we have taken the upper limit of ∞ for the age;

b(a) of course will tend to zero for large a, as in Figure 1.21(a), for example, and so we

could replace ∞ by a

say, where b(a) = 0fora > a

. Note that the birth rate b(a)

only appears in the integral equation expression (1.56) and not in the differential equa-

tion. Equation (1.54) is often referred to in the ecological literature as the Vo n Fo e rs te r

equation; the equation arises in a variety of different disciplines and theoretical biology

areas, such as cell proliferation models, for example, where cell age is important. The

main question we wish to answer with the model here is how the birth and death rates

b(a) and µ(a) affect the growth of the population after a long time.

One way to solve (1.54) is by characteristics (see, for example, the book by Kevor-

kian, 2000) which are given by

= 1onwhich

=−µn. (1.57)

These are the straight lines

a =



t +a

, a > t

t −t

, a < t

(1.58)

as shown in Figure 1.22. Here a

, t

are respectively the initial age of an individual at

time t = 0 in the original population and the time of birth of an individual. The second

of (1.57), which holds along each characteristic, has a different solution accordingly as

a > t and a < t, that is, one form for the population that was present at t = 0, namely,

a > t, and the other for those born after t = 0, that is, a < t. On integrating the second

of (1.57), using da/dt = 1 and (1.58), the solutions are

38 1. Continuous Population Models for Single Species

a = t −t

a < t

a > t

a = t +a

Figure 1.22. Characteristics for the Von Foerster equation (1.54).

n(t, a) = n(0, a

) exp



−



µ(s) ds



, a > t,

where n(0, a

) = n(0, a −t) = f (a −t) from (1.55), and so

n(t, a) = f (a −t) exp



−



a−t

µ(s) ds



, a > t. (1.59)

For a < t,

n(t, a) = n(t

, 0) exp



−



µ(s) ds



and so, since n(t

, 0) = n(t − a, 0)

n(t, a) = n(t −a, 0) exp



−



µ(s) ds



, a < t. (1.60)

In the last equation n(t − a, 0) is determined by solving the integral equation (1.56),

using (1.59) and (1.60) to get

n(t, 0) =



b(a)n(t −a, 0) exp



−



µ(s) ds





∞

b(a) f (a − t) exp



−



a−t

µ(s) ds



da.

(1.61)

Although this is a linear equation it is not easy to solve; it can be done by iteration

however.

We are mainly interested in the long time behaviour of the population and in partic-

ular whether it will increase or decline. If t is large so that for practical purposes t > a

for all a then f (a − t) = 0 and all we require in (1.61) is the ﬁrst integral term on the

right-hand side. The solution is then approximated by n(t, a) in (1.60), although it does

not satisfy the boundary condition on a in (1.56). It is still not trivial to solve so let us

1.7 Population Model with Age Distribution 39

return to the original partial differential equation (1.54) and see if other solution forms

are possible.

We can look for a similarity solution (see, for example, Kevorkian 2000) of (1.54)

in the form

n(t, a) = e

γ t

r(a). (1.62)

That is, the age distribution is simply changed by a factor which either grows or decays

with time according to whether γ>0orγ<0. Substitution of (1.62) into (1.54) gives

=−

[

µ(a) +λ

]

and so

r(a) = r(0) exp



−γ a −



µ(s) ds



. (1.63)

With this r(a) in (1.62) the resulting n(t, a) when inserted into the boundary condition

(1.56), gives

γ t

r(0) =



∞

b(a)e

γ t

r(0) exp



−γ a −



µ(s) ds



and hence, on cancelling e

γ t

r(0),

1 =



∞

b(a) exp



−γ a −



µ(s) ds



da = φ(γ), (1.64)

which deﬁnes φ(γ). This equation determines a unique γ,γ

say, since φ(γ) is a mono-

tonic decreasing function of γ . The sign of γ is determined by the size of φ(0);see

Figure 1.23. That is, γ

is determined solely by the birth, b(a), and death, µ(a), rates.

The critical threshold S for population growth is thus

S = φ(0) =



∞

b(a) exp



−



µ(s) ds



da, (1.65)

where S > 1 implies growth and S < 1 implies decay. In (1.65) we can think of

exp [−



µ(s) ds] almost like the probability that an individual survives to age a, only

the integral over all a is not 1.

The solution (1.62) with (1.63) cannot satisfy the initial condition (1.55). The ques-

tion arises as to whether it approximates the solution of (1.54)–(1.56), the original prob-

lem, after a long time. If t is large so that for all practical purposes n(t, 0) in (1.61)

requires only the ﬁrst integral on the right-hand side, then

n(t, 0) ≈



b(a)n(t −a, 0) exp



−



µ(s) ds



da, t →∞. (1.66)

40 1. Continuous Population Models for Single Species

φ(γ)

φ(0)

Figure 1.23. The growth factor γ

is determined by the intersection of φ(γ) = 1:γ

> 0ifφ(0)>1, γ

< 0

if φ(0)<1.

If we now look for a solution of this equation in the similarity form (1.62), substitution

of it into (1.66) then gives (1.64) again as the equation for γ . We thus conjecture that

the solution (1.62) with r(a) from (1.63) and γ from (1.64) is the solution for large

time t of equation (1.54), with initial and boundary conditions (1.55) and (1.56). It is of

course undetermined to the extent of a constant r(0) but, since our main question is one

of growth or decay, it is not important to know r(0) since it does not affect either. The

important parameter is the threshold parameter S in (1.65) from which long time effects

of alterations in the birth and death rates can be assessed.

Exercises

1 A model for the spruce budworm population u(t) is governed by

= ru



1 −



−

1 +u

where r and q are positive dimensionless parameters. The nonzero steady states are

thus given by the intersection of the two curves

U(u) = r



1 −



, V (u) =

1 +u

Show, using the conditions for a double root, that the curve in r, q space which

divides it into regions where there are 1 or 3 positive steady states is given paramet-

rically by

r =

(1 +a

)

, q =

−1

Exercises 41

Show that the two curves meet at a cusp, that is, where dr/da = dq/da = 0, at

a =

√

3. Sketch the curves in r, q space noting the limiting behaviour of r(a) and

q(a) as a →∞and a → 1.

2 In Section 1.3 we showed that it was not possible to have a limit cycle (periodic)

solution for a simple ﬁrst order (nondelay) equation

= f (N).

A seeming counterexample is N(t) = 2 +sin t (and any number of similar forms).

Determine the f (N) for which this is a solution of the differential equation and

explain why it is not a counterexample.

3 If the per capita birth rate of a population is given by r[1 − a(N − b)

] where r, a

and b are positive parameters, write down a population model equation of the form

dN/dt = f (N). Nondimensionalise the equation so that the dynamics depend on a

single dimensionless parameter k = b(a/r)

1/2

.Ifu is your nondimensional popula-

tion, sketch f (u) for k > 1andk < 1 and discuss how the qualitative behaviour of

the solution changes with k and the initial condition.

4 The predation P(N) on a population N(t) is very fast and a model for the prey N(t)

satisﬁes

= RN



1 −



− P



1 −exp



−

ε A



, 0 <ε 1,

where R, K, P and A are positive constants. By an appropriate nondimensionalisa-

tion show that the equation is equivalent to

dτ

= ru



1 −



−



1 −exp



−



where r and q are positive parameters. Demonstrate that there are three possible

nonzero steady states if r and q lie in a domain in r, q space given approximately by

rq > 4. Could this model exhibit hysteresis?

5 A continuous time model for the baleen whale (a slightly more complicated model

of which the International Whaling Commission used) is the delay equation

=−µN(t) +µN(t − T )[1 + q{1 −[N(t − T )/K ]

}].

Here µ(> 0) is a measure of the mortality, q(> 0) is the maximum increase in

fecundity the population is capable of, K is the unharvested carrying capacity, T is

the time to sexual maturity and z > 0 is a measure of the intensity of the density-

dependent response as the population drops.

Determine the steady state populations. Show that the equation governing small

perturbations n(t) about the positive equilibrium is

42 1. Continuous Population Models for Single Species

dn(t)

≈−µn(t) − µ(qz − 1)n(t − T ),

and hence that the stability of the equilibrium is determined by Re λ where

λ =−µ −µ(qz −1)e

−λT

Deduce that the steady state is stable (by considering the limiting case Re λ = 0) if

µT <µT

π −cos

−1

−1)

1/2

, b = qz − 1 > 1

and stable for all T if b < 1.

For T = T

+ ε, 0 <ε 1 show that to O(1) the period of the growing

oscillation is 2π/[µ(b

−1)

1/2

], b > 1.

6 The concentration of carbon dioxide in the blood is believed to control breathing lev-

els through a delay feedback mechanism. A simple delay model for the concentration

in dimensionless form is

dx(t)

= 1 −axV(x(t − T )), V (x) =

1 + x

where a and m are positive constants and T is the delay. For given a and m a critical

delay T

exists such that for T > T

the steady state solution becomes linearly

unstable. For T = T

+ ε,where0<ε 1, carry out a perturbation analysis and

show that the period of the exponentially growing solution is approximately 4T

7 A model for the concentration c(t) of arterial carbon dioxide, which controls the

production of certain blood elements, is given by

dc(t)

= p − V (c(t − T ))c(t) = p −

max

c(t)c

(t − T )

where p, b, a, T and V

max

are positive constants. (This model is brieﬂy discussed

in Section 1.5.) Nondimensionalise the equation and examine the linear stability of

the steady state. Obtain a relation between the parameters such that the steady state

is stable and hence establish the existence of a bifurcation value T

for the delay.

Obtain an estimate for the period of the periodic solution which bifurcates off the

steady state when T = T

+ε for small ε.

8 A similarity solution of the form n(t, a) = e

γ t

r(a) of the age distribution model

equation

∂n

∂t

∂n

∂a

=−µ(a)n

satisﬁes the age boundary equation

n(t, 0) =



∞

b(a)n(t, a) da

Exercises 43



∞

b(a) exp



−γ a −



µ(s) ds



da = 1.

Show that if the birth rate b(a) is essentially zero except over a very narrow range

about a

> 0 the population will die out whatever the mortality rate µ(a).Ifthere

is a high, linear in age, mortality rate, say what you can about the birth rate if the

population is not to die out.

2. Discrete Population Models

for a Single Species

2.1 Introduction: Simple Models

Differential equation models, whether ordinary, delay, partial or stochastic, imply a

continuous overlap of generations. Many species have no overlap whatsoever between

successive generations and so population growth is in discrete steps. For primitive or-

ganisms these can be quite short in which case a continuous (in time) model may be a

reasonable approximation. However, depending on the species the step lengths can vary

widely. A year is common. With fruit ﬂy emergence from pupae it is a day, for cells it

can be a number of hours while for bacteria and viruses it can be considerably less. In

the models we discuss in this chapter and later in Chapter 5 we have scaled the time-step

to be 1. Models must thus relate the population at time t +1, denoted by N

t+1

,interms

of the population N

at time t. This leads us to study difference equations, or discrete

models, of the form

t+1

= N

F(N

) = f (N

), (2.1)

where f (N

) is in general a nonlinear function of N

. The ﬁrst form is often used to

emphasise the existence of a zero steady state. Such equations are usually impossible to

solve analytically but again we can extract a considerable amount of information about

the population dynamics without an analytical solution. The mathematics of difference

equations is now being studied in depth and applied in diverse ﬁelds: it is a fascinating

subject having given rise to some totally unexpected phenomena some of which we

discuss later. Difference equation models are also proving of use in a surprisingly wide

spectrum of biomedical areas such as cancer growth (see, for example, the article by

Cross and Cotton 1994), aging (see, for example, the article by Lipsitz and Goldberger

1992), cell proliferation (see, for example, the article by Hall and Levinson 1990) and

genetics (see, for example, the chapter on inheritance in the book by Hoppensteadtand

Peskin 1992 and the book by Roughgarden 1996.) It has recently been shown to be of

astonishing use in dynamic modelling of marital interaction and divorce prediction; we

discuss this application in Chapter 5. The largest use to date is probably in ecology; the

book by Hassell (1978) gives numerous examples, see also the more recent excellent

book by Kot (2001).

From a practical point of view, if we know the form of f (N

) it is a straightforward

matter to evaluate N

t+1

and subsequent generations by simply using (2.1) recursively.

2.1 Simple Models 45

Of course, whatever the form of f (N

), we are only interested in nonnegative popula-

tions.

The skill in modelling a speciﬁc population’s growth dynamics lies in determining

the appropriate form of f (N

) to reﬂect known observations or facts about the species

in question. To do this with any conﬁdence we must understand the major effects on the

solutions of changes in the form of f (N

) and its parameters, and also what solutions

of (2.1) look like for a few specimen examples of practical interest. The mathematical

problem is a mapping one, namely, that of ﬁnding the orbits, or trajectories, of nonlinear

maps given a starting value N

> 0. It should be noted here that there is no simple con-

nection between difference equation models and what might appear to be the continuous

differential equation analogue, even though a ﬁnite difference approximation results in

a discrete equation. This becomes clear below.

Suppose the function F(N

) = r > 0; that is, the population one step later is simply

proportional to the current population. Then from (2.1),

t+1

= rN

⇒ N

= r

. (2.2)

So the population grows or decays geometrically according to whether r > 1orr < 1

respectively; here r is the net reproductive rate. This particularly simple model is not

very realistic for most populations nor for long times but, even so, it has been used with

some justiﬁcation for the early stages of growth of certain bacteria. It is the discrete

version of Malthus’ model in Chapter 1. A slight modiﬁcation to bring in crowding

effects could be

t+1

= rN

, N

= N

1−b

, b constant,

where N

is the population that survives to breed. There must be restrictions on b of

course, so that N

≤ N

otherwise those surviving to breed would be more than the

population of which they form a part.

Fibonacci Sequence

Leonardo of Pisa, who was only given the nickname Fibonacci in the 18th century, in

his arithmetic book of 1202 set a modelling exercise involving an hypothetical growing

rabbit population. It consists of starting at the beginning of the breeding season with

a pair of immature rabbits, male and female, which after one reproductive season pro-

duce two pairs of male and female immature rabbits after which the parents then stop

reproducing. Their offspring pairs then do exactly the same and so on. The question is

to determine the number of pairs of rabbits at each reproductive period. If we denote the

number of pairs of (male and female) rabbits by N

then normalising the reproductive

period to 1 we have at the tth reproductive stage

t+1

= N

+ N

t−1

, t = 2, 3,... . (2.3)

This gives, with N

= 1, what is known as the Fibonacci sequence, namely

1, 1, 2, 3, 5, 8, 13,... .