Predator Prey Systems

In Predator-Prey systems of differential equations, we have two populations that are inter-related.

We will look for plots of population numbers both over time, and populations in the phase plane (which means the axes in the plots are the population functions that satisfy the differential equations).

In this introductory look at systems of differential equations, we will look at numerical solutions although some parts can be solved by hand.

The function H(t) represents a host that is preyed upon, such as tuna.

The function P(t) represents a predator, such as sharks.

This is an example of an autonomous DE, since the derivatives to not have an explicit time dependence.

(dH/dt)=α H - β H P,
(dP/dt)=δ H P-γ P.

Equilibrium Points

Solve a system of equations to determine the equilibrium points. This can easily be done by hand.

(dH/dt)=α H - β H P=0
(dP/dt)=δ H P-γ P=0

H(α -β P)=0
P(δ H-γ)=0

H=P=0 (no populations)

α -β P=0
δ H-γ=0

H=(γ/δ)
P=(α/β)

Clear[α,β,δ,γ,P,H,c,eqpts,equil,sol]

equation1=α H - β H P==0
equation2=δ H P-γ P==0
eqpts=Solve[{equation1,equation2},{H,P}]

equil={H,P}/.eqpts[[1]]

P as a function of H

Since the system of DE is autonomous, we can take the ratio of the DEs and eliminate the parameter t.

((((dP/dt)))/(((dH/dt))))=((δ H P-γ P)/(α H - β H P)),
(dP/dH)=((δ H P-γ P)/(α H - β H P)),
(dP/dH)=((P(δ H -γ ))/(H(α  - β P))),
(((α  - β P))/P)dP=(((δ H -γ ))/H)dH,
∫(((α  - β P))/P)dP=∫(((δ H -γ ))/H)dH
α ln |P| - β P = δ H - γ ln|H|+ C

Here is the implicit solution which relates P and H.

Notice as a family of curves we have a constant in the solution.

sol[H_,P_,c_]=α*Log[Abs[P]]-β*P==δ*H-γ*Log[Abs[H]] + c

To plot using a computer, we need to pick values for the parameters in the model.

α=10.7;
β=0.667;
δ=0.9;
γ=3.3;

equil
sol[H,P,c]

Here is a plot that shows the solution curves for a variety of c values.

equilpt=Show[Graphics[{Purple,Disk[equil,{0.75*15/40,0.75}]}]];
plot[c_]:=ContourPlot[Evaluate[sol[H,P,c]],{H,0,15},{P,0,40},BaseStyle->20,FrameLabel->{"Host","Parasite"},PlotPoints->50,AspectRatio->1,ContourStyle->{{Red,Thick}}]

Manipulate[
Show[plot[c],equilpt]
,{c,10,20,0.1}]

Show[plot[15],plot[16],plot[17],plot[18],plot[19],equilpt]

What’s missing from our analysis is the time dependence of the solution, since we eliminated t to get the above solution.

To get the time dependence, we need to use numerical solutions, either

direction fields

stream plots

numerical solutions with Euler’s method

numerical solutions using NDSolve

Streamplot or Direction Field to Show Time Dependence

Streamplots are like direction fields, but rather than showing the direction of the solution at specific points in the phase plane, they essentially follow a trajectory from a given starting point.

Streamplots are still numerical approximations to the solution.

fieldplot=VectorPlot[{α H-β H P,δ H P-γ P},{H,0,15},{P,0,40},AspectRatio->1,VectorScale->{Tiny,Tiny,None},BaseStyle->20,FrameLabel->{"Host","Parasite"},VectorStyle->Orange]

streamplot=StreamPlot[{α H-β H P,δ H P-γ P},{H,0,15},{P,0,40},StreamPoints->20,BaseStyle->20,FrameLabel->{"Host","Parasite"}]

Here we compare the two, and see they show the same information, just from a different perspective.

Show[fieldplot,streamplot]

Manipulate[
Show[plot[c],equilpt,fieldplot]
,{c,10,20,0.1}]

Numerical Solution of IVP Using Euler's Method

Here is one way to implement Euler's method for systems in Mathematica, using recursion relations:

Clear[f1,f2,x,y,t]
α=10.7;
β=0.667;
δ=0.9;
γ=3.3;

f1[H_,P_]=α H - β H P
f2[H_,P_]=δ H P-γ P

