Predator-Prey

from Red Blob Games
1 Jan 2021

During the covid pandemic I’ve been learning about dynamic systems in biology from a textbook, Modeling Life[1]. To go deeper into the math I’ve been watching Robert Ghrist’s lectures[2]. To go deeper into the biology I’ve been watching the MIT Systems Biology class[3]. One of the early 2-variable examples is the predator-prey model (“Lotka Volterra”). I’m not doing all the homeworks but instead I wanted to pick a handful of projects to dive into, and I decided this would be the first.

These were originally rough notes that I polished over time as I learned the material. See my blog post[4] to see how earlier versions of this page looked.

 1  Lotka-Volterra#

The main idea is that there’s a prey animal \(x\) 🐇 and a predator animal \(y\) 🦊. When the prey population increases, it provides more food for predators to multiply.

Starting prey: 🐇 predators: 🦊
α β γ δ

The Lotka-Volterra rules are:

  1. Birth: the prey increases on its own: \(+α \cdot x\)
  2. Malnourishment: the predators decrease on their own: \(-γ \cdot y\)
  3. Predation: the prey decreases due to predators: \(-β \cdot x \cdot y\)
  4. Predation: the predators increase due to prey: \(+δ \cdot x \cdot y\)

These rules aren’t a particularly good model but they’re simple. This model is widely studied and seems like a good starting point for me to learn this type of system. It also has some surprising behavior that I want to explore in the simulations. The rules give us these differential equations which could be implemented in code:

\[\begin{aligned} \frac{dx}{dt} & = \alpha x - \beta x y \\ \frac{dy}{dt} & = \delta x y - \gamma y \end{aligned}\]

dx_dt = alpha*x - beta*x*y
dy_dt = delta*x*y - gamma*y

x += dx_dt * dt
y += dy_dt * dt

The chart above shows the populations over time. You can see they’re repeating. There’s a different type of chart, the phase diagram[5], which plots the two populations together but without time:


Starting prey: 🐇 predators: 🦊
α β γ δ

This shows how the two variables related, but not how fast they’re moving. Both of these 2d diagrams are views of a 3d diagram (time, predator, prey), which might be fun to make but I didn’t get around to doing so. It’d look like this:

Sketch of a spiral diagram showing predator, prey, and time all in one 3d chart

Viewed from the side or top, it’s a plot of predators or prey over time. But viewed from the end, it’s a phase plot. It’d be a cool animation but it’s a bit of work so I didn’t make it. I did read a little bit about how to draw vector (this paper[6] and this blog post)[7] and I experimented a little bit[8].

What else can I explore?

 1.1 Period

The period of the cycle is how long it takes before it repeats. The Wikipedia page[9] says the period is \(2π / \sqrt{\color{#00a}α \color{#a44}γ}\). For my default values of \(\textcolor{#00a}{α} = \frac{2}{3}, \textcolor{#a44}{γ} = 1\), that means the period is around 7.695 seconds. But my numerical integration doesn’t produce that! It’s not even close. It turns out that the formula only applies at the equilibrium point, and I didn’t find a formula at any other point.

Would it be interesting to display the period for every set of parameters?

{ need to try it and see }

{My implementation is fast enough that I can calculate period, min, max, mean and then do something with it, like plot them for many parameter values at a time. “Ladder of abstraction”. Since there are cycles, I need to map each point to a cycle, and then I can plot the cycle vs the period/min/max/mean. I could use the max as the cycle number, and then I’d be plotting period/min/mean vs max? Alternatively I could use a 2d input with starting }

 1.2 Equilibrium

The equilibrium point is when the population stays the same.

When I compare this to the equations

\[\begin{aligned} \frac{dx}{dt} & = \textcolor{#00a}{\alpha} x - \textcolor{#00a}{\beta} x y \\ \frac{dy}{dt} & = \textcolor{#a44}{\delta} x y - \textcolor{#a44}{\gamma} y \end{aligned}\]

I was a bit surprised that the equilibrium for predators \(y\) depends only on the parameters in the prey equation \(\frac{dx}{dt}\), and the equilibrium for prey \(x\) depends only on the parameters in the predator equation \(\frac{dy}{dt}\).

I would expect that the equilbrium for prey would depend on the prey birth rate, \(\color{#00a} α\), but it doesn’t! And the equilibrium for predators does not depend on the death rate of predators. This both means the model can be surprising (good for math!), and that the model likely does not match reality (bad for biology!).

There’s also an equilibrium point when both predator and prey are at 0, but that’s not so interesting.

 1.3 Mean, min, max

The mean should be the average value over one cycle. I don’t see this mentioned on the Wikipedia page or the textbook.

But I discovered by calculating and plotting the mean that it’s equal to the equilibrium!

Similarly, I can calculate the min and max numerically but I don’t see an analytical way to calculate this.

 1.4 Intervention

One of the more fascinating things I read in the Modeling Life textbook is that there’s counterintuitive behavior when intervening in the population:

Modeling Life textbook, describing that reducing predator populations can sometimes cause the predator population to increase
page 7, Modeling Life, with my highlights+annotation

Reducing the predator population can cause the predator population to go up higher than it had been. The textbook says that it depends on which phase of the cycle you’re in when you reduce the population. I wanted to see for myself, so I wrote the simulation:

Starting prey: 🐇 predators: 🦊
α β γ δ

Try this: click/drag in the middle and set the predator (red) population. Try it near the bottom and the top of the cycle. See what happens to the predators and also the prey (blue). Fascinating!

Reducing the predator population slightly can either increase or decrease the maximum predator population, but reducing it by a large amount almost always increases the amount of predators! Is there a way to know when the population will decrease vs increase? At the time I wrote this page, I hadn’t gotten to that part of the book yet, but the answer is yes:

Modeling Life textbook, describing the use of phase plots to predict the effect of intervention
page 60, Modeling Life

 1.5 Implementation

I used numerical integration and placed the results into Javascript typed arrays. This calculation turned out to be just as fast in the browser (Javascript, no type declarations) as in C++ or Rust (native code, type declarations, optimizing compiler). That was a bit of a surprise to me.

I had thought I should put the calculation in Web Workers so that it doesn’t interrupt the main flow of the page but I never got around to it.

I wrote some notes about accuracy of numerical integration and about the performance of various implementations I tried.

 2  Other models#

other stuff?

 3  More reading#

Email me , or tweet @redblobgames, or comment: