Predator–prey cycles in a chemostat¶

C1 Research Computing and Software Development · Lecture 1 demonstration

Rotifers (Brachionus calyciflorus, the predator) feeding on green algae (Monoraphidium minutum, the prey) in a continuously diluted culture vessel (a chemostat). Experiment C9, a run with a pulsed nutrient supply (the medium alternates between 160 and 0 µmol N/l every four days), of

B. Blasius, L. Rudolf, G. Weithoff, U. Gaedke, G. F. Fussmann, Long-term cyclic persistence in an experimental predator–prey system, Nature 577, 226–230 (2020), www.nature.com/articles/s41586-019-1857-0.

Data: figshare, doi:10.6084/m9.figshare.10045976 (CC BY 4.0); the file C9.csv is figshare file 18105215.

In [1]:

import sys

# Which interpreter runs this notebook, and is it a virtual environment?
print(sys.executable)
print("inside a venv:", sys.prefix != sys.base_prefix)

import sys # Which interpreter runs this notebook, and is it a virtual environment? print(sys.executable) print("inside a venv:", sys.prefix != sys.base_prefix)

/Users/boris/GitHub/c1-research-computing-2026/.venv/bin/python3.12
inside a venv: True

In [2]:

import urllib.request
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

DATA = Path("C9.csv")
URL = "https://ndownloader.figshare.com/files/18105215"  # C9.csv, doi:10.6084/m9.figshare.10045976
if not DATA.exists():
    urllib.request.urlretrieve(URL, DATA)

# columns: time (days), algae (10^6 cells/ml), rotifers (animals/ml), ...
data = np.genfromtxt(DATA, delimiter=",", skip_header=1)
t, algae, rotifers = data[:, 0], data[:, 1], data[:, 2]
print(f"{len(t)} samples over {t[-1]:.0f} days")

import urllib.request from pathlib import Path import matplotlib.pyplot as plt import numpy as np DATA = Path("C9.csv") URL = "https://ndownloader.figshare.com/files/18105215" # C9.csv, doi:10.6084/m9.figshare.10045976 if not DATA.exists(): urllib.request.urlretrieve(URL, DATA) # columns: time (days), algae (10^6 cells/ml), rotifers (animals/ml), ... data = np.genfromtxt(DATA, delimiter=",", skip_header=1) t, algae, rotifers = data[:, 0], data[:, 1], data[:, 2] print(f"{len(t)} samples over {t[-1]:.0f} days")

141 samples over 140 days

Time series¶

In [3]:

fig, ax = plt.subplots(figsize=(10, 3.5))
ax.plot(t, algae, "o-", ms=2, lw=0.8, color="C0", label="algae (prey)")
ax.set_xlabel("time (days)")
ax.set_ylabel("algae ($10^6$ cells/ml)", color="C0")
ax2 = ax.twinx()
ax2.plot(t, rotifers, "o-", ms=2, lw=0.8, color="C3", label="rotifers (predator)")
ax2.set_ylabel("rotifers (animals/ml)", color="C3")
ax.set_title("Experiment C9")
plt.show()

fig, ax = plt.subplots(figsize=(10, 3.5)) ax.plot(t, algae, "o-", ms=2, lw=0.8, color="C0", label="algae (prey)") ax.set_xlabel("time (days)") ax.set_ylabel("algae ($10^6$ cells/ml)", color="C0") ax2 = ax.twinx() ax2.plot(t, rotifers, "o-", ms=2, lw=0.8, color="C3", label="rotifers (predator)") ax2.set_ylabel("rotifers (animals/ml)", color="C3") ax.set_title("Experiment C9") plt.show()

Phase plane¶

A predator–prey cycle traces a loop in the plane of the two populations: prey grow, predators follow, prey collapse, predators starve.

In [4]:

fig, ax = plt.subplots(figsize=(5, 4))
ax.plot(algae, rotifers, lw=0.7)
ax.set_xlabel("algae ($10^6$ cells/ml)")
ax.set_ylabel("rotifers (animals/ml)")
plt.show()

fig, ax = plt.subplots(figsize=(5, 4)) ax.plot(algae, rotifers, lw=0.7) ax.set_xlabel("algae ($10^6$ cells/ml)") ax.set_ylabel("rotifers (animals/ml)") plt.show()

Period and phase lag¶

The samples are almost daily; interpolating onto a daily grid allows a Fourier transform (the dominant period; in this run it is the period of the nutrient pulses) and a cross-correlation (the delay of the predator behind the prey).

In [5]:

days = np.arange(t[0], t[-1])
a = np.interp(days, t, algae) - algae.mean()
r = np.interp(days, t, rotifers) - rotifers.mean()

freq = np.fft.rfftfreq(len(days), d=1.0)
power = np.abs(np.fft.rfft(a)) ** 2
period = 1 / freq[1:][np.argmax(power[1:])]

lags = np.arange(-6, 7)
corr = [np.corrcoef(a[max(0, -k):len(a) - max(0, k)], r[max(0, k):len(r) - max(0, -k)])[0, 1] for k in lags]
lag = lags[int(np.argmax(corr))]

print(f"dominant period: {period:.1f} days (nutrient pulses every 4 days)")
print(f"rotifers lag algae by about {lag} days (a quarter period: {period / 4:.1f} days)")

days = np.arange(t[0], t[-1]) a = np.interp(days, t, algae) - algae.mean() r = np.interp(days, t, rotifers) - rotifers.mean() freq = np.fft.rfftfreq(len(days), d=1.0) power = np.abs(np.fft.rfft(a)) ** 2 period = 1 / freq[1:][np.argmax(power[1:])] lags = np.arange(-6, 7) corr = [np.corrcoef(a[max(0, -k):len(a) - max(0, k)], r[max(0, k):len(r) - max(0, -k)])[0, 1] for k in lags] lag = lags[int(np.argmax(corr))] print(f"dominant period: {period:.1f} days (nutrient pulses every 4 days)") print(f"rotifers lag algae by about {lag} days (a quarter period: {period / 4:.1f} days)")

dominant period: 7.8 days (nutrient pulses every 4 days)
rotifers lag algae by about 2 days (a quarter period: 2.0 days)

Next steps¶

The cycles call for a model: prey growth, predation, predator mortality and the nutrient supply of the chemostat. Building, testing and scaling that model is the running project of the course, the attractor package.