Population Dynamics

Mathematical Modeling

Richard Hitt
Department of Mathematics and Statistics
University of South Alabama
Mobile, AL 36688

Fibonacci Growth

Perhaps one of the first mathematical models of a population growth was due to Leonardo Pisano (Fibonacci) in 1202.  In this setting, we assume we have a newborn pair of rabbits who will eventually mate.  We assume that rabbits are born in oppositely sexed pairs which become sexually mature at age 1 month and that the gestational term is 1 month.  

Let F[n] denote the number of pairs of rabbits at the end of each month.  Lets look at different ways we can compute the number of rabbits at the end of the nth month.

Simple but Inneficient Mathematica Code for the Fibonacci Recurrence

The most natural was to describe the Fibonacci relation to Mathematica is as follows.

F[1] = F[2] = 1;
F[n_] := F[n-1] + F[n-2]

Note that if we want to make a table of Fibonacci numbers, this formula for F[n] is horribly inefficient.  Explain why.  You can use the Timing function in Mathematica to measure how long a calculation takes (in CPU seconds, not wall clock time).  Of course, the value of the Timing output will vary from machine to machine.

Timing[Table[F[n],{n,1,15}]]

Improved Code

Mathematica can be convinced to cache its previously computed information about a sequence as follows.

Clear[F2];
F2[1] = F2[2] = 1;
F2[n_] := F2[n] = F2[n-1] + F2[n-2]
Timing[Table[F2[n], {n,1,100}]]

This approach is great if n is not too large.  However, it is easy to exceed Mathematica's recursion depth if n is large.  You can experiment with this approach, but you should clear Mathematica's memory of the cached F2 values if you run Timing tests so you are not using previously stored information.

Another Approach

The disadvantage in the above improvement is that in order to compute a Fibonacci number, all the intermediate numbers are stored as rules.  This will eventually require a large amount of storage for a term in the sequence if the index n is very large.  In fact, for sufficiently large n, the process will fail due to exceeded capacity.  Also, the lookup procedure for rules is slow.  These problems can be overcome with an iterative procedure.

F3[n_] :=
        Block[{f = 1, f1 = 0},
          Do[ {f, f1} = {f + f1, f}, {n - 1} ];
          f
        ]

Timing[Table[F3[n],{n,1,100}]][[1]]

Timing[F3[1000]]

This procedure really flies.  It does not have the recursion depth restriction of the previous procedure.

Graphs

To determine the growth rate of the Fibonacci sequence, we can plot the points.

ListPlot[Table[F3[n],{n,1,10}],AxesOrigin->{0,0}]

Is this exponential?  To see, we can plot the log of the data.  In fact, we can let the data go further out.

ListPlot[Log[Table[F3[n],{n,1,100}]]]

Since the log plot looks linear, the data are probably exponential in growth, at least over the range we are examining.  This observation is made more precise in the next subsection.

A Closed Form Method

In this subsection, we will derive the standard closed form formula for the Fibonacci sequence to show that the terms grow exponentially.

The method used to compute F3 can be rewritten in matirx form.

A = {{1,1},{1,0}};

Now the Fibonacci sequence can be described by a matrix equation with appropriate initial conditions.

X[n_] := A . X[n-1]; X[1] = {{1},{0}};

Since the eigenvectors of the matrix A are distinct, the basis can be changed to an eigenbasis (a basis consisting of distinct eigenvectors) in order to solve the system.

T = Transpose[Eigenvectors[A]];

Under this new basis, the coupled recurrence equations become uncoupled since the coefficient matrix A is transformed into a disgonal matrix of eigenvalues.

MatrixForm[Simplify[Inverse[T] . A . T]]

This shows that the formula, X[n], for the nth term of the Fibonacci sequence is a linear combination of powers, say nth, of the two eigenvalues.  A formula for X[n] can now be determined as follows.

Clear[c1, c2, e1, e2];
{e1,e2} = Eigenvalues[A];
sol=
  Solve[ 
    {c1 + c2 == 0, c1 e1 + c2 e2 == 1},
    {c1, c2} 
  ]

{c1, c2} = {c1, c2} /. sol[[1]]

This shows that a closed form formula for the nth term of the Fibonacci sequence is given by

