Lab 10: The Volterra Predator-Prey Model

In this lab, you will plot the predator-prey model in both discrete and continuous time.

Let H[t] represent the numbers of prey and P[t] represent the numbers of predators.  The prey are assumed to grow exponentially, but they are
captured by predators at a rate that depends on both the availability of prey and predators.

The discrete time recursions are thus:

Prey:  		H[t+1] = H[t] + r H[t] - b H[t] P[t]
Predator:	P[t+1] = P[t] + c H[t] P[t] - d P[t]

which uses the following parameters:
	the growth rate of the prey in the absence of the predator (r),
	the capture rate at which predators contact and kill prey (b),
	the rate at which eaten prey are turned into predator babies (c),
	and the death rate of predators (d).

The analog in continuous time is:

Prey:  		dH/dt = r H[t] - b H[t] P[t]
Predator:	dP/dt =  c H[t] P[t] - d P[t]

There are two main ways of graphing these equations.  You can look at how each species (the predator or the prey) changes over time or you can look at the relationship of predators to prey.  This second type of graph is called a phase-plane analysis.

Graphs of discrete system

Here is a Mathematica simulation giving the number of
prey (popH) and the number of predators (popP) over time:

Clear[popH,popP]
popH[r_,b_,c_,d_,t_]:=popH[r,b,c,d,t]=
	popH[r,b,c,d,t-1]*(1+r-b*popP[r,b,c,d,t-1])
popP[r_,b_,c_,d_,t_]:=popP[r,b,c,d,t]=
	popP[r,b,c,d,t-1]*(1+c*popH[r,b,c,d,t-1]-d)
popH[r_,b_,c_,d_,0]:=500
popP[r_,b_,c_,d_,0]:=45

(The initial numbers are chosen aribtrarily, try using
different numbers to see if they affect the dynamics.)

Try plotting the number of prey over time using a species-time graph:

ListPlot[Table[popH[0.5,0.01,0.0002,0.1,t],{t,0,50}]]

Now try plotting the number of predators as a function
of the number of prey using a phase-plane graph
(each dot will represent a different time point):

ListPlot[Table[{popH[0.5,0.01,0.0002,0.1,t],
	popP[0.5,0.01,0.0002,0.1,t]},{t,0,50}],
	PlotRange->{{0,1000},{0,100}}]

Which is the x-axis and which is the y-axis?
[Hint:  Look at the initial numbers]

Graphs of continuous system

Here is a Mathematica program to numerically
evaluate the differential equations giving the number of
prey (nH) and the number of predators (nP) over time
in the continuous model:

Clear[system,nH,nP]
system[r_,b_,c_,d_,nH0_,nP0_,tmin_,tmax_] := 
	system[r,b,c,d,nH0,nP0,tmin,tmax] =
	NDSolve[{	nH'[t] == nH[t]*(r-b*nP[t]),
				nP'[t] == nP[t]*(c*nH[t]-d),
				nH[0]  == nH0,
				nP[0]  == nP0},
			{nH[t],nP[t]}, {t,tmin,tmax}]

Try plotting the number of prey over time (a species-time graph):

Plot[Evaluate[ nH[t] /. 
	system[0.5,0.01,0.0002,0.1,500,45,0,100] ],{t,0,100},
	PlotRange->{{0,100},{0,1000}}]

Now try plotting the number of predators as a function
of the number of prey using a phase-plane graph.
To do this, we'll use the Mathematica command 
ParametricPlot[{fx, fy}, {t, tmin, tmax}], which plots 
the function fy on the y-axis against fx on the x-axis
from time tmin to tmax.

ParametricPlot[Evaluate[{nH[t],nP[t]} /. 
		system[0.5,0.01,0.0002,0.1,500,45,0,100] ],
	{t,0,100},
	PlotRange->{{0,1000},{0,100}}]

Questions for the lab

Sketch one of the above phase-plane diagrams, making sure
to specify which axis corresponds to prey and which to predator.

In which direction do the spirals go?
[Try plotting the graphs for one, two, five, and ten generations.]

Can you make sense of this biologically?

In the discrete model, what is the effect of halving both
the intrinsic growth rate of prey (r) and the capture rate (b),
but keeping r/b constant?

In the continuous model, what is the effect of changing the initial
numbers of prey and predators?

Can you make sense of the fact that if the spirals form closed loops in
the continuous-time model that they must spiral out in the discrete-time
model?  [Think about the difference in time step.]

If the spirals formed closed loops in a discrete-time model, what
must they do in the analogous continuous-time model?