h=0.01;
t[0]=0;
x[0]=5.0;
y[0]=5.0;
t[i_]:=t[i]=t[i-1]+h
x[i_]:=x[i]=x[i-1]+h*f1[x[i-1],y[i-1]]
y[i_]:=y[i]=y[i-1]+h*f2[x[i-1],y[i-1]]

euler=Table[{x[k],y[k]},{k,0,200}];

eulerplot=ListPlot[euler,PlotStyle->{Green,Thick},Joined->True,AxesLabel->{"x","y"},PlotRange->{All,All},BaseStyle->20]

Notice the solution suffers from serious numerical error, since if we decrease h we get a solution that at least to the eye appears to be periodic.

There is still error in this solution, but the error is much smaller.

Clear[f1,f2,x,y,t]
α=10.7;
β=0.667;
δ=0.9;
γ=3.3;

f1[H_,P_]=α H - β H P
f2[H_,P_]=δ H P-γ P

h=0.0001;
t[0]=0;
x[0]=5.0;
y[0]=5.0;
t[i_]:=t[i]=t[i-1]+h
x[i_]:=x[i]=x[i-1]+h*f1[x[i-1],y[i-1]]
y[i_]:=y[i]=y[i-1]+h*f2[x[i-1],y[i-1]]

euler=Table[{x[k],y[k]},{k,0,20000}];

eulerplot=ListPlot[euler,PlotStyle->{Green,Thick},Joined->True,AxesLabel->{"x","y"},PlotRange->{All,All},BaseStyle->20]

Numerical Solution of IVP Using NDSolve

Here is a numerical solution to the system of differential equations, in terms of interpolating functions.

Notice we need to include initial conditions to get this numerical solution, here we chose (5,5).

The last line of code is just the Mathematica syntax that is required to store the solution from NDSolve in a way that we can easily use.

DE1=H'[t]==α H[t]-β H[t] P[t]
DE2=P'[t]==δ H[t]P[t]-γ P[t]

solDE=First[{H,P}/.NDSolve[{DE1,DE2,H[0]==5,P[0]==5},{H,P},{t,0,5}]]

Here is a plot of the two populations as a function of time.

Plot[{solDE[[1]][t],solDE[[2]][t]},{t,0,5},PlotRange->All,PlotStyle->{Orange,Purple},PlotLegends->{"Host","Parasite"},BaseStyle->20,AxesLabel->{"time","Population"}]

Here is a plot of the two populations in the phase plane.

ParametricPlot[{solDE[[1]][t],solDE[[2]][t]},{t,0,5},PlotRange->All,PlotStyle->{Thick},BaseStyle->20,AxesLabel->{"Host","Parasite"},AxesOrigin->{0,0}]

Here we determine the value of c for the solution that passes through (5,5), so we can show an animation of the populations changing.

Solve[sol[5,5,c],c]

max=5000;
point[i_]:={solDE[[1]][5/4*i/max],solDE[[2]][5/4*i/max]};

Manipulate[
Show[fieldplot,equilpt,plot[14.6971],Graphics[{Black,Disk[point[i],{1.0*15/40,1.0}]}],AspectRatio->1,PlotRange->{{0,15},{0,40}}]
,{i,0,max,0.05}]

Notice the ball speeds up and slows down! This was not evident from any of the other solutions we had.

Modifications to Predator-Prey Models

As with any population model, the predator-prey model is just a model of the situation and we can modify the model to be more realistic. For example, consider the following model, which includes a term that deccreases the rate of change of the Hosts when the host population becomes larger than K (a logistic type growth)

(dH/dt)=α H (1-(H/K))- β H P,
(dP/dt)=δ H P-γ P.

There are now three equilibrium points:

Clear[α,β,δ,γ,P,H,c,K,eqpts,equil,sol]

equation1=α H (1-(H/K))- β H P==0
equation2=δ H P-γ P==0
eqpts=Solve[{equation1,equation2},{H,P}]

The Streamplot shows us the behaviour. Notice the predators die out in this model.

α=10.7;
β=0.667;
δ=0.9;
γ=3.3;
K=50;

eqpts

streamplot=StreamPlot[{α H(1-(H/30))-β H P,δ H P-γ P},{H,0,15},{P,0,40},StreamPoints->20,BaseStyle->20,FrameLabel->{"Host","Parasite"}]

Notice we do not get periodic solutions anymore!

The solutions spiral towards an equilibrium point.