Introduction to Ordinary Differential Equations
APAC-ANU Teaching Module
Stephen Roberts, Daniel MacKinley, and Peter Humburg

Contents

1  Introduction
2  A Population Problem
    2.1  A Very Simple ODE
    2.2  Testing the solution against the data
    2.3  Numerical solutions in SCILAB
    2.4  Numerical versus Analytic methods
    2.5  Compartmental Models
3  A more complex ODE - Limited Population growth
    3.1  Formulating the Differential Equation
    3.2  Equilibrium Solutions
    3.3  Logistic Equation
4  Systems of ODEs - Modelling Predators
    4.1  Getting a handle on coupled ODEs
    4.2  Numerical solutions to the fox and rabbit problem
    4.3  Phase Plane Analysis
    4.4  Equilibrium Populations and Null Clines
    4.5  Explaining the oscillations
    4.6  Adding Logistic Growth
5  Extended Predator-Prey Model
    5.1  Functional Response
    5.2  Numerical Response
    5.3  Implementing the Extended Model
6  Case Study: Competition, Predation and Diversity
Version of Document for Printing

1  Introduction

This is a brief introduction to simple Ordinary Differential Equations (ODEs), using them to model physical process and to the considerations of approximating their solution computationally. Ordinary Differential Equations are equations that relate the rate of change of some quantity (such as position, velocity, concentration) with respect to a quantity such as time, with the value of the quantity. In particular the derivative of a quantity such as position with respect to time may be proportional to the position. Unlike simple algebraic equations, the solution to a differential equation involves not just finding a particular number which satisfies the equation, but a particular function which satisfies the equation. ODEs arise naturally in a huge number of areas, including physics, biology, chemistry and ecosystems. They crop up in any area that the rate of change of something is related to how much of it there already is. In the language of calculus, differential equations relate the first or higher derivatives of an unknown function to the function itself. We examine what ODEs mean in one very particular class of problems: modelling the populations of various species of living things in an ecosystem. As we all know, the more animals there are making babies the faster the population of animals grows. We can thus expect differential equations to be useful at describing this. ODE population models are hopefully easy to understand yet behave in a fairly interesting way. More importantly, the techniques that we develop are generally applicable to very large classes of differential equations. In this module we will first see how differential equations can be used to model a range of population models. We will introduce the idea of a compartmental model. A number of useful theoretical tools will be introduced, such as the phase plane diagrams. We will use SCILAB, an open source numerical computing environment to approximate the solution of our differential model, and also see how to use SCILAB's more sophisticated ODE solvers. SCILAB's syntax is very similar to MATLAB. If you know either of these software packages you should have no problems to follow the examples in this module. The latest SCILAB version is available for free download from http://www.scilab.org. If you want to learn more about SCILAB, try the SCILAB primer and the numerical methods tutorials at http://comptlsci.anu.edu.au/tutorials.html.

Throughout this module you will find inserts like this. They contain questions and and small exercises about the material covered in the section above them. Try to answer them. They should help you to get some practical experience with ODEs.

 

2  A Population Problem

Let's look at some (fake) population figures. Say we wanted to predict what would happen with an accidental release of rabbits into a previously rabbit-free National Park. Fortunately for us, there are records available of what happened last time that rabbits were released into a rabbit free property, and we can use these records to try and find the rule that they obey. From a nearby national park, we get the following population figures:

Month	Rabbit population	Monthly increase	Percentage increase
January	200	-	-
February	220	20	10%
March	243	23	10.5%
April	268	25	10.3%
May	296	28	10.4%
June	326	30	10.1%
July	358	32	9.8%
August	395	37	10.3%
September	434	39	9.9%
October	478	44	10.1%
November	527	49	10.3%
December	577	50	9.5%
January	634	57	9.9%

What is the pattern here? How can we simulate it computationally? How do we begin to tackle this? One way might be to look at how the data changes each month. We look at the percentage rate of change of population: It seems that the rabbit population increases by approximately 10% each month.

2.1  A Very Simple ODE

