22 Lotka-Volterra predator-prey dynamics

This practical runs over two sessions, after the two predation lectures. In the first session you build the classic Lotka-Volterra predator-prey model in Excel and explore its cycles. In the second session you add two things that Neal (Predator-Prey Interactions) and the lectures identify as stabilising: a carrying capacity for the prey, and a refuge where some prey cannot be caught.

Learning outcomes:

• Build a predator-prey model in Excel from its equations, and check it against known values.
• Use the zero isoclines to explain the direction of the cycles and to predict how each parameter moves the equilibrium.
• Explain why the spreadsheet’s cycles grow, and show that this comes from the size of the time step rather than from the biology.
• Explain how a prey carrying capacity and a prey refuge stabilise the interaction.

22.1 Background: the model

We use Neal’s notation (Predator-Prey Interactions). The lecture slides show the same symbols next to Gotelli’s (V for H, α for a, β for c and q for d):

\[\frac{dH}{dt} = rH - aHP \qquad\qquad \frac{dP}{dt} = cHP - dP\]

Symbol	Meaning	Value in the worksheet
\(H\)	number of prey (“herbivores”)	starts at 1000
\(P\)	number of predators	starts at 20
\(r\)	prey growth rate (per prey, per unit time)	0.25
\(a\)	attack rate: prey caught per prey, per predator, per unit time	0.01
\(c\)	conversion efficiency: new predators per prey, per predator, per unit time	0.00008
\(d\)	predator death rate (per predator, per unit time)	0.1

The prey grow exponentially and are removed at a rate proportional to how often prey and predators meet (\(aHP\)). The predators are born at a rate proportional to the prey they catch (\(cHP\); \(c\) is smaller than \(a\) because many prey are needed to make one new predator) and die at a constant rate \(d\).

Zero isoclines. The prey stop changing when \(rH = aHP\), i.e. when \(P = r/a\): a horizontal line on a graph of \(P\) against \(H\). Below it (few predators) the prey increase; above it they decrease. The predators stop changing when \(cHP = dP\), i.e. when \(H = d/c\): a vertical line. To the right of it (many prey) the predators increase; to the left they decrease. The lines cross at the equilibrium \(H^* = d/c\), \(P^* = r/a\). Following the arrows round the four regions gives cycles, with the predator peaks a quarter of a cycle behind the prey peaks. For small cycles the period is about \(2\pi/\sqrt{rd}\).

The spreadsheet version. Excel cannot solve the differential equations directly. Instead it takes small steps of length \(\Delta t\) (cell F8, “dt”), each time adding the change over that step:

\[H_{t+\Delta t} = H_t + \Delta t\,(rH_t - aH_tP_t) \qquad\qquad P_{t+\Delta t} = P_t + \Delta t\,(cH_tP_t - dP_t)\]

Each step uses the numbers at the start of the step, so the populations always react a little late. With \(\Delta t = 1\) this lag is large enough to change the behaviour of the model, as you will see.

22.2 Your task

This practical runs over two sessions. In the first you build the basic model in Excel and explore its cycles; in the second you add a prey carrying capacity and a prey refuge, and then answer four questions.

22.2.1 Part 1 (session 1, about 10 minutes): build the model

Download the Excel worksheet here. It has one sheet for each part (Part 1 Basic model, Part 3 Prey K and Part 4 Prey refuge). Most of it is ready for you. The colours show where you work: yellow cells are where you type formulas, and orange cells hold the parameters you change.

Cells	What they hold
F4 to F8 (orange)	The parameters: \(r\), \(a\), \(c\), \(d\) and the time step \(\Delta t\) (F9 holds \(K\) or \(H_r\) in Parts 3 and 4)
B11 and C11	The starting numbers, \(H = 1000\) and \(P = 20\)
A11 to A3011	Time, going up by \(\Delta t\) each row
B4 to C8	The ends of the two zero isoclines and the equilibrium \((H^*, P^*)\), calculated from the parameters
Two graphs	\(H\) and \(P\) against time (first 300 time units), and \(P\) against \(H\) with both isoclines

