F
neither
Of indix.
as its b
native t.
“Story 1
have beer

Thi
individua
“99 profi:
fuflction ‘

Solu
“mam

 

logical ac

 

distribut‘.
hrtalllty I

1

W1 t3

 

ABSTRACT
THE THEORY OF PHYSIOLOGICAL AGE
DISTRIBUTIONS IN BIOLOGICAL POPULATIONS
By

John Van Sickle

For many species of plants and animals, chronological age is
neither convenient nor useful to the population biologist as a measure
of individual maturity. A physiological feature of an organism, such
as its body size or chemical composition, provides an attractive alter-
native to chronological age as a frame of reference for discussing life
history phenomena. Unfortunately, most demographic population models
have been structured on the basis of chronological age.

This thesis is a study of a dynamic population model in which
individuals are classified by physiological age. The physiological
age profile of the population is represented by a number density
function which obeys a first-order partial differential equation.

Solutions of the differential equation are presented. The
asymptotic behavior of the model is discussed and well-known stability
results for chronological age distributions are generalized to physio-
logical age distributions. The shape of the stable physiological age
distribution is shown to be a function of individual growth as well as
mortality and fecundity.

With body size as a specific measure of physiological maturity,

 

 

 

 

Equatio
equatio:
distrib:
previous
tion she
be used
freqlean
Pi:
context c
available
analysis

regulate

 

 

John Van Sickle

the model is applied to some problems in fish population dynamics.
Equations for estimating biomass production are derived from the dynamical
equations for the size distribution. A steady-state pOpulation size
distribution is defined and applied to an exploited fish population
previously modeled on the basis of chronological age. A third applica-
tion shows how the stable size distribution predicted by the model can

be used to compute mortality rates from empirical population size-
frequency curves and a knowledge of individual growth rates.

Finally, time-varying growth rates are discussed, still in the
context of fish populations. Fish are known to respond to changes in
available resources by altering their growth rates. A stability
analysis and simulations demonstrate how individual growth in fish can

regulate pOpulation growth in a limited environment.

THE THEDRY OF PHYSIOLOGICAL AGE

DISTRIBUTIONS IN BIOLOGICAL POPULATIONS

BY
Q“ 0‘“
John van Sickle

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR.OF PHILOSOPHY
Department of Electrical Engineering and Systems Science

1975

DEDICATION

To Jane, with love.

 

 

I at
friendship
dynamics.
of my com:
K. Y. Lee,
K- Y- Lee
Systems, J
other ecol.
were Of 9:.

I we
and the Na

during the

 

ACKNOWLEDGMENTS

I am deeply grateful to my advisor Erik D. Goodman for his
friendship, constant encouragement, and unerring insight into system
dynamics. For their interest and advice, I thank the other members
of my committee: Drs. R. H. Boling, Jr., R. O. Barr, D. J. Hall,

K. Y. Lee, and H. E. Koenig. I would especially like to thank Dr.
K. Y. Lee for his guidance in my early study of distributed parameter
systems. Thanks, too, to Drs. Patricia and Earl Werner and my many
other ecologist friends and colleagues whose long discussions with me
were of great benefit to this thesis.

I would also like to thank Drs. H. E. Koenig and W. E. Cooper
and the National Science Foundation (Grant GI-ZO) for financial support

during the past three years.

iii

Chapter

TABLE OF CONTENTS

INRODUCTIW O O O O O O O . O O O O O O O O O O O O 0

THE MODEL: DEFINITION AND BACKGROUND . . . . . . .

Physiological Age . . . . . . . . . . . . .
The POpulation Balance Model . . . . . . . .
Relations to Other Models . . . . . . . . .
Solutions for Time-Invariant Parameters . .
Interpretation of Solutions . . . . . . . .

ASYMPTOTIC PROPERTIES--THE STABLE P-AGE DISTRIBUTION

Asymptotic Properties of C-Age Models . . .
A Renewal Equation for B(t), the Birth Rate
Relations in an Exponential Population . . .
The Shape of the Stable P-Age Distribution .
The Stable P-Age Distribution--Discussion .

APPLICATIONS TO FISH POPULATION DYNAMICS . . . . . .

Introduction . . . . . . . . . . . . . . . .
Biomass and Production . . . . . . . . . . .
The Growth Function, 9(2), for Fish

Populations . . . . . . . . . . . . . . . .
Recruitment and Seasonal Periodicity . . . .
Application to a Specific Model . . . . . .
Catch Curves and Mortality Estimation . . .

THE REGULATION OF POPULATION GROWTH BY INDIVIDUAL

GROWTH

BIBLIOGRAPHY . .

APPENDIX

Density-Dependence in Fish Populations . . .
Stability Theory . . . . . . . . . . . . . .
Simulation-~Goals and Methods . . . . . . .
Simulation Results . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . .

Page

ll
14
18

26
26
28
31
34
38
43

43
47

52
57
63
67
73

73
77

89
99

103

108

Table

5.1

5.2

 

 

 

Table

5.1

5.2

LIST OF TABLES

Model Population Features as a Function of Life
History Parameters . . . . . . . . . . . . . . . . . . . 96

Settling Times for Representative Model Populations.
Mortality Rate of .116 Results in 50% Loss/Yr., While

11 = .231 Results in 75% Loss/Yr. Settling Times are
MeasuredinGenerations................ 98

Figure

2.1

2.2

2.3a

 

2.31)

2.4

 

3.1

4.1

4.2

4.3,
4.31.

4.4

4.5

4.5

 

Figure

2.1

2.3a
2.3b

2.4

3.1

4.1

4.2

4.3a

4.3b

4.4

4.5

4.6

LIST OF FIGURES

Physiological age (p-age) profile at a point in time,
t1, showing reproductive feedback and mortality loss . .

Trajectories for the von Foerster equation . . . . . . .
Characteristic curves for sigmoidal growth . . . . . . .
Growth rate curve for sigmoidal growth . . . . . . . . .
Trajectories of idealized cohort: a) Case 1 - Concave
growth function and high mortality at all p-ages. b)
Case 2 - Concave growth function and high mortality up

to 21, low mortality beyond 21 . . . . . . . . . . . . .

Hypothetical p-age frequency distribution for stationary
mpulation O O O C O O O O O O O O O O O C O O O O O O O

Idealized age distribution at fixed point in time,
showing year classes arising from seasonal reproduction

Idealized size distribution of fish showing overlapped
year classes due to decreasing growth rate . . . . . . .

von Bertalanffy growth curve . . . . . . . . . . . . . .
Growth rate for von Bertalanffy growth curve . . . . . .
Monthly series of length frequencies for the New South
Wales Tiger Flathead population in the east Australian
trawl fiSherY O O O O I O O O O O O O O I O O O O O O O
Steady-state size distributions associated with a

theoretical fishery model. Parameters from Beverton
and Halt [1957] O O O O O O O O O O O O O O I O O O O 0

Catch curve for the Pacific Halibut population, 1925-1926.

Abscissa--mean length of successive 5 centimeter length

groups, in millimeters. Ordinate--logarithm of the number

of fish taken at each length interval . . . . . . . . .

vi

Page

19
20

20

22

36

46

54

56

56

60

65

68

 

 

Figure

4.7

5.1

5.2

5.3
5.4

L1

58.!

Sa:

Re!

Figure

4.7

5.1

5.2

5.3

Length-frequency curve for whitefish in Shakespeare
Island Lake, Ontario, showing mode at 17-18 inches

due to decreasing growth rate and low mortality in
that length range . . . . . . . . . . . . . . . . . . .

Growth of fish population in a limited environment.
The population initially has only a few individuals . .

Net reproductive rate versus density for reproductive
melg. O O O I 0-. O O 0' O 0‘. O O O O O O 0 CI. 0 0

Sample trajectories of model population growth . . . .
Sample trajectories of model populations . . . . . . .

Rectangular grid for centered difference scheme . . . .

vii

Page

71

78

90

92

109

CHAPTERl

INTRODUCTION

Quantitative models with a single dynamic variable are inadequate
fbr a realistic description of most populations. If functional or time
lag terms are included, ordinary differential or difference equations for
total numbers can reproduce many observed population trajectories.
However, such complex equations usually can not be interpreted in a
biologically useful manner.

Much of the problem lies in the assumption of these models that
all members of the population are identical. Numerous field and labora-
tory populations have demonstrated how the heterogeneous composition of
a population.profbundly affects its dynamics. Population models must
begin to account for this hetereogeneity.

Oster and Takahashi [1974, p. 483] express this in a different
way by stating that ". . .populations are intrinsically distributed
parameter systems.“ For populations whose members are differentiated
by chronological age, a highly developed set of distributed parameter
models is available which includes integral equation, partial differ-
ential equation, and discrete matrix forms. Using this set of models
as a theoretical base, demographers have developed a comprehensive
methodology for age-specific parameter estimation and prediction of
future age distributions [Keyfitz, 1968] .

There are many ways to classify individuals other than age. Such

1

features as individual body size, dietary requirements, or chemical
composition might strongly influence how individual birth and death
phenomena "sun up" to give population birth and death rates.

Alternate methods of differentiating among individuals become
especially important when describing non-human populations. Chrono-
logical age is difficult to assess for a great many species of plants
and animals, and it is usually more convenient to keep track of other
physiological characteristics of individuals.

Even if chronological age can be determined, it is often an
unreliable indicator of an organism's sexual maturity, fertility, or its
chances for future survival. For most species of invertebrates, plants,
and fish, these vital parameters of population growth are strongly
dependent on environmental factors such as available food or ambient
temperature. For such populations, a physiological, rather than
chronological, measure of age or maturity may lead to more realistic
and useful models of population dynamics.

This thesis is a study of a distributed parameter model in which
individuals are classified by physiological age. It is assumed that
these individuals live in a spatially homogeneous population and that
imigration and emigration are negligible. Sex differences among
individuals are also ignored, since sex ratios tend to remain fairly
constant in most populations.

The study has two principal Objectives. The first objective is
to define and analyze a continuous-time dynamical system for the
physiological age distribution of a population. Various forms of the
system equations which will be used have already been applied to specific

physiological characteristics of specific populations, but the system

will be analyzed in the more general context of an arbitrary measure of
physiological age. Chapter 2 of the thesis defines the model and gives
a brief review of its previous applications to biological populations.
Solutions of the system equations in their autonomous form are presented
and interpreted.

Chapter 3 examines the asymptotic (in time) behavior of the
autonomous model. A. J. Lotka's familiar stability theory of chrono-
logical age distributions is generalized to apply to physiological age
structures.

The second main objective is to thoroughly understand the crucial
role played by growth in population dynamics. In this thesis, "growth"
refers to the dynamic relationship between an individual's physiological
age and its chronological age. Individual growth rates become vital
parameters along with mortality and fecundity rates. In addition,
individual growth is the mechanism through which many population systems
are coupled to their environment.

In Chapters 2 and 3, we see how time-invariant growth rates
affect the rate of increase of a population and the shape of its
physiological age distribution.

Chapter 4 is an application of the theory of Chapters 2 and 3 to
fish populations. Individual growth in length or weight has been long
recognized as a key factor in fish population dynamics, and body size
is employed as a particular measure of physiological age in applying
the model to fish. The chapter briefly discusses growth in fish and
then goes on to attack three quite different problems of observing and
modeling fish populations which are assumed to have time-invariant

individual growth rates.

Chapter 5 retains the general context of fish populations classi-
fied by body size, but the emphasis shifts to time-varying growth rates.
Specifically, we study the population in relation to environment by
assuming individual growth is a function of population density. This
growth response to limited food or space is felt by many fishery
biologists to be a powerful regulator of fish population densities. A
perturbation analysis and computer simulations of the model show the
regulatory properties of a flexible growth response.

The arguments of this thesis will show that physiological age
models provide the added dimension of individual growth to population
dynamics at little or no cost in mathematical tractability or ease of

interpretation.

CHAPTER 2

THE MODEL: DEFINITION AND BACKGROUND

2.1 Physiological Age
Suppose 2 represents some measure of maturity or physiological
development of individuals in a population. The variable 2 could be
chronological age, body size, chemical composition, or any other physio-
logical feature having a bearing on individual reproduction or mortality.
The physiological age, z, of an individual is dynamically described

by a differential equation:

A
3—3.: g‘z’t). (2.1-1)

Here, g(z,t) is the growth rate of the physiological age variable, 2.
The term "physiological age" will be abbreviated by ”p-age", and here-
afterz will be referred to as a p-age variable.

Entomologists, notably Hughes [1962] and Lefkovitch [1965], were
among the first to use stages of maturation, rather than chronological
age, in structuring their models. This approach was partly motivated,
no doubt, by the discrete, easily recognizable, life stages seen in
many insect species.

Other researchers have been attracted to the concept of p—age by
the relative ease with which p-age, as opposed to chronological age,
can be determined in many organisms. Examples are the early fish
pOpulation model of Baranov [1918] and the model of Usher [1966] for

forest management.

More recently, the need for generalized age-structured models,
together with the rise of computer simulation as a tool of the population
biologist, has fostered the theoretical development and implementation
of continuous [Striefer and Istock, 1973; Streifer, 1974; Weiss, 1968]
and discrete [Boling, 1973; Coulman, Reice, and Tummala, 1972] models

employing "critical" state variables based on p—age.

2.2 The ngulation Balance Model

 

We now introduce a number density function, N(z,t) with the meaning
2
1 ~ ~ I
that g N(z,t)dz is the total number of individuals in a population

2
between the p-ages of 21 and 22 at time ta N(z,t) is thus a dynamic,

p-age number distribution.

Its equation of motion obeys the conservation law:

8N ,t a
5141—)- + 3; [g(z,t)N(z,t) ]= -u(z.t)N(z,t) . (2.2-l)

Here, g(z,t) is defined by (2.1-1), and u(z,t) is the instantaneous
loss rate of individuals from the population due to predation and
natural mortality.

Equation (2.2-1) is a "balance" equation. It says that the rate
of growth of individuals through p-age 2, given by ‘%;'[§(z,t)N(z,ti],
and the loss rate of individuals of p-age 2, given by [u(z,t)N(z,t)],

are balancedfiby the net rate of change of numbers at p-age 2, expressed

3N(z,t)
3t '

To describe the dynamics of a reproducing population, Equation

by

(2.2-l) must be accompanied by a reproductive feedback which will appear

as a boundary condition. Let the birth rate of individuals of p—age z

at time t be denoted by B(z,t). Then, the birth rate can be defined

as the functional:

zb+a

B(z,t) = f b(z,2,t) N(z,t) d2. (2.2-2)

In this equation, b(z,z,t) is the rate at which individuals of p—age 5
give birth to neonates of p-age 2. 2b is the p—age at which individuals
first reproduce, and a is the width of the breeding interval.

Equation (2.2-2) is the most general form for a reproductive
feedback, but when it is coupled with (2.2-l) it appears as a distributed
boundary condition, and the resulting system equations yield little to
analysis. For most biological populations, with the notable exception
of cell populations, a restricted form of (2.2-2) will entail only a
minor loss of realism. We will hereafter assume that all neonates have
the same p-age, denoted by z . Then (2.2-2) becomes

0

zb+a
é b(E,t) N(E,t) dE. (2.2-3)
b

HD'

B(t)

Now B(t) is defined as the total instantaneous birth rate, and b(z,t) is
the p-age-specific individual fecundity.

The connection between the reproductive functional, B(t), and a
boundary condition for (2.2-l) is found by integrating (2.2-1) over all
p-ages to get

a Q ~ ~ m cm -
'5: £0 N(z,t)dz] + [9‘ ot)N( pt) - 9(zo,t)N(zo,t):l _

- f” u(E,t)n(§,t)di. (2.2—4)
20

8

. A ‘w ~ ~
Define Nt(t) = g N(z,t)dz as the total nunber in the population. Also,
0

assume that no organisms survive past some upper limit in p—age so that

9(°°,t) N(°°,t) = 0. Then (2.2-4) becomes

dNt(t) m
_...._dt = g(zo,t)N(zo,t) - £0 u(z,t)N(z,t)dz. (2.2-5)
Since 6"t

dt = B(t) - D(t), where B(t) and B(t) are population birth and

death rates, respectively, we must have
g(zo,t) N(zo,t) = B(t). (2.2-6)

Using (2.2-3) the boundary condition is also written as
+0

g(zo,t) N(zo,t) = éb b(2,t) N(E,t)d2. (2.2-7)

With reasonable restrictions on g(z,t), to be discussed shortly, and an
initial condition, N(z,0) = NI(z), the system (2.2-l), (2.2-7) prescribes
a unique solution surface above the (p-age)x(time) plane for the number
density N(z,t)-

Figure 2.1 is a sketch of a p-age distribution, N(z,t1), at an
instant in time, t1, showing the relationships implied by (2.2-1) and
(2.2-7). Similar sketches are drawn by Oster and Takahashi [1974].
Derivations of various forms of the balance Equation (2.2-1) in the
context of pOpulations can be found in several papers [Von Foerster,
1957; Sinko and Streifer, 1967: Weiss, 1968: Levin and Paine, 1975;
Lee, et al., 1975].

For the sake of tractability, an important restriction must be

made on the growth rate in (2.2-1). It is necessary that g(z,t) be

A A

 

 

 

 

 

 

 

 

 

 

 

 