We write this using calculus: It seems that perhaps we can model this rabbit population using some function which increases by 10% of its current value each unit of time. Let us represent the population of rabbits at a given time with a function X(t). This would be represented by

dX dt = 0.1X,     X(0) = 200.

Writing our problem like this, relating a single unknown function of one variable, X(t) and its derivatives, gives us what we call an Ordinary Differential Equation - an ODE. The solution to this ODE is:

X(t) = 200 e0.1t

How do we get that? It's funny you should ask that: We notice that our equation can be written:

1 X dX dt = 0.1

Integrate

⌠ ⌡ 1 X dX dt dt = ⌠ ⌡ 0.1 dt   so    ⌠ ⌡ 1 X dX = 0.1 ⌠ ⌡ dt

So lnX = 0.1t + C for some arbitrary constant C. Taking exponentials X(t) = A e^0.1t, where A=e^C. Using the fact that X(0) = 200 we get 200 = Ae^−kt₀ and A=200e^−(0)(0.1) Thus

X(t) = 200e0.1(t)

Now we have solved the equations we set ourself. Does it model the data?

2.2  Testing the solution against the data

To see if our solution really fits the data we use SCILAB to produce a plot of the data as well as of the function X(t) = 200e^0.1(t). Take a look at this piece of SCILAB code:

function x=f(t)
  x=200*exp(0.1*t)
endfunction

tData = 0:12;
tFunction=0:0.1:12;

yData = [200;...
         220;...
         243;...
         268;...
         296;...
         326;...
         358;...
         395;...
         434;...
         478;...
         527;...
         577;...
         634];
         
yFunction = feval(tFunction,f);

xbasc();

plot2d(tData,yData,style=[-1,2]);
plot2d(tFunction,yFunction,style=[3,3]);
legend("Measured data","Analytical solution",2);

