Elements of Mathematical Ecology

Elements of Mathematical Ecology provides an introduction to classical and

modern mathematical models, methods, and issues in population ecology.
The first part of the book is devoted to simple, unstructured population
models that, for the sake of tractability, ignore much of the variability
found in natural populations. Topics covered include density dependence,
bifurcations, demographic stochasticity, time delays, population interactions
(predation, competition, and mutualism), and the application of optimal
control theory to the management of renewable resources. The second part
of this book is devoted to structured population models, covering spatially
structured population models (with a focus on reaction-diffusion models),
age-structured models, and two-sex models. Suitable for upper level students
and beginning researchers in ecology, mathematical biology and applied
mathematics, the volume includes numerous line diagrams that clarify the
mathematics, relevant problems throughout the text that aid understanding,
and supplementary mathematical and historical material that enrich the main
text.

MARK KOT is Associate Professor in the Department of Applied Mathe-

matics at the University of Washington, Seattle, USA.
Elements of Mathematical Ecology

MARK KOT
Department of Applied Mathematics,
University of Washington

CAMBRIDGE
UNIVERSITY PRESS
CAMBRIDGE UNIVERSITY PRESS
Cambridge, New York, Melbourne, Madrid, Cape Town, Singapore, Sao Paulo

Cambridge University Press

The Edinburgh Building, Cambridge CB2 8RU, UK

Published in the United States of America by Cambridge University Press, New York

www.cambridge.org
Information on this title: www.cambridge.org/9780521802130
© Mark Kot 2001

This publication is in copyright. Subject to statutory exception

and to the provisions of relevant collective licensing agreements,
no reproduction of any part may take place without the written
permission of Cambridge University Press.

First published 2001

Reprinted 2003 with corrections

A catalogue record for this publication is available from the British Library

Library of Congress Cataloguing in Publication data

Kot, Mark, 1956-
Elements of mathematical ecology / Mark Kot.
p. cm.
Includes bibliographical references (p.).
ISBN 0 521 80213 X - ISBN 0 521 00150 1 (pb.)
1. Population biology - Mathematical models. 2. Ecology - Mathematical models.
3. Population biology — Mathematics. 4. Ecology — Mathematics. I. Title.
QH352.K66 2001
577'.01'51 - d c 2 1 00-065165

ISBN 978-0-521-80213-0 hardback

ISBN 978-0-521-00150-2 paperback

Transferred to digital printing 2008

Contents

Preface page vii

Acknowledgments ix

I. UNSTRUCTURED POPULATION MODELS 1

A. SINGLE-SPECIES MODELS 3
1 Exponential, logistic, and Gompertz growth 3
2 Harvest models: bifurcations and breakpoints 13
3 Stochastic birth and death processes 25
4 Discrete-time models 43
5 Delay models 70
6 Branching processes 93

B. INTERACTING POPULATIONS 107

7 A classical predator-prey model 107
8 To cycle or not to cycle 116
9 Global bifurcations in predator-prey models 140
10 Chemostat models 161
11 Discrete-time predator-prey models 181
12 Competition models 198
13 M utualism models 220
C. DYNAMICS OF EXPLOITED POPULATIONS 237
14 Harvest models and optimal control theory 237

II. STRUCTURED POPULATION MODELS 265

D. SPATIALLY STRUCTURED MODELS 267

15 Formulating spatially structured models 267

v
v i Contents

16 Spatial steady states: linear problems 276

17 Spatial steady states: nonlinear problems 294
18 Models of spread 311
E. AGE-STRUCTURED MODELS 345
19 An overview of linear age-structured models 345
20 The Lotka integral equation 353
21 The difference equation 365
22 The Leslie matrix 377
23 The M cKendrick-von Foerster PDE 391
24 Some simple nonlinear models 401

F. SEX-STRUCTURED MODELS 413

25 Two-sex models 413
References 425
Author index 443
Subject index 447
Preface

Ecology is an old discipline. The discipline was christened in 1866 by

Ernst Haeckel, a well-known German evolutionary biologist. Haeckel was a
neologist - he loved to invent new scientific terms. His most famous gems
are phylogeny and oecologie. Oecologie and ecology take their derivation
from the Greek oikos, house or dwelling place. Ecology, as envisioned by
Haeckel, is the study of the houses and the housekeeping functions of plants
and animals. It is the scientific study of the interrelationships of organisms,
with each other, and with their physical environment. The idea of ecology
is even older (Worster, 1994). It is closely related to 18th century notions of
the balance or economy of nature reflected, most clearly, in Linnaeus's 1749
essay Oeconomia Naturae (Stauffer, 1960).
Ecology is also a diverse discipline. After all, it has all of life to account for.
In the old days, it was common to divide ecology into two subdisciplines:
autecology, the ecology of individual organisms and of populations, and
synecology, the study of plant and animal communities. Ecology is now
divided into many subdisciplines (see Table 1).
Several subdisciplines use mathematics. For example, behavioral ecology
makes extensive use of game theory and of other brands of optimization. It
is impossible to cover all of these subdisciplines in one short book. Instead, I
focus on population ecology and engage in occasional forays into community
ecology and evolutionary ecology. This book could, and perhaps should, have
been entitled The Dynamics of Biological Populations.
The material in this book has been used to teach a two-semester course.
There is, therefore, a dichotomy in these notes. The first semester of the
course is devoted to unstructured population models, models that, in ef-
fect, treat organisms as 'homogeneous green gunk'. Unstructured population
models have the advantage, at first, of simplicity. As one adds extra bits

vn
viii Preface

Table 1. Branches of ecology

Synecology Landscape ecology

Systems ecology
Community ecology
Autecology Population ecology
Evolutionary ecology
Behavioral ecology
Physiological ecology
Chemical ecology

