Population Dynamics

We are going to look at two models for population growth of a species. The first is a standard exponential growth/decay model that describes quite well the population of a species becoming extinct, or the short-term behavior of a population growing in an unchecked fashion. The second, more realistic model, describes the growth of a species subject to constraints of space, food supply and competitors/predators.

Exponential Growth/Decay

We assume that the species starts with an initial population (Subscript[P, 0]). The population after n time units is denoted  (Subscript[P, n]). Suppose that in each time interval, the population increases or decreases by a fixed proportion of its value at the beginning of the interval. Thus (Subscript[P, n]) = (Subscript[P, n - 1]) + r (Subscript[P, n - 1]), n ≥ 1. The constant r represents the difference between the birth rate and the death rate. The population increases if r is positive, decreases if r is negative, and remains fixed if r = 0. We can express this in Mathematica via

P[n_] := (1 + r) P[n - 1]

Let's compute two populations at five year intervals for different values of r:

r = 0.1; P[0] = 100; Table[P[n], {n, 0, 100, 5}]

r = -0.1; P[0] = 100; Table[P[n], {n, 0, 100, 5}]

In the first case, the population is growing rapidly; in the secondd, decaying rapidly. In fact, it is clear from the model that, for any n, the quotient (Subscript[P, n + 1])/(Subscript[P, n]) = (1+ r), and therefore it follows that (Subscript[P, n]) = (Subscript[P, 0])(((1 + r))^n), n ≥ 0. This accounts for the expression exponential growth/decay. The model predicts a population growth without bound (for growing populations), and is therefore not realistic. Our next model allows for a check on the population caused by limited space, limited food supply, competitors and predators.

Logistic Growth

The previous model assumes that the relative change in population is constant, that is
		(((Subscript[P, n + 1]) - (Subscript[P, n]))/(Subscript[P, n])) =  r.
Now let's build in a term that holds down the growth, namely
		(((Subscript[P, n + 1]) - (Subscript[P, n]))/(Subscript[P, n])) =  r - u (Subscript[P, n]).
We shall simplify matters by assuming that u = 1 +  r , so that our recursion relation becomes
		(Subscript[P, n + 1])= u (Subscript[P, n])(1 - (Subscript[P, n])),
where u is a positive constant. In this model, the population P is constrained to lie between 0 and 1, and should be interpreted as a percentage of a maximum possible population in the environment in question.

This time we will use the Mathematica commands Nest and NestList to automate the computations. Here are explanations of the commands

?Nest

?NestList