This script produces a plot that shows the measured data as black crosses and the function X(t) as a green line. The result is shown in Figure 1. Figure 1: Comparison of the data with the analytical solution of our equation X(t) = 200e^0.1(t). Success! Or at least, close enough to it. We appear to have produced a model for a rabbit population, and solved a differential equation along the way. Examining Figure 1 we notice that our analytical solution is diverging slightly from the measured data. But the result is still accurate enough for our purpose. Note that our `solution' in this case still has a variable in it - t. We haven't found a single specific number, but a function that describes the behaviour of the populations of rabbits for all times t > 0.

Can you explain the difference between the analytical solution and the data? Try to think of any assumptions we made on the way. Can you find a better choice of rate parameter which provides a better fit to the data.

 

2.3  Numerical solutions in SCILAB

There is another way to get a solution to the rabbit problem. Instead of solving the ODE explicitly we can use SCILAB to solve the differential equation for us using the following code:

function xDot=g(t,x)
  xDot=0.1*x;
endfunction

tData = 0:12;
tFunction = 0:0.1:12;

xODE = ode(200,0,tFunction,g);

xbasc();

plot2d(tData,yData,style=[-1,2]);
plot2d(tFunction,xODE,style=[3,3]);
legend("Measured data","Numerical solution",2);

This produces the graph in Figure 2. Figure 2: Comparison of the data with the numerical solution of the differential equation [dX/dt] = 0.1X. It looks identical to Figure 1, and, in fact, it is to 5 decimal places! The trick here is: SCILAB can estimate numerical solutions to the ordinary differential equation without us having to go to the trouble of solving them. How? Can SCILAB tell us explicitly that it has found that the solution to our ODE is X(t) = 200e^0.1(t)? Well, no. What SCILAB can do is tell us how many rabbits there are at any time t using its internal numerical ODE solver, without ever knowing what the explicit formula is. Examining our program a little closer, this little snippet

function xDot=g(t,x)
  xDot=0.1*x
endfunction

is SCILAB-ese for

dX dt = 0.1X

and this command

xODE=ode(200,0,tFunction,g)

is enough to find the solution numerically for a given set of time values.

The numerical and analytical solutions to the rabbit problem are very similar but not quite identical. Can you produce a plot in SCILAB that reveals the difference?

 

2.4  Numerical versus Analytic methods

We have just found a way of telling SCILAB: "Give me an estimate of what will happen to the population of rabbits that is constantly increasing at a rate of 10% per month, if they start with a population of 200." Without ever mentioning the explicit solution to our problem, SCILAB has done a very passable job. What is the difference between the explicit formula solution and this computed numerical solution then? The thing is that the SCILAB version is just as easy for more complicated differential equations, and even for ones that are impossible to find numerical solutions for. However, the numerical solution is not exactly the same as the analytic solution. As we pointed out before, in this case they agree to 5 decimal places- this is enough to be identical to the human eye, but very far from identical! However, this sort of accuracy might be enough - in the case of our rabbit population model, where the population figures are a little imprecise anyway, we would probably never notice the difference. How do we know if we have an extra 0.01 of one rabbit? In any case, working out how accurate this ODE solver is, and choosing how accurate to try to make it is a big subject - and one of the major goals of this course.

2.5  Compartmental Models

How can we relate these to the differential equation that we have constructed? Let's pick apart our assumptions. In particular, let's look at what we assume is happening when we blithely model 'growth' in a rabbit population. For a start, rabbits are not immortal. Obviously, rabbits are dying and being born all the time. The population growth that we have described will only occur if the birth rate outweighs the death rate. If we had a very high death rate for some reason, say disease or starvation, we might expect to find a decreasing population of rabbits. Our rabbit population model is one of the class of models using ordinary differential equations known as compartmental models.

     net rate of change of rabbit population      =      rate of births      −      death rate     

     rate of births      =βX(t).

     death rate      =αX(t).

Adding these together, we get the following differential equation

dX dt = (β−α) X

If we set r=β−α then we may write this as

dX dt = r X

Our initial problem was precisely an equation describing a model with r=0.1. This technique of considering the inflow and outflows to the rabbit population will be used later on when we consider more complex sets of equations.

inflow       −−−−→    world        outflow −−−−→

Staying on our hypothetical rabbit population, inflow becomes rabbit births, and outflow becomes rabbit deaths. Considering the model in this fashion is useful if we want to extend it by having multiple compartments with complex interaction - for example, by looking at the interaction of predators and prey, or even more complicated systems - say rabbits, pasture, cows, foxes, native marsupials...

     rate of change of population      =      rate of births      −      rate of deaths     

3  A more complex ODE - Limited Population growth

It doesn't take much thought to realise that with our current model the rabbit population will increase forever, until eventually there are more rabbits than there can possibly fit in one little national park. If our model were correct then we should expect there to be exponentially expanding masses of rabbits to be forcing their way into the room right now. This is clearly (hopefully?) not the case. Let's say that, in the wild, the rabbit population exhibits a different sort of behaviour: one where they start out growing exponentially, but as they approach the capacity of the land to bear them, they begin to slow down their rate of population change - say, as they start to consume all the vegetation, rabbits have a greater risk of staving to death, litters are less likely to survive, rabbits are more likely to fight... If we extended our original rabbit population data, that sort of pattern of growth might gradually approach some finite number - the `carrying capacity' - which denotes the largest number of rabbits that can live on a piece of land indefinitely. But actually, this is an ecological fallacy - There is in general no set carrying capacity for a given ecosystem, and assuming there is one usually has disastrous consequences! But we're here for the maths, not the wildlife management, right? In short, we're looking for something like the S-curve in Figure 3. Figure 3: Development of rabbit population with carrying capacity. The population is gradually approaching the maximum number of individuals. This graph reflects the following assumptions:

• When rabbit populations are small, they will grow approximately exponentially.
• As an (initially small) rabbit population grows, individuals eventually will compete for limited resources.
• A given environment can support only a limited number of individuals.

This limit number is called the carrying capacity (denoted K). In our logistic equation, similarly, we are witnessing not a sudden freeze in the population of rabbits, where a bunch of immortal bunnies hang about producing no offspring. Rather, we are assuming that at the carrying capacity the birth rate has equalled the death rate, and so the growth of the population stops.

3.1  Formulating the Differential Equation