Your job on the sheet Part 1 Basic model is to type the two formulas that move the populations forward one time step.

1. In cell B12 (yellow), type the formula for the prey:

   • English Excel: =MAX(0,B11+$F$8*($F$4*B11-$F$5*B11*C11))
   • Danish Excel: =MAKS(0;B11+$F$8*($F$4*B11-$F$5*B11*C11))

   Compare it with the prey equation above: B11 is \(H_t\), C11 is \(P_t\), $F$4 is \(r\), $F$5 is \(a\) and $F$8 is \(\Delta t\). The $ signs keep the formula pointing at the parameter cells when you copy it down. MAX(0, ...) stops the population going below zero.

2. In cell C12 (yellow), write the matching formula for the predators yourself. Use the predator equation: you need \(P_t\) (C11), \(H_t\) (B11), \(c\) ($F$6), \(d\) ($F$7) and \(\Delta t\) ($F$8).

3. Select B12 and C12 together and copy them down to row 3011 (double-click the small square at the bottom right of the selection). The graphs fill in as you go.

4. Check your sheet. Time 1 should give \(H = 1050\) and \(P = 19.6\), and time 2 should give \(H = 1106.7\) and \(P = 19.286\). If not, check the $ signs and which parameter each cell reference points to.

22.2.2 Part 2 (session 1, about 25 minutes): explore the model

For each task, write down your prediction first (a sentence or a quick sketch of the phase plane), then change the values and check.

1. Cycles. With the starting values, which peaks first, prey or predators? How long is it from one prey peak to the next? In the phase plane, which way does the trajectory go round (clockwise or anticlockwise), and why? Use the four regions marked out by the isoclines.
2. Equilibrium. Calculate \(H^* = d/c\) and \(P^* = r/a\). Start the model exactly there (B11 and C11). What happens? What happens if you start a little way away, and further away?
3. The time step. With \(\Delta t = 1\) the cycles grow until the prey are wiped out (when?). Now set \(\Delta t\) (F8) to 0.5 and then to 0.1, with everything else unchanged. The time-series graph always shows the first 300 time units. What happens to the growth of the cycles? What does this tell you about where the growing cycles come from? What would you expect from the continuous model in the lecture?
4. Moving the isoclines. With \(\Delta t = 0.1\), predict which isocline moves, and in which direction, if you double each of \(r\), \(a\), \(c\) and \(d\) in turn (one at a time, putting each back afterwards). Check with the phase plane.
5. (Optional) Volterra’s principle. A pesticide kills both prey and predators: it lowers \(r\) to 0.2 and raises \(d\) to 0.12. What happens to the equilibrium numbers of prey and predators? Is that what you would have guessed?

22.2.3 Part 3 (session 2, about 15 minutes): a carrying capacity for the prey

Real predator-prey systems rarely show the neutral cycles of the Lotka-Volterra model. Neal (Predator-Prey Interactions) uses the graphical model of Rosenzweig and MacArthur (1963) to show how other features of the biology change the outcome. Here we add two of them to the spreadsheet, one at a time. Keep \(\Delta t = 1\) unless a task says otherwise, and for each value you try, predict first (sketch the isoclines), then run the model to check.

Without predators, real prey populations do not grow for ever: they grow logistically towards a carrying capacity \(K\). The prey equation becomes

\[\frac{dH}{dt} = rH\left(1 - \frac{H}{K}\right) - aHP\]

and the prey isocline is no longer horizontal: it slopes down from \(P = r/a\) at \(H = 0\) to \(P = 0\) at \(H = K\). The predator equation does not change.

Use the sheet Part 3 Prey K (\(K = 2000\) in F9). The predator formula is already filled in.