of biology, these models become more realistic and more challenging. The
topics in the first half of the book include density dependence, bifurcations,
demographic stochasticity, time delays, population interactions (predation,
competition, and mutualism), and the application of optimal control theory
to the management of renewable resources.
Variety, and variability, are the spice of life. We frequently ascribe dif-
ferences in the success of individuals to differences in age, space (spatial
location), or sex. The second half of this book is devoted to structured pop-
ulation models that take these variables into account. I begin with spatially-
structured population models and focus on reaction-diffusion models. There
is also tremendous interest in metapopulation models, coupled lattice maps,
integrodiflference equations, and interacting particle systems (Turchin, 1998;
Hanski, 1999). However, my colleagues and I tend to leave this material for
our advanced course. I follow with an overview of age-structured population
models in which I compare integral equations, discrete renewal equations,
matrix population models, and partial differential equations. I conclude with
a brief introduction to two-sex models.
The emphasis in these notes is on strategic, not tactical, models (Pielou,
1981). I am interested in simple mechanistic models that generate interesting
hypotheses or explanations rather than in detailed and complex models
that provide detailed forecasts. You will also find many equations, but few
formal theorems and proofs. Applied scientists and pure mathematicians
both have reason to be offended ! Because of the interdisciplinary nature
of my class and because of my own preference for solving problems over
proving theorems, I have tried to hold to a middle course that should appear
natural to applied mathematicians and to theoretical biologists. I hope that
this middle course will appeal to a broad range of (present and future)
scientists. Failing that, I hope that you, gentle reader, can use this book as
a springboard for more detailed applied and theoretic investigations.
Acknowledgments

I have been blessed with excellent teachers and students. I wish to thank all
my teachers, but especially William K. Smith, W. Tyler Estler, Richard H.
Rand, Simon A. Levin, William M. Schaffer, Paul Fife, Jim Cushing, Stephen
B. Russell, and Hanno Rund.
Stephane Rey coauthored Chapter 9. Other former students, Michael G.
Neubert, Emily D. Silverman, and Eric T. Funasaki, will recognize work that
we published together. Michael Neubert used this material in a class and
provided a number of useful comments and criticisms.
Several cohorts of students studied this material as Mathematics or Ecol-
ogy 581 and 582 at the University of Tennessee or as Applied Mathematics
521 at the University of Washington. I thank these students for their enthu-
siasm and hard work. I am grateful to the University of Tennessee and the
University of Washington and to my colleagues at these institutions for the
chance to teach these courses.
The early drafts of this book could not have been written without several
valuable pieces of software. I thank Joseph Osanna for troff, Brian W.
Kernighan and Lorinda L. Cherry for eqn, Jon L. Bentley and Brian W.
Kernighan for grap, Michael Lesk for ms, tbl, and refer, James J. Clark
for groff, Bruce R. Musicus for numeqn, Nicholas B. Tufillaro for ode, and
Ralph E. Griswold, Madge T. Griswold, and the Icon Project for the Icon
programming language.
It has been a pleasure working with Cambridge University Press. I wish
to thank Alan Crowden, Maria Murphy, Jayne Aldhouse, Zoe Naylor, and
especially Sandi Irvine for all their efforts.
Finally, I want to thank my parents for their encouragement and interest
and my wife, Celeste, for her encouragement, support, and desire to purchase
the movie rights.

Knoxville, Tennessee Mark Kot

IX
Part I Unstructured population models
Section A
SINGLE-SPECIES MODELS

1 Exponential, logistic, and Gompertz growth

Tradition dictates that we begin with a simple homogeneous population.

This population is that 'homogeneous green gunk' that I referred to in the
preface. I will represent the number (or sometimes the density) of individuals
in this population by N(t). I will also make frequent reference to the rate
of change, dN/dt, and to the per capita rate of change, (l/N)dN/dt9 of this
population.
Let us assume that all changes in this population result from births and
deaths and that the per capita birth rate b and per capita death rate d are
constant:

The difference between the per capita birth and death rates, r = b — d, plays
a particularly important role and is known as the intrinsic rate of growth.
Equation (1.1) is commonly rewritten, in terms of r, as

One must also add an initial condition, such as

iV(0) = No, (1.3)

that specifies the number of individuals at the start of the process.

Equation (1.2) is a linear, first-order differential equation. It is easily
integrated, either as a separable equation or with an integrating factor, and
it possesses the solution
N(t) = Noert. (1.4)

This solution grows exponentially for positive intrinsic rates of growth and
4 A. Single-species models

1 dN
N dt

N
Fig. 1.1. Per capita growth rate.

dN
~dt

rN

N
Fig. 1.2. Population growth rate.

decays exponentially for negative intrinsic rates of growth. It remains con-

stant when births balance deaths.
Three different graphs capture the behavior of this system. In Figure 1.1,
I have plotted the per capita growth rate as a function of the population
size. The per capita growth rate remains constant for all population sizes:
crowding has no effect on individuals. However, the growth rate for the
entire population (Figure 1.2) increases with number as each new individual
adds its own undiminished contribution to the total growth rate. The result
(Figure 1.3) is a population that grows at ever-increasing rates.
The population size N* = 0 is an equilibrium point. Since there is no
 Exponential, logistic, and Gompertz growth 5

r >0

Fig. 1.3. Population trajectory.

immigration or emigration in this model, populations that start at zero stay

at zero. For positive r, this equilibrium is unstable. After small perturbations,
the population moves away from zero. For negative r, this equilibrium is
asymptotically stable. Small perturbations now decay back to zero. I will say
more about equilibria and stability later.

Problem 1.1 Monod's^ nightmare

Escherichia coli is a bacterium that has been used extensively in microbio-
logical studies. Escherichia coli cells are rod shaped; they are 0.75 pirn wide
and 2 jum long. Under ideal conditions, a population of E. coli doubles in
just over 20 minutes.
(1) What is r for E. coli?
(2) If No = 1, how long would it take for an exponentially growing population of
E. coli experiencing ideal conditions to fill your classroom?

There are several defects with this simple exponential model:

(1) The model has constant per capita birth and death rates and generates limitless
growth. This is patently unrealistic.
(2) The model is deterministic; we have ignored chance or stochastic effects.
Stochastic effects are particularly important at small population sizes.