With that in mind, we might try to construct an improved model, using differential equations again: As before we assume

     rate of births      =βX(t).

Now assume that the per capita death rate depends on population

     per capita death rate      =α+γX(t).

A larger population results in a higher per capita death rate, for all the reasons we just listed. So

     rate of deaths      =(α+γX(t))X(t).

Now, for convenience we choose r=β−α. Then our final differential equation looks like this:

dX dt =rX−γX2    ( = βX−αX−γX2)

Does this give us a limited population size? Has it produced a model of rabbit populations on land with a limited carrying capacity? We answer that by looking at the equilibrium solutions.

3.2  Equilibrium Solutions

If we have formulated the equations we set out to, we would like to know that at some number of rabbits, the population will not be growing, i.e. [dX/dt]=0. Finding static population levels like this is looking for equilibrium solutions to the differential equation. We look for solutions to the equation

dX dt = 0 i.e. rX−γX2 = 0

Solving

rX−γX2 = 0 ⇒ X(r−γX) = 0 ⇒ X=0    or     X=r/γ

There are two equilibrium solutions. The first, X(t)=0 correspond to no rabbits at all - and since it takes rabbits to make rabbits, it makes sense that one possible way of having a steady population of rabbits is to have no rabbits at all. The second equilibrium solution, X=r/γ is more interesting. This is precisely the carrying capacity that we sought earlier. That is K=r/γ. Hence our differential equation becomes:

dX dt = rX   1− γ r X   dX dt = rX ( 1−X/K )

These fixed points are clearly depending on the initial populations of rabbits. However, the solution as a whole does not. That is, we can have the same differential equation describing our rabbits' population for a different starting population. This is important! One of the common characteristics of ODEs is that solutions depend on both the differential equations, but also on the constraints that we set on the problem - the initial conditions and the boundary values. When we find a particular function that satisfies both the initial conditions of our problem and the differential equation that defines it, we say that we have solved the initial value problem (the IVP) for that ODE. In general certain characteristics of ODE solutions are given by the ODE itself, and certain ones are given by the initial conditions. For example, in this case, the equilibrium points of our solutions are given by the ODE that defines it. But whether the rabbit population ever approaches this equilibrium value requires us to know the initial value.

To see this, ignore for a moment the fact that rabbits are supposed to have non-negative populations and look at what happens to the function if we set X(0)=−10.

 

3.3  Logistic Equation

If

dX dt =rX(1−X/K)

Then we can solve the equation by separation of variables, much as before! The solution to this ODE is given by

X(t)= K 1+m e−rt

where m=[K/(x₀)]−1

4  Systems of ODEs - Modelling Predators

We can take this technique further, solving several related functions simultaneously. The quintessential example of this are the Predator-Prey ODEs. We will consider a particular example of these, the Lotka-Volterra equations. Let's examine the consequences of introducing a new species whose growth is determined by how many rabbits there are to eat...

     Rate of change of rabbitpopulation      =    Rabbitbirth rate    −    Rateof fox rabbit kills   

     Rate of change of foxpopulation      =    Foxbirth rate    −    Rateof fox deaths   

A little thought reveals that we will need not one but two differential equations to reflect this. Since there are two populations, there are two related functions, one for each. Let X(t) be the number of prey per unit area and Y(t) be the number of predator per unit area. Let β₁ > 0 be the prey per-capita birth rate and α₂ be the predator per-capita death rate. Let the prey deaths per-capita be c₁ Y(t) (i.e. dependent on the number of predators). Let the rate of predator births per capita be c₂ X(t). So, rate of prey births = β₁ X(t), rate of predator deaths = α₂Y(t). Rate of prey killed by predator = c₁ X(t)Y(t) and rate of predator births = c₂X(t)Y(t). The term X(t)Y(t) can be thought of s representing the chance of an encounter between prey and predators at a given moment. Each 'encounter' entails a fox eating a rabbit, and a decrease in the rabbit population, and a smaller increase in the fox population (because a fox surely eats more rabbits than it has babies in its life...)

dX dt = β1 X − c1 X Y. (1) dY dt = c2 X Y − α2 Y. (2)

