Rosenzweig-MacArthur and rate tipping
Fast prey and slow predators whose carrying capacity shrinks at a rate, which tips them when it is fast.
- Name:
rosenzweig_macarthur - Field:
biology/ecology - Formalism: ordinary differential equations
- References: Rosenzweig and MacArthur 1963; Datseris and Parlitz 2022
Factory
@register_model(
field='biology/ecology',
tags=('predator-prey', 'ode', 'tipping', 'rate-induced tipping', 'slow-fast'),
references=('Rosenzweig and MacArthur 1963', 'Datseris and Parlitz 2022'),
defaults={'until': 150.0, 'dt': 0.5, 'solver': ScipyODESolver(rtol=1e-8, atol=1e-300)},
)
def rosenzweig_macarthur(
prey: float | None = None,
predators: float | None = None,
*,
kappa: float = 0.01,
eta: float = 0.8,
phi: float | None = None,
phi_end: float | None = None,
rate: float = 0.002,
) -> System:
"""Return fast prey and slow predators whose carrying capacity shrinks at a rate, which tips them when it is fast.
By default the crowding rises from `eta (1 - eta) / (1 + eta)` to `1 - eta`, the range where the coexistence point
at `x = 1 / (1 - eta)` is stable, and the run starts there. A slow rise lets the populations track the point, and a
rise faster than about 0.0004 makes the prey collapse for a while before they return.
In a collapse the prey fall to about `1e-30`, far below any absolute tolerance, and an error of either sign there
decides whether they return or, below zero, run into the pole of the feeding at `x = -1 / eta`. The registered
solver therefore controls the relative error only, with `ScipyODESolver(rtol=1e-8, atol=1e-300)`.
Args:
prey:
The initial prey `x`, or `None` for the coexistence point.
predators:
The initial predators `y`, or `None` for the coexistence point.
kappa:
The ratio of the time scales of the prey and the predators.
eta:
The handling time of the predators.
phi:
The initial crowding, or `None` for `eta (1 - eta) / (1 + eta)`.
phi_end:
The final crowding, or `None` for `1 - eta`.
rate:
The rate at which the crowding rises.
Returns:
system:
The model, with the populations `prey` and `predators`.
Examples:
>>> from skcomplex.models.biology.ecology import rosenzweig_macarthur
>>> round(rosenzweig_macarthur().elements[0].n, 6)
5.0
"""
phi = eta * (1 - eta) / (1 + eta) if phi is None else phi
phi_end = 1 - eta if phi_end is None else phi_end
prey = 1 / (1 - eta) if prey is None else prey
predators = (1 - phi * prey) * (1 + eta * prey) if predators is None else predators
return System(
elements=[PopulationElement('prey', n=prey), PopulationElement('predators', n=predators)],
interactions=[
RosenzweigMacArthurInteraction(
'predation',
prey='prey',
predators='predators',
kappa=kappa,
eta=eta,
phi=phi,
phi_end=phi_end,
rate=rate,
),
],
)
Components
class PopulationElement(CompartmentElement):
"""Population of one species, a compartment whose state `n` holds its size.
Args:
label:
The label of the population.
visual:
How the element is drawn, or `None` for `default_visual_`.
n:
The initial size.
Examples:
>>> from skcomplex.models.biology.common.populations import PopulationElement
>>> PopulationElement('rabbits', n=40).n
40
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
n: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual, n=n)
class RosenzweigMacArthurInteraction(Interaction):
"""Rosenzweig-MacArthur predation of fast prey by slow predators, whose prey crowding grows at a constant rate.
In units of the predators' time scale, `kappa dx/dt = x (1 - phi x) - x y / (1 + eta x)` and
`dy/dt = x y / (1 + eta x) - y`. The crowding `phi`, the inverse of the prey's carrying capacity, rises at the rate
`rate` from `phi` to `phi_end` and stays there.
Args:
label:
The label of the interaction.
prey:
The path of the prey population `x`.
predators:
The path of the predator population `y`.
kappa:
The ratio of the time scales of the prey and the predators.
eta:
The handling time of the predators.
phi:
The initial crowding of the prey.
phi_end:
The final crowding of the prey.
rate:
The rate at which the crowding rises.
Examples:
>>> from skcomplex.models.biology.ecology import RosenzweigMacArthurInteraction
>>> RosenzweigMacArthurInteraction('predation', prey='prey', predators='predators', rate=0.01).rate
0.01
"""
def __init__(
self,
label: str,
*,
prey: Annotated[str, RoleSpec()],
predators: Annotated[str, RoleSpec()],
kappa: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.8,
phi: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
phi_end: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
rate: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label)
self.prey = prey
self.predators = predators
self.kappa = kappa
self.eta = eta
self.phi = phi
self.phi_end = phi_end
self.rate = rate
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the rates of the prey and the predators under the ramp of the crowding.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
x, y = s.prey.n, s.predators.n
crowding = clip(s.phi + s.rate * s.t, high=s.phi_end)
feeding = x * y / (1 + s.eta * x)
return [RateLaw(x, (x * (1 - crowding * x) - feeding) / s.kappa), RateLaw(y, feeding - y)]
Equations
prey.n: diff(n_prey(t), t) == -((n_prey(t)*where(rate_predation*t + phi_predation > phi_end_predation, phi_end_predation, rate_predation*t + phi_predation) - 1)*n_prey(t) + n_predators(t)*n_prey(t)/(eta_predation*n_prey(t) + 1))/kappa_predation
predators.n: diff(n_predators(t), t) == n_predators(t)*n_prey(t)/(eta_predation*n_prey(t) + 1) - n_predators(t)