f Jacques Monod (1910-1976) was the recipient of a 1965 Nobel Prize for Medicine for his work on
gene regulation. He also conducted innovative experimental studies on the kinetics and stoichiometry
of microbial growth (Panikov, 1995).
6 A. Single-species models

1 dN

r—

I
N
K
Fig. 1.4. Decreasing per capita growth rate.

(3) The model ignores lags. The growth rate does not depend on the past. More-
over, the population responds instantaneously to changes in the current popu-
lation size.
(4) We have ignored temporal and spatial variability.

Let us start with the first defect.

What are the factors that regulate the growth of populations? There
have been two schools of thought. In 1933, A. J. Nicholson, an Australian
entomologist, published a seminal paper in which he stressed the importance
of density-dependent population regulation. Nicholson (1933), the British
ornithologist David Lack (1954), and others argued that populations are
regulated by biotic factors such as competition and disease that have a
disproportionately large effect on high-density populations. The opposing
view, promulgated by the Australian entomologists H. G. Andrewartha
and L. C. Birch (1954), is that populations are kept in check by abiotic,
density-independent factors, such as vagaries in the weather, that have as
adverse an effect on low-density populations as they do on high-density
populations.
The dispute between these two schools occupied ecology for most of the
1950s (Tamarin, 1978; Kingsland, 1985; Sinclair, 1989). Density-dependent
and density-independent factors may both be important in regulating pop-
ulations. From a modeling perspective, however, it is easier to start with
density-dependent regulation.
Consider a per capita growth rate that decreases linearly with population
 Exponential, logistic, and Gompertz growth 7

dN
~dt

Fig. 1.5. Parabolic population growth rate.

size,
1 dN
(1.5)
N dt
(see Figure 1.4). This decrease in the per capita growth rate may be thought
of as an extremely simple form of density-dependent regulation. Note that
the per capita growth rate falls to zero at the carrying capacity K.
The population's growth rate,

dN N
(1.6)

is now a quadratic function of population size (see Figure 1.5). Equation (1.6)
is known as the logistic equation or, more rarely, as the Pearl-Verhulst
equation. It has an exact analytical solution. Figure 1.6 illustrates this solution
for two different initial conditions. You are asked to find this closed-form
solution in Problem 1.2. Since few nonlinear differential equations can be
solved so easily, I will concentrate on a general method of analysis that
emphasizes the qualitative features of the solution.
Equation (1.6) has two equilibria, AT* = 0 and N* = K; at each of these
two values, the growth rate for the population is equal to zero. Near N* = 0,
N2/K is small compared to N so that

dN
rN. (1.7)
It
8 A. Single-species models

K —

Fig. 1.6. Logistic growth.

For r > 0, small perturbations about N* = 0 grow exponentially; the

equilibrium N* = 0 is unstable.

Problem 1.2 Exact solution of the logistic equation

Show that the logistic equation has the solution
K
N(t) = (1.8)
1 e n
+ ft - 0 ~
(1) treating the logistic equation as a separable equation, and
(2) treating the logistic equation as a Bernoulli equation.

Close to AT* = X, we instead introduce a new variable that measures the

deviation of N from K:
x = N - K. (1.9)
Substituting N = K + x into equation (1.6) gives us
dx _ r 2
rX (1.10)
It ~ ~ ~ KX '
and since x is small for N close to K, we have that
dx
——
tt T X. (1.11)
dt
For r > 0, small perturbations about N* = K decay exponentially; the
equilibrium AT* = K is asymptotically stable. For positive r, solutions to the
 Exponential, logistic, and Gompertz growth 9

logistic equation (see Figure 1.6) are essentially a combination of exponential

growth, close to zero, and of exponential decay, close to the carrying capacity.
Equations (1.7) and (1.11) imply that the solution of the logistic equation
is concave up just above the origin and concave down just below the carrying
capacity. It stands to reason that an inflection point lies between the origin
and the carrying capacity. This inflection point can be found by setting the
derivative of both sides of logistic equation (1.6) equal to zero:

It follows that the inflection point is at N = K/2.

Mathematical meanderings
Consider the differential equation

This equation is autonomous in that / does not contain any explicit dependence
on t. I have introduced several concepts that are useful, not only for the logistic
differential equation, but for many autonomous, first-order, ordinary differential
equations. Let's formalize these concepts.

Definition We say that N = N* is an equilibrium point (also known as & fixed

point, critical point, rest point) if
/(AT) = 0. (1.14)

Definition An equilibrium point N* is Lyapunov stable if, for any (arbitrarily

small) e > 0, there exists a 3 > 0 (depending on e) such that, for all initial
conditions N(to) = No satisfying | No — N* | < d, we have | N(t) — N* \ < e for
all t > to. In other words, an equilibrium is stable if starting close (enough)
guarantees that you stay close.

Definition An equilibrium point N* is asymptotically stable (in the sense of

Lyapunov) if it is stable and if there exists a p > 0 such that
lim IJV(t) - ATI = 0 (1.15)
t->co
for all No satisfying
| No - A T | < p . (1.16)
Thus an equilibrium is asymptotically stable if all sufficiently small initial
deviations produce small excursions that eventually return to the equilibrium.
At this point, the only interesting question is a practical one: is a given
equilibrium point stable or unstable?
10 A. Single-species models

Theorem Suppose that N* is an equilibrium point and that f(N) is a con-

tinuously differentiable function. Suppose also that ff(N*) =/= 0. Then the
equilibrium point N* is asymptotically stable if f'(N*) < 0, and unstable if
fW) > 0.
Proof Consider an equilibrium for which /'(N*) < Oandletx(t) = N(t) — N*.
If we expand f(N) about AT, we obtain

J = f(N*) + f'(N*)x + g(x). (1.17)

Since N* is an equilibrium, equation (1.17) reduces to

J = f'{Nm)x + g(x), (1.18)