1. In B12, type the new prey formula and copy it down to row 3011:

   • English Excel: =MAX(0,B11+$F$8*($F$4*B11*(1-B11/$F$9)-$F$5*B11*C11))
   • Danish Excel: =MAKS(0;B11+$F$8*($F$4*B11*(1-B11/$F$9)-$F$5*B11*C11))

   Check: time 1 should give \(H = 925\), \(P = 19.6\).

2. With \(K = 2000\), what happens to the cycles? Where do the populations end up? Compare this with the point where the isoclines cross.

3. Try \(K = 5000\) and \(K = 10000\) (a richer environment for the prey). What happens to the cycles, and to the final number of predators?

4. Try \(K = 1000\), which is smaller than \(d/c\). What happens to the predators, and why? Use the isoclines.

22.2.4 Part 4 (session 2, about 10 minutes): a refuge for the prey

Some prey can hide where the predators cannot reach them (burrows, crevices, dense vegetation). In the sheet Part 4 Prey refuge, the first \(H_r\) prey (cell F9) cannot be caught: the predators only meet the \(H - H_r\) prey outside the refuge. Both formulas are filled in, so you only need to change \(H_r\).

1. Set \(H_r = 0\). Check that you get the same results as Part 1.
2. Try \(H_r = 50\), 100, 200 and 400. What happens to the size of the cycles, and to the risk of extinction?
3. How does the refuge move the predator isocline? (Look at the equilibrium in cell B8.) Why do the predators need more prey at equilibrium?

22.2.5 Part 5 (session 2, about 10 minutes): the questions

Use what you found in Parts 2 to 4 to answer these questions. Write a few sentences for each, and use a sketch of the isoclines where it helps.

1. In the basic model, which parameters set the equilibrium numbers of prey and of predators? Explain in words why the equilibrium number of prey depends only on predator parameters, and the equilibrium number of predators only on prey and attack parameters.
2. Why do the cycles in your spreadsheet grow (with \(\Delta t = 1\)), when the continuous Lotka-Volterra model predicts cycles that keep the same size? What does this tell you about using discrete time steps in models?
3. How do a prey carrying capacity and a prey refuge each change the outcome? Explain each using the isoclines.
4. The Lotka-Volterra model assumes, among other things, that predators can eat without limit (no satiation). Name two other assumptions of the model that are unlikely to hold in real populations, and say how each might change the dynamics.
5. (Optional) Neal (Predator-Prey Interactions) describes Gause’s experiments with Paramecium (prey) and Didinium (predator). When Gause added a sediment where Paramecium could hide, Didinium ate all the prey in the open water and then starved, and the prey took over the microcosm. Your refuge model in Part 4 keeps both species. Which assumptions of the model could explain the difference? (Hints: in your model, what happens to the prey born in the refuge once there are more than \(H_r\)? And look at how low the predator numbers go in Part 1: can a real population have that many?)

22.3 More: the same model in R

Two optional R chapters (extra homework) solve the same model, with the same parameter values and symbols: Continuous-time Lotka-Volterra predator-prey model in R (the model in the lecture, solved accurately) and Discrete-time Lotka-Volterra predator-prey model in R (the spreadsheet model, with any time step).

22.4 Takeaways

• Predator-prey feedbacks can generate cycles without any outside forcing, with the predators a quarter of a cycle behind the prey.
• The parameters set the isoclines and the equilibrium (\(H^* = d/c\), \(P^* = r/a\)); in the continuous model the starting numbers set the size of the cycles.
• The growing cycles of the spreadsheet come from its time steps: with small steps they almost disappear.
• A prey carrying capacity and a prey refuge both damp the cycles towards a stable equilibrium.

22.5 References

Lotka, A. J. 1925. Elements of Physical Biology. Williams & Wilkins.

Rosenzweig, M. L. and R. H. MacArthur. 1963. Graphical representation and stability conditions of predator-prey interactions. American Naturalist 97: 209–223.

Volterra, V. 1926. Variazioni e fluttuazioni del numero d’individui in specie animali conviventi. Memorie della Reale Accademia Nazionale dei Lincei 2: 31–113.