Introduction

In this chapter we will explore how to do some basic programming in Mathematica.  Our vehicle for doing this will be exploring the phenomenon of extinction in a discrete predator-prey model for the behavior of coyotes and jackrabbits in Northern Utah/Southern Idaho.  The basic equations we will examine are

(Subscript[x, n+1])=(1+a) (Subscript[x, n]) (1-b (Subscript[y, n])),

(Subscript[y, n+1])=(1+c) (Subscript[y, n])(1-((Subscript[y, n])/(k+d (Subscript[x, n])))).

The parameter and variables in the equations are:

({{Variable, Description}, {(Subscript[x, n]), Population density of bunnies, year n}, {(Subscript[y, n]), Population density of coyotes, year n}, {a, Birth rate of uneaten bunnies}, {b, Predation rate}, {c, Small - population growth rate of coyotes}, {d, Carrying capacity of coyotes per bunny}, {k, Carrying capacity of coyotes without bunnies}})

Given an initial set of population values, ((Subscript[x, 0]), (Subscript[y, 0])), the two equations recursively generate a series of population values, ((Subscript[x, n]), (Subscript[y, n])), in successive years.

As we will see below, a very common occurrence for solutions to these equations is for one or both species to go extinct, or hit zero population values.  Largely this is because the model above always allows predation to occur; the term b (Subscript[x, n]) (Subscript[y, n]) reflects the expected value of predation based on "mass action" assumptions: the idea that the rate of individual encounters is proportional to the population densities.  In other words, this model gives the bunnies no place to hide, and extinction in this model is analogous to the impact of habitat destruction on a predator-prey system.  The prevalence of extinction in this sort of model reflects the fact that prey species depend very strongly on heterogeneous environments - having places to hide themselves or their young, the presence of refuges, and having some way to avoid contact with predators.

We will start our analysis by defining the predator and prey equations into Mathematica.

prey[xx_,yy_]:=(1+a)*xx*(1-b*yy)
pred[xx_,yy_]:=(1+c)*yy*(1-yy/(k+d*xx))

Steady States and Stability

We will begin by determining equilibrium solutions for the population, and then try to determine if these steady state solutions are stable in the above model.  In a discrete model like the one we are using, a steady state is reached when there is no year-to-year change in the population values.  The equilibrium solution will satisfy

x = prey(x, y)

y = pred(x, y)

where the prey and pred functions are the ones defined above.  We can solve these equations using Mathematica in the following manner:

sstates=Solve[{x==prey[x,y],y==pred[x,y]},{x,y}]

Here we see that there are three steady state solutions.  The only solution that does not have the extinction of one or both species is,

ss2=sstates[[2]]

The stability analysis of a steady state in a discrete model proceeds much as the stability analysis of fixed points in ordinary differential equations.  The behavior of solutions in the vicinity of a fixed point is controlled by the linearization of the system near the fixed point.  Write the system of equations in vector form as

(Subscript[v, n+1])=F((Subscript[v, n]))

where (Subscript[v, n])=(([(Subscript[x, n]),(Subscript[y, n])])^T) is a 2 by 1 vector and F is a vector valued function with components given by prey and pred, the right hand side of the model, above.  If a steady state solution is given by a vector w, and we write the solution corresponding to a (small) perturbation as

(Subscript[v, n])=w+ϵ (Subscript[u, n]),

then the behavior of the perturbation is controlled by the linear system

(Subscript[u, n+1])=J(w)(Subscript[u, n]).

Here J is the Jacobian of F.  We say the equilibrium solution is stable if the size of the perturbation, (Subscript[u, n]), shrinks over time.  That is, for stability we require that

|(Subscript[u, n+1])|⩽|(Subscript[u, n])|

for all perturbations.  From linear algebra, since

(Subscript[u, n+1])=J(w)(Subscript[u, n]),

we know that the stability condition is satisfied only if all of the eigenvalues of J have (complex) magnitude less than one.  If any eigenvalue has modulus bigger than one, the fixed point w is unstable (or the perturbation will grow).

Let's see about the stability of the fixed point determined above.  First let's define the Jacobian,

jm=D[{prey[x,y],pred[x,y]},{{x,y}}]

If we evaluate this matrix at the relavant steady state solution we get,

jmss=Simplify[jm/.ss2];
MatrixForm[jmss]

We can use Mathematica to find the eigenvalues of this matrix as follows:

evals=Eigenvalues[jmss]

Here we see a pair of eigenvalues.  From the form of the eigenvalues it is not clear whether the fixed point is stable or not.  Let's assign values to all of the parameters except one, then examine how the stability is affected by variations in that one parameter.  Since birth rates for the bunnies are the most susceptible to climatic variations, let's assign values to all of the parameter except a.

b=0.18
c=1.53
d=1.25
k=1.03

With these values the eigenvalues now have the form:

evals

To analyze the stability as a function of a, we could plot the modulus of the eigenvalues vs a.

p1=Plot[Abs[evals[[1]]],{a,0,2},PlotStyle->Directive[Red,PointSize[0.015]],Mesh->45];
p2=Plot[Abs[evals[[2]]],{a,0,2},PlotStyle->Directive[Blue,Thick]];
p3=Plot[1,{a,0,2},PlotStyle->Directive[Black,Dashed]];
Show[p1,p2,p3]

Note that for values of a between about 0.2 and 1.25 the fixed point is stable, and that the fixed points for values of a outside that range are unstable.  That is, if the birth rate is either too large or too small, then the steady state solution is unstable.

Generating Solutions

In this section we will examine how to use the recursion relations to generate solutions.  We will do this by writing a Maple procedure that will except as input initial population values and return as output the yearly populations and the plot structures for the graphs of these population values vs time.  The following is one such way of writing said procedure:

pp[x0_,y0_]:=Module[{xplot,yplot},
xd=Array[0,51];
yd=Array[0,51];
xd[[1]]=x0;
yd[[1]]=y0;
Do[
xd[[jj+1]]=Max[prey[xd[[jj]],yd[[jj]]],0];
yd[[jj+1]]=Max[pred[xd[[jj]],yd[[jj]]],0],
{jj,1,50}];
xplot=ListPlot[xd,PlotStyle->Directive[Blue,Thick],Joined->True];
yplot=ListPlot[yd,PlotStyle->Directive[Red,Thick],Joined->True];
{xd,yd,xplot,yplot}
]

Ok.  Note that this module uses the pred and prey functions and generates outputs: xd, yd, xplot, and yplot.  Also note that this procedure will calculate populations for 50 years and that the right hand side of the recursion equations have been modified so that negative population can not be generated by the numerical calculations.

To use this procedure we will first need to assign values to all of the parameter.  Since we have assigned values to all of the parameters except a in a command line above, we will now choose a values for a.  First we will look at a choice of a for which the fixed point is unstable.

a=1.4

For this value of a, the corresponding steady state is

ss2

If we choose initial population values near the steady state and call the procedure we get,

sols=pp[3.2,3.0];

We have suppressed the output, but we now have a list whose first and second components are the populations and the third and fourth components are graphs of those po;ulations.If we wish to see the graphs we could,

Show[sols[[3]],sols[[4]],AxesLabel->{"Years","Pops"},PlotRange->All]

This shows that the bunny populations grows without bounds and the coyotes appear to be insignificant.  Lets just look at the coyote graph

Show[sols[[4]],AxesLabel->{"Years","Pops"},PlotRange->All]

It appears that the coyotes die out at about year 24.  Let's look at the populations for years 0 to 30.

xpop=Table[{k-1,sols[[1,k]]},{k,1,31}];
ypop=Table[{k-1,sols[[2,k]]},{k,1,31}];
xp=ListPlot[xpop,PlotStyle->Directive[Blue,Thick],Joined->True];
yp=ListPlot[ypop,PlotStyle->Directive[Red,Thick],Joined->True];
Show[xp,yp,AxesLabel->{"Years","Pops"},PlotRange->{0,10}]

Based on the above graph, we see that for these initial populations and value of a, the coyote population (in red) goes to zero at about year 24.  That is, the coyotes go extinct, leaving no predators for the bunnies, who proceeded to grow without bound.  This is a common occurrence when the coexistence fixed point is unstable and the bunny birth rates are high.  Additional simulations, which we won't pursue here, shows that when coexistence is unstable and bunny birth rates are low, then the most common catastrophe is that the bunny population goes extinct, while the coyotes subsist by eating cactus (or in the case of Idaho coyotes, potatoes).

If we try a value of a that generates a stable fixed point then we get the following:

a=1.0

Now the steady state is

ss2

and solutions to the model are

sols=pp[3.0,2.0];

Show[sols[[3]],sols[[4]],AxesLabel->{"Years","Pops"},PlotRange->{0,4}]

In this case the populations oscillate about the fixed point with the amplitude of the oscillation decreasing until  we have almost constant populations.  Since we know that bunny and coyote populations do not remain constant from year to year, we would conclude that either the birth rate for bunnies is not in a stable range, or our model is not quite sophisticated enough to describe the problem.  However, as with most real world models, the goal is not to reproduce every detail in the computer, but to be able to analyze general concepts (in this case about bunnies and coyotes).

The Spectrum of Extinction

