4 - Predator-Prey models#
A first model#
Of course, species in the natural world do not just engage in friendly competition for resources. Nature is ‘red in tooth and claw’, with many species relying on others as a primary food source: foxes and rabbits; lions and zebras; ladybirds and aphids. Such predator-prey interactions will have a similar structure to the competition models we considered above, but now one species benefits from the interaction. Calling the prey density \(N\) and the predator density \(P\), the simplest model would look like this:
The prey has an independent birth rate \(a\), but its death rate depends on the predator density (with parameter \(b\)) with more predators leading to higher prey death. In contrast, the predator has an independent death rate \(c\), but its birth rate depends on the prey density (with parameter \(d\)), with predators needing to eat enough prey to get the energy to reproduce.
Let’s find the equilibria of the system. There is an extinction equilibrium at \((N,P)=(0,0)\), and a coexistence equilibrium at \((N,P)=(c/d,a/b)\). The stability of these steady states is calculated, as before, by finding the Jacobian of the system:
Have a go
Write out the Jacobian matrix for the predator-prey system.
Click for solution
We can see that the trivial equilibrium, \((N,P)=(0,0)\), is unstable (since \(a>0\)). What about the coexistence equilibrium? This evaluates to:
Here, \(T=0\) and \(D=ac>0\). These conditions mean that the eigenvalues at this equilibrium have zero Real part and are purely Imaginary. The behaviour near the steady state is then for ‘centres’, a family of neutrally-stable closed orbits around the equilibrium (note, these are not limit cycles). We therefore predict that these two populations will constantly cycle: with low predation the prey density increases; this produces more food for predators, and so the predator numbers begin to rise; this in turn starts to push the prey numbers back down; finally, with their food supplies falling, the predator numbers also drop down. A phase portrait of these dynamics would look like this,
Figure: Sketch of phase portrait for predator-prey model showing a centre.”
A more realistic model#
While it is nice to have demonstrated the existence of cycles, these are not structurally stable, which suggests something may not be quite right. From our previous models we might already think of two additions to make the model more realistic. One is that we ignored any density-dependence, meaning prey in particular could grow to infinite levels. Another is that we assumed predation was linearly dependent on prey density, whereas we previously argued it should be saturating. Including these assumptions modifies the model to,
We can see we still have the trivial equilibrium, and some coexistence equilibrium \((N^*,P^*)\). There is now also a new ‘prey-only’ equilibrium at \((N,P)=(a/\alpha,0)\) . The Jacobian of this modified system is,
Prey-only equilibrium
In this case the Jacobian reduces to,
Therefore stability depends on the growth rate of predators. If the growth rate is negative the prey-only equilibrium is stable, but if it is positive the prey-only equilibrium is unstable.
Coexistence equilibrium
Using what we know about this equilibrium to our advantage, the Jacobian becomes,
From here we can find that the determinant is given by, \(D=bdAN^*P^*/(N+A)^3>0\), and the trace by, \(T=-\alpha N^*+bN^*P^*/(N+A)^2\). Stability therefore depends on whether the trace is positive or negative, which we can see depends on the balance of intraspecific competition in the prey and predation: if competition is high relative to predation, we will have a negative trace and the equilibrium is stable; if predation is high relative to competition, we will have a positive trace and the equilibrium is unstable.
Notice that if we vary one parameter, \(\alpha\), say, we could move the system from being a stable spiral to being an unstable spiral, since we move from a negative to a positive trace. This is a change in stability, and therefore a bifurcation, but one that we have not yet seen in practice and which cannot occur in one-dimensional systems.
Hopf bifurcation#
The previous three bifurcations are the standard types that occur in one and higher dimensions. There are further bifurcations that only occur in two or more dimensions. One of particular importance is the Hopf bifurcation. In this case changing a parameter turns a stable spiral in to an unstable spiral. At this point a unique, stable closed orbit emerges from the equilibrium. This results in the population cycling/oscillating. It is not very easy to draw a bifurcation diagram in this case, but we can see the process by looking at the phase portraits as we vary a parameter.
The example in the figure below comes from the predator-prey model above. Here we see we initially have a stable spiral, but that as we increase \(\alpha\) the equilibrium loses stability and starts to spiral away. However, this outward trajectory does not continue to spiral out but tends towards a closed orbit (which can be seen clearly in the right-hand figure). An initial condition from further away starts to spiral in but also approaches this closed orbit. This is a much stronger result for population cycles than the centres we saw in the simple model. Hopf bifurcations have important consequences as populations which are fluctuating may be much harder to measure or control.
Figure: Example of a Hopf bifurcation in the predator-prey model, equations (1.4.3)-(1.4.4) with a type II functional response. Parameter values: \(a=5\), \(b=1\), \(c=0.5\), \(d=2\), \(A=1\).
Have a go: coding
To use the code, click the rocket icon at the top of the page, and then click Colab. This will take you to a live version of this page where you can run the code.
Scroll down to the code cell below and either click in the cell and hit Shift-Enter, or click on the ‘Play’ icon in the top-left of the cell. You may see a warning about loading content from Github - click ‘Run anyway’.
This code will produce a time-course and phase portrait for the dynamics of a predator and its prey. The dynamics are governed by the following ordinary differential equations:
After you run the code, velow the code cell you will see a set of sliders to choose the values of the 6 parameters, and boxes for you to enter two lots of initial conditions for \(N\) and \(P\). You will see two plots appear - a time course for the Prey and the Predator, and a phase portrait (including nullclines (red/black) and 2 trajectories (blue/green)). You can then change the parameter values or initial conditions to update the plots.
Note that the sliders/plot in this webpage are static and do not change - you need to load the Colab version for them to work.
Show code cell source
################
# HAVE A GO: CODING
################
# Import the necessary libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import ipywidgets as widgets
from ipywidgets import interactive
aw=widgets.FloatSlider(min=0.5, max=5, step=0.25, value=4, description='a')
bw=widgets.FloatSlider(min=0.5, max=5, step=0.25, value=1, description ='b')
cw=widgets.FloatSlider(min=0.5, max=5, step=0.25, value=0.5,description='c')
dw=widgets.FloatSlider(min=0.5, max=5, step=0.25, value=2,description='d')
alphaw=widgets.FloatSlider(min=0, max=5, step=0.25, value=4,description='alpha')
Aw=widgets.FloatSlider(min=0, max=5, step=0.25, value=1,description='A')
N0_1w=widgets.BoundedFloatText(value=1.5,min=0,max=5,description='1st N(0)')
P0_1w=widgets.BoundedFloatText(value=6,min=0,max=10,description='1st P(0)')
N0_2w=widgets.BoundedFloatText(value=0.5,min=0,max=5,description='2nd N(0)')
P0_2w=widgets.BoundedFloatText(value=4,min=0,max=10,description='2nd P(0)')
u1 = widgets.HBox([N0_1w, P0_1w])
u2 = widgets.HBox([N0_2w, P0_2w])
#Lotka-Volterra dynamics
def predprey(x,t,a,b,c,d,alpha,A):
N=x[0]
P=x[1]
dN = N*(a-b*P/(N+A)-alpha*N)
dP = P*(-c+d*N/(N+A))
return [dN,dP]
def runpredprey(a,b,c,d,alpha,A,N0_1,P0_1,N0_2,P0_2):
N0=[N0_1,P0_1]
N1=[N0_2,P0_2]
tc = np.linspace(0, 50, 1000)
Nc = odeint(predprey, N0, tc,args=(a,b,c,d,alpha,A))
Nd = odeint(predprey, N1, tc,args=(a,b,c,d,alpha,A))
plt.rcParams['figure.figsize'] = [12, 4]#
plt.rcParams.update({'font.size': 16})
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(tc, Nc[:,0], "r", label="Prey")
ax1.plot(tc, Nc[:,1], "k", label="Predators")
ax1.plot(tc, Nd[:,0], "r:")
ax1.plot(tc, Nd[:,1], "k:")
ax1.set(xlabel='Time', ylabel='Densities')
ax1.legend()
ax1.axis([0,50,0,10])
nn=np.linspace(0,10,100)
nnull=(a-alpha*nn)*(nn+A)/b
pnull=A*c/(d-c)
ax2.plot(Nc[:,0],Nc[:,1],'b')
ax2.plot(Nd[:,0],Nd[:,1],'g')
ax2.plot(nn,nnull,'r')
ax2.axvline(x=pnull)
ax2.axis([0, 5, 0, 5])
ax2.set(xlabel='Prey', ylabel='Predators')
ax2.axis([0,2,0,10])
plt.show()
return()
out = widgets.interactive_output(runpredprey, {'a': aw, 'b': bw, 'c': cw, 'd': dw, 'alpha': alphaw, 'A': Aw, 'N0_1': N0_1w, 'P0_1': P0_1w, 'N0_2': N0_2w, 'P0_2': P0_2w})
display(alphaw,Aw)
display(aw,bw,cw,dw)
display(u1,u2, out)
3 key points#
The simplest predator-prey model produces ‘centres’ - structuarally unstable cycles.
A more realistic predator-prey model can produce stable ‘limit cycles’.
We cann the transition from an equilibrium to a limit cycle a Hopf bifurcation.