Food webs

A food web where basal species grow logistically and consumers live on Lotka-Volterra predation.

• Name: food_webs
• Field: biology/ecology
• Formalism: ordinary differential equations
• References: Argun et al. 2021; Pimm 1982

Factory

@register_model(
    field='biology/ecology',
    tags=('food web', 'predator-prey', 'ode'),
    references=('Argun et al. 2021', 'Pimm 1982'),
    defaults={'until': 100.0, 'dt': 0.5},
)
def food_webs(
    populations: Mapping[str, float] | None = None,
    *,
    growth: Mapping[str, tuple[float, float]] | None = None,
    mortality: Mapping[str, float] | None = None,
    links: Sequence[tuple[str, str, float, float]] | None = None,
) -> System:
    """Return a food web where basal species grow logistically and consumers live on Lotka-Volterra predation.

    The default web has frogs and mice at its base, snakes that eat frogs and mice, and hawks that eat snakes and
    mice. Its four species coexist at a stable equilibrium.

    Args:
        populations:
            The initial size of every species by name, or `None` for the default web.

        growth:
            The growth rate and the carrying capacity of every basal species, or `None` for the default web.

        mortality:
            The death rate of every consumer, or `None` for the default web.

        links:
            The links `(predator, prey, a, b)`, with the attack rate `a` and the conversion rate `b`, or `None` for the
            default web.

    Returns:
        system:
            The model, with one population per species.

    Examples:
        >>> from skcomplex.models.biology.ecology import food_webs
        >>> [element.label for element in food_webs().elements]
        ['frogs', 'mice', 'snakes', 'hawks']
    """
    populations = WEB_POPULATIONS if populations is None else populations
    growth = WEB_GROWTH if growth is None else growth
    mortality = WEB_MORTALITY if mortality is None else mortality
    links = WEB_LINKS if links is None else links
    interactions: list[BaseInteraction] = [
        LogisticGrowthInteraction(f'growth_{name}', target=name, r=r, K=K) for name, (r, K) in growth.items()
    ]
    interactions += [MortalityInteraction(f'death_{name}', target=name, d=d) for name, d in mortality.items()]
    interactions += [
        PredationInteraction(f'{predator}_eat_{prey}', predator=predator, prey=prey, a=a, b=b)
        for predator, prey, a, b in links
    ]
    return System(
        elements=[PopulationElement(name, n=n) for name, n in populations.items()],
        interactions=interactions,
    )

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 LogisticGrowthInteraction(GrowthInteraction):
    """Verhulst's growth limited by a carrying capacity, `dn/dt = r n (1 - n(t - tau) / K)`, without harvest.

    A delay `tau` above 0 gives Hutchinson's delayed logistic equation. The generalized logistic of Richards, with an
    exponent on the crowding, is `GrowthInteraction(form='logistic', p=...)`.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the element or group whose population grows.

        r:
            The per capita growth rate at low density.

        K:
            The carrying capacity.

        tau:
            The delay of the limiting factor, 0 for none.

        state:
            The name of the state of the target that holds the population.

    Examples:
        >>> from skcomplex.core.elements import CompartmentElement
        >>> from skcomplex.core.system import System
        >>> from skcomplex.engine.runtime import Simulation
        >>> from skcomplex.models.biology.common.populations import LogisticGrowthInteraction
        >>> system = System(
        ...     elements=[CompartmentElement('rabbits', n=1.0)],
        ...     interactions=[LogisticGrowthInteraction('growth', target='rabbits', r=1.0, K=50.0)],
        ... )
        >>> trajectory = Simulation(system, dt=0.1).check().prepare().run(until=30.0).trajectory_
        >>> round(float(trajectory.values('rabbits.n')[-1]), 3)
        50.0
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
        K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
        tau: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
        state: str = 'n',
    ) -> None:
        super().__init__(label, target=target, r=r, K=K, tau=tau, form='logistic', state=state)

    def _get_crowding(self, namespace: SimpleNamespace) -> Expression:
        """Return the crowding `n / K` of Verhulst's logistic.

        Args:
            namespace:
                The namespace of the symbols the rate can use, with the namespace of the target.

        Returns:
            crowding:
                The ratio of the population the limiting factor reads to the capacity.
        """
        crowding: Expression = self._get_limited(namespace) / self._get_capacity(namespace)
        return crowding

    def _get_rate(self, namespace: SimpleNamespace) -> Expression:
        """Return the logistic growth, with no harvest.

        Args:
            namespace:
                The namespace of the symbols the rate can use, with the namespace of the target.

        Returns:
            rate:
                The rate of the population.
        """
        return self._get_growth(namespace)

class MortalityInteraction(RateInteraction):
    """Death or loss at a constant per capita rate, `dn/dt = -d n`.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the element or group whose population dies.

        d:
            The per capita death rate.

        state:
            The name of the state of the target that holds the population.

    Examples:
        >>> from skcomplex.core.elements import CompartmentElement
        >>> from skcomplex.core.system import System
        >>> from skcomplex.engine.runtime import Simulation
        >>> from skcomplex.models.biology.common.populations import MortalityInteraction
        >>> system = System(
        ...     elements=[CompartmentElement('foxes', n=8.0)],
        ...     interactions=[MortalityInteraction('death', target='foxes', d=0.5)],
        ... )
        >>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
        >>> round(float(trajectory.values('foxes.n')[-1]), 3)
        2.943
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
        state: str = 'n',
    ) -> None:
        super().__init__(label, target=target, state=state)
        self.d = d

    def _get_rate(self, namespace: SimpleNamespace) -> Expression:
        """Return the loss `-d n`.

        Args:
            namespace:
                The namespace of the symbols the rate can use, with the namespace of the target.

        Returns:
            rate:
                The rate of the population.
        """
        loss: Expression = -namespace.d * getattr(namespace.target, self.state)
        return loss