which may be viewed as a perturbation of a linear, constant-coefficient differen-
tial equation. Note that g(x) consists of higher-order terms. In particular, g(x)
satisfies g(0) = 0 and also g'(0) = 0. This, along with the continuity of g'(x)
(which follows from the continuity of f'(N)), guarantees us that for each e > 0
there is a small d neighborhood about zero wherein | g\x) | < e. As a result,

g(x) = f g'(s)ds < e\x\. (1.19)

Jo
It follows that

J < /'(N')x + e|x|. (1.20)

For small enough d and e and f\N*) ^ 0, the higher-order terms cannot
change the sign of dx/dt. As a result, small enough perturbations will decay;
the equilibrium is asymptotically stable. A similar argument can be made to
show that f'(N*) > 0 implies instability. •
Take another look at Figure 1.5. You should be able to ascertain the stability
of the two equilibria by inspection with this theorem. What do you think
happens when /r(iV*) = 0?

Historical hiatus
The concepts of exponential and logistic growth arose gradually. A few people
played especially important roles in the development of these concepts.
John Graunt (1662) was a 'collector and classifier of facts' (Hutchinson, 1980).
He was also the inventor of modern scientific demography. Graunt tabulated
the Weekly Bills of Mortality for London. These bills listed births and deaths;
they were used as an early warning system for the plague. Using these bills,
Graunt estimated a doubling time for London of 64 years. This is an extremely
short period of time. Graunt posited that if the descendants of Adam and Eve
 Exponential, logistic, and Gompertz growth 11

(created in 3948 BC, according to Scaliger's chronology) doubled in number

every 64 years, the world should be filled with 'far more People, than are now
in it.' By my calculation, this would amount to
2(i662 +3948)/64 ^ 287.7 ^ 15 x 1Q26 ^ 2 00 million people/cm 2 . (1.21)

Graunt was clearly aware of the power of exponential growth.

Sir William Petty (1683) faulted Graunt for ignoring the biblical flood. He
started the clock with Noah (t0 = 2700 BC with No = 8). Petty also felt that
Graunt's estimate for a doubling time was misleading, since much of London's
increase was due to immigration. Petty estimated the doubling time for England
to be between 360 and 1200 years. However, a doubling time of 360 years left
Petty with a population projection,
8 x 2 12175 » 36994, (1.22)
that was far too small. Petty therefore proposed that the human growth rate
had fallen steadily in postdiluvian times. He produced a table 'shewing how the
people might have doubled in the several ages of the world'; the table exhibited
a steady increase in the doubling time, much as one would expect for logistic
growth.
The Reverend Thomas Robert Malthus (1798) is famous for having written
An Essay on the Principle of Populations. The essence of this book may be
represented with a simple quasi-equation:

a geometrically an arithmetically much

growing Hh growing = human (1.23)
population food supply misery
Many of Malthus's conclusions had already been anticipated by Graunt, Petty,
and others. However, Malthus is justly famous for stating the case so clearly.
Malthus's book had tremendous influence on Charles Darwin and Alfred Rus-
sel Wallace and, in effect, provided them with the material basis for natural
selection.
Pierre-Francois Verhulst (1845) was a Belgian who presented an entirely
modern derivation of the logistic equation. His work went unappreciated during
his own lifetime and he died in relative obscurity.
Raymond Pearl and Lowell Reed (1920) rediscovered the logistic equation and
launched a crusade to make the logistic equation a law of nature' (Kingsland,
1985). They published over a dozen papers between 1920 and 1927 promulgating
this law. Some of their extrapolations and ad hoc pastings of logistic curves
were questionable, but they did make the logistic equation famous.

The logistic equation allows us to handle limited growth in a natural way.

This equation also has a long history. However, there is nothing sacred about
this equation. Other models, derived differently, exhibit many of the same
properties.
12 A. Single-species models

Problem 1.3 The Gompertz (1825) equation

Consider a population that grows with an intrinsic rate of growth that decays
exponentially:
= r e atjV (L24)
IF °" -
Solve this nonautonomous ordinary differential equation. Sketch typical solu-
tion curves. What is the 'carrying capacity' and how does it differ from the
carrying capacity of the logistic equation? Where is the inflection point for
your solution? How does this differ from the logistic equation? Show that
the Gompertz equation can be rewritten as

£-•»-(£)•
where K is the carrying capacity.

The Gompertz equation was originally formulated as a law of decreasing

survivorship (Gompertz, 1825; see also Easton, 1995). It was popularized as
a growth curve by Winsor (1932). It has been used to model the growth of
plants (Causton and Venus, 1981) and of tumors (Wheldon, 1988). It also
appears in fisheries ecology in the Fox (1970) surplus yield model.

Recommended readings

Banks (1994) is encyclopedic in his treatment of logistic-like growth models.

Arrowsmith and Place (1992), Drazin (1992), and Hale and Kocak (1991)
provide useful introductions to the qualitative analysis of ordinary differential
equations. Cole (1957), Hutchinson (1980), and Kingsland (1985) recount the
fascinating histories of demography and of the logistic differential equation.
2 Harvest models: bifurcations and breakpoints

In Chapter 1, I emphasized the equilibrial properties of some simple popu-

lation models. First-order autonomous differential equations are dominated
by their equilibria. Solutions tend to stay at equilibria, move monotonically
towards or away from equilibria, or go off monotonically to infinity. This
emphasis on equilibria may seem boring. It isn't. If a model has a parameter,
and you change that parameter, equilibria may collide; they may be de-
stroyed or they may survive and exchange stability. The resulting qualitative
changes in the behavior of a dynamical system are called bifurcations.
Consider a population of fish that is growing logistically and that is being
harvested. Let us assume that the catch of fish per unit effort is proportional
to the stock level JV,

dN
(-%) -qEN (2.1)
dt V

(Schaefer, 1954; Clark, 1990). The harvest rate is the product of three terms:
the fishing effort £, a proportionality constant q that measures catchability,
and the stock level N. The product of the catchability and of the effort, qE,
is the fishing mortality; it has the same dimensions as r and will play an
important role in what follows.
We have equilibria, AT, whenever the growth rate of the fish population
equals the harvest rate:

