1. Introduction and Background
1.1
Continuous in Time vs. Discrete in Time Models
1.2 A Discrete in Time Model
2. The Mathematical Models
2.1 Model
1: Malthusian Growth
2.2 Model 2:
Growth with Carrying Capacity
2.3
Model 3: Growth with Carrying Capacity and Harvesting
2.4 Model 4: A Predator-Prey
Model
3. Fortran 90 Code implementing Models 1-4
4.
Computational Experiments and Model Refinement
4.1 Computational
Experiments
4.2 Model
Refinements
Mathematical models play an important role in describing both biological and physical phenomena and are frequently used to describe the state of a biological or physical system as it evolves in time. Generally, the biological sciences are more difficult to model than the physical sciences, partly because they are more complex and partly because there are no fundamental biological laws analogous to say, Newton' law in the physical sciences. We begin our study of biological systems by studying four simple population models: (1) Malthusian Growth, (2) Growth with Carrying Capacity, (3) Growth with Carrying Capacity and Harvesting, and (4) A Lotka-Volterra Predator-Prey Model
Mathematical models of time dependent processes can be split into two categories depending on how the time variable is to be treated. A continous in time mathematical model is based on a set of equations that are valid for any value of the time variable. A discrete in time mathematical model is designed to provide information about the state of the physical system only at a selected set of distinct times.
The solution of a continuous in time mathematical model provides information about the physical phenomenon over a continuum of time values. The solution of a discrete in time mathematical model provides information about the physical system at a finite number of time values. Continuous in time mathematical models have two advantages over discrete in time models: (1) they provide information at all times and (2) they more clearly show the qualitative effects that can be expected when a parameter or an input variable is changed. Discrete in time models have two advantages over continuous in time models: (1) they are less demanding with respect to skill level in algebra, trigonometry, calculus, differential equations, etc., and (2) they are better suited for implementation on a computer.
We will present both discrete in time and continuous in time models for each of our four population examples. However, all of our development, algorithms, and code will be presented for the discrete case.
To develop a discrete in time mathematical model, we mark off a time line as indicated in Figure 1.
<-- h -->
--x------------x----------x----------x----------> (t-axis)
t0 t1 t2 . . .
The points on the t-axis marked with the symbol " x " are called
grid points or mesh points. The grid we have defined is called
a uniform grid or a regular grid since each pair of adjacent
grid points are separated by a constant time interval. This time interval
between grid points, h, is called the time step, the grid spacing
or the mesh spacing. Following the notation introduced in Figure
1, we can represent an arbitrary grid point on the t-axis as follows.
ti = h*i for i = 0,1,2,...
Note that ti+1 - ti = h, which is just a restatement of the fact that successive grid points are separated by a time interval of h.
In 1798 the Englishman Thomas R. Malthus proposed a mathematical model for population growth. His model, although simple, has become a basis for the modeling of biological populations. Malthus observed that unchecked by environmental or social constraints, it appeared that human populations doubled every 25 years, regardless of the initial population size. Or, said another way, he observed that populations increase by a fixed proportion over a given period of time, and that, in the absence of constraints, this proportion was not affected by the size of the population. For example, this means that if a population of 100 individuals increased to a population of 135 over the course of, say, five years, then a population of 1000 would increase to 1350 over the same period of time.
The Malthusian population model has one variable and one parameter. A variable is the quantity that we are interested in observing and in this case is the population which we denote by P; P is the state variable for the model and characterizes the behavior of the population over time. A parameter, on the other hand, is usually a constant although it is possible for a parameter to change over time. The parameter in the Malthusian model is the population growth rate, which we denote by r. It is customary to define r by r = b - d, where b represents the population birth rate and d represents the population death rate. Can you see that the model predicts either population growth without bound or inevitable extinction, depending whether the growth rate is positive or negative.
To write down a mathematical equation for the Malthusian model, we assume that additions to the population are made at the descrete times t1, t2, ..., rather than continuously in time and that the time intervals are all the same: ti+1 - ti = h for i = 1, 2, ... . Let Pi denote the population size during the time period ti to ti+1, and let r = b - d denote the population growth rate. Further, we assume that our model has population P0 at time t0. The Malthusian model can then be represented by the following mathematical equation:
(*) Pi+1
= Pi + r*h* Pi, i = 0, 1, 2, 3, ... ,
P0
given
Eqn. (*) is called a difference equation; it defines the change in the population P between two different time periods: i and i+1. The difference equation is of first order because the indices, i+1 and i, differ by 1.
If Eqn. (*) is rearranged to read r = {[Pi+1 - Pi]/Pi}/h, we see that r is the per capita growth rate for the time interval ti to ti+1. To simplify things somewhat, we will assume that h = 1 and note that in this case r represents the population growth rate per unit time. So, if we are tabulating a population at monthly time intervals, r would represent growth rate per month; if the tabulation were at yearly time intervals, r would represent the growth rate per year. The Malthusiam model then becomes
(1) Pi+1
= Pi + r*Pi, i = 0, 1, 2, 3, ...
,
P0
given
Example: Suppose that a human population begins at 260 million (US population in 1994) and has a growth rate of 0.05 per year. Determine the population at the end of 1, 2, 3, and 4 years. Expressing all populations are expressed in millions, we then have from Eqn. (1):
So, the population after 4 years is 316.03 million. What is the trend of the population as i becomes large, say P100? Is this a realistic model of human population growth?
Problem 1: The computer program pop.f90 (Model 1) inputs an initial population, P0, a yearly growth rate r, and produces a population history over a 100 year period using Eqn. (1). Run the program using P0 = 10 and r = 0.05. Run the program a second time using P0 = 100 and r = -0.05. Graph Pi vs. i for i = 0, 1, 2, ... for both populations. Discuss your results.
The population size simulated using the Malthusian model may be in favorable agreement with real-world population data during the early years of a population's growth. What is often observed in biological communities after a relatively long period of time is a steady state where significant changes in population are not observed unless a major environmental change is encountered. There are several limiting mechanisms that could cause the population size to approach a steady state or carrying capacity. Can you name a few?
Carrying capacity can be described as the population level at which the birth and death rates of a species match, resulting in a stable population over time. The logistic growth model explicitly incorporates the idea of carrying capacity. Letting Pi represent the population during time period i, r represent the growth rate per unit time, and K represent the carrying capacity, the logistic model is
(2) Pi+1
= Pi + r*Pi*(1 - Pi / K), i = 0, 1, 2,
3, ... ,
P0
given
Eqn. (2) is the same as Eqn. (1) except for the term -r * Pi * Pi / K; this term tends to modify the growth of the population P. Also note that
Examples of growth with a carrying capacity are yeast cells, fruit flies, and the height of of sunflower plants.
Problem 2: The population of Sweden from 1800 to 1920 can be modeled by the equations
where all populations are expressed in millions. (a) Run the computer program pop.f90 (Model 2) to produce a list of populations from 1800 to 1920 in yearly increments. Graph your results. (b) Run your program again for the years 1800 to 1970. Compare your results with the actual populations: 1930 - 6.142; 1940 - 6.372; 1950 - 7.041; 1960 - 7.494; and 1970 - 8.040.
An interesting variation of Model 2 is obtained by adding harvesting to the model. Such harvesting in an animal population could occur as the result of hunting, fishing, or disease, and can be modeled by
(3) Pi+1
= Pi + r*Pi*(1- Pi / K) - H, i = 0,1,2,3,
...,
P0
given
where H is a harvesting constant.
Probem 3: Consider a population of sandhill cranes who are hunted because of their destructive foraging in grain fields. Let P0 = 100, r = 0.0987/year, and K = 194.6 and let all populations be expressed in thousands. Run the computer program pop.f90 (Model 3) using the harvesting constants H = 0, 2, 4 and 6. For each value of H describe the trend of the sandhill crane population. Does it appear to approach a stable value? If so, what is that value?
The models discussed thus far are single species models. One of the most interesting dynamics between populations has to do with interactions between species, such as the interaction of a predator species and its food species, its prey. The American biophysicist Alfred Lotka and the Italian mathematician Vito Volterra were the first (in 1925) to propose a model for interactions between predator and prey; their model is called the Lotka-Volterra model.
Consider a population consisting of two species: rabbits and foxes. The rabbit population tends to exhibit exponential growth unless held in check by a predator, such as a fox. Foxes on the other hand are territorial and have a relatively stable carrying capacity. Foxes also tend to breed better if they are well fed. In our model foxes eat rabbits at a rate proportional to the number of encounters between foxes and rabbits. This number of encounters is, in turn proportional to the number of foxes and the number of rabbits in a given region.
Denoting the rabbit population by R and the fox population by F, the Lotka-Volterra model describing the interaction between these two species is given by the following system of two first order difference equations:
(4) Ri+1
= Ri + gr*Ri*(1-Ri/K) -
dr*Ri*Fi ,
Fi+1
= Fi +gf*Fi*Ri - df*Fi , i
= 0, 1, 2, ...
R0
and F0 given
where
Problem 4: Run the computer program pop.f90 (Model 4) to determine the rabbit-fox populations defined by the following model constants: R0 = 1.05, F0 = 0.9, gr = 0.012, dr = 0.019, gf = 0.081, df = 0.0048, k = 1.0.. Populations are expressed in thousands. Describe the behavior of the two populations graphically. Graph Fi vs. Ri on a separate plot. Discuss this graph.
Here is a Fortran 90 source code implementing the all four models described in Sec. 2; it was developed using the NAG compiler. You are asked to specify the model that you want to solve as input.
To become familiar with the behavior of the various population models described in Sec. 2, perform the following experiments.
Experiment 1: An important thing to determine about any population model is the set of equilibrium points for the model. An equilibrium point is a population level from which the population doesn't change. Another name for an equilibrium point is steady state. Finding an equilibrium point is simply a matter of setting the population change to zero and then solving for the equilibrium population. For example, in Model 1, we write Eqn. (1) in the form
Pi+1 - Pi = r*Pi
and set Pi+1 - Pi = 0 to obtain
r*PE= 0
where PE denotes the equilibrium population. Solving this equation for PE, we see that PE must be zero. Of course, r could also be zero; in this case we would be modeling a population that never grows, and so any population size is an equilibrium size.
Equilibrium points are classified as stable or unstable. A stable equilibrium is one that is attractive; nearby population levels eventually settle down to a stable equilibrium. If a population level is near (but not at) an unstable equilibrium, the population levels will move quickly away.
Verify by computation that PE = 0 is an unstable equilibrium for the Malthusian model. Choose r = .01 and P0 to be a small positive number. What happens as the solution evolves? Does the situation change if r is negative, say r = -0.01?
What are the equilibrium points in Model 2? Verify by computation.
Experiment 2: Run Model 2 with input values r =0.07 (per year), K = 1000. and P0 = 150. Plot the resulting population for 200 years. Note the population's behavior. Is it what you expected? Try other initial populations, some less than K and some greater than K. What happens to the population for values of P0 greater than K?
Model Refinement 1: It is well documented that in simple cell growth and in insect population growth, the growth rate of a population is temperature dependent. On a daily or seasonal basis, temperature can be approximated as a periodic variation about some average temperature. For example, this behavior could be approximated by the mathematical expression
k*(T + A*sin(w*i))
where the w specifies the frequency (period) of the temperature variation, A specifies the amplititude of the temperature variation, T specifies the average temperature, and k would be a new growth constant. Note in the case of no variation in temperature (A = 0), k*T would play the same role as r did in models 1 and 2.
a) Write a program to compute and plot the expression k*(T + A*sin(w*i)) vs. i for i = 0, 1, 2, ...; select your own values for k, T, and A.
b) Replace r in the Malthusian model by the expression k*(T + A*sin(w*i)) and compute and plot the resulting population for i = 1, 2, 3, ... .. Note that the new model becomes
Pi+1 = Pi +k*(T + A*sin(w*i))*Pi, i
= 0, 1, 2, 3, ... ,
P0 given
Select your own model constants.
Model Refinement 2: The descrete in time mathematical model described for Malthusian growth (Model 1) can also be formulated as a continuous in time model:
P'(t) = r*P(t)
P(0) = P0
where prime (') denotes differentiation with respect to t. Those who have studied calculus will recognize that the exact solution to this model is given by
P(t) = P0*exp(r*t)
where exp denotes the exponential function. Write a program to evaluate P(t) above and compare it with your solution to Model 1 obtained from the program pop_f90. How well do they agree. Which do you think is "more accurate"?