class PredationInteraction(Interaction):
    """Lotka-Volterra predation, where encounters remove prey and feed predators.

    The prey loses `a * predator * prey` per unit time and the predator gains `b * predator * prey`.

    Args:
        label:
            The label of the interaction.

        predator:
            The path of the predator population.

        prey:
            The path of the prey population.

        a:
            The attack rate.

        b:
            The conversion rate of eaten prey into predators.

    Examples:
        >>> from skcomplex.models.biology.ecology import PredationInteraction
        >>> PredationInteraction('hunt', predator='foxes', prey='rabbits', a=0.4, b=0.1).prey
        'rabbits'
    """

    def __init__(
        self,
        label: str,
        *,
        predator: Annotated[str, RoleSpec()],
        prey: Annotated[str, RoleSpec()],
        a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
        b: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    ) -> None:
        super().__init__(label)
        self.predator = predator
        self.prey = prey
        self.a = a
        self.b = b

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the rates `-a * predator * prey` and `b * predator * prey`.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        encounters = s.predator.n * s.prey.n
        return [RateLaw(s.prey.n, -s.a * encounters), RateLaw(s.predator.n, s.b * encounters)]

Equations

frogs.n: diff(n_frogs(t), t) == -r_growth_frogs*(n_frogs(t)/K_growth_frogs - 1)*n_frogs(t) - a_snakes_eat_frogs*n_frogs(t)*n_snakes(t)
mice.n: diff(n_mice(t), t) == -r_growth_mice*(n_mice(t)/K_growth_mice - 1)*n_mice(t) - a_hawks_eat_mice*n_hawks(t)*n_mice(t) - a_snakes_eat_mice*n_mice(t)*n_snakes(t)
snakes.n: diff(n_snakes(t), t) == b_snakes_eat_frogs*n_frogs(t)*n_snakes(t) - a_hawks_eat_snakes*n_hawks(t)*n_snakes(t) + b_snakes_eat_mice*n_mice(t)*n_snakes(t) - d_death_snakes*n_snakes(t)
hawks.n: diff(n_hawks(t), t) == b_hawks_eat_mice*n_hawks(t)*n_mice(t) + b_hawks_eat_snakes*n_hawks(t)*n_snakes(t) - d_death_hawks*n_hawks(t)