rN* 1 = qEN" (2.2)

(see Figure 2.1). There are two equilibria,

]V* = K (1 - — ^ (2.3)
V r )
13
14 A. Single-species models

dN
~dt

qEN

N* K
Fig. 2.1. Low-mortality harvesting.

and
N* = 0. (2.4)

The second equilibrium corresponds to extirpation of the stock.

For small values of the fishing mortality, equilibrium (2.3) is asymptotically
stable. This is clear from Figure 2.1. For perturbations to the right of this
equilibrium, the harvest rate is greater than the growth rate and the stock level
decreases. For small perturbations to the left of this equilibrium, the growth
rate exceeds the harvest rate and the stock increases. In both instances,
perturbations decay and the stock level returns to AT, as suggested by the
arrows along the abscissa. Equilibrium (2.4), in turn, is unstable.
As we begin to increase the fishing mortality (Figure 2.2), equilibrium (2.3)
shifts to the left but maintains its stability. However, as we continue to
increase the fishing mortality, we eventually reach a point (Figure 2.3) where
the harvest exceeds the growth rate for all positive stock levels. Suddenly,
the equilibrium at the origin, corresponding to extinction, is stable.
I have drawn the location of the equilibria as a function of the fishing
mortality in Figure 2.4. In doing so, I have continued the first equilibrium
into the fourth quadrant. This does not make any biological sense, but it
does make eminent mathematical sense. I have represented stable branches
of equilibria with a thin line, and unstable branches with a thick line.
Branching diagrams that illustrate how the location and stability of solutions
depend on a parameter are called bifurcation diagrams. Figure 2.4 has a
transcritical bifurcation at qE = r: two branches 'collide' and exchange
 Harvest models: bifurcations and breakpoints 15

qEN

N
N* K
Fig. 2.2. High-mortality harvesting.

dN

qEN

Fig. 2.3. Severe overfishing.

stability. At this bifurcation there is a qualitative change in the behavior of

the system.
Figure 2.4 depicts changes in the equihbrial stock level with changes in the
fishing mortality. The equihbrial harvest rate or sustainable yield is equally
important. The sustainable yield at equilibrium (2.3) is

Y = qEN* = qEKll- — ). (2.5)

The graph of this function is a parabola (see Figure 2.5). Increasing fishing
16 A. Single-species models

N*

K-

Fig. 2.4. Transcritical bifurcation.

rK
T

qE

Fig. 2.5. Yield-effort curve.

effort increases sustainable yield - up to a point. Further increases in effort

lower the yield as the stock becomes increasingly overexploited and depleted.
For fixed catchabihty, the maximum sustainable yield (MSY) occurs when

(2.6)
r J
The corresponding optimal level of effort,

(2.7)
 Harvest models: bifurcations and breakpoints 17

dN
~~dt

N
Ko K
Fig. 2.6. Critical depensation.

produces the maximum sustainable yield

MSY _ f. (2.8)

Problem 2.1 Fox surplus yield model

Find the equilibria, the yield curve, and the maximum sustainable yield for
a fish population that is growing according to the Gompertz equation and
that is being harvested so that the catch per unit effort is proportional to
the stock size:
A AT / T<T \
-qEN. (2.9)

Our entire discussion has been premised on the fact that fish populations
grow logistically. However, some populations possess a threshold to growth.
Consider the simple differential equation
dN
ATfN (2.10)
—r=rN[—
dt side of\Kequation
By plotting the right-hand 0 (2.10) as a function of N (Fig-
ure 2.6), we see that there is a threshold at Ko. Solutions with initial con-
ditions above Ko approach the carrying capacity K, while those with initial
conditions below this threshold decay to zero. Since the net growth rate is
18 A. Single-species models

1 dN
~N~dt

N
K
Fig. 2.7. Allee effect.

dN

qEN

N
Ko K
Fig. 2.8. Critical depensation with harvesting.

negative at low population levels, this model exhibits critical depensation. In

addition, the per capita growth rate is no longer a monotonically decreasing
function of density but instead shows an Allee effect, or an increase in the
per capita growth rate, over certain ranges of density (see Figure 2.7).
If we add harvesting to model (2.10), we obtain

dN (N
— = rN — 1 - - ) - qEN (2.11)
dt \K0
 Harvest models: bifurcations and breakpoints 19

Stable
Unstable
qE
Fig. 2.9. Saddle-node bifurcation.

(see Figure 2.8). Equilibria occur when

(2.12)

There is a trivial equilibrium at N* = 0. The remaining equilibria satisfy the

quadratic equation

(2.13)

We can solve for N*(qE) with the quadratic formula. Alternatively, we

may plot the fishing mortality qE as a function of N* and reflect our graph
about the 45° line. Either way, we quickly generate an interesting bifurcation
diagram (see Figure 2.9).
There are three branches of equilibria (including N* = 0). The middle
branch consists of unstable equilibria that act as breakpoints (May, 1977)
between the equilibria on the two stable branches. As we increase the fishing
mortality, equilibria on the middle and upper branches approach each other.
Eventually they collide and disappear in a fold, tangent, or saddle-node
bifurcation. A fishery poised close to this bifurcation point may undergo
a catastrophic collapse after a small increase in fishing mortality. This
precipitous state of affairs is also reflected in the yield curve for this model
(see Figure 2.10). This yield curve shows a catastrophic collapse at high
enough mortality levels. Bifurcations and breakpoints pervade mathematical
ecology: small changes often have large effects.
20 A. Single-species models

qE
Fig. 2.10. Catastrophic yield curve.

Problem 2.2 Constant-rate harvesting

Consider the logistic differential equation with constant-rate harvesting:

Treat the harvest rate h as a bifurcation parameter. Find the critical value of
h for a bifurcation to occur. Which bifurcation is this? Sketch the bifurcation
diagram. What are the effects of overharvesting this system?

Problem 2.3 Some basic bifurcations

Let x be the state variable and let \i be a bifurcation parameter. Sketch
bifurcations diagrams for the following differential equations:
(1) saddle-node bifurcation

| = „ - x\ (2.15)

(2) transcritical bifurcation

- x\ (2.16)

(3) supercritical pitchfork bifurcation

- x\ (2.17)
it
 Harvest models: bifurcations and breakpoints 21
(4) subcritical pitchfork bifurcation

^ = j^x + x3. (2.18)

dt
Be sure to show the stability of each branch.

Problem 2.4 More than one bifurcation

Sketch the bifurcation diagram and characterize the bifurcations for the
following differential equation in which x is the state variable and \i is the
bifurcation parameter:

Y = fix + x3 - x5. (2.19)

Problem 2.5 A bunch of bifurcations

Sketch the bifurcation diagram for the following differential equation in
which x is the state variable and \i is the bifurcation parameter:

^ = x(9 - IIX){II + 2x - x2) \hi - 10)2 + (x - 3)2 - ll . (2.20)