This is the Lotka-Volterra predator-prey system, and was proposed independently in the 1920s by two ecologists. This equation has, in fact, no analytic solution - as much as we might play with calculus, we will not find an equation which encapsulates this without an ordinary differential term.

4.1  Getting a handle on coupled ODEs

What will the solution look like? If X = 0 (meaning, there is no prey), [dY/dt] = −α₂ Y, which implies exponential decay for the predators. That is, the predators die out. And if Y=0 (no predators) then [dX/dt] = β₁ X, which implies that the prey grow exponentially - exactly as in our very first rabbit population model. In fact, in the case that Y(t) = 0 and β₁ = 0.1 we have found exactly the same model. In the general case, though, with arbitrary values for the constants, if we wish to observe the behaviour of the equations over time, we require SCILAB to come to our aid.

4.2  Numerical solutions to the fox and rabbit problem

By now it should be fairly routine to produce SCILAB code to calculate and display a numerical solution for some arbitrarily chosen constants and initial conditions.

//define our constants

beta1 = 1;
alpha2 = 0.5;

c1 = 0.01;
c2 = 0.005;

//give the initial values
y0 = [200;80];

//define the predator-prey system of equations
function PopsDot=predprey(t,y)
  PopsDot(1) = beta1*y(1)-c1*y(1)*y(2);
  PopsDot(2) = c2*y(1)*y(2)-alpha2*y(2);
endfunction

//evaluate at the given time steps
t = [1.d-5:0.01:60];

y = ode(y0,0,t,predprey);

//plot the results (prey = blue curve)
xbasc();

plot2d(t,y(1,:),style = 4);
plot2d(t,y(2,:),style = 6);
legend("Rabbit Population","Fox Population");

Figure 4: Lotka-Volterra predator-pray system. The coupled rabbit and fox populations show oscillating behaviour. This produces the graph shown in Figure 4. Notice the oscillations in both populations. This makes some sense intuitively and reflects the behaviour observed in some natural systems.

This model has four parameters. Instead of the values we used above try α2 = 1.5. What happened? Can you explain this behaviour? Try other sets of parameters. How do the different parameters influence the result?

 

Can we quantify these oscillations analytically?

4.3  Phase Plane Analysis

We get an idea of how X and Y behave with respect to one another by plotting the `phase plane', which shows the direction of motion of the system at a given moment. This is perhaps best demonstrated by example. Figure 5 shows the phase plane for the Lotka-Volterra system discussed above. Figure 5: Phase plane for the Lotka-Volterra system. The `arrows' in the graph are produced by the following command in SCILAB:

xbasc();
plot2d(y(1,:),y(2,:),style=4);
fchamp(predprey,0,[0:20:240],[0:20:200]);

a = get("current_axes");
a.x_label.text = "Rabbits";
a.y_label.text = "Foxes";