OutputForm[c1 e1^n + c2 e2^n]

        1 - Sqrt[5] n     1 + Sqrt[5] n
       (-----------)     (-----------)
             2                 2
     -(--------------) + --------------
          Sqrt[5]           Sqrt[5]

Since e1 is less than one in absolute value, the first term will go to zero as n gets large.  Thus for large n, the Fibonacci sequence is approximately c2*e2^n -- exponential!

To implement this closed form description as a Mathematica definition is easy.

F4[n_] := Expand[c1 e1^n + c2 e2^n]

Notice that there is a considerable amount of algebra being done by Mathematica when using this formula to compute large terms in the sequence.

Malthusian Growth

Thomas Malthus (1766-1834) is generally credited with the idea that populations tend to grow exponentially.  Exponential growth (or decay) occurs whenever the rate of change over time of a variable is directly proportional to the value of the variable.  To see this in Mathematica, we use will K as the constant of proportionality.

solution1=
  DSolve[{P'[t] == K P[t], P[0] == a}, P[t], t]

We can plot several solution curves at once to see the behaviour of the solutions.

K = 1;
Plot[
  Evaluate[
    Table[
      P[t] /. solution1,
      {a,.4,4,.2}
    ]
  ],
  {t,0,4},
  PlotRange -> {0,20}
]
    

Logistic Growth

According to Malthus, although populations tend to increase exponentially, resources which support populations (e.g. food) tend to grow only linearly.  As the exponential growth of the population outstrips the linear growth of the resources, the resulting constraints would damp the population growth.

We can build in this dampening effect of excessive growth into the differential equation of the previous section.  Instead of P' being linear in P, we will make it quadratic in P.  This approach is called the logistic equation.  It was introduced by the Belgian statistician Pierre-Francois Verhulst (1804-1849).

P'[t_] = (1 - P[t]/K)*P[t]

With this model, the rate of change of P will be 0 when P=0 and when P=K.  Before we actually solve this DE, we will analyze it graphically using vector fields.

<<Graphics`PlotField`

Clear[K,f,P];
K=1;
f[t_,P_] := (1-P/K)*P

PlotVectorField[
  {1,f[t,P]},
  {t,0,4},
  {P,0,2},
  Axes->True,
  AxesLabel->{"t","P"}
]

Note that P' does not depend at all on t.  It only depends on P.  So all the columns of vectors in the above vector field are identical.

Clear[K,a,P]
K=1;

solution = 
  Table[
    DSolve[
      { P'[t] == (1 - P[t]/K)*P[t], P[0] == a }, 
      P[t],
      t
    ],
    {a,.1,1.7,.2}
  ]

We have to get rid of the outer set of braces to provide suitable input for the ParametricPlot function.  We do this with the Flatten command shown below.

Plot[
  Evaluate[P[t] /. Flatten[solution,1]],
  {t,0,4}
]

If, instead of having P' quadratic in P, we make P' cubic in P.  We can build a model which has three stationary points:  P = 0, T, K  (T<K).  The equilibrium point P=T represents the threshhold of a population - a number below which the population becomes extinct.  For 0<P<T, we want P'<0.  For T<P<K, we want P'>0.  And for P>K, we want P'<0.  A cubic function with these properties is easily constructed.

K=5; T=1;
g[P_] := (P/T - 1)*(1 - P/K)*P
Plot[g[P],{P,0,6},AxesLabel->{"P","P'"}]

We use this cubic function to assign values to P' in a differential equation.  The we plot solution curves for various values of the intial condition.

solution = 
  Table[
    NDSolve[
      { P'[t] == (P[t]/T - 1)*(1 - P[t]/K)*P[t], 
        P[0] == a 
      }, 
      P[t],
      {t,0,2}
    ],
    {a,.3,7,.3}
  ];

Plot[
  Evaluate[P[t] /. Flatten[solution,1]],
  {t,0,2},
  PlotRange->{0,7}
]

Lotka-Volterra Predator-Prey Populations

Vito Volterra was an Italian mathematician who lived from 1860-1940.  His son-in-law, Humberto D'Ancona was an Italian biologist who, in 1926, completed a statistical study of fish populations in the Adriadic Sea.  D'Ancona asked Volterra if there was a mathematical model which could explain the increase in predator fish (and corresponding decrease in prey fish) which he observed during the World War I period.  Within a couple of months, Volterra produced a series of models for the interaction of two or more species.  Alfred J. Lotka was an American biologist and actuary who independently produced many of the same models.

The simplest of their models can be described by letting

  x(t) = predator population at time t;
  y(t) = prey population at time t;
  a     =  the natural rate of decay of the predator population if there are no prey;
  b     =  the natural rate of growth of the prey population if there are no prey;
  c     =  the efficiency with which predators convert prey encounters into offspring,
  d    =  the rate of decrease of prey due to encounter with predators.

The simplest Lotka-Volterra system for modeling predator-prey interactions is

  x'[t] = -a x[t] + c x[t] y[t]
  y'[t] =  b y[t] - d x[t] y[t]

For a specific example, we consider the following.

eqn={x'[t]== -x[t] + x[t] y[t],
     y'[t]== 3y[t] - x[t] y[t],
     x[0]==x0,
     y[0]==1} ;

In order to graph several solution curves, we make a table of solutions.

tbl =
  Table[
    NDSolve[
      eqn, 
      {x,y},
      {t,0,6}
    ],
    {x0,.2,3,.3}
  ] ;

Then we can plot the solutions together.

ParametricPlot[
  Evaluate[
    { x[t],y[t] } /. Flatten[tbl,1]
  ],
  {t,0,6}
]

There are many variations on the Lotka-Volterra equations.  Here is one example which introduces an additional quadratic term in the rate of change of y.

eqn={x'[t]== -x[t] + x[t] y[t],
     y'[t]== 3y[t] - x[t] y[t] - y[t]^2,
     x[0]==x0,
     y[0]==1} ;

tbl =
  Table[
    NDSolve[
      eqn, 
      {x,y},
      {t,-3,3}
    ],
    {x0,.2,1.7,.3}
  ] ;

ParametricPlot[
  Evaluate[
    { x[t],y[t] } /. Flatten[tbl,1]
  ],
  {t,-3,3}
]

Competitive Exclusion

Instead of a predator-prey relationship, two species may compete with each other for the same niche in an ecosystem.  How can this relationship be modeled?  As you would expect, there are many ways to do this.  Perhaps the most direct way is to modify the predator-prey system as follows:

  x'[t] =  a x[t] - b x[t] y[t]
  y'[t] =  c y[t] - d x[t] y[t]

You get an opportunity to explain the differential equations in the exercises below!

First, we get a feel for the solutions to the system by plotting the vector field.

For a specific example, we consider the following.

eqn={x'[t]== 0.1 x[t] - 0.01 x[t] y[t],
     y'[t]== 0.2 y[t] - 0.05 x[t] y[t],
     x[0]==x0,
     y[0]==y0} ;

In order to graph several solution curves, we make a table of solutions.

tbl =
  Table[
    NDSolve[
      eqn, 
      {x,y},
      {t,0,20}
    ],
    {x0,1,9,1},
    {y0,1,19,2}
  ] ;

Then we can plot the solutions together.

ParametricPlot[
  Evaluate[
    { x[t],y[t] } /. Flatten[tbl,2]
  ],
  {t,0,20},PlotRange->{{0,10},{0,20}}
]

Exercises

Exercise 1

Discuss the different methods described in this notebook for computing the terms of the Fibonacci sequence.  Analyze the run times of the methods.  You can use the Timing function to gather data.  You can also count the number of operations required for each of the methods.

Exercise 2

Starting in 1790, the U. S. population has been measured every 10 years.  The following data have been collected (see http://www.census.gov/dmd/www/resapport/states/unitedstates.xls)

We want to find a model for the U. S. population growth.  To begin, we can build the data set from the above table and look at a point graph.

years   = Table[1790 + 10n, {n,0,21}];
pop     = {3.929214,5.308483,7.239881,9.638453,
            12.860702,17.063353,23.191876,31.443321,
            38.558371,50.189209,62.979766,76.212168,
            92.228496,106.021537,123.202624,132.164569,
            151.325798,179.323175,203.302031,226.542199,
            248.709873,281.421906};
data    = Table[{years[[i]], pop[[i]]}, {i,1,22
}];

ListPlot[
  data,
  AxesOrigin -> {1780,0},
  Prolog -> PointSize[.016]
]

Part (a)

Find an exponential model for the above data.  First, use only the data from the period 1790 through 1860.  Once you have this exponential model built, examine how well it conforms to the data after 1860.  Comment.  Graph the data points and the exponential curve on the same set of axes for comparison.

Next, build an exponential model which gives a best fit for all the data.  Again, graph the data and the curve on one set of axes for comparison.  Comment.

You may find the NonlinearFit function in the package Statistics`NonlinearFit` helpful.  You will need to consult the Help section to get information about using the NonlinearFit function.  You should carefully read the information about the function to get an understanding what it is doing.  

To get started with this process, you need to know the general form of the function you are trying to fit.  Although this is easy for the exponential case, you can let Mathematica do if for you anyway.  This will be useful in more complicated examples (like the logistic fit problem later in this exercise).

So we have one parameter, r, to vary in order to find the best fit to the dataset.

You can also leave the initial population unspecified and let the power of the NonlinearFit function use that extra variability to improve the exponential fit to the data.

Part (b)

Using the same population data, build a logistic model,  P' = r (1-P/K)P.  Use the same method as illustrated in Part (a).  Graph the logistic model together with the data points on one set of axes and discuss the suitability of the model.  Then extent the time axis forward in time and illustrate the carrying capacity (the K parameter in the model) for the model.  What is the value of K?  Can you offer explanations (based, say, in historical facts) for any deviations you see between the model and the data?

Exercise 3

Give a thorough explanation of the equations used in the basic predator-prey model.  You should fully explain the motivation of the terms used in the system of differential equations.

Exercise 4

In the variation of the predator-prey model discussed in the lab material above, discuss the effect of the term on the prey population.  Try to relate your discussion to topics we have already discussed.  HINT:  Try to find a way to use the word logistic.

Exercise 5

The equation dP/dt = r*P*(1-P/K) - H gives a logistic type model of a population with natural rate r, carrying capacity K and constant harvesting rate H (r, K, H positive).

Part A

The value H = rK/4 is called the critical harvesting rate.  Show that the population beocomes extinct if H exceeds this critical harvesting rate.

Part B

Plot portraits of solution curves in the quadrant P>0, t>0 for three cases: no harvesting (H=0); subcritical harvesting (use H=rK/8); and supercritical harvesting (use H=rK/2).  Let r=.001, K=1000, and choose various values for the initial population.

Exercise 6

Chapters 3 and 4 of the Herod, Shonkwiler, Yeargers text deal with populations.  Look over these chapters.  Carefully read section 4.4 on predator-prey problems, especially the discussion on stability around equilibrium points.

Part A

Now consider the differential equations

     x'[t] == -x[t] + x[t] y[t]
     y'[t] == 3y[t] - x[t] y[t] - y[t]^2
      x[0] == x0
      y[0] == y0

which were discussed briefly in the lab material.  How many equilibrium points does this system have?  Starting at the point (1,2), graph the solution curve as t ranges from 0 to 12.  Try a few other starting points and do the same thing (modify the upper bound on t if needed).  What do you conclude?  Can you think of a real-world situation which would make the y^2 term plausible (perhaps not!)?  Explain.

Part B

Now consider a system in which the prey species is limited with a logistic growth model, say

     x'[t] == -x[t] + x[t] y[t]
     y'[t] == y[t] (3-y[t]) - x[t] y[t]
      x[0] == x0
      y[0] == y0

Draw solution curves beginning at several different starting points and explain what is going on in the system.

Exercise 7

Consider the pair of differential equations given by

  x'[t] =  0.5 (1 - 0.25 x[t]) - 0.025 x[t] y[t]
  y'[t] =  0.1 (1 - 0.125 y[t]) - 0.006 x[t] y[t]

If x[t] and y[t] represent populations at time t, what happens to the two populations as time goes by?  Draw a set of solutions in the xy-plane.  Try to superimpose the nullclines on the solutions so you can see how they relate.  Explain why the results are different than in the competitive exclusion example given in the section.