at L -I
Be sure to show the stability of each branch.

All the bifurcations that we have considered are codimension 1. (The

codimension is the minimum number of control parameters that are needed
to characterize the bifurcation.) However, not all codimension-1 bifurcations
are equally robust. Small perturbations to the structure of the underlying
differential equations may affect different bifurcations in different ways.
To analyze the structural stability of a bifurcation that arises in a differ-
ential equation of the form

^ = /(AT, ,*), (2.21)

we consider the perturbed differential equation

/JN
(2.22)

where e is the small amplitude of the structural perturbation and g(N)

is some arbitrary function. We expand g(N) in a Taylor series about the
bifurcation point and keep only the lower-order terms (since we are only
interested in changes close to the bifurcation point). The function g(JV)
22 A. Single-species models

= -0.02K

Stable
Unstable

2r
Fishing mortality (qE)

AT

1=0

Stable
Unstable
I
2r
Fishing mortality (qE)

= + 0.02K

\ I
0 r 2r
Fishing mortality (qE)
Fig. 2.11. Migration as an imperfection.
 Harvest models: bifurcations and breakpoints 23
may be thought of as an imperfection. We would like to know whether this
imperfection causes qualitative changes in the bifurcation diagram.
The structural instability of transcritical bifurcations can be highlighted
by adding a small constant to harvest model (2.1),