From a management perspective, we might like to know the probable rate at which a species is going to go extinct when its refuges are removed.  Once land has been clear-cut, how rapidly must recovery occur to keep native species from going extinct?  The model above describes just such a refuge-less ecology, but the question is still not easy to answer.  Parameters in such a model (particularly birth rates) will vary from year to year with climatic changes, and measurements of population density for mobile species are error-prone to say the least.  To account for all this variability, we will do simulations in which the parameters and the initial values of the species vary randomly, but with known probabilities, and examine the probability that the ecosystem will collapse in a given year.  This probability we refer to as the 'spectrum of extinction.'

The most basic piece of information we need to know is not the multi-year history of a predator-prey system, but the year in which one or the other of the species dies off.  We will accomplish this using the   the Mathematica While loop and a simple counting algorithm.  First let's design a procedure that returns the year extinction occurs given initial values for the populations and a given value of a, the birth rate.  One such procedure is as follows.

Clear[a]

extinct[x0_,y0_,a0_]:=Module[{bn=x0,cn=y0,years,newbn,newcn},
a=a0;
years=0;
While[(bn>0 && cn>0 && years<50),
years++;
newbn=prey[bn,cn];
newcn=pred[bn,cn];
bn=newbn;
cn=newcn;
];
years
]

Notice that the While statement will increment the model until one of the populations goes extinct or the number of years reaches 50.  The output to this routine is the year where extinction occurs.

To test to see if this procedure works, try it with the values from the above analysis.  We should get a year of extinction of about 24.

extinct[3.2,3.0,1.4]

Good.  Next define and initialize an array to hold the information we will collect.  If we run the above procedure many times and count up how many times extinction occurred in year j, then we could build a relative frequency diagram of the likelihood that extinction will occur in a given time period.

num=Table[0,{j,1,50}];

Now that we have initialized our accumulator, we need only run the above procedure multiple times with random inputs (for the initial populations for now) and keep track of the year that extinction occurs.  Please note that Mathematica's  command returns a 12 digit positive integer.  The following loop runs the procedure 1000 times with initial bunny populations between 2.5 and 3.5, and initial coyote populations between 2.0 and 4.0.

Do[
icb=RandomReal[{2.5,3.5}];
icc=RandomReal[{2.0,4.0}];
year=extinct[icb,icc,1.4];
num[[year]]++,
{j,1,1000}
]

Now the variable, num, holds the frequency distribution of the spectrum of extinction.  For instance, out of 1000 simulations, the number of times that extinction occurred in year 24 is given by,

num[[24]]

Note that this value may (and probably will) change each time the loop is executed, so don't worry if you are following along and get a different value.  To plot the relative frequency of the year of extinction vs. the year we would do the following,

ListPlot[Table[{j,num[[j]]/1000},{j,1,50}],PlotStyle->Red,PlotMarkers->{"\[Square]",Small},PlotRange->All,AxesLabel->{"Year","Frequency"}]

While this graph brings up many interesting questions, the main question we wish to answer at this point is the average time to extinction.  This questions is answered by calculating the mean, or center of mass:

mass=Sum[num[[j]],{j,1,50}]
moment=Sum[j*num[[j]],{j,1,50}]

averagetimetoextinction=N[moment/mass]

From this we see that the mean time to extinction for this given value of a is about 14 years.  Given the multi-modal nature of the graph, this is a number to take with a grain of salt; a large fraction of the populations went extinct well before 14 years.

The Spectrum with Parametric Stochasticity

An effect which could potentially have much greater impact than uncertainty of population measurement would be variations in the parameters of the model, or parametric stochasticity.  In an ecological model like the one we are examining, many of the parameters are "known" only in the loosest sense, and may vary a good deal from year to year.  Take, for example, the birth rate of the bunnies.  In a very cold or dry year, the birth rate will be much lower because offspring have a much more difficult time surviving in the more difficult environment.  We might, therefore, view the value of a we examine as being the average of some yearly fluctuation due to climatic changes.  A completely sensible question, then, is how does fluctuations from year to year in the birth rates affect the spectrum of extinction?

First we need to create a basic solution procedure, like 'extinct' above.  In this procedure, however, we will build in the capacity for the birth rate, a, and the external carrying capacity, k, to fluctuate in some range.  The domain of the function we create will include initial population levels, x0 and y0, the center of the parameter ranges, a0 and k0, and the range of variation (or amount of +/- about a0 and k0), called arange and krange.  The range of the function would be the year in which either of the species goes extinct.  We will call this function 'extinct2', and define it below:

Clear[x0,y0,a,b,c,d,k,a0,k0,arange,krange,extinct2]
b=0.18;
c=1.53;
d=1.25;