The first command, plot2d we have seen before. In this case it plots X versus Y rather than X versus t. fchamp is a new one. It plots, for a selection of points, the direction that the system of populations of creatures is evolving. For each particular pair of possible values of the populations, X and Y, the rate of change of the system along each axis is given by [dX/dt] and [dY/dt] respectively. Each of the arrowed vectors in our plot is hence [dX/dt] [dY/dt]. If both [dX/dt] and [dY/dt] are `smooth' enough, this plot will show us how these values evolve in time. Indeed, the blue path depicts a system of populations whose curve follows those vectors. Note that the periodic behaviour we observed in Figure 4 corresponds to a closed curve in the phase plane diagram. Such a curve is called a limit cycle. For each different set of initial conditions, we will get a (possibly different) path through the phase plane diagram - for example, Figure 6 is a similar plot, but with starting fox populations of 1, 4, 16, and 64. Figure 6: Phase diagram for starting fox populations of 1, 4, 16, and 64.

Can you work out which starting population was used for each curve in Figure 6?

 

4.4  Equilibrium Populations and Null Clines

Now lets take a closer look at the phase plane diagram. You can see that some of the arrows are parallel to one axes. This means that either [dX/dt]=0 or [dY/dt]=0. The points where at least one of our differential equations is 0 all occur along a few lines. These are called null clines. In Figure 7 the null clines were added to a schematic representation of the phase plane diagram for our predator-prey model. Figure 7: The phase plane diagram is partitioned into four regions, each representing a different general behaviour of the two differential equations. They separate the plane into four regions. In each region the two differential equations show a different behaviour. This is indicated by the arrows in Figure 7.

Region I:

X decreases and Y increases.

Region II:

X decreases and Y decreases.

Region III:

X increases and Y decreases.

Region IV:

X increases and Y increases.

If both differential equations are 0 there is no change in the system. These points are equilibrium solutions, just like the one we encountered in section 3.2. These equilibrium solutions lie at the intersection of two null clines.

Before you continue to the next section take another look at Figure 7. There are four intersections between horizontal and vertical lines. Only two of these intersections are equilibrium solutions. Can you tell which ones? Why is there no equilibrium at the other two points?

 

4.5  Explaining the oscillations

If initially the populations are x₀ and y₀ and we travel along the curve we eventually return to x₀, y₀. To find the equilibrium points set [dX/dt] = [dY/dt] = 0. Then

X(β1 − c1 Y) = 0 (3) Y(c2X −α2) = 0 (4)

If X = 0 , then from equation 4, Y = 0. If β₁− c₁ Y = 0, Y = β₁/c₁ and X = α₂/c₂. The equilibrium points are (0, 0) and (α₂/c₂,β₁/c₁). Figure 8 shows a phase plane diagram with added null clines and marked equilibrium points. Figure 8: Phase plane diagram for the Lotka-Volterra predator-prey model with added null clines. Dark blue lines indicate points where [dY/dt] = 0 and red lines indicate points where [dX/dt] = 0. Equilibrium points are marked with crosses. Lets look at the equilibrium point (α₂/c₂,β₁/c₁). At this point the amount of rabbits that is eaten by foxes during a time interval is equal to the amount of rabbits born during that time. And the number of foxes born is equal the the number of foxes that die. Both populations are in an equilibrium. But what happens when we move away from this equilibrium point? Consider region II We know that X < [(α₂)/(c₂)] and Y > [(β₁)/(c₁)]. We can write this as X = [(α₂)/(c₂)]− δ₁, Y = [(β₁)/(c₁)] + δ₂ for δ_1/2 > 0. Using equations 1 and 2 we get

dX dt = X(β1 − c1( β1 c1 + δ2)) = −X c1 δ2 < 0 dY dt = Y(c2( α2 c2 − δ1) −α2) = −Y c2 δ1 < 0

Thus, both the predator and the prey population decrease in region II. Using similar calculations it is easy to see what happens in the other regions:

X > α₂/c₂, Y > β₁/c₁:

From Equation 2 Y increases, From Equation 1 X decreases.

X < α₂/c₂, Y > β₁/c₁:

From Equation 2 Y decreases, From Equation 1 X decreases.

X < α₂/c₂, Y < β₁/c₁:

From Equation 2 Y decreases, From Equation 1 X increases.

X > α₂/c₂, Y < β₁/c₁:

From Equation 2 Yincreases, From Equation 1 X increases.

Examine the phase plane diagram and describe the behaviour of the predator-prey system moving along trajectories in anti-clockwise direction. What happens when the system enters a new region? Try to identify the four regions from the phase plane diagram in Figure 4.

 

4.6  Adding Logistic Growth

Recall

dX dt = β1 X − c1 X Y. = X(β1 − c1 Y) dY dt = c2 X Y − α2 Y = Y(c2 X − α2) .

Previously we considered the case where β₁ = 1, α₂ = 0.5, c₁ = 0.01, c₂ = 0.005. These equations produced oscillations in time for various conditions and parameter combinations. The model assumes the prey grows exponentially in the absence of the predator. We could replace the growth term for the prey population with a density dependent growth with carrying capacity K. Then

dX dt = β1 X   1− X K   − c1 X Y. dY dt = c2 X Y − α2 Y.

We represent this in SCILAB as

K = 800;

function PopsDot=predprey(t,y)
  PopsDot(1)=beta1*y(1)*(1-y(1)/K)-c1*y(1)*y(2);
  PopsDot(2)=c2*y(1)*y(2)-alpha2*y(2);
endfunction

Which gives us a slightly different set of results. Figure 9 shows the development of the predator and prey populations over time. We see a damped oscillation where the changes in the rabbit and fox populations get less extreme as they approach an equilibrium. Figure 9: Introducing a carrying capacity for the prey population results in a damped oscillation. Figure 10: Phase plane diagram for the Lotka-Volterra model with carrying capacity for the prey population. The system approaches an equilibrium. This kind of behaviour is easy to observe on our phase plane diagram (Figure 10). We can see how the rabbit-fox system `spirals' towards an equilibrium point. It's a little more difficult to show that the solutions approach an equilibrium point, rather than circling it. We need a little more experience with partial derivatives, and a neat tool called the Jacobian matrix. An introduction to that is a little beyond the scope of this document!