u--Mortality
NUMBERS ‘,""”"""'”""""""" ""“ """""""’\\‘
N(z,t1) 1’ \\\
‘ D I
I
l
I
i
I
I
l
zO 2b zb+a P-AGE,z
B(tl) fb'Ndz
Reproduction

 

 

 

Figure 2.1. Physiological age (p—age) profile at a point in time, t

. 1'
showing reproductive feedback and mortality loss.

10

everywhere positive, so that an individual's p-age is continually
increasing. This restriction, together with the restricted description
of reproduction seen in (2.2-3), results in a one-to-one correspondence
between any individual‘s p-age and its chronological age. This cor-
respondence is necessary for the validity of the analytical solutions
to (2.2-1) which are discussed later in this chapter. The restriction
that individual growth must be positive results in a less realistic
model, especially if 2 represents body weight. Fortunately, the
restriction is necessary only for analytical solutions. Most detailed
applications of the model for predicting population trajectories would
require simulation, in which case negative growth could be easily
handled.

The vital statistics 9, b, and u are in general dependent upon
both the internal p-age structure of the population and external environ-
mental factors. The effects of intraspecific competition can be expressed
by allowing 9, b, and p to vary with the state function, N(z,t). The
inclusion of density-dependence in (2.2-l) and (2.2-7) is discussed
in Chapter 5.

It is also possible, in special cases, to factor out the explicit
time dependence of g, b, and D. For example, developmental rates in
insects and other poikilotherms are closely tied to temperature, and
this has prompted entomologists to discuss growth relative to a physio-
logical time scale based on degree days or "accumulated heat units."

Hughes [1962] used this concept to model the instar or life—stage
structure of an aphid population. Lee,et a1. [1975] model the cereal

leaf beatle and its parasites using an equation of the same form as

11

(2.2-l). They base their growth, mortality, and reproductive parameters
on a physiological time scale, so that, for example, the equation for
growth,(2.1—1L appears in their model as 9§-= g(z), where E represents
degree-days. dt
2.3 Relations to Other Models

The balance model is easily generalized to handle characterization

of individuals by more than one feature. In this case, 2 becomes a

vector, 2, with
z =Ezl,age; 22,1ength; Z3,weight;. . .zmj .

Equation (2.2-l) is expanded to the form

3N(E,t) m _a_ - - _ - _ _
3t 4' 1:1 321 Eli‘mfimmtfl - -u(z,t)N(z,t). (2.3 1)

Here, gi(z,t) is the growth rate of characteristic 21 [Oster and
Takahashi, 1974].

Sinko and Streifer [1967; 1969; 1971] discuss the model (2.3-1)
with z = chronological age and z

1 2
tions of Daphnia and the planarian worm Dugesia tigrina. In a model of

= body mass, and apply it to popula-

 

shrimp in an estuary, Billups, Wilson, and Pike [1971] apply (2.3-1)
to shrimp classified by body length and three spatial coordinates,
resulting in a four-dimensional state vector, N.

An interesting application of (2.3-1) to a quite different ecolog-
ical system is proposed by Levin and Paine [1974; 1975]. They use the

density function N(;,t) to represent component islands, or "patches"

12

within a heterogeneous natural community. Each patch is assumed to be
homogeneous with respect to biological properties, e.g., species com-
position, and the growth rate of patches is a function of intra- and
inter-patch interactions and localized environmental disturbances.

In this fashion, Levin and Paine hope to examine the long-ignored role
of spatial heterogeneity in structuring communities.

If, in (2.1-l), we let p—age, 2, be equivalent to chronological
age, a, then g%-= 1 = g. Hereafter, the abbreviation "c-age" will be
used for chronological age, a. Now (2.2-l) reduces to a simpler
equation:

3N(a,t) + 3N(a,t)

3t 3a = -u(a,t)N(a,t). (2.3-2)

This form of (2.2-1) is known as von Foerster's equation, after H.

von Foerster [1959] who used it to model cell populations. Since von
Foerster's work, balance models have been studied extensively in the
context of cell population dynamics. Trucco [1965] and Nooney [1967]
investigate (2.2-l) in some detail. Weiss [1968] and Bell and Anderson
[1967] present more generalized models of the form (2.3-l) for cell
populations. Bell and Anderson, for example, classify cells according
to the c-age and volume. For mitotic cell populations with individuals
classified by size, distributed boundary conditions must be used with
balance models since daughter cells at birth are assumed to be one-half
the size of their progenitors. This fact, together with the common
assumption by cell modelers that cells experience no natural mortality
before division, results in partial differential equation models whose

dynamics are quite different from those of the system (2.2-1), (2.2-7).

13

An excellent application of von Foerster's equation to host-parasite
systems has been presented by Oster and his coworkers [Oster and Takahashi,
1974; Auslander, Oster and Huffaker, 1974]. A variety of system dynamics
are discussed in relation to such parametric pr0perties as synchronous
maturation delays in host and parasite populations and periodic time
functions for mortality and fecundity.

A more complete review of applications of distributed parameter
models to various populations is found in Streifer [1974]. Because
of the complexity of the generalized model (2.3—l), a detailed analysis
of its dynamics will not be attempted here.

If the c-age/time continuum is discretized, then the well-known
matrix methods developed by Lewis [1942] and Leslie [1945] are appli-
cable. Sinko [1969] Shows the equivalence, in the limit as At + 0, of
the matrix methods and the von Foerster equation.

Discrete models based on p-age are a more recent development.
Lefkovitch [1965] subdivided an insect population into life stages and
defined a vector fi(t) whose component fii(t) gives the number of individ-
uals in the ith stage at time t. He then described a projection matrix
M; the (i,j) element of M is the proportional contribution of the numbers
in the ith stage at time t to the numbers in the jth stage at time
(t + At). Thus, n(t + At) = M 0 B(t). Lefkovitch recognized that the
numbers in the matrix M represented the effects on n of mortality during
the time interval [t, t + At], confounded with the effects of growth
during the same period. Because of this complication, he was unable
to give a clear biological interpretation of the matrix entries.

A more complex discrete p-age model was develOped by Coulman,

14

Reice and Tummala [1972]. In this model, stage-specific growth and
mortality rates are considered separately, as is done continuously in
(2.2-l). The equivalence in the limit as At +'O is shown in Kharkar,
Cbulman and Barr [1974].

A stochastic version of (2.2-1) permits individual differences in
growth rate due to genetic variability or to microenvironmental variations
arising from the spatial distribution of organisms. A term is added
to (2.2-l) to reflect the dispersion in p-age of individuals, so that
the model takes the form

3N( z,t) + 3[g(z,t)N(z, t)]
3t 32

 

2
%”%;z [v(z,t)N(z,t)] =
-u(zlt)N(zlt)o (2.3-3)

In (2.3-3) the individual growth rate is a random variable with mean
g(z,t) and variance v(z,t). Higher moments of the probability distri-
bution associated with growth are assumed negligible in their effects

on population dynamics. The diffusion equation (2.3-3) and its applica-
tion to populations is discussed by Bailey [1964]. A derivation can be

found in Kharkar [1973] or Lee, et a1. [1975].

2.4 Solutions for Time-Invariant Parameters

As a first step, assume that (2.2-l) and (2.2-7) describe a popu-
lation whose growth, mortality and fecundity rates are independent of
time. The system (2.2-l), (2.2-7) is then autonomous and, for specific
functions 9, u, and b, its solution depends only on an initial p-age

distribution, NI(z).

Equation (2.2-l), upon application of the product rule for

15

differentiation, becomes

3N(z,t)
at

+ g(z) 2E§§LEL = - [g'(z) + u(z)] N(z,t), (2.4-l)

where:

, A dg(z)
g (2) dz .

By the method of characteristics [Courant and Hilbert, 1966],

(2.4-l) is equivalent to the set of differential equations

dt = 9(2) -[g'(z)+u(z)]

 

‘ (2.4-2)

To facilitate writing solutions to (2.4—2), define

A Z2 d2
= ——~—-° 2.4-3
Then T(zz,zl) is the time required to grow from p—age 21 to p-age 22.

Since we assumed that all neonates are of p-age z

the quantity T(z,z )

0'

is simply the chronological age, a, of individuals of p-age z.

0

It was required that g(z) > 0 for all values of 2, so T(zz,zl)
is a monotonically decreasing function in 21 for any fixed 22, and it

thus has an inverse function in 21, defined as
A A '1
z(zz,T) — T (22,21). (2.4 4)

This definition means that 2(22,T) is the p-age at which an organism
must be at time (t - T) in order to grow to p-age 22 at time t.
Now, using the functions T,§, and the initial condition N(z,O) =

NI(z), we have the solution to (2.4-2):

16

f

= A , z 9'(E)+u(2) ~ _
N(z,t) NI(z(z,t)) expl: g(z,t)l 9(2) ]dz (2.4 4a)

 

for all p-ages, 2, such that T(z,zo) 2.t' and

 

2 ~ ~
= _ _ 9' (z)+}1(z) ~ _
N(z,t) N(zo,t T(z,zo)) exp]: £0 [ g(z) ] dz](2.4 4b)

for all p-ages, 2, such that T(z,zo) < t.

The solution thus has two parts. The first part, (2.4-4a),
describes that portion of the p-age distribution, N(z,t), whose members
are survivors of the initial population, NI(z). The second equation,
(2.4-4b), characterizes the portion of the distribution whose members
were born after time t = O, and entered the distribution via the
boundary,N(zO,-).

The exponential terms in Equations (2.4-4) can be simplified as

follows:
2 , ~ 2 ~
exp — f 2 SFigl-dfi - 2 Elél-di =
z1
9(2 ) 2 ~ ‘
l xp -f 2 11(2) d2 (2.4-5)
g(zz) 21 9(2) °

Now, the exponential factor on the right-hand side of (2.4—5) will be

redefined as the p:age-specific survivorship. That is,

 

2 ~
_ 2 v(z) ~ A _

where 2(22,zl) is the proportion of organisms alive at p-age 21 which

survive to reach p—age 22 after being exposed to the mortality rate ”(2).

Using this simpler notation, the solution (2.4-4) becomes

= ,. g(fi) A _
N(z,t) NI(z(z,t)) g(z) £(z,z), (2.4 7a)

for T(z,zo)‘z_t, and

0)

N(z,t) = N(zo,t - T(z,zo)) 37:7-

£(z,zo) (2.4-7b)

for T(z,zo) < t.

To illustrate the solution equations, take the simple case of

constant growth rate, i.e., g(z) = G. Then

 

2 1 ~
T(22,zl) —£ '6' dz "" A t

and g(z,t) = z - gt, since fit is the increment in p-age during the time
interval (0,t).

The solution (2.4-7) then becomes

 

 

N(z,t) = NI(z—§t) £(z,z-§t), (2.4-8a)
2-2
for ‘3 t, and
= - l -
N(z,t) N(zo,t [29 :l) (2,20), (2.4 8b)
z-z
for < t.

 

If 20 a 0, this case of constant growth rate is equivalent to a c-age
model which uses time units of fit instead of t, and the Equations (2.4-8)
become identical to those given for the von Foerster equations by Trucco

[1965]. If u a O, the Equations (2.4-8) are identical to the results

in Lee, et a1. [1975].

18

2.5 Interpretation of Solutions
To better understand the nature of the solutions (2.4-4) we will
sketch the characteristic curves, defined by (2.4-2), for two cases.

Figure 2.2, reproduced from Skellam [1967], shows a typical family of

trajectories for the von Foerster equation where z a = c-age. In

1 and 20 = 0.

this case, the solution is given by (2.4-8) with g

The Equation (2.4-2) reduces to

 

dt = da = dN (2.5-1)
-u(a)

In Figure 2.2, the family of oblique lines in the age-time plane
of slope-g%-= l are life-history trajectories, and they are also the
characteristic curves of the partial differential Equation (2.3-2).
This family of 45° lines in the age-time plane is known to demographers
as a Lexis diagram [Keyfitz, 1968]. Along a life history trajectory,

the number density N decreases in accordance with the equation

 

= da. A pulseéwave type initial condition is shown along the
(t = O)-age axis, and a sinusoidal birth rate, B(t), provides the
boundary condition along the (a = O)-time axis.

Let us contrast the curves in Figure 2.2 with a common, but more
complex, function for 9. Suppose the p-age variable 2 represents body
weight. Growth in body weight is described, for many organisms, by a
sigmoidal or logistic weight-for-age function, where body weight 2 tends
toward an upper limit, 2m, as t +-m [Thompson, 1942]. Figure 2.3a shows
a family of sigmoidal growth curves. The growth function g(z) = 3%-

gives the slope of any curve at a particular body weight, 2. It can be

seen that g(a) must have the Shape shown in Figure 2.3b. It is generally

19

3(9)

 

 

 

 

Figure 2.2. Trajectories for the von Foerster equation. From
Skellam [1967].

20

   
 

 

 

 

A
z
m .__.._.
Body
Weight,.g
= 9(2)
2 :=
0 Time
a)
Growth Rate
dz
g(z) -- dt
1 1
U I
20 21 Z2 200
b) Body Weight, 2

Figure 2.3. a) Characteristic curves for sigmoidal growth.
b) Growth rate curve for sigmoidal growth.

21

a concave function of 2, first increasing, then decreasing to zero
as 2 + z”.

The family of curves in Figure 2.3a is analogous to the set of
Oblique straight lines in the age—time plane of Figure 2.2. The curves
in Figure 2.3a are the characteristic curves defined by Equation (2.4-2)
employing a concave growth function g(z).

In order to solve (2.4-l) when it contains a growth function g(z)
which approaches zero for a finite value 2w, it is necessary that
[g'(z) + 11(2)] > 0 as 2 + zoo so that the exponential terms in the solu-
tion (2.4-4) do not become unbounded as 2 + a”.

Figure 2.2 illustrates the dynamics of N(z,t) for the simple case
_g = 1 and z = a. What dynamics might be seen for N(z,t) if a concave
function like that in Figure 2.3b is used for g(z)? Let us examine the
progress of a simple idealized population density, N(z,t), through
a p-age distribution. Assume that 2 represents body weight, g(z) is
described by Figure 2.3b, and N(z) = fl a constant.

Imagine starting with a single cohort of individuals in the
weight range 2

to (20 + s) at time t (Figure 2.4a). The total number

0
of individuals in the population is equivalent to the area under the

O

rectangular number density curve. Assume that s is small enough so that
the number density of the cohort can be approximated by an average value
over the width of the cohort. Two cases for the fate of the cohort

will be described.

CASE 1 - High mortality, u, over the full weight range (Figure 2.4a).
As individuals increase in weight, up to 21, the larger individuals

in the cohort grow at a faster rate than the smaller ones. The result

22

NUMBERS, N(z,to)
"(20.) 4

/N(z,tl) N(ZItz)

 

P_—————————_———-

 

width > s I I ‘(”
H’ I "
l L - iii 1
a) z 2 +5 if ' '
O O 1 22 2m

P-AGE (Body Weight)

 

 

 

 

NUMBERS, N("l'tb‘o’
N(zr')
l
I
|
I N(ZItZ)
I N(z,t ) i i
i 1 l |
| / i I
| I m I I
l l I l l
I 1 1 l l 1 i 1
I I I
b) 20 zo+s 21 22 EB

P-AGE (Body Weight)

Figure 2.4. Trajectories of idealized cohort: a) Case 1 - Concave growth
function and high mortality at all p-ages. b) Case 2 - Con-
cave growth function and high mortality up to 21, low
mortality beyond 21.

23

is that the cohort spreads out, or "expands", to encompass a broader

weight range. At a later time, t , when the cohort has reached a mean

1

weight for individuals of z , the weight range of the cohort is broader

1

than 3, its initial width. See N(z,tl) in Figure 2.4a. The total area

under the number density curve however, has been reduced by a factor

e-“(tl-t )

of 0 due to mortality losses.

Once individuals reach the mean weight 2 the expansion process

1'

in the cohort changes to a compression process. From 21 onward, the

growth rate is continually decreasing (Figure 2.3b), and larger organ-
isms in the cohort are growing more slowly than smaller ones. The
result is that the cohort becomes compressed into a smaller weight range

as smaller individuals catch up to larger ones in size. At a time t2 > t1

and a mean weight of 2 near 2 the cohort looks like N(z,tz) in

ml

2
Figure 2.4a. Its width has been reduced from the width seen at time t1,
and the area under the number density is again decreased, this time by
a factor of e-U(t2-t1).
From this qualitative picture, it is clear that the magnitude, N,
of a number density at any point (z,t) is the result of the force of

mortality combined with "compression" or "expansion" of the distribution

in individual weight differences.

CASE 2 - High mortality up to weight Z low mortality for larger

1,
individuals (Figure 2.4b).

Suppose now that the dynamics of the cohort are identical to those
of Case 1, up to mean weight 21. After this point the growth rate

again decreases, as in Figure 2.3b, but this time assume individuals

are subject to little or no mortality in the weight range (21,22).

 

24

If this is the case, individuals growing beyond 2 will begin to

l
"stack up", i.e., the cohort will be compressed into an increasingly
narrow weight range, as before, but very few individuals will be lost

due to mortality. The result is that, at time t and mean weight 2

2 2’
the cohort will look like N(z,tz) in Figure 2.4b. The area under the
number density curve at t2 has decreased only slightly from the area at
t1, due to the small mortality loss during the interval [t1,t2], but
the cohort is now squeezed into a much narrower weight range that the

range seen at t Thus, the number density at time t2 has a greater

1.
magnitude than it did at t1.

These cohort growth dynamics provide a mechanistic basis for
describing the phenomenon of "stunting" in populations, which will be
discussed again in Section 5.2.

The quantitative tradeoff between mortality and growth rates in
changing the magnitude of the number density can be determined by
examining the solution (2.4-4b). Suppose in some small time increment

At, an organism of p-age z grows to p-age (z + 62).

Then, according to (2.4-4b),

2+62 '(2)+ (2) ~
N(z+62,t+At) = N(z,t)exp f - 3-§7§%-—- dz .(2.5-2)
2

It is clear that the relative magnitudes of N(z + 52, t + At) and
N(z,t) depend on the sign of the integral exponent in (2.5-2) which in

('~+~
turn depends on the sign of the integrand, - l (z;(§;z)] , over the

range [2, z + 62]. Since 9(2) > 0 over this interval, we can focus on

the numerator of this fraction, - [912) + u(2)] .

 

 

25

In both Case 1 and Case 2, up until p—age 21, the growth rate was
increasing, so over each small p-age interval, [2, 2 +<52], the inequality
-[g'(2) + IK2)] < 0 was true. Thus, N(z f 62,t + At) < N(z,t). The
number density of the cohort decreased in magnitude over each small
increment (52,At).

When p-age/thme increments for p-ages greater than 21 are considered,
a difference appears between Cases 1 and 2. In both cases, the growth
rate is decreasing so that g'(2) < 0. Case 1, with a high mortality
rate, evidentally has Ig'(2)l < u(2) so that -[g'(2) + u(2)] < 0. Hence,
each small p-age/time increment results in a decrease in the number
density, according to (2.5-2). The sum of these decreases is seen in
the magnitude difference between N(z,tl) and N(z,tz) in Figure 2.4a.
< 2 < 2 ,

1 2

must have Ig'(2)| > u for 21‘< 2 < 22. In this case, -[g'(2) +p(2)] > 0,

and each p-age/time increment increases the nunber density. The sum of

On the other hand, Case 2, with very low mortality for 2

these increases is seen in the difference in height between the number
density N(z,tl) and the number density N(z,tz) in Figure 2.4b.

In summary, as a cohort grows through the p—age profile, the
magnitude of the nunber density, N, will rise or fall in accordance
‘with the sign of the quantity, -[g'(2) + u(2)].

From this analysis we begin to see the relationship between an
observed p-age distribution of a population and the growth-mortality
interaction which gives rise to the p-age profile. Chapter 3 will
examine this relationship for a steady-state population in which

individuals of every p-age are represented.

CHAPTER 3
ASYMPTOTIC PROPERTIES--THE STABLE P-AGE DISTRIBUTION

3.1 Asymptotic Properties of C-Age Models

What is the fate of the p-age distribution and the total population
numbers as t gets "large", that is, on the order of very many generations?
This question can be answered for the system (2.2-1), (2.2-7) restricted
to time-invariant fecundity, mortality and growth rates. To provide
the answers, essentially the same methods employed by A. J. Lotka in
his pioneering work on the dynamics of c—age distributions will be used.
Lotka's well-known results will be compared with the p—age model's
behavior as t + ”.

Let us briefly review Lotka's three central conclusions about the
limiting (in time) dynamics of populations whose members have time-
invariant, c-age-specific fecundity and mortality rates.

A. As t gets large, the population c-age structure approaches a

fixed c—age distribution. That is, if da is a small c-age increment, .

N(z,t)da
Nt(t)

only on individual mortality and fecundity rates [Sharpe and Lotka, 1911].

then the ratio ‘+ c(a)da, where c(a) is a function depending
The quantity c(a)da is the prOportion of individuals in the population
aged a to (a+da).

B. For large t, the total population birth rate, B(t), defined
by Equation (2.2-3), grows or decays exponentially in time, as does the
total number, Nt(t). Both B(t) and Nt(t) grow or decay at a constant

rate, r, called the ”intrinsic rate of natural increase" [Lotka, 1925].

26

27

C. If the stable c-age distribution of an exponentially-growing
population is perturbed, it will eventually return to a stable config-
uration and the population will again grow exponentially. Lotka showed
this using an ingenious graphical method [Lotka, 1922], but the stability
of the c-age distribution is more easily seen by considering a pertur-
bation at some time 2 as defining a new initial condition, N(a,%), for l
the c-age distribution. Then, as the population and its c-age structure
develop from time 2 onward, according to results A and B, a stable c-age
distribution and exponential growth will again emerge for t >> I.

Lotka's methods and the results A - C form the basis for modern
population theory in the age-time continuum, and his work has been
refined and expanded by many researchers. The derivation of these
results can be found in a number of books and papers, e.g., Keyfitz
[1968], Langhaar [1972], Lotka [1925], and will not be repeated here.

In this chapter, results A - C will be derived for populations
characterized by physiological age. The methods used are similar to
Lotka's and the usual terminology and notational conventions of con-
tinuous c-age theory will be used wherever possible.

In anticipation of the stability results, consider the relationship
Ibetween p-age and c-age which was discussed in Section 2.2. There it
was pointed out that there is a one-to-one correspondence in the
:restricted p-age model between an individual's p-age and its c-age.
Thus, as a population acquires a stable c-age distribution, one would
expect it to also show a stable p-age distribution.

The advantage of pursuing a proof of A - C for p-age which is

sinfilar to Lotka's proof for c-age is twofold. First, the method helps

28

to illustrate how the p-age model generalizes the classical c-age theory.
Secondly, the exact shape of the stable p-age distribution will emerge

as a function of growth, mortality, and fecundity rates.

3.2 A Renewal Equation for B(t), the Birth Rate

 

Define G(t) as the instantaneous birth rate due to survivors of
the initial population, NI(2). As t increases, G(t) + 0, since all
members of the initial population eventually die or grow past the p-age
breeding interval, [2b, 2b + a]. The birth rate now obeys the following

integral equation:

zb+a
B(t) = G(t) + f b(2)N(2,t)d2. (3.2-1)

Since G(t) accounts for births from.members of the initial population,
the density function N(2,t) in the integrand must refer to those
individuals in the population born after time t = O. N(2,t) in (3.2-l)
is therefore given by the solution (2.4-5b), and substitution into

(3.2-l) yields

zb+0' ~ ~ g(zo) ~ ~
B(t) = G(t) + g b(2) N(20,t-T(2,2O)) 3737—, 2(2,20)d2. (3.2-2)

b

Using (2.2—6), this becomes

I+a £(2,zo)
= “' - “‘, —-—-..— d“. .2-
B(t) G(t) + £1) b(2) B(t T(2 20)) gm 2 (3 3)

This renewal equation is the equivalent, for p-age, of the well-known
Lotka equation for c-age. In fact, (3.2—3) can be transformed to the

Lotka equation by a simple change of variables in the integral term

 

to permit integration over the time domain. Since d2 = dt, it follows that

g(z)

 

=_l

 

 

'1‘:

29

T(2 +a,2 )

B(t) = G(t) + {(zbb O

b(2(E)) B(t-E) £(z(E),2 )dE. (3.2-4)
,20) 0

Also, since b(2(t)) E O for t > T(zb + a, 20) andt:< T(zb,20), equation
(3.2-4) can be written as
t

B(t) = G(t) + £ B(t—t)b(z(2)) £(z(2),zo)d2. (3.2-5)
Equation (3.2-5) is the familiar Lotka renewal equation. Theorems for
the existence and uniqueness of solutions to this equation are given
in Keyfitz [1968] and Bellman and Cooke [1963]. The solutions to
(3.2-5) are presented by Lotka [1925] and are also found in Keyfitz
[1968] and Langhaar [1972].

The unique solution to (3.2-5) can be written as
B(t) = Qle + E le . (3.2-6)

The constants Qk' k = 1, 2,. . .,mq depend upon G(t) [Keyfit2, 1968].
Keyfitz also shows that there is a unique, real root, r1, for (3.2-6)
and that the complex roots, rk, k = 2,. . .ém: all have moduli which are
strictly less than Irll. Hence, B(t) 2 Qle 1 for large t--that is, the

birth rate becomes a simple exponential, growing at the rate r The

1.
root 1:1 is equal to r, the intrinsic rate of increase, and r will replace

r1 in subsequent equations.

Since (3.2-5) and (3.2-3) are equivalent, the solution (3.2-6)

*with.its asymptotic properties must also be the unique solution of (3.2-3).

That.is, the birth rate associated with the model (2.2—1) using time—

invariant parameters becomes exponential as t +-w.

rt
Let us assume that t is large enough so that B(t) = 916 ,

30

and G(t) = 0. Then by Equation (3.2-3),

Zb+a 242,20)
= ~ — ~' —--—~—— d~. . ..
B(t) f b(2) B(t 1(2 20)) 9(2) 2 (3 2 7)
That is,

 

zb+a 'rT(z’zO) ~ d2. (3.2-8)
1 = f e g(z)
2b

Equation (3.2-8) is the characteristic equation for the root, r,
and it can be solved numerically by a variety of methods to give r.
[Keyfit2, 1968]. It was claimed above that r is the intrinsic rate of
increase for the total population numbers. To see this, use the solution

(2.4-7b), which is appropriate for large t:

9(20)
N(2,t) = N(zo,t-T(2,20)) g(z) 2(z,20). (3.2-9)

 

Stated in terms of a birth rate, this becomes

2(2,20)

N(z,t) = B(t-T(z,zo)) (3.2-10)

g(z) '

If B(t) is exponential, then

r[t"T(Z:Z )] 2(2,2 ) 3_2-11
“‘2'“ “he 0 [fl ( )

Now, the total numbers are given by Nt(t):

m m r[t-I(2,z )] (2,2 ) ~ _
N (t) = f N(2,t)d2 = f Q e 0 l 0 dé]. (3.2 12)
t Z0 Z0 1 9(2)

Differentiation of (3.2-12) confirms that

dt ’ rNt° (3.2-13)

31 j

Hence, both the birth rate, B(t), and the total numbers, Nt(t) are

exponential in time for large t.

3.3 Relations in an Exponential Population
The exact shape of the stable p-age distribution can now be easily

derived. As before, let

N(z,t)dz

c(2,t)dz =
Nt(t)

(3.3-1) ‘-

be the proportion of individuals in the population whose p-age is between
2 and (2 + d2) at time t.
Using (3.2-10) and the definition of Nt(t), we have

B(t-T(z,zo)) 2(z,zo)

 

-——7;T—-d2
c(z,t)dz = ~ git; z ) . (3.3-2)
[m B(t-T(szo)) ' 0 ~
20 9(2)

Assuming t is large, we use (3.2-ll) to write this as

 

-r‘l'(2,2o)mz'z )
e 0 d
_ 9(2) 2
c(2)d2 - -rT(2}2 ) ~ . (3.3-3)
m e 0 2(2,2o)
f -——-=—-d2
20 9(2)

In (3.3-3) c(2,t) has been replaced by c(z) since it becomes independent
of time. The function c(z) in (3.3-3) gives the shape of the stable
,p-age distribution.

We will return to examine the shape function, c, in more detail,
but first let us generalize the definitions of some vital population

statistics which are widely used in c-age theory. Define the function

V by

32

2 +01 -rT(2.z )lb(2)2(2,z )
‘I‘(r) 9 lb e O 0
Z. L

 

g(2) d2. (3.3-4)
Then W(r) = l is the characteristic equation, (3.2-8). The quantity
W(0) is known as the net reproductive rate, R0, of the population
[Keyfit2, 1968]. That is, R0 9 W(O) is the average number of offspring
which will be produced by an arbitrary newborn individual during its

lifetime. To see this, note that b(2) is an instantaneous rate, and

b(2)dz
9(2)

of p-age 2(t) during the time interval [t,t+dt]. Since £(2,2 ) is the

= b(2(t))dt is the number of offspring produced by an individual

 

b(2)£(2,20)
proportion of individuals which survive to p-age 2, 9(2) d2

is the number of offspring expected from each newborn as it grows

through the p-age interval [2, 2 + d2]. Therefore,

L+a b(2)£(2,2 )
91(0) =21” [— °]d2

zb 1: 9(2)
is the average offspring production during an individual's lifespan.

In a stable population, r, the intrinsic rate of increase, can

be written
r = S - d (3.3-5)
where S and d are the per capita birth and death rates, respectively.

The birth rate 5 is also defined as

— _B(t)
b(t) — N (t)

t

(3.3-6)

 

That is,

33

- _ B(t)
b(t) — 2(2, 20)
zfm8(t-T(2,20 )) Wdz
ZO (

r
But B(t) = Qle t, so b(t) becomes a constant:

 

b = m-rT(2, z0 0)[?(2,:0 I] (3.3-7)
f:e d2
g() n

Once b has been determined from r, 2, and 9 via (3.3-7), the death rate,

d, can be obtained from (3.3-5).

 

If r = 0, i.e., the population birth and death rates balance, the kfi
dN
population is said to be stationary, and EE£-= 0.

 

Examination of (3.3-4) and the definition of R0 confirms the
following intuitive relationships between R0 and r: If R.0 < 1, then

r < 0, and the population numbers are decreasing. If RO > 1, then r > O
and Nt(t) is increasing. If R.0 = 1, then r = O and the population
is stationary.

We have shown that Lotka's asymptotic results A - C and many of
the vital statistics useful to demographers can be rigorously obtained in
the more general context of physiological age structures.

The asymptotic results A - C were recognized for discrete p—age
:models by Lefkovitch [1965]. In his discrete matrix model, described
Ibriefly in Section 2.3, the shape of the p-age distribution is given
by the eigenvector associated with the real eigenvalue of the projection
.matrix, M, which has the largest absolute value. The Equation (3.3-3)

has the advantage of expressing the stable p-age distribution directly

in terms of r and individual growth and survivorship functions.

34

3.4 The Shape of the Stable P-Age Distribution

Let us take a closer look at the shape of the stable p-age dis-
tribution, c(2), as specified in (3.3-3). How does it compare with
the stable c-age distribution derived by Lotka?

The stable c-age distribution is given [Keyfit2, 1968] by:

e’ra2(§,0) . (3.4—1)

5” e'ra2(§,0)da

c(a) =

As would be expected, this is a special case of (3.3-3) with g(z) = l
and z = a, chronological age.

For the sake of simplicity, assume that we are discussing the
p—age and c-age structures, cs(z) and cs(a), of stationary populations.
Then r = 0, and

£(2,20)

_.1(Z) I (3.4-2)
2(2,20)

c (2)
S 00

~

———w—dz
20 9(2)

while

2(a,O) ¥__.
fm2(§,0)d§
O

 

(3.4-3)

c (a)

Incidentally, the Equations (3.4-2) and (3.4-3) may also be
«obtained by setting %%-= O in partial differential Equations (2.2-l)
auni (2.3-2), respectively, and solving the resulting ordinary differé
ential equations for the stationary number density, N.

In (3.4-3), the shape function cs(a) must be a nonincreasing

:function of c-age, since it is a normalized, c-age—specific survivorship

curve .

 

II

 

 

 

 

35

The function cs(2) is more complex, however. Its most interesting
feature is that it need not be a nonincreasing function of p-age, 2.
This is more clearly seen by using the definition of survivorship,

(2.4-6), along with Equation (2.4-5), to rewrite c (2) as
s

_ z b'(£)+u<‘é)]~:l
expl: £0 5 9“?) dz ( )' n
. 3.4-4

_ Fg‘(£)+n(£) ~
zexPEé‘o [ g(i‘) )d’2:ld2

 

 

 

C(z)=
S

 

The nature of this shape function becomes apparent if it is

 

differentiated to get

 

dc (2)
s a __ g'(z)+n(2) _
-—-——-—dz 1: g(2) cs(2). (3.4 5)
dc (2)

Since 9(2) and cs(2) are both positive functions, the sign of d2
depends on the quantity -[g' (2)-111(2)]. If, over some p-age range
[21,22], the quantity -[g' (2)-111(2)] is greater than zero, then the shape
function cs(2). will actually increase in 2 over the interval [21,22] .
One would see a stationary p-age distribution which looked like
Figure 3.1.

The curve in Figure 3.1 illustrates the steady-state culmination

of the "compression" and "expansion" phenomena described for a cohort

density function in Section 2.5. The positive slope of the shape curve

do (2) .
{—32—- >01 seen between 21 and 22 is the steady-state version of Case
2 in Section 2.5. A rapidly decreasing growth rate and low mortality,
i.e., - [g'(z)+u(2)] > 0, force a "compression" in the number density,

N, and result in a ”stacking up" of individuals over the interval

36

cs(z)

   
 

-[g'(z)+u(2)] < 0

I
I
-[g'(z)+u(z)]>0:
i-[g'(z)+u(z)l < O

 

_—— ———

 

 

 

20 21 Z

P-AGE ' Z

N

Figure 3.1. Hypothetical p-age frequency distribution for stationary
pOpulation.

37

[21,22]. For the population.with a full spectrum of p-aged individuals
represented, the number density increases in 2 over the interval [21,22],
as it did for the cohort in Section 2.5.

A third factor, besides individual growth and mortality, shapes
the stable p-age distribution if the population is not stationary. In
this case, the total population numbers are increasing or decreasing
at a rate r, and the shape function c(2) is given by (3.3-3). 'To put
(3.3-3) in forms similar to (3.4-4) and (3.4-5), it is first necessary
to recall the definition of T(2,2o), Equation (2.4-3). The exponential
term in (3.3-3) can then be written as

-rT(z,2

e 0) = eprEEéz -—£—-d2].
O

92)

As in the derivation of(3.4-4),Equations (2.4-5) and (2.4-6) are used

£(2,20)
9(2)

terms of (3.3-3) are reassembled, they appear as

exp [£2 [ti-2' (EH um] d5]
O

to express the term of (3.3-3) in exponential form. When the

 

9(2) (
C(Z) . ’ ~ ‘— . A A _ 3.4-6)
{0 exp [- 4‘2 [r 9(2); (2)-J d2] d2
0 0
Upon differentiation, this becomes
____d°(z’ .. .. L—E—r "z” ‘2’] c(2). (3.4-7)
d2 9(2)

Taxis equation shows that c(2) increases or decreases depending upon the

sign of the quantity [r + g"(z) 4» 11(2)].

 

38

3.5 The Stable P-Age Distribution--Discussion

Equation (3.4-6) describing the stable p-age distribution lays
bare some prdblems encountered by an experimenter observing a p-age
distribution in a laboratory or wild population which is assumed to have
a stable p-age structure. Individual growth rate, 9, mortality,I1,
and the intrinsic rate of increase, r, interact in a complex manner to ’3
determine the shape of the observed p-age frequency distribution. I
Determination of any one of these vital statistics from an observed

p-age distribution requires knowledge of the other two.

 

This problem has been appreciated by many modelers. This section hi
will briefly discuss three different models from the literature, and
will concentrate on the efforts of each modeler to use experimental
p-age frequencies to determine growth and mortality rates.

Bell and Anderson [1967] seem to be the only other authors to
recognize the relationships between growth and.mortality expressed by
(3.4-6). They employ a partial differential equation like (2.3-l) to
describe cell populations whose members are classified by cell volume
and.c-age. The c-age variable is eliminated through simplifying assump-
‘tions, and a model of the form (2.2-l) results, with the p-age variable
(3f cell volume. Bell and Anderson assume a population is growing
exponentially, so that a stable cell volume distribution is achieved.
'Their equation (17) for the stable volume distribution is equivalent,
except for a normalizing factor, to Equation (3.4-6) with 11(2) = 0.

They discuss the problems of determining the intrinsic rate of
increase and the volume-specific rates of cell growth and division when

the only available data are volume frequency distributions of experimental

39

cultures. Due to the nature of cell life histories, Bell and Anderson's
discussion can not be easily generalized to other populations. For
example, although one can safely assume that few cells die before
division, the very process of cell division removes parent cells from
the population. The mortality rate can be assumed = 0 in cell popula-
tions, but a loss rate term due to cell division must be added. The .
shape of the stable volume distribution expressed by their equation F}
(24) is a function of cell division and growth rates and the intrinsic

rate of increase. Bell and Anderson conclude that it is impossible to

 

uniquely determine cell division and growth rates from volume frequency L}
data, but they suggest methods of estimating "reasonable" rates in
an interative fashion.

This provides an example of the possible value of (3.4-6).

Because Bell and Anderson could derive the exact dependence of the cell
volume distribution on growth and cell division rates, expressed in
equations analogous to (3.4-6), they were able to devise procedures

for estimating these vital statistics from experimental volume distribu-
tions.

The prOblem of estimating vital statistics from experimental p—age
distributions also arises with discrete p-age models. As we have
already mentioned, Lefkovitch [1965] originally could not quantitatively
interpret the elements of his projection matrix, M, in terms of growth
and mortality. More recently [1971], he has clarified the relationship

lmetween growth and mortality by writing the matrix M as:

M = w'lAw (3.5-1)

40

Here, A is the classical Leslie matrix containing c-age-specific
survivorships and fecundities as entries. W is a diagonal matrix with
each entry wii = mean p-age of individuals in the ith c-age class.
Lefkovitch shows that the matrix Equation (3.5-l) transforms a c-age
matrix, A, into a p-age matrix, M, via the (p-age):(c-age) growth data
contained in W. F1
In the context of Lefkovitch's model, the procedure for estimating
vital statistics can be summarized as follows: a sequence of observa-

tions from p—age distributions which are assuned to be stable are used

 

to estimate the entries in M for a specific population [Lefkovitch,
1965]. Given M, Equation (3.5-l) indicates that A and W are not uniquely
determined. Independent estimates of growth can be used to estimate W.
Then A, with its fecundity and mortality statistics, is computed from

A = WMW-l. This method is a discretized alternative to the use of
(3.4-6). The most obvious disadvantage of the discrete method is that
it is not suitable for analytical purposes. The quantitative tradeoffs
among mortality, growth, and fecundity in determining the stable p-age
distribution cannot be traced through the discretization process and

its attendant matrix manipulations.

It remains to be seen whether Lefkovitch's methods for parameter
estimation in the discrete model will satisfy Equation (3.4-6) as time
intervals and p-age classes become infinitesimal in length.

A discrete approach similar to that of Lefkovitch is taken by
Usher [1966] in his model of selection forests. A vector of size

classes for trees based on trunk girth is determined and a transition

matrix like that of Lefkovitch is used to update the size-distribution

41

vector. In Usher's model, the projection matrix takes a particularly
simple form. The update interval is chosen so that trees can advance

at most one size class during the interval. This assumption results in
a projection matrix, M, with non-zero entries only on the diagonal. In
addition, mortality is assumed to be zero during the update interval,

so the entries in M are a function solely of growth rate. These assump-
tions make it possible to compute size-specific growth rates directly
from two sequential size-frequency observations.

The three examples given above have a common theme. In each case,
the authors attempt to compute a vital population rate-—growth,fecundity,
or mortality rate, or the intrinsic rate of increase. The values of
other vital statistics are assumed to be known. The authors then use
empirical p-age distributions, which are assumed stable, to determine
the unknown vital rate, which may be a function of p-age.

The methods of Bell and Anderson and of Usher seem to produce
valuable results. However, each method is applicable only to that
small class of populations with life-history phenomena similar to that
of the population for which it was designed.

Although the method of Lefkovitch has a greater generality than
those of Bell and Anderson or Usher, it must be applied to data on a
case-by-case basis and cannot make predictions based on vital rates
which are functions, not data values.

We claim that the equations in Sections 3.3 and 3.4 for stable
or stationary p-age distributions are a basis for procedures to compute
'vital rates of any stable population classified by p-age which has

Observable p-age distributions.

 

42

R0 and r can be determined from (3.3-4) with a knowledge of b(2),
£(2,20), and 9(2). The crude birth rate, b, is given by (3.3-7),
which uses growth and survivorship functions, and r. Equations (3.4-4)
and (3.4-6) give the relationships among cs(2) , c(2) , g(z) , 11(2), and r.
Chapter 4 will give an example of the use of (3.4-4) in determining
mortality rates for fish populations.

As a caution, it should be noted that application of these formulae
will not be as simple as it may appear. The equations of Chapter 3 are
expressed in terms of continuously varying functions of physiological
age. There will be a considerable amount of work involved in developing
a sampling theory for estimation of these vital parameters. Even the
use of discrete-age and discrete-time data in the parameterization of the
continuous (c-age)x(time) models of Lotka and von Foerster still presents
statistical problems, as a reading of Keyfitz [1968] shows. Nonetheless,
it is clear that the equations for vital statistics in this chapter
can increase the utility of observing and modeling populations on a

demographic basis other than chronological age. Future work on this

model will proceed toward this goal.

 

CHAPTER 4

APPLICATIONS TO FISH POPULATION DYNAMICS

 

4.1 Introduction
Partly because of their importance as a food source for man, fish r1
populations have received a great deal of attention by ecologists.

I
Fishery biologists have long understood that a detailed knowledge of j

 

fish population dynamics is essential for the efficient management of bJ
exploited fisheries.

For these reasons, there is a large amount of data available on the
vital statistics of fish populations, and a good deal of work has been
done in the develOpment and practical use of age-structured models of
fisheries. A good sumary of such models can be found in Gulland [1962] .

A key factor in fish population models is individual growth in
body size. Fish are exploited by man on a size-specific basis, and the
total biomass yield of a fishery cannot be determined without knowledge
of population size distributions.

' However, the dynamics of size distributions in fish populations
have been difficult to examine, for several reasons. In the first place,
individual fish typically show very labile growth. Individuals which
survive the very high probability of death in the egg and larval stages
'will respond to changes in ambient temperature or available resources
by changing their growth rates [Weatherley, 1972] . This indeterminate
growth is generally felt to be the principal form of self-regulation

43

 

I

44

fbr many populations of marine and freshwater fish [Weatherly, 1972;
Backiel and LeCren, 1967].

The ages of fish in wild populations are often difficult to
determine. A great deal of effort has been expended in trying to age
individuals from various species, and the methods which are currently
used, such as finding aging marks on scales or various bones, are quite
tedious [deBont, 1967]. These aging methods decrease in accuracy “1‘
for older fish due to generally decreased growth rates and individual

differences in growth for fish of the same age [deBont, 1967;

 

Weatherley, 1972]. EJ

A further compligation in the attempt to study dynamics of fish
populations is that individual growth, fecundity, and mortality are
largely size-dependent rather than age-dependent. It has already been
‘mentioned that losses due to human harvesting are size-specific, and
there is evidence that natural predators also select fish on the basis
of size [Earl Werner, pers. comm.]. Several authors discuss the
correspondence of fecundity [Weatherley, 1972; Bagenal, 1967] and
growth rate [Weatherley, 1972; Parker and Larkin, 1959] with Size as
opposed to age.

For these reasons it seems clear that dynamic population models
based on size distributions might be easier to work with than those
based on population age structure. However, most of the theoretical
models which have been developed to estimate yields, recruitment, and
survivorship in commercial fisheries are age-structured. This is due
in part, no doubt, to the marked seasonal reproduction of most

temperate zone fish. Females in most populations spawn during only a

 

II

45

few weeks in the spring. Thus, a fish population is most often visualized
as comprised of a set of "year classes", which form the peaks of

traveling waves moving through an age distribution. Figure 4.1 shows

a snapshot in time of an idealized age distribution arising from seasonal

reproduction. The peaks in the waves are separated by approximately

one year . kl
The most widely-used models [Beverton and Holt, 1957; Ricker,
1958] employ this concept of year classes with simplified assumptions

concerning mortality and recruitment to describe population age structure.

 

The age distributions are then combined with size-at-age data or fitted u]
growth curves to give a picture of the population size distribution.

The remainder of this chapter will try to show how the p-age
model developed in Chapters 2 and 3 can be a useful alternative to
models based on year classes in studying fish populations. Rather than
discuss a single application of the p-age model in detail, the chapter
aim is to illustrate how a variety of new perspectives on population
dynamics emerge when individuals are classified by p-age instead of
c-age.

Section 4.2 shows the relationships of the p-age model to comonly-
used methods of computing production in fish populations. Sections
4.3 and 4.4 examine fish growth and recruitment, respectively, in the
context of the p-age model of Chapters 2 and 3. The two sections discuss
growth and recruitment models used by Beverton and Holt [1957] in their
comprehensive commercial fishery model, and Section 4.5 examines the
p-age distributions associated with their model populations. In Section

4.6 we turn to a different tOpic. This section examines the problem

46

[A one year

.JL.

NUMBERS

 

 

 

 

AGE

Figure 4.1. Idealized age distribution at fixed point in time, showing
year classes arising from seasonal reproduction.

 

47

of estimating fish mortality rates from sampled size or age frequency
distributions.

In this chapter, the p-age variable, 2, will represent body size,
either length or weight, of individual fish. The state variable of the
model is the function N(2,t), the number density of fish of size 2 at

time t. The dynamics of a fish population are then described by the

 

system (2.2-l), (2.2-7). From a modeling standpoint, one advantage
of using this system to describe a fish population is immediately

apparent. All three vital rates of fish population dynamics--growth,

 

fecundity, and mortality--are explicitly contained in the dynamical bJ

 

equations.

4.2 Biomass and Production

Computation of biomass production for a population is of primary
interest to fishery biologists as well as other resource managers.
Production over a time interval [tl'tZ] is defined as the total biomass
elaborated by the fish population, regardless of the ultimate fate of
the biomass [Chapman, 1967]. Thus, production during [t1,t2] is
equivalent to the total net change in pOpulation biomass, plus the sum
of the weights of all fish that die during [t1,t2] with weight being
assessed at the moment of death [Chapman, 1967].

From the model (2.2-l), (2.2-7), a differential equation can be
derived for population biomass. If the size variable, 2, is chosen to
be body weight, then

gm 2 N(2,t)d2 Q NB(t)é population biomass. (4.2-1)
0

48

When differentiated, this takes the form

dt 2 3t '

0

by Liebniz rule. Using (2.2-l), this can be written

dNB(t)
dt

 

= {m 2 _ g(z,t)N(Zpt)] d2 _ fa) Eu(2't)N(§'t)dE. Ht
0 32 20

If the first term is integrated by parts, the equation becomes

 

I
dN (t) a, m E”
'?fi?"= zog(zo,t)N(zo,t) + £0 g(z,t)N(2,t)d2 - £0 2u(2,t)N(2,t)d2,

assuming no individuals survive to size 2 = m, so that N(m,t) = 0.

This can be also written

dNB(t)

00 ~ ~ ~ oo~ ~ ~ -.
dt 2 203(t) + { g(zvtIN‘zvtIGZ - g zu(z.t)N(z,t)dz.

O O
(4.2-2)

Equation n.2-2) is an equation for net population biomass change
which shows clearly the three possible sources of biomass gain or loss
in any closed population. The first term in (4.2-2) is the rate of
addition of new biomass due to reproduction, the second term is the
addition of new biomass due to growth, and the third term is the total
loss rate due to mortality.

Assume now that we wish to compute production in a fish pOpulation
between time t and t . Also assume that no reproduction occurs between

1 2

t1 and t2 so that B(t) = O in that time interval. This is a common

49

assumption, since production is often computed over the summer months
beginning just after the spring spawning.
Then, production = [net biomass change] + [total biomass lost to

mortality]. That is,

t an (E) t
p = £ 2 i.— dt + £ 2 I” 211(2,2)N(2,2)d2 at.
1 at 1 2O H

Using (4.2-2) this equation can be rewritten as:

 

t .
p = £ 2 [F 9(2,E)N(2,'t)d2] at. (4.2-3) J
1 25O p

The use of (4.2-3) requires full knowledge of the number density
N(z,t) over the interval [t1,t2]. This is rarely practical for wild
populations, so (4.2-3) must be simplified. Suppose the production of
a year class, or cohort, is being computed as it grows during [tl,t21.

Then N(2,t) at any fixed time t, is non-zero over only a
restricted size range centered roughly on the mean size of the similarly-
aged individuals in the cohort. Thus, the inner integral in (4.2-3)
can be approximated by

g” g(2,t)N(2,t)d2 = g(2(t))Nc(t). (4.2-4)

0

In this approximation, 2(t) is the mean size of individuals in the
cohort at time t, and Nc(t) is defined as the number of indiviudals in

the cohort. Since g(2(t)) - £15-, production can be written as

dt
p a»;t2 93—95-1- u (eat (4 2-5)
t dt c ° °

1

50

Equation (4.2-5), with slightly different notation, is equivalent
to the basic production formula (9.19) of Beverton and Holt [1957].
Equation (4.2-5) is also the mathematical equivalent of the well-known
”Allen curve" graphical method for computing production [Chapman, 1967;
Weatherley, 1972], as is noted by Beverton and Holt [1957].

A slightly different approach is taken by Ricker [1958]. He
uses a simple equation for net biomass gain in a pOpulation, and then
integrates over [t1,t2], as above, to get a formula for production.

Ricker's formula can also be derived as follows: once again,
assume that biomass input from reproduction is negligible during

[tl,t2]. Then Equation (4.2-2) becomes

dNB = [0° A. ~ d~ foo ~ ~ ~ d~
at Z g(ZOt)N(Z,t) Z - Z zu(z't)N(z't) Z

0 0

If G(z,t) = ELELEL> is defined as the instantaneous growth rate per

unit of body weight, then

dN
—-§ = f" [G(2,t) - u(2,t)]§ N(2,t)d2. (4-2-6)

dt 0
run», suppose that growth per unit body weight and mortality are
constants across all sizes during [t1,t2] . Let G(z,t) = G and

(z,t) = 8. Then (4.2-6) becomes

dNB A
dz— = (G - fi)NB.

mexn this simplified equation for biomass, Ricker derives the production

formula P = GB where B is the mean biomass during [t1,t2]. Details can

be found in Ricker [1958] or Chapman [1967] .

 

Hui-'5 -.

51

The preceding derivations are valuable for several reasons. They
show that both of these commonly used production estimates are derived
via approximations to the same detailed expression for population
biomass (4.2-2). It appears that,this demonstration has not been
made before. The derivations also clearly show the mathematical
nature of the approximations inherent in the simplified production l‘
estimates. The formulae of Beverton and Holt and Allen are accurate ; a.
to the extent that Equation (4.2-4) is a good approximation. This

approximation, in turn, depends on the extent to which individuals in I

 

the same age cohort have similar body weights. L]
The accuracy of Ricker's formula, on the other hand, depends on
the constancy, in time and size, of growth and mortality rates for the
population size distribution being considered. Thus, Ricker's model
is most suitable for single year classes with the added restriction
that production be computed over short time intervals, as has been
noted by Chapman [1967]. Both production estimates assume that
reproduction does not contribute immediately to biomass change rates.
The apprOximations outlined above should be kept in mind by those workers
‘who attempt to modify fishery production equations for use in other
populations whose size distributions may show quite different dynamics
from those of fish.
In turn, the derivations provide a partial validation of the
Ibiomass Equation (4.2-2) and thus of the physiological age model

(2.2-1) 1 (2. 2-7) .

52

4.3 The Growth Function, g(z), for Fish Populations

Chapters 2 and 3 showed how both individual growth, g(z), and
its derivative, g'(2) = 33" influence the dynamics of size distribu-
tions. Growth in fish has been extensively studied, and there is a
wide variety of models, both theoretical and empirical, from.which
to choose.

Most of the models of fish growth have one thing in common. They
present growth in body weight as being a sigmoidal function of c-age.
That is, fish grow at an accelerating rate up to some weight, after
which growth rate decreases monotonically so that older fish appear
to approach an asymptotic upper weight limit. This "logistic" type
growth is shown by Figure 2.3. A single curve from the family shown
in Figure 2.3a is a sigmoidal growth curve, and Figure 2.3b shows a
typical growth rate curve which generates sigmoidal growth.

Curves of the type shown in Figure 2.3 are felt by D'Arcy
Thompson to describe the general course of growth in many organisms, and
in his treatise, "On Growth and Form" [1942], he gives numerous
empirical examples of this curve, e.g., growth in humans, bamboo,
whales, tortoises, and beanstalks.

There are a number of functions which describe sigmoidal growth
curves, and the particular function used depends on whether one wishes
to assume an upper size asymptote, and on the steepness and symmetry
properties of the sigmoid curve. Comparisons of various growth
functions are found in Weatherley [1972], Holt [1962], and in
Beverton and Holt [1957].

In passing, it should be noted that this chapter assumes that

53

9(2) does not depend on time. Growth of fish in temperate-zone waters
certainly has a seasonal component, with growth rates being reduced
in the winter months due to lower water temperatures [Weatherley,
1972], but we will follow the lead of other fish population modelers
and temporarily ignore this environmental component of growth.

The most striking effect of the sigmoidal growth pattern on the F1
dynamics of a size distribution is due to the decrease in growth rate
for larger sized fish. This effect is illustrated by Figure 4.2,

which is a size-frequency distribution shown at one point in time. As

 

a result of decreasing growth rates, the modes in the size-frequency

M
1"...)

distribution due to year classes begin to overlap for older and larger
fish [deBont, 1967; Holt, 1962].

The overlapping is further aggravated in many cases by differ-
ences in growth among individuals of the same age [Weatherley, 1972].
For many populations it becomes impossible to identify age-size
relationships for more than the first or second year-class by examining
size-frequency distributions.

This blurred modality in size-frequency curves is the major
obstacle encountered in using such curves to determine ages of larger-
sized fish, and the difficult aging methods mentioned in Section 4.1
must be attempted. In the language of Section 2.5, the overlapping
is caused by the "stacking up" or "compression", along the body
size dimension, of different-aged individuals due to decreasing.
growth rates.

A well known, sigmoidal growth formula, the von Bertalanffy

growth equation, will be used to illustrate how a specific function,

54

NUMBERS

Year classes

 

I \ / \\ I, ‘\
/ \ \
I { i \
,I [I \ 1 \ \

 

 

SIZE (length or weight)

Figure 4.2. Idealized size distribution of fish showing overlapped
year classes due to decreasing growth rate.

55

g(z), can be used in the population model. The von Bertalanffy growth
equation describes growth rate in body weight as a concave function
of weight [von Bertalanffy, 1957, 1968].

The equation for growth can be written

d2 1 3 2 3
mm = 5‘; = 3K [200/ 2 / - 2]. (4.3-1)
Integration of this curve gives:
-k(t-t0) 3
2(t) = 20 + (20° - 20) [l - e ] . (4.3-2)

In these equations, 20° is the asymptotic weight limit for the
particular species and ecosystem being considered, and K is a parameter
expressing the rate of approach to the upper size limit. The functions
are plotted in Figure 4.3. Notice that the growth curve is slightly
asymmetrical, with the maximum growth rate seen about 1/3 of 2”. The
family of characteristic curves in Figure 2.3 is composed of von
Bertalanffy growth curves. Since the equation has only been applied
in fish to post-larval growth, we let 20 be the size of young fry at
time to.

This growth function has been widely used and discussed in fish
population models. It is advocated by Beverton and Holt.who use it
as a basis for their theoretical age-structured.models [Beverton,

1954; Beverton and Holt, 1957; Holt, 1960]. Beverton and Holt [1957)
show methods of estimating 2m} and K from growth data, and they also
discuss the dependence of the parameters on population density.

Criticism of the von Bertalanffy equation as applied to fish growth

can be found in Parker and Larkin [1959].

56

 

 

z —— _ — — — ——————_—
CD
BODY
WEIGHT
.32 ———
CD
20 TIME
a) von Bertalanffy growth curve.
GROWTH
RATE,
g(z)

 

 

BODY WEIGHT
b) Growth rate for von Bertalanffy growth curve.

Figure 4.3.

57

4.4 Recruitment and Seasonal Periodicipy

The reproductive feedback in fish populations is poorly under-
stood. Since mortality in the egg and larval stages is very high
(>99.9+ % loss), it has been very difficult to relate the number of
surviving fry of a given size to the size of the breeding population
which produced them [Bagenal, 1967]. Thus, the basic models of
Ricker [1958] and Beverton and Holt [1957] simply assume that a

constant number, R, of recruits of a certain size, are added to the

" I11 k..- 3

population each year. For the moment, let 2R = the size of recruits. ,

 

We will also assume, along with the modelers mentioned above, E}
that recruitment, R, is constant from year to year and that yearly
deaths balance yearly births. The hypothetical fish population is in
a steady-state condition, since the net increase in population numbers
'is zero from year to year, but such a population is still dynamic in
that a train of wavelike modes due to year classes moves through the
age and size distributions.

This produces a population which is in a periodic type of
steady-state condition. Weatherley [1972, p. 155] defines the broadened
steady-state concept as follows:

To say that a population is in 'steady state' means that

such changes in growth, biomass, recruitment and mortality

as occur, do so with a regular pattern and magnitude so

that at a particular point in each phase of a population

'cycle', age and size structures will tend to be similar.
Total numbers, Nt(t), in this steady-state population would show a
"limit cycle" behavior, as described by May [1973]. POpulation numbers
oscillate up and down in a stable, periodic fashion between fixed maxima

and minima.

58

In the p-age model, Equation(3.4-2)gives the shape of the size
distribution for a stationary population, but it can be adapted to this
fish population which has yearly, pulsed reproduction.

Define a time-averaged size distribution,

i {HT N(z,t)dt. (4.4-1)

N(z,t) = T

Here, t is arbitrary and T is the period of the reproduction cycle,
about 1 year in this case.

Assume that the population has been in steady-state for several
generations so that N(z,t) can be replaced by the expression for the
solution (2.4-7b):

g(zR)
g(z)

N(z,t) =

 

HIF‘

ft+T N(z ,t-T(z,z )) 2(z,z )dE. (4.4-2)
t R R R

Notice that z replaces z in this usage of Equation‘2-4-7b), because

R 0

2R is the size at which recruits are assumed to enter the population.

Now assume that a total of R recruits enter the population at a

constant rate over a time interval TR. The interval TR would be on

the order of a few weeks to a few months for most temperate zone fish.

z
Then,

§}-, for E in [t,t+T] and in TR

N(zR.E~-T(2.2R))g(zR) = R .
, for E in [t,t+T] but not in TR

Here, the interval TR occurs sometime within the yearly period [t,t+T].
Equation (4.4-2) can now be written

l(z,z )

- R R
= — *- . 4.4-3
N(z) T g(z) ( )

 

59

In (4.4—3) we have replaced N(z,t) by N(z) since N becomes time inde-
pendent. Notice that N(z) is equivalent, except for a normalization
factor, to the shape function cs(z) defined in Equation (3.4-2) for
stationary pOpulations.

This time-averaging procedure smooths out any periodic "square
wave" reproductive input. Beverton and Holt [1957] use an instantaneous
yearly pulse for recruitment in their models, and this is handled in E}
(4.4-2) simply by making N(zR,t-T(z,zR)) a Kronecker-delta impulse

function in the variable t. Ricker [1958], on the other hand, assumes

 

that recruitment occurs at a constant rate over the whole year in his J
model. In this case TR = T and the time average in (4.4-3) still 9
applies. The remainder of this chapter will discuss the time—averaged
size distribution, N(z), of steady-state fish populations with pulsed
reproduction.
We ask the reader to look again at the idealized size-frequency
distribution for fish in Figure 4.2 and imagine its fluctuations during
the course of a year. The train of year-class modes will move to the
right in the distribution, and sometime during the year, a new wave of
recruits will appear at the left side of the distribution. However,
because of the small growth rate in larger-sized fish the position of the
blurred year-class modes comprising the right-hand section of the curve
should change little, and the relative fluctuations of numbers of larger-
sized fish should be comparatively small during the year.
This observation is illustrated by Figure 4.4, a year-long series

of monthly length-frequency curves drawn by Fairbridge for the tiger

flathead, an important species of the east Australian trawl fishery.

60

 

~= m
jll\j\

'°' "v \ ”\i,
°" \\\

3m 19“

 

 

 

 

 

 

v

A“ I“.

1mm. 1 \L

 

 

III! m0

 

/
I].

 

v

I" I“.

{\R
i a“, -
m M F\\\_
J
f
\/

 

 

7*

III I94. I

'7

\xii,

 

"1]
m... j \

Nu~_
so it 5:\"=1?"""”