extinct2[x0_,y0_,a0_,arange_,k0_,krange_]:=
Module[
{years=0,bn=x0,cn=y0,newbn,newcn},
While[(bn>0 && cn>0 && years<50),
years++;
a=a0+arange*RandomReal[{-1.0,1.0}];
k=k0+krange*RandomReal[{-1.0,1.0}];
newbn=prey[bn,cn];
newcn=pred[bn,cn];
bn=newbn;
cn=newcn;
];
years
]

Now that we have the function defined, it would be wise to test it and make sure it works.  If we set the range of variation for a and k to zero, we should get the same answer as before given we use the other inputs as before:

extinct2[3.2,3.0,1.4,0.0,1.03,0.0]

And the result is 24, just as it should be.

Now we can do some simulations, just as before, to determine the spectrum of extinction.  This time, however, we will allow for a fifty percent variation in birth rates and a ninety-nine percent variation in external carrying capacity.  Just to make sure it works we will only do a few (100) simulations first.

num=Table[0,{j,1,50}];
Do[
icb=RandomReal[{2.5,3.5}];
icc=RandomReal[{2.0,4.0}];
year=extinct2[icb,icc,1.4,0.65,1.03,1.029];
num[[year]]++,
{j,1,100}
]

A graph of the relative frequency of this run is as follows:

ListPlot[Table[{j,num[[j]]/100},{j,1,50}],PlotStyle->Red,PlotMarkers->{"\[Square]",Small},PlotRange->All,AxesLabel->{"Year","Frequency"}]

A run of 1000 iterations gives the following

num=Table[0,{j,1,50}];
Do[
icb=RandomReal[{2.5,3.5}];
icc=RandomReal[{2.0,4.0}];
year=extinct2[icb,icc,1.4,0.65,1.03,1.029];
num[[year]]++,
{j,1,1000}
]

ListPlot[Table[{j,num[[j]]/1000},{j,1,50}],PlotStyle->Red,PlotMarkers->{"\[Square]",Small},PlotRange->All,AxesLabel->{"Year","Frequency"}]

mass=Sum[num[[j]],{j,1,50}]
moment=Sum[j*num[[j]],{j,1,50}]

averagetimetoextinction=N[moment/mass]

Reviewing the results depicted in the above graphs and center of mass calculations, first notice there are differences between 100 and 1000 samples.  The more samples that are taken, the truer a representation of the actual distribution of collapse times.  The spectrum also displays a lot less natural periodicity, which is sensible since varying parameters smears out the natural period of the model.  The final thing to notice is that if we do the basic statistics to calculate the mean time to extinction, we find that the expected time to extinction has dropped to 7.7 years.  In this case, the mean value reflects the location of a true peak in the spectrum, representing the most likely time at which species have become extinct.  The fact that this number has dropped significantly from 14, the value without parametric stochasticity, illustrates a common-sense fact from ecology: environmental fluctuations exacerbate the effects of habitat destruction.

Exercise 10

One model for predator-prey interaction between bunnies and coyotes in the intermountain region of North America is

(Subscript[x, n+1])=(1+a)((Subscript[x, n])-b (Subscript[x, n]) (Subscript[y, n]))
(Subscript[y, n+1])=c (Subscript[x, n]) (Subscript[y, n])+(1-d) (Subscript[y, n])

where:

({{(Subscript[x, n]), The density of bunnies in year n}, {(Subscript[y, n]), The density of coyotes in year n}, {a, The birth rate of bunnies (about 0.8)}, {b, The rate of predation of bunnies by coyotes (about 0.0735)}, {c, The birth rate of coyotes in presence of bunnies (about 0.03)}, {d, Starvation rate of coyotes in absence of bunnies (about 0.25)}})

The parameter most susceptible to seasonal variations is the birth rate of bunnies, a.  Typical population density values are about 10 for the bunnies and 5 for the coyote populations.

1.

Write a Mathematica program using a procedure that will allow you to experiment easily with this system.  Inputs should be the population values, the birth rate of bunnies, and the death rate of coyotes.  The output from your program should be a single graph showing the two populations plotted simultaneously.  Try a couple of different death rates.  How does altering the death rate seem to change the dynamics?

2.

Write a program to examine the spectrum of extinction in this problem.  Your program should include the following features:

(a)  Some random fluctuation in initial measurements of population levels by 50% around some base measurement.
(b)  Random yearly fluctuations in birth and death rates by as much as 80% around some base level.
(c)  A relative frequency diagram which illustrates the number of extinctions per year over however many simulations you choose to do.
(d)  A calculation of the expected time to extinction in this model, using the data from your simulation.

Discuss the behavior of the spectrum of extinction and how it relates to the expected value you calculated.