Find the null clines and equilibrium points for this new model. You can check your results by comparing them to Figure 12

 

Figure 11: Phase plane diagram for different starting populations. We can plot multiple starting populations with comparative ease within SCILAB.

K=400

xbasc();
y0 = [2;200];
y = ode(y0,0,t,predprey2);
plot2d(y(1,:),y(2,:),style=1);

y0 = [10;200];
y = ode(y0,0,t,predprey2);
plot2d(y(1,:),y(2,:),style=4);

y0 = [50;200];
y = ode(y0,0,t,predprey2);
plot2d(y(1,:),y(2,:),style=3);

fchamp(predprey2,0,[0:20:500],[0:20:240]);

You can examine the output of this script in Figure 11. Note that we changed the carrying capacity for the prey population to 400 to get a nicer plot. SCILAB can't easily compute null-clines for us, but we have already done enough algebra to persuade it to at least plot them for us. The lines marked in red are null-clines, and the crosses mark equilibrium points, on the intersection of an X and a Y null cline. Figure 12: Phase plane diagram for different starting population. Null clines are shown by red lines, the equilibrium points are marked with crosses. The result is quite revealing. We see that the behaviour of solutions is, in some sense, conditioned by the equilibrium solutions and null-clines.

5  Extended Predator-Prey Model

In section 4 we modelled a predator-prey system with a simple Lotka-Volterra model using the two differential equations

dX dt = r1X−c1XY dY dt = c2XY−α2Y

This model resulted in oscillatory behaviour of the system. We then refined our model by adding logistic growth for the prey population:

dX dt = r1   1− X K   −c1XY dY dt = c2XY−α2Y

This led to a loss of the oscillatory behaviour. Instead we observed damped oscillations while the system approached an equilibrium. But in nature we observe oscillatory behaviour! Now we want to extend our model even further. Maybe we can get the oscillatory behaviour back. In practice, field ecologists measure

• functional response     F(X)
• numerical response     N(X,Y)

Functional response is the rate at which a single predator kills prey, as a function of prey population. Numerical response is the per capita growth rate of the predators which may depend on prey and predator populations.

dX dt = r1X   1− X K   −F(X)Y dY dt = N(X,Y)Y

5.1  Functional Response

Assume a single predator takes time t_h to handle each prey. Consider a time interval [0,T] with

T=Th+Ts

where T_h is the time taken to handle all prey encountered before time T and T_s is the time left to search for new prey.

     number of prey eaten in timeT by one predator      =c1TsX(t)=Ne

     rate of prey eaten byone predator      = c1TsX(t) T

Now T=T_h+T_s=N_et_h+T_s=c₁T_sX(t)t_h+T_S.

     rate of prey eaten byone predator      = c1X 1+thc1X             F(x)→ 0 as x→ 0 F(x)→ 1 th as x→∞     

This is known as Holling type II response function.

5.2  Numerical Response

Now we assume logistic growth for predator population

r2   1− Y K2   .

But assume the carrying capacity K₂ is proportional to prey density, ie K₂=c₃X. Hence