FISH-[EIIYH (CI)

 

 

 

 

 

 

 

 

Figure 4.4. Monthly series of length frequencies for the New South
Wales Tiger Flathead population in the east Australian
trawl fishery. Taken from Fairbridge [1951].

61

Since spawning in this species is protracted over a six-month period, the
year-class modes are not distinctly seen. However, one can see that the
changes in the length-frequency curve for sizes larger than about 45 cm.
are small compared with fluctuations in the curve for smaller fish.
Therefore, the number of different time points chosen for sampling
such a population should not be crucial in building an estimate of the
seasonal average, N(z). This observation may give some comfort to .1
experimenters who would like to sample trawl catches as infrequently
as possible.

To recapitulate: it has been shown that, for a steady—state J

 

population (annual births = annual deaths) with constant annual recruit-
ment, the periodically-varying size distribution can be replaced by a
time-averaged size distribution which is stationary in time. Furthermore,
for larger fish, the average distribution, N(z), gives a good approxi-
mation to the distribution N(z,t) at any instant t.

The average distribution N(z), or its normalized equivalent,
cs(z), can now be used to examine the relationships between the growth
and survivorship functions, g and 2, since these parameters are unaf-
fected by the time-averaging procedure.

The problems associated with time-varying parameters such as
reproduction in age-structured models merit a short digression. Recently,
several authors have begun to examine population models with time-varying
parameters and have obtained satisfying results. Lopez [1967] defines
ergodic prOperties of a population. The tendency for a closed population
to approach a stable c-age distribution when subjected to time-invariant,

age-specific vital rates is called strong ergodicity. The tendency for

62

a closed population to have a c-age distribution which is independent

of its shape in the distant past and is determined exclusively by the
history of fertility and mortality is referred to as weak ergodicity.
Lotka proved the theorem for strong ergodicity of c-age distributions,
and Chapters 2 and 3 of this thesis prove the strong ergodicity property
of p-age distributions.

Lapez [1967] proves a weak ergodic theorem for the continuous
c-age model with some restrictions. Special cases of weak ergodicity
dealing with periodically varying vital rates have recently received
deserved attention. As a result, a good deal of biological realism
can be added to analytically tractable population models. Oster and
Takahashi [1973], using the von Foerster model (2.3—2), discuss the
results of synchrony or asynchrony between an environmentally-cued
reproductive pulse and the individual maturation delay.

In an excellent paper on seasonal periodicity, Skellam [1967]
shows the weak ergodic property for a Leslie matrix model whose parameters
vary periodically. The normalized age-frequency vector is shown to
depend asymptotically only on the season of the yearly cycle at which
it is computed.

Cull and ngt [1974] also prove a weak ergodic theorem for the
Leslie mode1--name1y, that a projection matrix with periodic entries
results in a stable c-age distribution for the average value of the
asymptotic population vectors over one period.

We conjecture that weak ergodic properties shown for the discrete
age—time models will also be seen in continuous age-time models with

periodic vital rates. Further, it is felt that these same ergodic

63

properties will be seen in p—age distributions. The time-averaging
results of this chapter for the steady-state size distribution with
periodic reproduction are a special case of the weak ergodic property

for the continuous p-age model.

4.5 ‘Application to a Specific Model

A well-known theoretical fishery model can now be examined from
the p-age point of view. In their basic model for a fishery, Beverton
and Holt [1957] assume a steady-state population of the type described
in Section 4.4.

The model of Beverton and Holt was specifically derived for
computing yields of exploited populations, and its many results will
not be reviewed here. However, the fundamental components of the model
will be outlined.

In addition to the steady-state assumption, Beverton and Holt
assume that the number of fish surviving to age a, out of R fish

-(F+u)(a-aR)

recruited into the population at age a , is given by Re .

R
The total instantaneous mortality rate is the sum of a natural mortality
rate, U, and a fishing mortality rate, F, both constants. The expression
is assumed to be valid for fish up to an age ax, the maximum age of
individuals sampled for the fishery.

The weight of an individual fish aged a is given by the von
ZBertalanffy formula (4.3-2). Using these expressions for survivorship
and growth, Beverton and Holt derive equations for biomass yield, ‘.

;production, and other statistics of interest for a given pOpulation.

Let us discuss the size distribution, N(z), associated with the

 

64

growth and survivorship functions described above. Combining (4.4-3)

and (2. 4-5) yields

 

" ____R____ _ Z g'(EHME) ~ _
N(z) T9(ZR) exp]: £1; 9(2) :ldz. (4.5 1)