^ = rN (l - *P\ + / - qEN. (2.23)

at \ A /
For / positive, we have immigration. For / negative, we have emigra-
tion. Figure 2.11 shows the bifurcation diagram for this equation for small
emigration (/ = — 0.02 K), no migration (/ = 0.0), and small immigra-
tion (/ = +0.02K). For small emigration, the transcritical bifurcation is
replaced by two fold bifurcations. The first fold bifurcation leads to a catas-
trophic collapse of the fishery. For small immigration, in contrast, there
is no bifurcation. Since the transcritical bifurcation undergoes such a pro-
found qualitative change in its behavior when we subject the underlying
model to a small perturbation, we consider this bifurcation to be structurally
unstable.
We have looked only at small, constant perturbations. One might also
consider small linear perturbations or, for the pitchfork bifurcation, small
quadratic perturbations. Also, a structurally unstable bifurcation may often
be made structurally stable by embedding it in a higher-dimensional extended
control parameter and phase space. However, this takes us to topics such as
the versal unfolding of bifurcations and catastrophe theory that are beyond
the scope of this book.

Problem 2.6 Imperfection theory

Consider the three differential equations

(1) x = f(x, fi,e) = e + fix - x2,

(2) x = /(x, fi9e) = e + fix - x3,
(3) x = /(x, ft, e) = e + [i — x2.

For each equation, plot the bifurcation diagram for e = —0.25, e = 0.0, and
e = 0.25. Show that transcritical and pitchfork bifurcations are structurally
unstable to small constant perturbations.
24 A. Single-species models

Recommended readings

This section borrows heavily from chapter 1 of Clark's (1990) classic book
on the optimal management of renewable resources. May (1977) wrote a
review article on breakpoints that preceded Ludwig et a/.'s (1978) extremely
important qualitative analysis of bifurcations and outbreaks in the spruce
budworm system. Wiggins (1990) may be consulted for a more thorough
introduction to bifurcation theory, codimension, and structural stability.
3 Stochastic birth and death processes

The deterministic equations

f - "<
and
N
(3.2)

ignore elements of chance that may be important in determining the growth

of populations, especially at small population sizes. Chance may enter envi-
ronmentally, due to storms, freezes, etc., or demographically, due to natural
variation in the birth and death rates about their averages. We will con-
sider some simple models that include demographic stochasticity. Stochastic
models are more informative, but harder to analyze, than deterministic
models.
Linear stochastic processes pose fewer problems than their nonlinear coun-
terparts. I will begin with a simple linear birth process that was introduced
by Yule (1924) to model the evolution of new species and by Furry (1937) to
model particle creation in cosmic ray showers. The Yule-Furry process will
act as a stepping-stone to an even more realistic (Feller-Arley) process in
which deaths are also possible. We will see that extinction may occur even
if the average birth rate exceeds the average death rate (/? > fi). We will
also derive the odds of extinction for the simple birth and death process in
terms of the average birth rate, /?, the average death rate, ju, and the initial
population size no.
I will concentrate on linear stochastic models. However, nonlinear stochas-
tic models are also important. In analyzing nonlinear deterministic models,
I focused on simple equilibria. The stochastic analog of an equilibrium is a
statistically stationary state, a limiting equilibrium distribution of population

25
26 A. Single-species models
sizes. Unfortunately, most closed (no immigration or emigration) density-
dependent birth and death processes only have the trivial stationary state. If
you cannot continue to grow, it is just a matter of time before a string of bad
luck knocks you to extinction. Even so, many nonlinear processes possess
two time scales. Over the short term, populations approach a statistically
quasistationary state. This is a stationary state for a conditional probability
distribution - conditional on no extinction. Over the long term, this quasi-
stationary state 'leaks' to extinction. In certain cases, one can determine the
mean time to extinction. I will briefly touch on one or two topics from the
theory of nonlinear stochastic processes towards the end of this chapter.

Linear birth process

The Yule-Furry process is a Markov process - the future is independent of

the past, given the present - with a continuous parameter (time). The state
space for this process is the potential number of individuals in the population
at any instant of time and, since this space is countable, it is customary to
refer to this Markov process as a continuous-time Markov chain.
As usual, I will let N(t) be the number of individuals at time t. N(t) is
now a random variable. Accordingly, I will let

pn(t) = P[N(t) = n], n = 0, 1, 2 , . . . , (3.3)

represent the probability that the population size, N(t)9 takes the value n.
I will begin by allowing births, but not deaths. I will assume that each
individual can give birth to new individuals and that each individual acts
independently of all others. For a single individual, I will also assume that

P {1 birth in (t,t + At] | N(t) = 1} = j8 At + o{At\ (3.4a)

P {>1 birth in (t,t + At] \ N(t) = 1} = o{At\ (3.4b)
P {0 births in(M + At] \ N(t) = 1} = 1 - p At + o(At). (3.4c)

Thus the probability of a birth to an individual in some small time is (to a

first approximation) proportional to that time, and the probability of more
than one birth is higher order in that time. The order symbol o(At) denotes
quantities for which

lim ^ = 0. (3.5)

For exactly one birth amongst n individuals, we must have one individual
that gives birth and n — 1 individuals who do not give birth. This can happen
 Stochastic birth and death processes 27
in n ways. Thus

P{1 birth in (t, t + At] \N(t) = n) = n [p At + o(At)] [1 - P At + o{At)]n~l

= nPAt + o(At). (3.6)

Similarly,
P {>1 birth in (t,t + At] |AT(t) = w } = o(At) (3.7)

and

P {0 births in (t,t + At] \N(t) = n) = 1 - n/3 At + o(At). (3.8)

We are now ready to write down an equation for the pn(t) of equation (3.3).
There can be n individuals at time t + At if there were n — 1 individuals at
time t and one birth occurred, or if there were already n individuals at time
t and no births occurred,

pn(t + At) = pn-i(t)P {1 birth in(t,t + At] \N(t) = n- 1}

+ pn(t)P {0 births in (t,t + At] \N(t) = n }. (3.9)

We thus have that

pn(t + At) = (n-WAtpn-iit) + (1 - nPAt)pn(t) + o(At). (3.10)

After some simple algebra, we obtain
+ A
»" ^ " ""{t) - -ntm + (. - l)^,-,«) + °^>. (3.11)
In the limit as At goes to zero this reduces to
dpn
= -nppn + (n- l)Ppn-i. (3.12)
dt
This chain of ordinary differential equations must be augmented with the
initial condition

M0) = (3J3)
\0, n+ n0.
Also, since this is a pure birth process,

pn{t) = 0, n < n0. (3.14)

Since the rate of change of pn(t) depends only upon pn(t) and on the
preceding p n _i(t), we can work our way up this chain. Say that we wish to
28 A. Single-species models

£ = 0.01
n0 = 100.0

0
0 1 2 3 4 5
Fig. 3.1. Probability of no births.

find the probability of staying at HQ. By equations (3.12), (3.13), and (3.14),
we have that

(3.15a)
dt
= 1. (3.15b)

This has the straightforward solution

pno(t) = (3.16)

(see Figure 3.1). The probability of staying at the initial condition thus
decreases exponentially with time. This probability decreases more rapidly
for large initial populations and for large birth rates.
Okay, how about the probability of being at no + 1 ? In light of equa-
tion (3.16), equations (3.12) and (3.13) reduce to

(3.17a)

(3.17b)

This linear, nonhomogeneous, first-order differential equation can be solved

using an integrating factor, leading to

l -e-l") (3.17c)

(see Figure 3.2). The probability of no + 1 individuals first increases, but

 Stochastic birth and death processes 29

Pl(t)

0.8-
/? = 0.01
0.6- no = 100.0
0.4-

0.2-

0
0 1 2 3 4 5
Fig. 3.2. Probability of one birth.

then decreases (as added births occur). It peaks at

(3.18)

This process can be continued by induction:

_ (no "+" m ~ M (3.19)

Equivalently,

—e (3.20)

for each n > no (see Figure 3.3). This is a negative binomial distribution in
which the chance of success in a single trial, exp(—/?£), decreases exponen-
tially with time.
With the probabilities pn(t) in hand, it is easy to verify that the expected
value, variance, and coefficient of variation are given by

E[N(t)] = ^ np,(0 = n0eP\ (3.21)

00

Var[iV(O] = E ^ 2 P,(0-£ 2 [iV(0]

n= 0

= wo ( 1 — e~l (3.22)
30 A. Single-species models

0.8-
^ = 0.01
0.6- no = 100.0
0.4-

0.2-

0
0 2 4 6 8 10
Fig. 3.3. Probabilities of 0, 1, 2, 3, and 4 births.

V[N(t)] = (3.23)
E[N(t)] n0
so that
lim V[N(t)] = ,/-. (3.24)

(The expected value and the variance may be derived directly, or by deriving
differential equations for each of these quantities. I will soon present an
alternative method for computing the mean and variance that uses the
probability generating function.)

Problem 3.1 Mean and variance

Derive the formulae for the expected value and the variance.

For positive /?, the mean population size, E[N(t)]9 increases exponentially.
The variance also increases. However, the coefficient of variation tends
towards a constant that depends only on the initial population size. This
coefficient is small if the initial population is large. If the initial population
size is instead small, any single run of the stochastic process is likely to differ
greatly from the corresponding deterministic process.
Figures 3.4 and 3.5 show two sets of 10 simulations of the birth process
 Stochastic birth and death processes 31

0 100 200 300

Fig. 3.4. Ten realizations of a birth process with large

N(t)
50—|

40-

30-

20-

10-

0
0 100 200 300

Fig. 3.5. Ten realizations of a birth process with small UQ.

for the same parameter (/? = 0.01) but for two different initial conditions
(no = 50 and no = 5). The second set clearly has more variability.

Problem 3.2 The simple death process

Construct and analyze a continuous-time Markov chain in which individuals
persist until they die. Do not replace dead individuals. Assume that an
32 A. Single-species models

individual alive at time t dies with probability fiAt + o(At), in (£, t + At]
and that each individual acts independently of all others.

Problem 3.3 Extinction time

Assume that there are HQ individuals present at the start of a simple death
process. Let T be the extinction time for this process. What is the probability
density function for T? What is the mean time to extinction?

Linear birth and death process

The obvious drawback to the simple birth process is that individuals never
die. Life, in contrast, is risky: some individuals die before siring young.
Accordingly, we will assume that for each individual
birth in (t, t + At] \N(t) = 1} = j8 At + o(At), (3.25a)
P{1 death in(t,t + At] \N(t) = 1} = \iAt + o{At\ (3.25b)
F{no change in (t,t + At] \N(t) = 1} = 1 - (j8 + p) At + o(At). (3.25c)
There now is some probability of dying as well as of giving birth. The
probability of several events (births and/or deaths) is taken to be o(At).
Individuals are still assumed to be independent and so, amongst n individuals,
P{lbirth in(t,t + At] \N(t) = n} = n\fiAt + o(At)]
x[l-(P + n)At + o(At)] n~l. (3.26)
As a result,
P {1 birth ]n(t,t + At] \N(t) = n) = n/3 At + o(At).
Similarly,
P {1 death in(t,t + At] \N(t) = n) = njuAt + o(At) (3.27)
and
P {no change in(t,t + At] \N(t) = n) = 1 - n(j8 + fi)At + o(At). (3.28)
With these probabilities in hand, we are now ready to derive an equation
for the probability, pn(t), that the population is of size n. There can be n
individuals at time t + At if there were n — 1 individuals at time t and
one birth occurred, if there were n + 1 individuals at time t and one death
 Stochastic birth and death processes 33

occurred, or if there were already n individuals at time t and no births or

deaths occurred,

pn(t + At) = (n - l)PAtpn-i(t) + (n +

+ [1 - n(P + n)At]pn(t) + o(At). (3.29)

After some simple algebra and after taking the limit as At goes to zero, we
obtain

*jj± = P ( n - l ) ^ - i + n(n + l ) p B + i - ( P + v)npn. (3.30)

We must augment this differential equation with the initial condition

(331)

and with the stipulation that pn(t) = 0 for n < 0.

You may wonder why we bothered with the simple birth process when the
(substantially more realistic) birth and death process is not much harder to
derive. It is not much harder to derive, but it is much harder to solve. For
the simple birth process, the rate of change of pn(t) depends only on pn(t)
and on the preceding pn-\{t). For the birth and death process, this same pn(t)
depends not only on pn(t) and pn-i(t), but also on the, as yet, unsolved-for
Pn+l(t).
The cure for this problem is to solve for all the pn(t) as one. The easiest
way do so is to introduce the probability generating function

F(t,x) = ^Pnit)^. (3.32)

n= 0

(Alternatively, one can introduce a moment generating function.) Since the

probabilities pn(t) are functions of time, F(t, x) is a bivariate function.
However, we will often think of F as a univariate function of x, with
the probabilities as the coefficients in the power series expansion of this
function. Since every pn(t) is less than or equal to unity, the generating
function converges for all | x | < 1 (by direct comparison with the geometric
series). Since the pn(t) sum to 1, F(t, x) also converges for x = 1. There are
a number of other properties of the probability generating function F(t, x)
that are worth highlighting:

(1) The probability of being extinct at time t, po{t\ is given by

po(t) = F(t, 0). (3.33)

34 A. Single-species models
(2) Taylor's theorem permits us to expand a function in terms of its derivatives at
zero. Thus, by calculating derivatives, we can find the sequence of probabilities
associated with a given probability generating function,
1 8nF
(3.34)
x=0
(3) The probability generating function allows us to compute the average or
expected value of N(t) without tedious calculations involving discrete sums:

E[N(t)] = (3.35)

(4) The probability generating function also allows us to compute the variance of
N(t) with only a little more effort. In particular,
d2F
(n2 -n)pn(t)
dx2 x=l

= E[N2(t)] - E[N(t)]. (3.36)

However,
Var[JV(0] = E[N2(t)] - E2[N(t)]9 (3.37)
so that
\d2F d_F_ 8FY
Var[JV(*)] = (3.38)
'dx2 ~dx dx)

We can, in other words, compute all the probabilities and statistics that we
need in a straightforward way with the probability generating function.
So, how do we find this wonderful function? Well, the probability gener-
ating function F(t, x) arises as the solution to a partial differential equation.
If we differentiate the generating function, equation (3.32), with respect to
time,
00
dF
(3.39)

and make use of differential equations (3.30) for the derivatives of the
probabilities,
^F oo oo oo

nx", (3.40)
n=0 n=0 n=0
d_F_

n=0 n=0

(3.41)
n=0
 Stochastic birth and death processes 35

we arrive at the partial differential equation

dt dx
Equivalently,

h (j6x - /i)(l - x ) — = 0. (3.43)

This partial differential equation has the initial condition

F(0,x) = xnQ (3.44)

since the population starts at no with probability 1.

Problem 3.4 Moment generating function

The moment generating function, M(t, 9), is defined as
00
ne
M{t,6) = T Pn(t)e . (3.45)

Show that
dM
(3.46)

Var [N(t)] = ^_m\2\ , (3 .47)

and that the moment generating function for the linear birth and death
process satisfies

=[«/_,) + „,,-._,„ ^ . (3.48)

Mathematical meanderings
To solve partial differential equation (3.43) we must use the method of char-
acteristics. Consider the change of variables (t, x) —> (u,v). We will use u as
an independent variable that measures distance along a characteristic curve
and v to parametrize the initial conditions. The value of v will specify which
characteristic curve we are on. Along each characteristic curve (for fixed v), F
will depend solely on u.