Custom Model — Lotka-Volterra¶

SeapoPym is not limited to marine ecosystem models. Any system of ODEs on a grid can be expressed as a Blueprint. This example builds the classic Lotka-Volterra predator-prey model from scratch.

$$\frac{dN}{dt} = \alpha N - \beta N P \qquad \frac{dP}{dt} = \delta \beta N P - \gamma P$$

where $\alpha$ is the prey growth rate, $\beta$ the predation rate, $\delta$ the conversion efficiency, and $\gamma$ the predator mortality rate.

In [1]:

import time

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from seapopym.blueprint import Blueprint, Config, functional
from seapopym.compiler import compile_model
from seapopym.engine import simulate

PALETTE = ["#1B4965", "#62B6CB", "#E8833A", "#5FA8D3"]
plt.rcParams.update({"figure.facecolor": "white", "axes.facecolor": "white", "axes.grid": True, "grid.alpha": 0.3})

import time import matplotlib.pyplot as plt import numpy as np import xarray as xr from seapopym.blueprint import Blueprint, Config, functional from seapopym.compiler import compile_model from seapopym.engine import simulate PALETTE = ["#1B4965", "#62B6CB", "#E8833A", "#5FA8D3"] plt.rcParams.update({"figure.facecolor": "white", "axes.facecolor": "white", "axes.grid": True, "grid.alpha": 0.3})

1. Define the Physics¶

Each term of the ODE system is a pure function decorated with @functional. The decorator registers the function in SeapoPym's registry and declares its unit signature for compile-time validation.

In [2]:

@functional(
    name="lv:prey_growth",
    units={"N": "dimensionless", "alpha": "1/s", "return": "1/s"},
)
def prey_growth(N, alpha):
    """Prey exponential growth: +alpha * N."""
    return alpha * N


@functional(
    name="lv:predation",
    units={
        "N": "dimensionless",
        "P": "dimensionless",
        "beta": "1/s",
        "delta": "dimensionless",
        "prey_loss": "1/s",
        "predator_gain": "1/s",
    },
    outputs=("prey_loss", "predator_gain"),
)
def predation(N, P, beta, delta):
    """Predation interaction: prey loses beta*N*P, predator gains delta*beta*N*P."""
    interaction = beta * N * P
    return -interaction, delta * interaction


@functional(
    name="lv:predator_death",
    units={"P": "dimensionless", "gamma": "1/s", "return": "1/s"},
)
def predator_death(P, gamma):
    """Predator natural mortality: -gamma * P."""
    return -gamma * P

@functional( name="lv:prey_growth", units={"N": "dimensionless", "alpha": "1/s", "return": "1/s"}, ) def prey_growth(N, alpha): """Prey exponential growth: +alpha * N.""" return alpha * N @functional( name="lv:predation", units={ "N": "dimensionless", "P": "dimensionless", "beta": "1/s", "delta": "dimensionless", "prey_loss": "1/s", "predator_gain": "1/s", }, outputs=("prey_loss", "predator_gain"), ) def predation(N, P, beta, delta): """Predation interaction: prey loses beta*N*P, predator gains delta*beta*N*P.""" interaction = beta * N * P return -interaction, delta * interaction @functional( name="lv:predator_death", units={"P": "dimensionless", "gamma": "1/s", "return": "1/s"}, ) def predator_death(P, gamma): """Predator natural mortality: -gamma * P.""" return -gamma * P

2. Build the Blueprint¶

The Blueprint declares state variables, parameters, and the process DAG — which functions to call, with which inputs, producing which outputs. Tendencies define how derived fluxes feed back into state variables via Euler integration.

In [3]:

blueprint = Blueprint.from_dict(
    {
        "id": "lotka-volterra",
        "version": "1.0",
        "declarations": {
            "state": {
                "prey": {"units": "dimensionless", "dims": ["Y", "X"]},
                "predator": {"units": "dimensionless", "dims": ["Y", "X"]},
            },
            "parameters": {
                "alpha": {"units": "1/s", "description": "Prey growth rate"},
                "beta": {"units": "1/s", "description": "Predation rate"},
                "delta": {"units": "dimensionless", "description": "Conversion efficiency"},
                "gamma": {"units": "1/s", "description": "Predator mortality rate"},
            },
            "forcings": {},
        },
        "process": [
            {
                "func": "lv:prey_growth",
                "inputs": {"N": "state.prey", "alpha": "parameters.alpha"},
                "outputs": {"return": "derived.prey_growth"},
            },
            {
                "func": "lv:predation",
                "inputs": {
                    "N": "state.prey",
                    "P": "state.predator",
                    "beta": "parameters.beta",
                    "delta": "parameters.delta",
                },
                "outputs": {
                    "prey_loss": "derived.prey_loss",
                    "predator_gain": "derived.predator_gain",
                },
            },
            {
                "func": "lv:predator_death",
                "inputs": {"P": "state.predator", "gamma": "parameters.gamma"},
                "outputs": {"return": "derived.predator_death"},
            },
        ],
        "tendencies": {
            "prey": [
                {"source": "derived.prey_growth"},
                {"source": "derived.prey_loss"},
            ],
            "predator": [
                {"source": "derived.predator_gain"},
                {"source": "derived.predator_death"},
            ],
        },
    }
)