For von Bertalanffy growth, 9' (z) is obtained by differentiating (4.3-1)

to get . F1

 

 

 

2 1/3
g'(2) = K 2 [:03] - 3 . (4.5-2)
Thus, letting T = 1 year, N(z) becomes 1
1/3 :3}
x [:2 [if] 4:] + (my)
N(z) = (R) exp - fz 2 1/3 2/3 dz . (4.5-3)
ZR zR 3K [20° 2 ~23

This equation gives the size distribution for the fish populations
described by Beverton and Holt's basic model. Figure 4.5 shows N(z) for
two different unexploited populations of North Sea plaice, using parameter
values estimated by Beverton and Holt [1957]. For these populations,
Beverton and Holt are interested in fish between the ages of aR = 3.7

years and aA = 18 years. For X = 0.1, 2R1 and le are the corresponding

upper and lower size limits, while 2 and z)‘ are upper and lower size

R2
limits, respectively, for K = 0.2. The curve for K = 0.1 describes a

2

”canonical" pOpulation whose statistics are felt to accurately reflect
pre-WWII conditions in the pOpulation. Beverton and Holt use this set
of parameters as a baseline for parameter perturbation analyses. This
particular size distribution could be compared with the original weight

frequency data to validate the derived values of 11, K, and zoo.