Thus Nest[f, 1, 10] computes f(f(...(f(1))...) with 10 evaluations of the function, and NestList[f, 1, 10] lists all the intervening iterates. Now let's compute a few examples.

Clear[q]

q[u_][x_] := u*x*(1 - x)

Nest[q[0.5], 0.5, 20]

NestList[q[1], 0.5, 20]

NestList[q[1.5], 0.5, 20]

NestList[q[3.4], 0.5, 20]

In the first computation, we have used the Nest command to compute the population density after 20 time intervals, assuming a logistic growth constant u = 0.5, and an initial population density of 50%. The population seems to be dying out. In the remaining examples, we kept the initial population density at 50%, and we used NestList to generate a list of densities for the first 20 time intervals. The only thing we varied was the logistic growth constant. In the second example, with a growth constant u = 1,once again the population is dying out. In the third example, with a growth constant of 1.5 the population seems to be stabilizing at 33.3...%. Finally, in the last example, with a constant of 3.4 the population seems to oscillate between densities of approximately 45% and 84%.

These examples illustrate the remarkable features of the logistic population dynamics model. This model has been studied for more than 150 years, its origins lying in an analysis by the Belgian mathematician Verhulst. Here are some of the facts associated with this model. We will corroborate some of them with Mathematica. In particular, we shall use the ListPlot and BarChart features in Mathematica to display some of the data.

(1) The logistic constant cannot be larger than 4.

In order for the model to work, the output at any point must be between 0 and 1. But the parabola u x (1-x), for 0 ≤ x ≤ 1, has its maximum height when x = 1/2, where its value is u/4. To keep that number between 0 and 1, we must restrict u ≤ 4. Here is what happens if u is bigger than 4:

Nest[q[4.5], 0.9, 10]

(2) If  0 ≤ u ≤ 1, the population density tends to zero for any initial configuration.

Nest[q[0.8], 0.99, 100]

lis1 = NestList[q[1], 0.75, 20]

Now let's use ListPlot and BarChart to display this data.

ListPlot[lis1, PlotJoined -> True];

<< Graphics`Graphics`

BarChart[lis1];

(3) If 1 < u ≤ 3, the population will stabilize at density 1 - 1/u for any initial density other than zero.

The third of the original four examples corroborates the assertion (with u = 1.5 and 1 - 1/u = 1/3). In the following examples, we set u = 2, 2.5 and 3 respectively, so that 1 - 1/u equals 0.5, 0.6 and 0.666..., respectively. The convergence in the last computation is rather slow (as one might expect from a boundary case -- or bifurcation point).

Nest[q[2],  0.25,  100]Nest[q[2],  0.75,  100]

lis2 = NestList[q[2.5], 0.5, 20];

ListPlot[lis2, PlotJoined -> True, PlotRange -> {0.5, 0.7}];

lis3 = NestList[q[3], 0.75, 100];

BarChart[lis3, Ticks -> {None, Automatic}];

(4) If 3 < u < 3.56994..., then there is a periodic cycle.

The theory is quite subtle. For a fuller explanation, the reader may consult Encounters with Chaos, by Denny Gulick (McGraw-Hill), section 1.5. In fact there is a sequence
	 (Subscript[u, 0]) = 3 < (Subscript[u, 1]) = 1 + (Sqrt[6]) < (Subscript[u, 2]) < (Subscript[u, 3]) < ... < 4,
such that: between (Subscript[u, 0]) and (Subscript[u, 1]) there is a cycle of period 2; between (Subscript[u, 1]) and (Subscript[u, 2]) there is cycle of period 4; and in general, between (Subscript[u, k]) and (Subscript[u, k + 1]) there is a cycle of period (2^(k + 1)). In fact one knows that, at least for small k, one has the approximation  (Subscript[u, k + 1]) ≈   1 + (Sqrt[3 + (Subscript[u, k])]). So

u1 = N[1 + Sqrt[6]]

u2approx = 1 + (Sqrt[3 + u1]) // N

This explains the oscillatory behavior we saw in the last of the original four examples (with (Subscript[u, 0]) < u = 3.4 < (Subscript[u, 1])). Here is the behavior for (Subscript[u, 1]) < u = 3.5 < (Subscript[u, 2]). BarChart is particularly effective here for spotting the cycle of order 4.

lis4 = NestList[q[3.5], 0.75, 100];

BarChart[lis4, Ticks -> {None, Automatic}];

(5) There is a value u < 4 beyond which -- chaos!

It is possible to prove that the sequence (Subscript[u, k]) tends to a limit (Subscript[u, ∞]). The value of (Subscript[u, ∞]) , sometimes called the Feigenbaum parameter, is aproximately 3.56994...  Let's see what happens if we use a value of u between the Feigenbaum parameter and 4.

lis5 = NestList[q[3.7], 0.75, 100];

ListPlot[lis5, PlotJoined -> True, PlotRange -> {0.2, 1}];

This is an example of what mathematicians call a chaotic phenomenon! It is not random -- the sequence was generated by a precise, fixed mathematical procedure, but the results manifest no discernible pattern. Chaotic phenomena are unpredictable, but with modern methods (including computer analysis), mathematicians have been able to identify certain patterns of behavior in chaotic phenomena. For example, the last figure suggests the possibility of unstable periodic cycles and other recurring phenomena. Indeed a great deal of information is known. The aforementioned book by Gulick is a fine reference, as well as the source of an excellent bibliography on the subject.