blueprint = Blueprint.from_dict( { "id": "lotka-volterra", "version": "1.0", "declarations": { "state": { "prey": {"units": "dimensionless", "dims": ["Y", "X"]}, "predator": {"units": "dimensionless", "dims": ["Y", "X"]}, }, "parameters": { "alpha": {"units": "1/s", "description": "Prey growth rate"}, "beta": {"units": "1/s", "description": "Predation rate"}, "delta": {"units": "dimensionless", "description": "Conversion efficiency"}, "gamma": {"units": "1/s", "description": "Predator mortality rate"}, }, "forcings": {}, }, "process": [ { "func": "lv:prey_growth", "inputs": {"N": "state.prey", "alpha": "parameters.alpha"}, "outputs": {"return": "derived.prey_growth"}, }, { "func": "lv:predation", "inputs": { "N": "state.prey", "P": "state.predator", "beta": "parameters.beta", "delta": "parameters.delta", }, "outputs": { "prey_loss": "derived.prey_loss", "predator_gain": "derived.predator_gain", }, }, { "func": "lv:predator_death", "inputs": {"P": "state.predator", "gamma": "parameters.gamma"}, "outputs": {"return": "derived.predator_death"}, }, ], "tendencies": { "prey": [ {"source": "derived.prey_growth"}, {"source": "derived.prey_loss"}, ], "predator": [ {"source": "derived.predator_gain"}, {"source": "derived.predator_death"}, ], }, } )

Process DAG¶

In [4]:

blueprint.to_graphviz()

blueprint.to_graphviz()

Out[4]:

3. Configure & Compile¶

The Config provides concrete parameter values, initial conditions, and execution settings. The compiler validates units, infers shapes, and produces a CompiledModel ready for JAX execution.

We use a single grid cell (1×1) — Lotka-Volterra is a purely temporal model, so no spatial structure is needed.

In [5]:

DAY = 86400.0  # seconds per day

config = Config.from_dict(
    {
        "parameters": {
            "alpha": xr.DataArray(0.04 / DAY),  # prey growth: ~4% per day
            "beta": xr.DataArray(0.005 / DAY),  # predation rate
            "delta": xr.DataArray(0.5),  # conversion efficiency: 50%
            "gamma": xr.DataArray(0.1 / DAY),  # predator mortality: ~10% per day
        },
        "forcings": {},
        "initial_state": {
            "prey": xr.DataArray(np.array([[42.0]]), dims=["Y", "X"]),
            "predator": xr.DataArray(np.array([[7.0]]), dims=["Y", "X"]),
        },
        "execution": {
            "time_start": "2000-01-01",
            "time_end": "2002-12-31",
            "dt": "1d",
        },
    }
)

model = compile_model(blueprint, config)
print(f"Compiled: {model.time_grid.n_timesteps} timesteps, dt = 1 day")

DAY = 86400.0 # seconds per day config = Config.from_dict( { "parameters": { "alpha": xr.DataArray(0.04 / DAY), # prey growth: ~4% per day "beta": xr.DataArray(0.005 / DAY), # predation rate "delta": xr.DataArray(0.5), # conversion efficiency: 50% "gamma": xr.DataArray(0.1 / DAY), # predator mortality: ~10% per day }, "forcings": {}, "initial_state": { "prey": xr.DataArray(np.array([[42.0]]), dims=["Y", "X"]), "predator": xr.DataArray(np.array([[7.0]]), dims=["Y", "X"]), }, "execution": { "time_start": "2000-01-01", "time_end": "2002-12-31", "dt": "1d", }, } ) model = compile_model(blueprint, config) print(f"Compiled: {model.time_grid.n_timesteps} timesteps, dt = 1 day")

Compiled: 1095 timesteps, dt = 1 day

4. Simulate¶

In [6]:

t0 = time.perf_counter()
state, outputs = simulate(model)
elapsed = time.perf_counter() - t0

prey = outputs["prey"].values[:, 0, 0]  # (T,) — single grid cell
predator = outputs["predator"].values[:, 0, 0]
time_axis = np.arange(len(prey))

print(f"Simulated {len(prey)} timesteps in {elapsed:.2f}s")
print(f"Prey     — min: {prey.min():.1f}, max: {prey.max():.1f}")
print(f"Predator — min: {predator.min():.1f}, max: {predator.max():.1f}")

t0 = time.perf_counter() state, outputs = simulate(model) elapsed = time.perf_counter() - t0 prey = outputs["prey"].values[:, 0, 0] # (T,) — single grid cell predator = outputs["predator"].values[:, 0, 0] time_axis = np.arange(len(prey)) print(f"Simulated {len(prey)} timesteps in {elapsed:.2f}s") print(f"Prey — min: {prey.min():.1f}, max: {prey.max():.1f}") print(f"Predator — min: {predator.min():.1f}, max: {predator.max():.1f}")

Simulated 1095 timesteps in 0.03s
Prey     — min: 15.4, max: 77.1
Predator — min: 2.1, max: 22.2

5. Results¶

Population Dynamics¶

Phase Portrait¶

The phase portrait reveals the characteristic closed orbits of the Lotka-Volterra system. The equilibrium point is at $N^* = \gamma / (\delta \beta)$, $P^* = \alpha / \beta$.

Summary¶

This example demonstrated SeapoPym's genericity: any ODE system can be expressed as a Blueprint, compiled, and simulated — with no code changes to the framework.

Step	Code
Define physics	3 @functional functions (~20 lines)
Declare model	1 Blueprint.from_dict() call
Configure	1 Config.from_dict() call
Compile & run	compile_model() → simulate()

The same pipeline applies to marine ecosystems (LMTL), biogeochemical models, or any dynamical system on a grid.

Next: Gradient-Based Optimization