65

NUMBERS x (R/65)

aR = 3.7 years

ax = 18 years
200 = 2867 grams
1.0 .1 u = 0.1 year-1
F = 0.0 year-1

.75.

.50.

 

. 25-1

 

 

 

 

1655 2672

1....
N
N
0‘
O
U)

SIZE (in grams)

Figure 4.5. Steady-state size distributions associated with a theoretical
fishery_mode1. Parameters from Beverton and Holt [1957].
(See text for explanation.)

 

66

As was pointed out in Section 3.4, a size-frequency curve can
actually increase in 2, if the condition [g'(z)+u(z)] < 0 holds over
some interval [21,22]. For the size distribution N(z) described by
(4.5-3), g'(z) decreases monotonically to (—K) as 2 increases to zm.
(See Figure 4.3b and Equation 4.5-2.)

Thus, if (-K + (U+ F)) < 0, the size distribution N(z)*'“>, as
z-+ z”. Since this can not be the case in an actual size distribution,
caution must be exercised in choosing values for K, H, and F for the
Beverton-Holt model.

Even if the maximum size of fish seen in the population is a
value 2A < zm, one may still get size distributions increasing in z.
This is illustrated by the graph in Figure 4.5 with K = 0.2.

The fish in this hypothetical pOpulation grow much faster than
those in the plaice population with K = 0.1. Notice that, with K = 0.2,
the number density of fish increases for fish larger than about 1300
grams.

The effects of altering K on the yield of a population are discussed
in Section 17.5 of Beverton and Holt's monograph. They conclude that a
population of faster-growing fish (K = 0.2) would produce a greater
biomass yield to fishing at low values of F than one with slower
growing fish.

This is corroborated by comparing the two curves in Figure 4.5.
The steady-state population with K = 0.2 shows a greater number of
larger-sized fish available in the exploitable size range, and one
'would expect greater yield from such a population if it were not fished

intensively enough to radically change its size structure. Figure 4.5

67

is also a striking illustration of the problems encountered in assuming,
along with Beverton and Holt, that fish growing at widely different
rates (K = 0.1 or K = 0.2) become liable to capture at a fixed age,
rather than a fixed size. Fishing gear must select fish largely on the
basis of size, and fish at the assumed age of recruitment, a , in the

R
two populations in Figure 4.5 differ in weight by a factor of 5.

4.6 Catch Curves and Mortality Estimation

In the preceding section, N(z) was used to predict the size
distribution associated with a given set of growth and mortality para-
meters for a steady-state population. This procedure can be reversed.
Empirical size distributions can be used to estimate growth and
mortality statistics for populations.

A catch curve is the size-frequency distribution of a sample of

 

the population at one point in time. The population is usually presumed
to be in steady-state for analysis of such distributions. The sample

is often taken from the catch of a commercial fishing vessel. Catch
curves are used primarily for assessing survivorship in fish populations,
whether mortality is due to natural causes or fishing. Chapter 2 of
Ricker's Handbook [1958] is a detailed account of the analysis of catch
curves for mortality estimation, and most of this section is based on
his work.

A typical catch curve is shown in Figure 4.6. Such curves are
Characteristically dome-shaped. The ascending left-hand limb reflects
increasing selectivity of fishing gear for larger fish and is not
representative of the actual size distribution. The descending right-

hand limb, however, is taken to be a true sample of the actual population

 

 

68

 

 

 

 

1 A J- A l J
425 525 615 723 m 925 '0” N 25 11 15

H . 4

Figure 4.6. Catch curve for the Pacific Halibut pOpulation, 1925-1926.
Abscissa--mean length of successive 5 centimeter length
groups, in millimeters. 0rdinate--logarithm of the number
of fish taken at each length interval. From Ricker [1958].

size distribution. In our notation, the right-hand limb of a catch

curve is an estimate of a segment of N(z), the steady-state size
distribution.

Fishery biologists, like the modelers mentioned in Section 3.5,
recognize that a size-frequency distribution of a steady-state population
shows the compounded effects of growth and mortality. To estimate
mortality in a population, growth must be "factored out" of the size-
frequency curve. For this reason, most recent authors use catch curves
showing chronological age, rather than body length or weight, along the
abscissa [Ricker, 1958]. Mortality is easily determined from catch

curves which are c-age-frequency distributions. A catch curve based on

c-age is equivalent to an empirical estimate of c-age—specific

69

survivorship, 2(a,0), except for a scale factor (see Equation (3.4-l)
with r = 0). The definition of survivorship, Equation (2.4-6), can
then be used to estimate mortality.
However, catch curves based on age are often difficult to obtain.
First, a large sample must be obtained and a size-frequency catch curve
drawn. Then, subsamples of the original sample are taken for age
determination. From the resulting size-at-age data, the size—frequency
curve is transformed to an age-frequency distribution [Seber, 1973:
Ricker, 1958].
There are situations when the direct use of a size-frequency
curve is more attractive. To quote Ricker [1958 , p. 72],
For example, when it is a question of assembling a repre-
sentative sample of the catch from a widely scattered fishery,
it may be necessary to sample so many fish that determina-
tion of the age of all of them becomes very tedious. . . .
In such a situation there would be two curves available:
(a) a curve of mean length against age, based on a relatively
limited body of data, and (b) a representative curve of the
logarithm of frequency against length, based on all the
samples available.
To correspond to Ricker's treatment, the logarithm of N(z) will

1x3 used as a description of catch curves. From Equation (4.5-l) we

get

R

2 g' (any (‘2')
T9(zR)

~ d; + 1n 0 (406-1)
R g(z)

lnN(z) = - f
2

Tina slope of the logarithmic size distribution takes the simple form

 

g_[_1m'uzn a _l:g'(z)+11(z)]. (4.6-2)
dz 9(2)

‘ "-='"j

 

w 115-.

70

An estimate of mortality is then given by

d[1nN(z)]

u(2) = -g(z) dz

-g'(z). (4.6-3)

Ricker is able to handle a size-frequency curve directly only if
he assumes growth, mortality, and the slope of the curve are all constant.

Letting growth and.mortality rates be 9 and u, respectively, and

d lnN . . - . ‘
—L_dz££)—]— =..1' , Ricker derives the estimate for mortality, ‘ .q
u = i'g. (4.6-4)

 

See Sections 28 and 2G of Ricker's Handbook. This is clearly a special

 

case of Equation (4.6-3) with g and u constant, and g'(z) = 0.

Ricker then states that Formula (4.6-4) is useful only for
segments of size frequency curves where growth is constant. To illustrate
this prOblem, he uses a model to generate size distributions (Ricker's
Figure 2.11) showing the "stacking up” phenomenon described in Chapters
2 and 3. “

Ricker states that ". . .when the mortality rate is small and
fairly steady, and rate of increase in size is decreasing at a moderate
rate. . .,” the size distribution increases in z[£bid,, p. 72]. This
is a qualitative description of our criterion that [g'(z)+u(z)] < 0.
for the distribution to be increasing in 2 (see Figure 3.1).

An example of "stacking up" in an actual fish population is
shown in Figure 4.7 which is a smoothed catch curve drawn by Hart [1932]
for whitefish in an Ontario lake. Since a fine-meshed net was used,

the curve is representative of the actual length distribution for fish

longer than about 6 inches. The slight mode in the distribution at

 

71

 

ll
10’
fl!

3

lENG'fll - MILLIHETEIS

3333.

i

6“
411

 

 

 

t
1

 

IUHIEIS
“a
w

 

 

 

 

 

'\

 

 

 

 

 

 

A
l V
\

 

 

 

 

 

 

 

 

 

 

 

 

§

lENG'flI - INCHES

1m

 

{b

§

1!”

 

 

Figure 4.7.

Length-frequency curve for whitefish in Shakespeare Island

Lake, Ontario, showing mode at 17-18 inches due to decreasing
growth rate and low mortality in that length range.

Hart [1932].

From

 

 

72

17-18 inches is the result of fish "piling up", in Hart's words, due to
low mortality and decreasing growth rates in this length range. This
length-frequency curve is also a nice example of the coalescence of
year classes in a size distribution due to decreased growth, as
described in Section 4.3 and Figure 4.2.

Ricker points out that Formula (4.6-4) clearly will not work for

a catch curve such as the one in Figure 4.7. For such a distribution,

dlnN(z)
dz

conclude that there is a negative mortality rate over that size interval.

is positive over a size interval. Using (4.6-4) one would

Equation (4.6-3) enables one to deal with changing rates of growth
and sinuous catch curves to estimate size-specific mortalities.

Ricker also notes that mortality rates estimated using (4.6-4)

". . .always tend to be too small, if absolute rate of increase in length
is decreasing with age." [Ibid,, p. 73] Equation (4.6-3) explains why
this is the case. If growth rate is decreasing with age/size, then

g'(z) is negative, and the term [-g'(z)] is a positive quantity not
accounted for in the approximation (4.6-4). Thus, mortality rates are
underestimated by an amount roughly equal to [-g'(z)], the absolute rate
of decrease in the individual growth rate. This observation can refine
mortality estimates based on size-frequency curves.

The various forms of N(z) given in this chapter provide a general
mathematical description of size distributions in terms of growth and
mortality. They should help in increasing the utility of size distri-
butions, which are more easily obtained than age distributions, in

esthmating vital statistics of fish populations.

 

'F _ .'

A,

 

CHAPTER 5

THE REGULATION OF POPULATION GROWTH BY INDIVIDUAL GROWTH

5.1 Density-Dependence in Fish Populations

 

In earlier chapters, pOpulations were modeled for which individual
growth was time-independent, fixed presumably by the physiology of each
organism as expressed by its physiological age. Many organisms, however,
have a highly plastic growth response to changes in their environment.

In particular, individual growth rates for some species are strongly
affected by intraspecific competition fOr resources such as food or

space. This phenomenon is perhaps most striking in plant populations.
Harper [1971] points out that individual plant growth responds to
competition so flexibly that biomass production per unit area remains
virtually constant for many species over a wide range of number densities.

In this chapter, the p-age model (2.2-l), (2.2-7) will be used
to study the dynamics of pOpulations whose members respond to their
limited environment by altering their growth rates. Both analysis and
stmulation.will be used to uncover the behavior of the model when growth
rate, 9, is assumed to depend on some measure of pOpulation density.

As in Chapter 4, the discussion will use fish pOpulations for
a paradigm. A central theme of Weatherley's [1972] book is that
individual growth response to competitive pressure plays a key role in
regulating fish populations, and we will draw heavily on his observations.

A number of empirically based relationships between individual

73

 

 

74

growth and population density in fish have been suggested, all of which
state a negative correlation between growth rate and total numbers or

. . ' gg_ 1
biomass. For example, Backiel and LeCren [1967] propose that dD d 6-
describes the effect of density, D, on growth rate, 9. This implies
that variations in density at a low population level result in much
greater changes in growth than do variations at a high density level.
In contrast, from experiments in warm-water fish ponds, Hepher [1967] h]

found a log-log relationship between 9 and D. Hepher's model states

that a relative decrease in individual growth rate is a constant

 

multiple of the relative increase in density. j
Another more detailed connection between growth and density is ;

proposed by Beverton and Holt [1957], who also give an excellent general

discussion with numerous examples of density-dependent growth. Using

the von Bertalanffy growth equation (see Section 4.3), they conclude,

on both theoretical and empirical grounds, that the asymptotic weight

limdt, zm, decreases with increasing D, while the parameter K remains

unaffected by density.
All of these models implicitly evaluate density relative to the

carrying capacity of the environment, which will be denoted by C. Another

advantage of discussing density-dependence in the context of fish

populations is that managers of hatcheries and farm ponds have, for some

time worked with a practical, quantitative definition of carrying

capacity. Hepher [1967], for example, defines the carrying capacity

as the maximum weight of fish which can be sustained by a pond. A

‘population would be considered to be at the carrying capacity if it

were composed of individuals having zero growth, and if it showed no

net biomass production over some time interval.

75

To construct a more realistic model of density-dependent growth,

a detailed picture of the dietary habits of a particular species would
be necessary. One might divide the observed size range, [zo,zm], of
fish into subintervals to account for the switching of preferred

_diets shown by individuals of many species as they increase in size.
Associated with each subinterval, j, would be a carrying capacity, Cj'
based on the availability of prey items for fish in that size range.
Growth of individuals in the subinterval, j, is then a function of the
total numbers, Dj' relative to Cj’ For example, the size range of
bluegills in the experimental ponds of Hall, et al. [1970] could be
divided into two intervals. Individuals less than about 35mm. in
length fed largely on zooplankton, while individuals greater than 35mm.
in length preyed predominantly on benthic insect larvae. Unfortunately,
we can not, as yet, deal analytically with growth models of this
complexity.

The relationships between density and the other vital rates, v(z)
and b(z), must also be considered. Mortality in adult fish is generally
felt to be density-independent for most species [Beverton and Holt,
1957: Weatherley, 1972]. However, some researchers see density-dependent
mortality during larval stages as a prime component of natural population
regulation [Backiel and LeCren, 1967]. The model used.bere deals with
postrlarval fish, so it.will be assumed that density does not influence
:mortality.

The influence of density upon reproduction is still an open

‘problem for fish. Most fish modelers do not attempt to follow the

dynamics of hatching and larval mortality. Instead, stock-recruitment

 

76

models are used to relate the current density of fish to future expected
recruitment of juveniles into the population at a given size or age
[Ricker, 1958; Beverton and Holt, 1957]. These models show recruitment
increasing with increasing population density, up to a certain point.
At higher densities, recruitment begins to decrease, due to increasingly
density-dependent egg and larval mortality.

In this chapter, it is assumed that recruitment of fish sized
20 does not depend on density. That is, b(z) will be the instantaneous
rate of production of fish sized 20, per breeding individual of size 2.
The size 20 will be the size ofithe smallest post-larval fish in the
population, since it is assumed that mortality in post-larval fish is
density-independent.

With this definition of 20, the rate b(z) is no longer a true
representation of fecundity, since fecundity is usually defined as the
number of eggs laid by an individual [Pielou, 1969]. Since b(z) stands

for the rate of production of live offspring sized 2 per individual

OI
of size 2, it is more accurately a fertility rate, and it will be
referred to as such for the remainder of the chapter. Keep in mind

that, with these definitions of z and b(z), the time lag between the

0
egg stage and the attainment of size 20 is being ignored.

The seasonal component of reproduction will also be temporarily
ignored in order to apply the model (2.2-l), (2.2-7) with time-
invariant mortality and fertility rates.

In stmmary, then, we are studying a hypothetical fish population

with the following characteristics:

i. The population is growing in a limited, but continuously

 

ll ' '

77.

renewable, environment with carrying capacity C.

ii. Individual mortality rates, u, and fertility rates, b, are
time invariant, i.e., density-independent, but they may be
functions of body size, 2.

iii. Individual growth, g, is negatively correlated with pOpulation

density.

5.2 Stability_Theory

Tb understand growth as a regulatory mechanism, we will outline
weatherley's [1972] description of the growth of a fish population in a
virgin environment, beginning with only a few colonizers. Figure 5.1,
reproduced from Weatherley [1972], illustrates this process.

Initially, with abundant food and space, growth rates will be
high. Fish mature early, and large numbers reproduce before being lost
to mortality. As numbers and biomass increase, individual growth rate
decreases in response to increasing competition. With decreasing growth
rates, fewer individuals survive to reach a size at which they can
reproduce, and the total population growth slows down. Eventually,
total numbers and biomass approach an upper limit. At this population
level, growth has become adjusted to permit a balance between replace-
ment and loss of individuals.

The fate of the population size structure is not quite so clear.
It would be expected that due to the integrative nature of the breeding
interval, any initial size distribution in this population would evolve
to a stable size distribution in the same fashion as the populations of
Chapter 3 with fixed growth rates. Weatherley [1972] claims that a

more complex size-structure will be seen as the pOpulation grows.

 

 

78

Individual Growth Rate Population Growth

-~“\\/ \ max
‘5

 

\
‘ '1
\. ,
\ «.1
\
DENS ITY \ GRW‘I‘H

(Numbers or RATE' 9 . r

Biomass) j

\ 1;"

\\
‘s
Ԥ
\ \ —--.
1 O

 

 

 

TIME —"

Figure 5.1. Growth of a fish population in a limited environment.
The population initially has only a few individuals.
From Weatherley [1972].

 

 

79

However, it is not immediately clear that the density-dependent growth
rate could not resonate with an unstable size structure in such a
manner as to produce persistent oscillations in population numbers
about the upper asymptote of Figure 5.1.

If the size distribution did stabilize after several generations,
the population would become stationary, with balanced birth and death
rates, fixed total numbers and biomass, and unchanging individual
life histories.

To model the population described above, a simple and realistic
relationship between growth and density is assumed--namely that, at any
size ’2‘, the growth rate g(§,D,C) is monotonically decreasing function of
density, 0. That is, growth rates are greatest when the individual has
no competition for space or food, and they decrease with increasing D,
approaching zero as D + C. All three of the growth-density relationships
proposed in Section 5.1 follow this general pattern.

To model the dynamics of a population using (2.2-1), (2.2-7)
with a growth function such as the one described above, it is necessary
to consider density as a state variable of the model. If, for example:
D is equated to total population biomass, then Equation (4.2-2) for
biomass must be solved along with (2.2-l). If D is measured by total
population numbers, Nt(t), then equations (2.2-5) and (2.2-7) can be

combined to give the differential equation:

-—-—— = f b(§)N(E,t)dE - f u(z)N(z.t)d2- (5-2-1)
dt 2b 20

A useful tool for analyzing the effect of density-dependent growth

on population dynamics is the net reproductive rate, R0, defined by

80

Eguation (3.3-4) with r = 0. Since the mortality and fertility para-
meters are assumed to be time-invariant, the net reproductive rate can
be expressed as a function of D, which varies with time.

Using the definition of survivorship, (2.4—6), in Equation (3.4-3)

the net reproductive rate can be written as

213”“ b(2) ‘ E 11(2) A ~ .,
R0(D,C) - £1) W exp [- fzo W dz] dz. (5.2-2)

Recall that R.() > 1 implies an increasing (in total numbers) pOpulation,

while RO < 1 implies a decreasing population.

An important feature of the net reproductive rate is that its
definition and implications for population growth over the long term
are completely independent of the status of the c-age or size distribur
tion. Equation (3.3-4) merely expresses the relationship between R0
and the characteristic equation of a stable population. Thus, Equation
(5.2-2) can be used to predict population growth or decay as a function
of growth, density, mortality and fertility, without regard to the form
Of the number density N(z,t).

Suppose for the moment that a population is stationary-~individual
growth and density have mutually adjusted so that RO(D,C) = 1. What
will be the result of a perturbation in D, either through the removal or
addition of individuals? Will the population return to an equilibrium
State? One would intuitively expect it to do so, since increased
density results in decreased growth, which in turn should depress the
birth rate, thus lowering the density. Likewise, a decrease in D

Should speed up growth--the resulting increase in reproduction should

force density to increase once again.

 

 

81

Upon differentiating (5.2-2) with respect to 0, however, it can

be seen that it is possible for R0 to actually increase with increasing

density. A closer examination of (5.2-2) and the definition of the
.fertility rate, b(z), reveals the problem. In addition to the negative
feedback 100p between R.o and density described in the previous paragraph,
there is a positive feedback loop Operating through the instantaneous

birth rate.

To see this, consider the definition of RO on page 32 , and set

£(z,zo) = l. The quantity

2 +0. ~
_ b b(2) ~ _
R6 - f g(E,D,C) (5.2 3)

must be the total number of offspring produced by each individual which
survives to grow all the way through the breeding interval. This

quantity is clearly a function of growth, g. Although the £32§_of
fertility, b(z), remains constant, a low growth rate means that individuals
remain in the breeding interval for a longer time, and the total
fertility, R3 , is increased. An increase in growth rate similarly
implies a decrease in total fertility, since individuals now reproduce

at a constant rate for a shorter period of time.

Apparently, reproduction, in the sense of total fertility, has
not really been uncoupled from density. This result also suggests that
the integral Equation (2.243) for total birth rate B(t) is, in most
cases, a poor model for reproduction when it is used for a population
*with time-varying growth rates but time-invariant fertility rates. Such

a.model implies that slower-growing individuals produce more offspring

 

82

than faster-growing ones. Species of fish are not likely to be seen
with this type of reproductive feedback.

This reproductive model provides a positive feedback on population
density. As density increases, growth decreases and total fertility
increases, thus increasing density still further. The opposite occurs
if density decreases. The pOpulation.will have a stable equilibrium
only if the negative feedback effects of altered survivorship up to the
breeding interval can more than compensate fer the positive feedback
effect of density on total fertility.

A.much more realistic model is obtained if we allow the fertility
rate, b(z), to vary with density, but assume that total fertility, R6,

is constant. This can be done by using a fertility function of the form
b(z.D.C) = 6(2) g(z.D.C). (5.2-4)

This model says that slower-growing individuals produce offspring
at a reduced rate compared to faster-growing ones. The net result is
that the total number of offspring expected from any individual which
survives through the breeding interval is constant, regardless of growth
rates. Now the total fertility is given by

+0

11* = fzb Saws,
0 2b

which does not depend on growth. This reproductive model, with constant
total fertility, will be called model g, The original model, with time-
invariant fertility rates, will be referred to as model A?

The net reproductive rate for model §_takes the form

83

+0 4 x
RO(D,C) = 12:13 8(2) exp [— £2 ll—(-3"c2———-c1’2r::[d'§.. (5.2-5)

0 g(zoDIC)
The derivative of (5.2-5) is

dR +c E A A

__g= 2b A.» 11(2) ag<z.D.C) A .

dD fzb b(2) é m'c) 30 dz

0 .
_ a (2) A ~ -
exp [2 W) dz d2. (5 . 2 6)

0

Since growth is assumed to be a decreasing function of density,

the quantity £2, is negative for all 2 and D. All other terms in

BD
dR
(5.2-6) are always positive, so the derivative ESQ-is strictly less

than zero for all D.

Thus, the function RO has the properties shown in Figure 5.2. The
value D* satisfying RO(D*,C) = l is a stable equilibrium point for
population density. If D is perturbed away from D*, the net reproductive
rate will automatically adjust via the growth rate to oppose the
perturbation. Since RO is a monotonic function of D, the equilibrium
density is unique.

The above analysis is independent of the pOpulation size structure,
so there is still no mathematical guarantee that some combinations of the
vital rate functions 5, u, and g can not produce persistent oscillations
in the size distribution. However, the negative feedback loop described
above will act to oppose any changes in D which occur as the size

structure fluctuates. It is difficult to imagine a set of vital rates

in (2.2-l), (2.2-7) which would allow a continuously changing size

 

84

    
 

Rapid Individual Growth

 

 

 

RO(D.C)
I,.St:.able E uilibrium
1 "" q
I Slow Individual Growth
0 L. i
o 13* 'c

DENSITY, D

Figure 5.2. Net reproductive rate versus density for reproductive
model B.

 

 

ll]

 

85

distribution while keeping total population numbers or biomass constant
at D*. This is a further indication that the size structure in this
feedback model will ultimately stabilize.

Two other features of this density-limited form of the model are
worth noting. Figure 5.2 and Equation (5.2-5) suggest that the model
population density may stabilize well below the carrying capacity, C.
The actual value of D* at which the pOpulation stabilizes is determined
by the vital parameters of the particular species being considered, as
expressed by the equation RO(D*,C) = 1. In contrast, the familiar

logistic equation,

5N— _H
dt mu K)'

predicts that a population will saturate its environment and will only
stabilize when density, N, is equal to carrying capacity, K.

A second feature of the density-dependent p-age model is its
high degree of structural stability. This means that the system is
stable with respect to perturbations of its parameters [Rosen, 1970].
The requirement for any initial population to eventually stabilize at a
'unique, non-zero density D*‘f_C is simply that the population can grow
successfully in an unlimited environment. In other words, a necessary
and.sufficient condition for the population density to be stable is that

the following inequality be satisfied:

