Generalized stochastic predator-prey
Prey and predators born and killed one at a time, with crowding of the prey and saturating predation.
- Name:
stochastic_predator_prey - Field:
biology/ecology - Formalism: jump process
- References: Carmichael and Hadžikadić 2019; Rosenzweig 1971
Factory
@register_model(
field='biology/ecology',
tags=('predator-prey', 'jump process', 'gillespie', 'paradox of enrichment'),
references=('Carmichael and Hadžikadić 2019', 'Rosenzweig 1971'),
defaults={'until': 100.0, 'dt': 0.5},
)
def stochastic_predator_prey(
prey: float = 200.0,
predators: float = 200.0,
*,
r: float = 1.0,
K: float = 600.0,
a: float = 0.005,
h: float = 0.5,
e: float = 0.5,
m: float = 0.3,
stochastic: bool = True,
) -> System:
"""Return prey and predators born and killed one at a time, with crowding of the prey and saturating predation.
The deterministic limit has the coexistence point `x = m / (a (e - m h))`. It is stable while
`K < 2 x + 1 / (a h)`, and enriching the prey beyond that turns it into large cycles, which demographic noise drives
to extinction.
Args:
prey:
The initial number of prey.
predators:
The initial number of predators.
r:
The birth rate of the prey.
K:
The carrying capacity of the prey.
a:
The attack rate.
h:
The handling time.
e:
The number of predators born per prey eaten.
m:
The death rate of the predators.
stochastic:
Whether births, deaths, and kills are random events.
Returns:
system:
The model, with the populations `prey` and `predators`.
Examples:
>>> from skcomplex.models.biology.ecology import stochastic_predator_prey
>>> stochastic_predator_prey(K=1000.0).interactions[0].K
1000.0
"""
return System(
elements=[PopulationElement('prey', n=prey), PopulationElement('predators', n=predators)],
interactions=[
PredatorPreyEventsInteraction(
'events',
prey='prey',
predators='predators',
r=r,
K=K,
a=a,
h=h,
e=e,
m=m,
stochastic=stochastic,
),
],
)
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 PredatorPreyEventsInteraction(Interaction):
"""Births, crowding, saturating predation, and deaths of prey and predators as events of single individuals.
Prey are born at the rate `r x` and die of crowding at the rate `r x^2 / K`. Predators kill prey at the Holling
type II rate `a x y / (1 + a h x)`, are born at `e` times that rate, and die at the rate `m y`. In deterministic
mode the rates drive the populations as differential equations.
Args:
label:
The label of the interaction.
prey:
The path of the prey population `x`.
predators:
The path of the predator population `y`.
r:
The birth rate of the prey.
K:
The carrying capacity of the prey.
a:
The attack rate.
h:
The handling time.
e:
The number of predators born per prey eaten.
m:
The death rate of the predators.
stochastic:
Whether births, deaths, and kills are random events.
Examples:
>>> from skcomplex.models.biology.ecology import PredatorPreyEventsInteraction
>>> PredatorPreyEventsInteraction('events', prey='voles', predators='weasels', K=200.0).K
200.0
"""
def __init__(
self,
label: str,
*,
prey: Annotated[str, RoleSpec()],
predators: Annotated[str, RoleSpec()],
r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 200.0,
a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
h: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
e: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
m: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
stochastic: bool = True,
) -> None:
super().__init__(label)
self.prey = prey
self.predators = predators
self.r = r
self.K = K
self.a = a
self.h = h
self.e = e
self.m = m
self.stochastic = stochastic
def _get_change(self, member: SimpleNamespace, step: int, propensity: Expression) -> Sequence[BaseLaw]:
"""Return the laws that change the population of a member by a step at a propensity.
In stochastic mode the change is a jump of one individual, and otherwise a rate of the step times the
propensity.
"""
change: Interaction = (
JumpInteraction(self.label, target='member', changes={'n': step}, rate=lambda _: propensity)
if self.stochastic
else RateInteraction(self.label, target='member', state='n', rate=lambda _: step * propensity)
)
return change.check().prepare().run(SimpleNamespace(target=member)).laws_
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the five channels of births and deaths.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
prey, predators = s.prey, s.predators
x, y = prey.n, predators.n
kills = s.a * x * y / (1 + s.a * s.h * x)
return [
*self._get_change(prey, 1, s.r * x),
*self._get_change(prey, -1, s.r * x**2 / s.K),
*self._get_change(prey, -1, kills),
*self._get_change(predators, 1, s.e * kills),
*self._get_change(predators, -1, s.m * y),
]