N(X,Y)=r2   1− Y c3X  

Putting everything together we get an improved model:

dX dt = r1X   1− X K   − c1XY 1+thc1X dY dt = r2Y   1− Y c3X  

This is known as Holling-Tanner predator-prey model. This extended model now has six parameters:

r₁ > 0:

Prey birth rate.

K > 0:

Prey carrying capacity.

c₁ > 0:

Prey per capita death rate.

0 < t_h < T:

Time a single predator takes to handle each prey.

r₂ > 0:

Predator birth rate.

K₂ = c₃X > 0:

Predator carrying capacity.

5.3  Implementing the Extended Model

Again we can implement this in SCILAB.

//define our constants
r1 = 0.16;
r2 = 0.05;

c1 = 0.006;
c3 = 0.6;

th = 2;

K = 800;


//give the initial values
y0 = [200;92];

//define the predator-prey system of equations
function PopsDot=predprey3(t,y)
  PopsDot(1) = r1*y(1)*(1-y(1)/K)-(c1*y(1)*y(2)/(1+th*c1*y(1)));
  PopsDot(2) = r2*y(2)*(1-y(2)/(c3*y(1)));
endfunction

//evaluate at the given time steps
t = [0d-5:0.1:600];

y = ode(y0,0,t,predprey3);

//plot the results (prey = blue curve)
xbasc();

plot2d(t,y(1,:),style = 4);
plot2d(t,y(2,:),style = 6);
legend("Rabbit Population","Fox Population");

The result of this is shown in Figure 13. Examining the SCILAB code above you can see that we used r₁ = 0.16, r₂ = 0.05, c₁ = 0.006, c₃ = 0.6, t_h = 2, K = 800. This set of parameters gives us a nice periodic behaviour. Unfortunately these parameter values are not making much sense biologically. It seems unlikely that a fox needs two month to hunt down and eat one rabbit. Also note that the oscillations are now much slower than before. One period takes approximately 100 months! Figure 13: Development of a Holling-Tanner predator-prey system over time. Figure 14: Phase plane diagram for a Holling-Tanner predator-prey model.

Try different sets of parameters. How does the system react to different parameter choices? How does this compare to the different parameter sets for Lotka-Volterra model? Can you find a more realistic parameter set that results in periodic behaviour?

 

Figure 14 shows the phase plane diagram for this system. Again you can see the limit cycle. Note that we have adjusted the starting population of foxes from the previous example to get a starting condition on the limit cycle.

Try different starting populations for this model. Examine the resulting phase plane diagrams. What do you notice?

 

6  Case Study: Competition, Predation and Diversity

By now you should be fairly comfortable with handling ODEs. We have considered different population models, implemented them in SCILAB, and analysed the results. Now it is time to try your skills on a more complex example. In the last few sections we have considered systems with one prey and one predator species. Now we will add a second prey species into the mix. Consider two prey species N₁, N₂ and one predator N₃. In our case, we will be looking at foxes and rabbits again. The rabbits are competing for the same resources and we should consider the total number of rabbits for the per capita growth rate of both prey species. The foxes are eating both rabbit species but they may have a preference for one. We allow for this possibility by using different parameters for the interaction between foxes and the two different rabbit species. The general form for a system like this is

dN1 dt = N1(ε1−α11N1−α12N2− α13N3) dN2 dt = N2(ε2−α21N1−α22N2− α23N3) dN3 dt = N3(−ε3+α31N1+α32N2− α33N3)

This has a simple competition model between N₁ and N₂, and a predator-prey model between N₁ and N₃, and N₂ and N₃.

Using the above equations formulate a predator-prey model using a combined carrying capacity K1 for both rabbit species, i.e. K1 ≥ N1+N2, and a carrying capacity K2 for the foxes that depends on the amount of available prey. Implement this in SCILAB. Choose a suitable set of parameters and plot the result. Produce a time as well as a phase plane diagram. Interpret the results. Try different parameter sets to see how they influence the result.

 

Two species competition models often lead to extinction of one species. Nature seems to be more stable. Stability is enhanced with more interaction of species.