zbmo ~ " z (’2‘) ,7 ~
R (0.0 = f b(z) exp - f L———o-A dz dz > 1 (5.2-7)
0 2b 20 g(z,0,C) --

A11 sets of vital functions 8, g, and p which satisfy (5.2-7)

 

86

as a strict equality comprise the region of parameter space for which

the model is structurally unstable. This is a relatively small region.
For most vital rate functions, a small perturbation will not result in

a change of the inequality sign in (5.2-7). Unlike a density perturba-
tion, however, a perturbation of a vital rate in a stable population

will result in a new equilibrium value, D*, toward which population density
will tend.

According to this theory, a stable population with very labile
individual growth can respond successfully to enormous increases in
total fertility and/or decreases in mortality rate. An example of this
is seen in Weatherley's [1972] discussion of H. S. Swingle's experi-
ments with bluegill sunfish populations in farm ponds. Weatherley
describes these pOpulations, which are subject to very low mortality,
as being comprised of individuals which have virtually stopped growing.
There are a large number of juveniles and very few breeding adults.
These "stunted" pOpulations remained stationary for several years.

The theoretical arguments of this section show how a labile
individual growth rate is a powerful regulator of population growth in a
model in which total individual fertility remains constant (model 2).
The remainder of this chapter will use simulations to examine the manner
in which a growth-regulated population approaches its equilibrium

density.

5.3 Simulation—~Goals and Methods
The stability analysis of the last section reveals the asymptotic

behavior of the system (2.2-l), (2.2-7) when growth is density-dependent.

 

87

There are a number of questions, however, which can be asked about the
model dynamics for smaller values of t.

The first question which comes to mind is whether the model
population has similar trajectories for density and growth to those
described at the beginning of Section 5.2 for fish populations. We
also would like to settle the question of the fate of the size structure.
Does it stabilize, and, if so, how rapidly? How doesjixsstabilization
depend on initial conditions and on vital parameters?

Another interesting problem has to do with the measure of density
to which growth responds in controlling the population. Will the pOpu-
lation respond differently if D is measured in units of biomass rather
than numbers? These are questions which can best be answered by simula-
tion based on the model.

Simulation experimentsiku:the system (2.2-l), (2.2-7) are difficult
to design. The system has several parameters, some of which are functions
of body size. As a first step in reducing the parameter space to a
manageable size, this functional dependence was removed. Growth,
fertility, and mortality were made constant with reSpect to z.

The growth-density relationship was assumed to be as simple as

possible. Growth was defined by
(D C) - (1— 9-) (5 3-1)
9 I "’ gmax C - -

With this growth function, individuals grow at the rate gmax when no
competitors are present, and growth decreases linearly to zero as the
carrying capacity is approached. This is also the relationship between

density and growth implied by Figure 5.1.

88

These simplifying assunptions result in a seven-dimensional
parameter space, not including the initial condition, NI(z). The

A

parameters are gmax' u, b, C, zb, a, and zo. Once the simulation
nbdelwas operating correctly, the parameters could be checked one at
a time for their influence on system dynamics.

To further restrict the parameter set, an effort was made to
select values in the same range as those generally used for fish popula-

tions. The variable 2 represented body weight, and 2 was set equal to

0
l in arbitrary weight units, while the upper limit for weight was taken
to be 21 units. The simulation time scale and gmax were mutually set
so that an individual growing in an unlimited environment would reach
the maximum size in about 4 years. The size of first breeding, zb,
was taken to be 1/2 to 3/4 of the maximum size reached, and breeding
was, for most runs, allowed to continue up to the maximum size.

The mortality rate, u, and fertility rate, 6, were the parameters
most frequently altered to produce different trajectories. Mortality
ranged between 50 percent and 90 percent loss per year, while total
fertilities spanned a range of 10 to 2,000 offspring per individual.
These values are all similar to the ones employed by Weatherley in his
models of fish population growth and maintenance [Weatherly, 1972,

Ch. 6].

Some runs of the simulation were made with the reproduction model
ALof Section 5.2, but this model proved to be stable only within a
narrow range of parameter values. This was due to the positive feed-

back effect of density on total fertility which was described in Section

5.2. All of the runs reported in Section 5.4 employ the reproduction

89

model 2 with its more realistic assumption of constant total fertility.
This model showed stable behavior for virtually any combination of
fertility and mortality rates satisfying the inequality (5.2-7).

There was one set of conditions, however, under which the simula-
tion could "blow up." As developed here, the system cannot deal with
the growth function (5.3-l) if the pOpulation density overshoots the
carrying capacity. If this happens, then D > C, and growth, according
to (5.3-l) is negative. To handle this problem, it was necessary to
artificially set both growth and reproduction equal to zero whenever
D > C. The simulation continued to run until mortality reduced the
density below the carrying capacity, at which time reproduction and
growth began again.

The numerical techniques used for solution of the system (2.2-l),

(2.2-7) are discussed in Appendix A.

5.4 Simulation Results

Figure 5.3 shows the results of a typical simulation run. Both
numbers and.biomass follow a generally sigmoidal path to steady-state
values. This model population initially had 10 individuals, all
weighing between one and two units.

Oscillations about the sigmoidal curves are due to the size
structure of the model population. The peaks in the oscillations are
separated by one "generation”, where a generation is defined as the
time required to reach the size, zb, of first reproduction. Thus, as
the initial cohort.moves through the breeding interval, it generates a

pulse of young, producing a peak in total numbers. A corresponding

 

9O

 

 

 

 

P5000
Fm
1000‘ ”4000
800‘ T3000
iomass
(Arbitrary
eight
600‘ u = 50% loss/year P2333)
R6 = 40 offspring
2b = 11 wt. units
Numbers
zb+a = 20 wt. units
400‘ '1000
C = 10,000 wt. units

 

 

2001 o

 

J l
U U

o 56 100 150 230 230

TIME (in bimonthly units)

Figure 5.3. Sample trajectories of model population growth.

 

91

peak in biomass occurs in few time units later, as individual size
increases while nunbers decrease due to mortality. Numbers and biomass
begin to decrease rapidly, as mortality depletes both the new cohort
and the remnant of the reproducing initial population, and they reach

a minimum just before the new cohort reaches size 2b. When the new
ochort grows to size zb, a new reproductive pulse is produced, a genera-
tion later than the previous pulse, and the cycle begins again;

Two other dynamics of the population can be observed. As time
increases, the integrative nature of the breeding interval broadens
successive cohorts, so that the size structure eventually stabilizes.
This action accounts for the damping of the oscillations in Figure 5.3.
At the same time, the total density increases and growth rate simul-
taneously decreases. By the time density, which is equivalent to total
biomass for the population in Figure 5.3, has stabilized at D*, the
individual growth rate is only half as great as in the small initial
population.

This general pattern of population growth was followed by nearly
all the simulation runs made. Two variations on this theme are sketched
in Figure 5.4. Figure 5.4a shows total population numbers for a run
‘with a low value for 3 relative to u and a long breeding interval.

The population size structure loses the traveling wave of the repro-

duced initial condition within two generations. Since the difference
Ibetween fertility and.mortality rates is small, the population takes

longer to reach its equilibriun density.

Figure 5.4b is a very different case. In this model population,

individuals have very high total fertility and high mortality rates

92

 

 

 

 

 

 

 

 

lflmBERS
6000 .
4000 .
50% loss/year
20 offspring
11 wt. units
2000 .
20 wt. units
10,000 individuals
0 : :
0 150 200
TIME
a)
NUMBERS

6000.

4000. u = 75% loss/year
R3 8 2,000 offspring
2b = 11 wt. units

2000 zb+a - 20 wt. units

C = 7,500 individuals
0 l ‘51 . n : =
0 5b 1o'o 150 200
TIME

b)

Figure 5.4. Sample trajectories for model populations.

 

93

compared with the populations of Figures 5.3 and 5.4a. Oscillations in
total numbers and biomass showed no sign of abatement when the simulation
was run for 14 generations. Note, however, that population density is
not out of control, in the sense that D does not approach either zero

or C.

As Figure 5.4a illustrates most clearly, the pattern of pOpulation
growth shown in the simulations is identical to the dynamics described
at the beginning of Section 5.2 and in Figure 5.1. The model also
corroborates other qualitative descriptions by Weatherley of fish
population dynamics.

For example, Weatherley [1972] states that, as the pOpulation

density increases towards its asymptotic value, the size structure of

the population should contain an increasingly larger proportion of smaller

individuals due to declining growth rates. This feature is clearly
seen in the model populations such as the one in Figure 5.4a where the
‘population comes to contain a full spectrum of individual sizes early
in its growth. As the population nunbers and biomass climb the sigmoid
curve in Figure 5.4a, the mean size of individuals in the population
decreases.

The theory of Section 5.2, along with supporting evidence from
simulation runs, also provides an answer to one of’Weatherley's
questions about carrying capacity. In describing the expected effects
of a sudden improvement in food or living space for a steady-state
[population, Weatherley [1972, p. 160] wonders whether the population,
after increasing to a new higher density, would regain the same c-age

structure it possessed at the old, lower density.

 

IL:

94

According to Section 5.2, the answer to this question is "yes",
provided that density does not affect mortality or total fertility, and

that growth depends on relative density, i.e.,
D
g(zIDIC) = g(zra) (5.4-1)

To see this, recall that the density of a stationary population
satisfies the equation RO(D*,C) = l, for the net reproductive rate.
When relative density is considered, this equation reads RO[(D/C)*] a 1.
For reproduction model 2, this solution is unique. By the stability
analysis of Section.5. 2, any increase in C for a stationary population
wound force D to adjust to maintain the ratio (D/C)*. The pOpulation
would again become stationary at the higher absolute density, and growth
rates would return to their previous steady-state values, g(z,(D/C)*).

According to Equation (3.4-4), the shape of the stationary
size distribution depends solely on the vital rates, so it too would
return to its original configuration, as would the c-age distribution.

In Equation (5.3-l), growth depends on relative density. Using
this growth function, the simulation runs showed the behavior described
above, regardless of whether reproduction.model A;or E was used. The
carrying capacity, C, could be set at any value, and, if the other
‘parameters remained unchanged, the population would always stabilize
‘with the same size structure and at the same relative density.

Weatherley observes these dynamics in fish. He writes:

It is difficult to avoid the impression that growth flexibility
is an adaptive mechanism.by means of which fish populations
tend always towards the same population age/size structure

and the same individual growth rate, no matter what the nutri-
tional status of their environment. [Weatherley, 1972, p. 172]

 

95

The theory of Section 5.2 and the simulations are evidence that
any perturbations in the stationary density and/or size structure of
these model populations will be damped out. The return of the popula-
tion to its original stationary state will be characterized by damped
oscillations in total density. If a perturbation of a vital parameter
occurs, the same stabilizing behavior will be seen, although the new
stationary state of the population may be quite different from the
old one.

These damped oscillations in density suggest a measure of the
"resilience" [May, 1973] of the population in response to perturbations.
Borrowing from the vocabulary of control systems theory, we will define
the settling time of a population as the number of generations required
for D to settle within 10 percent of its steady-state value, D*,
following a perturbation. Since the oscillations in D are due primarily
to a fluctuating size structure, the settling time is really a measure
of how long it takes the size structure to stabilize.

In the simulations, D was measured either by total numbers or
biomass. Using Equation (5.2-5) for net reproductive rate and the steady-
state versions of Equations (5.2-1) and (4.2-2) for numbers and biomass,
one can easily predict what the steady-state values of numbers, biomass,
birth rate, and growth rate must be for any given set of parameters,
zo, 2b, a, 11, 8, gm, and C. Each simulation run canputed settling
times by continuously comparing current system states with their pre-
dicted steady-state values.

Table 5.1 sunlnarizes three features of the model populations

with respect to changes in life history parameters. The table shows

 

96

how a perturbation of a single parameter in a stationary population alters
the ultimate relative density, (D/C)*, mean size of individuals, 2, and

settling time of the new stationary population.

TABLE 5.1

Model Population Features as a
Function of Life History Parameters

 

 

b a zb u

 

 

 

(D/C)* + + - -
E - - + +
Settling

Time + - + +

 

 

 

 

 

Plus signs indicate that the parameterwamdthe population feature
both increase or both decrease, while a minus sign means that the para-
meter and the feature have an inverse relationship.

Many of the entries in Table 5.1 can be predicted from the theory,
and all were corroborated by simulation runs. The relationships shown
in the table generally conform to our intuitions about the life-histories
of individuals in limited environments. It is to be expected, for example,

that an increase in the width of the breeding interval, a, would result

in an increased birth rate and thus an increased stationary density, (D/C)*.

From the point of view of growth, the entries for; may seem somewhat
redundant in that they each have an opposite sign with respect to cor-
responding entries for (D/C)*. Any increase in (D/C)* means a reduced
growth rate, and thus, a smaller mean individual size.

0f greater interest, perhaps, are the relationships between

 

97

various parameters and settling time. The table implies, for example,
that a stationary population with very high fertility and mortality
rates would be rarely seen in the real world. Any slight perturbation
in the size structure of such a population, due perhaps to its environ-
ment, would result in traveling waves in the size structure and oscile
lations in total density which would take a long time to damp out. This
is the situation seen in the model population of Figure 5.4b, which has
a settling time greater than 14 generations.

A representative sample of settling times for various model
populations is shown in Table 5.2. The runs are listed in the table
in pairs. The model populations of each pair are identical in their
parameter sets and initial conditions. They differ only in the measure
of density which was used in (5.3-l) to control growth rate. One
population of each pair has a growth response to biomass while growth
in the other population responds to total numbers.

Two characteristics of this table are of special interest. The
first is that, for the runs reported, the settling time for total
population biomass was less than that for total numbers. This was, in
fact, a feature of all the runs, and, as the table shows, settling
times for biomass were shorter regardless of whether individual growth
was controlled by numbers or biomass.

Secondly, a comparison of settling times within each pair reveals
that, with the exception of pair B, the settling times for both numbers
and biomass were shorter when D was measured by numbers as opposed to

biomass. This relationship was true in six out of the seven paired

runs which were made.

 

98

Settling Times for Representative Model Populations.

Mortality Rate of .116 Results in 50% Loss/Yr., While u = .231
Results in 75% Loss/Yr. Settling Times are Measured in Generations.

 

 

 

 

 

 

 

 

 

 

Model Population Measure Settling Time Settling Time
Parameters of D for Numbers for Biomass
5| b a 45, zb = 11, Numbers 3.8 '2.2 I 7"
zb+a = 20, u = .231 Biomass 4.2 3.0 J
§_| b = 10, 2b = 11, Numbers 1.3 1.2
zb+a = 20, u = .116 Biomass 1.8 1.1 ‘
_gl b -- 10, zb = 14 Numbers 7.3 3.4 a"
zb+a = 20, u = .116 Biomass 10.0 4.9
3| b = 5, .zb = 14 Numbers 7.0 4.0
zb+a = 20, u a .116 Biomass 10.4 4.2

 

 

 

 

Taken together, these results suggest some differences between
model populations in which individual growth, and thus population growth,
are controlled by biomass rather than numbers. Tbtal biomass responds
more sluggishly than numbers to fluctuations in the size structure which
are induced by reproductive pulses. A reproductive pulse causes an
immediate sharp change in numbers, but its effect on biomass is delayed,
and thus reduced by mortality. Consequently, biomass fluctuations damp
out more quickly than fluctuations in numbers as the size structure
stabilizes.

This relative insensitivity of biomass to fluctuating size distri-
butions also explains why it is a "slappier" controlling variable than
numbers, in the sense that populations controlled by biomass usually

have longer settling times.

99-

These conclusions suggest a broader generalization: populations
regulated by biomass, as opposed to numbers, will allow greater fluctua-
tions in numbers, especially of young individuals, before responding

with altered individual growth rates.

5.5 Discussion

The model populations of the last section follow the same trajec-
tories for growth, total numbers, and size structure as are seen in real
fiSh populations. They may also be compared with studies of other
populations having density-dependent vital parameters.

Previous simulation studies of age-structured populations have
considered mortality and/or fertility to be functions of population
'density. The p-age model discussed here allows for the third vital
parameter, growth rate, to be influenced by population size. In addition,
it compares density with a carrying capacity to regulate population
growth, a feature not found in other distributed parameter models.

A brief review of these models and comparisons with the trajec-
tories shown here reveal a common pattern of growth in total numbers,
regardless of which vital rates are made density-dependent. Leslie
[1948] and Pennycuick, et a1. [1968] show how total numbers in matrix
model populations grow sigmoidally, settling down to a steady-state
through damped oscillations.

Motivated by laboratory populations such as Nicholson's blowflies
[Nicholson, 1957] and Frank's Daphnia cultures [Frank, 1960], several
modelers have been especially interested in the oscillatory nature of
the nunbers of biomass trajectory. We will briefly outline their

efforts and discuss the origins of oscillations seen in p—age and c-age

 

100

simulations. The most important feature of density-dependent simulations,
however, is that oscillations have a bounded amplitude.

Several investigators have induced persistent oscillations in
age structured models by explicitly employing time lags in the negative
feedback loops between density and the regulated vital parameters
[Leslie, 1959: Frank, 1960; Pennycuick, et al., 1968].

Fluctuations can also be produced by assuming that density affects
mortality or fertility on a c-age-specific basis. If this assumption
is coupled with c-age-specific rates which are very sensitive to small
changes over a short range of densities, stable oscillations in total
density will result [Usher, 1972]. Williamson [1974] carries this
sensitivity to an extreme by using one of two different sets of vital
parameters depending on the relation of density to some threshold, and
he too obtains limit cycle behavior.

In these systems as well, time lags give rise to the oscillations.
These lags are not explicit, however. Because of their age structures,
the models have inherent time delays. For example, an increase in
mortality rates for juveniles in the pOpulation will affect the total
birth rate only after a maturation delay.

A third method of producing persistent oscillations is to use
extreme values for mortality and fertility in model populations. This
effect is seen in Figure 5.4b and is also found in a simulation from
Pennycuick, et a1. [1968] in which a very high fecundity rate was used.
There is probably a time lag operating here, as well, although its
exact nature is unclear. The density-regulating factors apparently

cannot respond quickly enough to large fluctuations in population

 

101

numbers to allow c-age or size structures to stabilize. More work on
this is needed.

As a first step, we intend in future work to define settling times
fbr pOpulations with time-invariant values of g, b, and p. In such
pOpulations, the numbers or biomass trajectory would "settle" down to
trace an exponential curve of increase or decrease, as predicted by J
the theory of Chapters 2 and 3. These model populations will be compared r1
with density-regulated populations. If there is a time delay hidden

in the density feedback mechanism, then settling times for the

 

regulated populations should be longer than those for unregulated

*fi

populations.

Fluctuations in size structure, density, and vital rates in real
populations are the result of environmental perturbations, as well as
time lags. Whatever their nature, these fluctuations are kept under
control in our model populations, as in many fish populations, by a
highly flexible individual growth response. Those fluctuations which
are not damped completely, as in Figure 5.4b, tend to be centered about
a density 0* which is considerably less than C. Hence, minor perturba-
tions of the population are not likely to result in an overshooting of
the carrying capacity. This is an attractive alternative to the steady-
state "saturation" of the environment implied by the density-carrying
capacity relationship of the logistic equation.

A good deal of territory remains to be explored with the simulation
model. For example, it would be useful to quantify the relationships
in Table 5.1. Lefkovitch [1971] began a study of this sort when he

ranked various forms of age-specific fecundity and mortality with respect

102

to settling times for c-age distributions in unregulated populations.

His simulations indicate the large number of runs which would be necessary
to quantify Table 5.1. Perhaps a formal sensitivity analysis using the
equation for net reproductive rate would be more fruitful. The results
would be well worth the effort. As May [1973] points out, parameters

such as mortality rates and carrying capacities are much more likely to

be perturbed in natural populations than are size distributions or F1
number densities.

The most interesting and realistic simulation dynamics, however,

 

are likely to be seen when growth is made a function of p-age as well L}
as population density. The logical starting point at this next level
of model complexity is Beverton and Holt's density-dependent form of
the von Bertalanffy growth equation for fish populations (cf. Section
5.1). This avenue of research would provide a natural culmination to
the work of Chapters 3, 4, and 5.

In the future, most of our work with the distributed parameter
system will continue to focus on the phenomenon of individual growth
and its influence on populations. This is where the strength and
novelty of the physiological age model lies. For many populations,
notably fish, growth is a life-history parameter of equal importance
to fertility and mortality in determining population trajectories. We
feel that this model will be a useful tool for exploring the dynamics

of these populations.

 

B IBLIOGRAPHY

BIBLIOGRAPHY

Ames, W. F. 1969. Numerical Methods for Partial Differential Equations.
New York: Barnes and Noble.

 

Auslander, D. M., Oster, G., and.Huffaker, C. B. 1974. "Dynamics of
Interacting Populations," J. Franklin Institute, 297, No. 5:
345-376.

Backiel, T. and LeCren, E. D. 1967. "Some Density Relationships for
Fish Population Parameters." In The Biolggical Basis of Freshwater
Fish Production. Edited by S. D. Gerking. New York: Wiley.

Bagenal, T. B. 1967. "A Short Review of Fish Fecundity.” In The
Biological Basis of Freshwater Fish Production. Edited by S. D.
Gerking. New York: Wiley.

Bailey, N. T. J. 1964. The Elements of Stochastic Processes with
Applications to the Natural Sciences. New York: Wiley.

 

Baranov, T. I. 1918. "On the Question of the Biological Basis of
Fisheries." Rep. Div. Fish Management and Scientific Study of
the Fishing Industry, I:81-128. Moscow.

Bell, G. I. and Anderson, E. C. 1967. "A.Mathematical Model with
Applications to Cell Vblume Distributions in Mammalian Suspension
Cultures," Biophysical Journal, 7:329-351.

 

Bellman, R. E. and Cooke, K. L. 1963. Differential-Difference
Equations. New York: Academic Press.

von Bertalanffy, L. 1957. "Quantitative Laws in Metabolism and Growth,"
Quarterly Rev. Biology, 32, No. 3:217-231.

 

. 1968. General System Theory. New York: George Braziller.

 

Beverton, R. J. H. 1954. "Notes on the Use of Theoretical Models in
the Study of the Dynamics of Fish Populations." Misc. Cont. #2,
U.S. Fishery Lab., Beaufort, North Carolina.

and Holt, S. J. 1957. On the Dynamics of Exploited Fish
Populations. London: Her Majesty's Stationery Office.

 

103

 

104

Billups, C., Wilkins, B. and Pike, R. W. 1971. "Application of
Population Balance Models in the Prediction of Shrimp Size
Distributions." Coastal Studies Bulletin #6, Special Sea
Grant Issue, Center for Wetland Resources, Louisiana State
University, Baton Rouge, Louisiana.

Boling, R. H., Jr. 1973. "Toward State-Space Models for Biological
Populations," Journal Theoretical Biology, 40:485-506.

Chapman, D. W. 1967. "Production in Fish Populations.” In The Bio-
lggical Basis of Freshwater Fish Production. Edited by _
S. D. Gerking. New York: Wiley. [1

Coulman, G. A., Reice, S. R. and Tummala, R. 1972. "A Systems
Analytic Approach to Single Species POpulation Dynamics,"
Science, 175:518-521.

Courant, R. and Hilbert, D. 1966. Methods of Mathematical Physics,
Vol. II. New York: Wiley.

 

Cull, P. and VOgt, A. 1974. "The Periodic Limit for the Leslie Model,"
Math. Biosciences, 21:39-54.

 

deBont, A. F. 1967. "Some Aspects of Age and Growth of Fish in
Temperate and Tropical Waters.” In The Biological Basis of
Freshwater Fish Production. Edited by S. D. Gerking. New
York: Wiley.

Frank, P. W. 1960. "Prediction of Population Growth Form in Daphnia
Pulex Cultures," American Naturalist, 44:357-372.

Gulland, J. A. 1962. "The Application of Mathematical Models to
Fish Populations." In The Exploitation of Natural Animal
ngulations. Edited by LeCren, E. D. and Holdgate, M. W.

Hall, D. J., Cooper, W. E. and Werner, E. E. 1970. "An Experimental
Approach to the Production Dynamics and Structure of Freshwater
Animal Conmunities," Lmnolngd Oceanogram, 15, No. 6:839-928.

Harper, J. L. 1971. "The Regulation of Numbers and Mass in Plant
Populations." In ngulation Biologygand Evolution. Edited by
R. C. Lewcntin. Syracuse, New York: Syracuse University Press.

Hart, J. L. 1932. "Statistics of the Whitefish (Corggonus clupeaformis)
Populaticn of Shakespeare Island Lake, Ontario." Universigy of
Toronto Studies, Biol. Ser. #36.

 

Hepher, B. 1967. "Some Biological Aspects of Warm-Water Fish Pond
Management." In The Biological Basis of Freshwater Fish Production.
Edited by S. D. Gerking. New York: Wiley.

105

Holt, S. J. 1962. "The Application of Comparative Population Studies
to Fisheries Biology--An Exploration." In The Exploitation of

Natural Animal Populations. Edited by E. D. LeCren and M. W.
Holdgate.

 

Hughes, R. D. 1962. "A Method of Estimating the Effects of Mortality
on Aphid Populations," J. Anim. Ecol., 31:389-396.

Keyfitz, N. 1968. Introduction to the Mathematics of Populations.
Reading, Mass.: Addison-Wesley.

Kharkar, A. N. 1973. "A Stochastic Approach to Population Balance
Models." Ph.D. dissertation, Department of Chemical Engineering,
Michigan State University, East Lansing, Michigan.

, Coulman, G. A. and Barr, R. O. 1974. ”Population Models in
Ecological Systems: A Systems Perspective.” (To be published).

Langhaar, H. L. 1972. "General Population Theory in the Age-Time
Continuum," J. Franklin Institute, 293:199-214.

Lee, K. Y., Barr, R. 0., Gage, S. and Kharkar, A. N. 1975. "Formulation
of a Mathematical Model for Insect Pest Ecosystems--The Cereal
Leaf Beetle Problem," J. Theor. Biol. (in press).

Lefkovitch, L. P. 1965. ”The Study of Population Growth in Organisms
Grouped by States," Biometrics, 21:1-18.

. 1971. ”Some Comments on the Invariants of Population Growth."
In Sampling and Modeling Biological Populations and Population
‘Qypamics. Edited by G. P. Patil, E. C. Pielou and W. E. Waters.
University Park, Penn.: Penn. State University Press.

Leslie, P. H. 1945. "The Use of Matrices in Certain Population
Mathematics," Biometrika, 33:183-212.

. 1948. "Some Further Notes on the Use of Matrices in Population
Mathematics," Biometrika, 35:213-245.

. 1959. "The Properties of a Certain Lag Type of Population
Growth and the Influence of an External Random Factor on a
Number of Such Populations," Physiol. Ecol., 32:151-159.

Levin, S. A. and Paine, R. T. 1974. "Disturbance, Patch Fommation
and Community Structure," Proc. Nat. Acad. Sci. U.S.A., Vol.
71, No. 7:2744-2747.

. 1975. ”The Role of Disturbance in Models of Community
Structure.” In Eco-System Analysis and Prediction. Edited
by S. A. Levin. SIAM Conference Proceedings, Alta, Utah.

 

106
Lewis, E. G._ 1942. "On the Generation and Growth of a POpulation,"
Sankhya, 6:93-96.

Iopez, A. 1967. "Asymptotic Properties of a Human Age Distribution
Under a Continuous Net Maternity Function," Eggpgpgppy, 4:68-687.

Lotka, A. J. 1922. "The Stability of the Normal Age Distribution,"
Proc. Nat. Acad. Sci.,VIII:339-345.

. 1925. Elements of Physical Biology, Baltimore: Williams
and Wilkins. Republished by Dover Publications, 1956. r3

May, R. M. 1973. Stability and Complexity in Model Ecosystems. .ej
Princeton, New Jersey: Princeton University Press.

Nicholson, A. J. 1957. "The Self-Adjustment of Populations to Change,"
Cold Spring Harbor Sypposia Quantitative Biology, 22:153-173.

 

Nboney, G. C. 1967. "Age Distributions in Dividing Populations," B3
Biophysical Journal, 7:69-75.

Oster, G. and Takahashi, Y. 1974. "Models for Age-Specific Interactions
in a Periodic Environment," Ecol. Mon., 44:483-501.

 

Parker, R. R. and Larkin, P. A. 1959. ”A Concept of Growth in Fishes,"
J. Fish. Res. Ed. Can., 16:721-745.

Pennycuick,: C. J., Compton, R. M., and Beckingham, Linda. 1968.
"A Computer Model for Simulating the Growth of a POpulation, or
of Two Interacting Populations," J. Theor. Biol., 18:316-329.

 

Pielou, E. C. 1969. An Introduction to Mathematical Ecology, New York:
Wiley.

Ricker, W. E. 1958. "Handbook of Computations for Biological Statistics
of Fish Populations," Fish. Res. Ed. Can. Bull. 119:300.

 

Rosen, R. 1970. Qynamical System Theory in Biology, Vol. I. New York:
Wiley-Interscience.

Sharpe, F. R. and Lotka, A. J. 1911. ”A Prcblem in Age Distribution,”
Philosophical Mgg,, Ser. 6, XXI:435-438.

Sinko, J. W. 1969. "A.New Mathematical Model for Describing the Age-
Size Structure of a Population of Simple Animals." Ph.D.
dissertation, Department of Electrical Engineering, University
of Rochester, Rochester, New York.

and Streifer, W. 1967. "A New Model for the Age-Size Structure
of a Population," Ecolpgy, 48:910-918.

107

. 1969. "Applying Models Incorporating the Age-Size Structure
of a Population to Daphnia," Ecolggy, 50:608-615.

. 1971. "A Model for Populations Reproducing by Fission,"
Ecolggy, 52:330-335.

Skellam, J. G. 1967. ”Seasonal Periodicity in Theoretical Population
Ecology.” In Fifth Berkeley Symp. on Math. Statistics and
Probability, V61. IV. Edited by L. M. Lecam and J. Neymann.

Streifer, W. 1974. "Realistic Models in Population Ecology," Advances
in Ecolgical Research, 8:199-266. 1

, and Istock, c. A. 1973. "A Critical Variable Formulation of '41
POpulation Dynamics,” Ecology, 54, No. 2:392-396.

Thompson, D'Arcy W. 1942. On Growth and FormL_V01. I. Cambridge,
England: Cambridge University Press. s

313;"

 

Trucco, E. 1965. ”Mathematical Models for Cellular Systems: The b}
von Foerster Equation, Parts I and II," Bull. of Math. Biophysics,
27:285-303 and 449-471.

Usher, M. B. 1966. ”A Matrix Approach to the Management of Renewable
Resources with Special Reference to Selection Forests,"
Journal Applied Ecol., 3:355-367.

. 1972. "Developments in the Leslie Matrix Model," Sypp.
Brit. Ecol. Soc.,12:29-60.

VOn Foerster, H. 1959. "Some Remarks on Changing Populations,"
In The Kinetics of Cellular Proliferation. Edited by F. Stohlman,
Jr. New York: Grune and Stratton.

VOn Rosenberg, D. U. 1969. Methods for the Numerical Solution of
Partial Differential Equations. New York: American Elsevier.

weatherley, A. J. 1972. Growth and Ecolpgy of Fish Populations.
London: Academic Press.

weiss, G. H. 1968. ”Equations for the Age Structure of Growing
Populations,” Bull.Math. Biophysics, 30:427-435.

Williamson, M. 1974. "The Analysis of Discrete Time Cycles.” In
Ecological Stability. Edited by M. 8. Usher and M. H. Williamson.
London: Chapman and Hall.

\4 u
I I .
h . .u... i. .. 1oufilflu

 

APPENDIX

APPENDIX

NUMERICAL SOLUTION TECHNIQUES

A variety of numerical techniques are available for simulation of
the dynamics of system (2.2-l), (2.2-7). For the simulations reported
in Chapter 5, an explicit finite differencing scheme was employed.

Because of its desirable stability and convergence properties, a

 

centered difference equation was chosen to approximate the continuous
system.

Let Az and At be discretization increments of p-age and time,
respectively. As in Figure A.1, a grid is constructed on the (p-age)
x (time) plane, with p_p-age values along any time line. Let the

subscript (i,j) refer to the grid point (zi,t ). Then the centered

1
difference scheme is derived by approximating the partial derivatives
of Equation (2.2-l) at the center of a grid rectangle by first-order,

centered differences. That is,

3

 

 

 

N. . - "
1+1/2:J+1/2 ; .1. N1+1, j+1 Ni+l:i + Ni'ifl Nij (11.1)
at 2 A" At

likewise,

 

3‘9") 1+1/2. j+1/2 ; l]: gi+1.j+1Ni+1.j+1 f 91 .1'+1Ni.j+1
32 2 A2

(A.2)

+ gi+1.j“i+1.j ' gijNii]

AZ

108

109

 

 

 

 

 

 

 

n 1)-

I I

. I
I I

I
l I

N. . . .
1+1,j 1+1,j+l
N 1. — - -
i+1,0
P-AGE +
. N
2; 1+1/2.j+1/2
N. . N. .
Ni'o up a. - JJJ l-J—J+l

I | I

I g I

% i I

N N .
0!) OIJ+1

Figure A.1. Rectangular grid for centered difference scheme.

110

A four-point average is also taken for the mortality term.

~
—

i+1/2,j+1/2 _ 4

“-

(HN) [u

. . N. . + . . . . + . . . .
i+l,j+l i+l,j+l ui+1,jNi+l,j ui,j+1Ni,j+1

+ u (A.3)

..N..]
l] 13

The terms in (2.2-1) are replaced by their approximation A.1-A.3,

and the resulting difference equation is solved for N. This

i+l,j+1°

provides an explicit, or recursive, method of computing Ni+1 j+1 from
I

values N. ., N.., and N.

1+1.) 1) i,j+l' The boundary condition prOVides the

values NO j' at any time tj, and the initial condition NI(z) is dis-
I

cretized to give initial values {Ni }, i = 1,. . .n.

,0

The centered difference method requires values of g at time (j+1)

to compute the discretized p-age distribution {Ni } , i = 1,- - 00n-
I

j+1
Since gj+1 depends on the total population numbers, or biomass (if z =

body weight), at t.+ , a predictor-corrector algorithm is used. At

3 l

the beginning of a time step, g is set equal to gj. Then a distribution

j+1
{Ni,j+l} , i = 1,. . .n, is predicted, and the value of gj+1 is reVised.

Using the new value of gj+1, the distribution {Ni } , i = 1,. . .n,
I

j+1
is corrected. This process only required one iteration to give successive
estimates of gj+1 which were identical to three decimal places.

The integral expressions for reproduction, total numbers, and total
biomass were approximated using the trapezoid rule. More elaborate
methods, such as Simpson's rule, were tried, but there was no significant
effect on the solution curves.

A finite differencing technique is said to be stable if truncation

and roundoff errors do not grow as the numerical solution trajectory is

 

111

traversed, and it is said to be convergent if the difference equation

 

converges to the differential equation as At, Az +'0 [Ames, 1969].
Ames states that a generally applicable condition for stability of a
finite difference scheme is that 9 ° At :_Az, where g is the slope of
the characteristic curves in the (p-age) x (time) plane. In this simula-
tion of population growth, this restriction means that an individual
does not grow through more than one increment in p-age, A2, in a single
time step of length At. If the coefficients of the hyperbolic equation
(2.2-1) are constants then the centered difference method is stable and
convergent for almost any ratio of A2 and At [von Rosenberg, 1969].

With a time-varying growth rate, 9, it was found that unstable
behavior resulted whenever the ratio At/AZ was considerably smaller
than l/g. Incorporating Ames’ stability condition, we obtained the best

results using a time increment At which could be adjusted at each

iteration, j, to satisfy the condition:

Az
—— <

It was also found that the p-age increment, Az, could independently
affect the numerical solution. If A2 is too large relative to the width
of the breeding interval, a, then spurious oscillations of total numbers
of biomass are seen. These oscillations disappear, for a given value of
a, as A2 is decreased.

As an independent check on the accuracy of the numerical solution,
the net reproductive rate, R0, was computed at the conclusion of every
simulation run. The net reproductive rate is a function of growth,

fecundity, and mortality rates, as specified by Equation (5.2-5), and

 

112

if a pOpulation is stationary, it must satisfy the equation R0 = 1.

In all simulations which were allowed to run long enough so that the
numerical solution became stationary, the equation R0 = l was satisfied
to within il%.

A centered difference technique was applied by Billups, et a1.
[1971] to a simulation of an equation of the form (2.3-l) describing a
shrimp population. These authors also discuss an implicit differencing
scheme which involves solution of a large set of simultaneous linear
equations.

Another investigation of numerical schemes for first order,
hyperbolic models was carried out by Sinko [1969]. He tried a variety
of methods, including a combined centered difference-forward difference
equation. He also examined a hybrid method in which the p-age derivative
was approximated by finite differences, while the resulting set of
ordinary differential equations were solved in the time domain using
Runge-Kutta methods. Sinko finally settled on the technique of integrat-
ing along the characteristic curves of the partial differential equation
[Ames, 1969]. In this method, the characteristic curves themselves are
used to construct the discretization mesh. The solution is found by
integrating between nodes along characteristics in the mesh. First-order
differences are used to approximate the differential equations describing
the characteristics.

Sinko decided against using differencing techniques, largely
because of their tendency to generate occasional negative values for
the discretized density function Ni j' This phenomenon is also seen

I

in our simulations, and it appears to be due to the manner in which the

 

113

numerical techniques distort the wave-front of a pulse-wave type of
initial condition as it moves through the p-age domain. As a square
wave moves along the p-age axis in the numerical solutions, waves of
small amplitude and high frequency are propagated ahead of the original
numerical pOpulation wave. When the population numbers are near zero,
the troughs of these waves have negative number densities of small
magnitudes. These negative number densities become insignificant as

A2 is decreased.

As long as they are recognized as numerical artifacts, the negative

 

number densities do not imply that we are modeling a population with
"negative" individuals. Unlike Sinko's computer runs, the simulations
reported in Chapter 5 are not intended to accurately reproduce the
trajectory of a particular population. Since we have been more interested
in the general behavior of the system with a density-dependent feedback

on the growth rate, the centered-differencing method was more than
adequate.

Sinko's comparison tests between differencing techniques and the
method of characteristics were conducted using a time-invariant growth
function. In future work, we intend to compare the performance of the
method of characteristics and the centered—difference method when both

use time-varying growth and mortality functions.