Lotka-Volterra

Prey and predators whose populations oscillate through predation.

• Name: lotka_volterra
• Field: biology/ecology
• Formalism: ordinary differential equations
• References: Lotka 1925; Volterra 1926

Factory

@register_model(
    field='biology/ecology',
    tags=('predator-prey', 'ode', 'oscillation'),
    references=('Lotka 1925', 'Volterra 1926'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def lotka_volterra(
    prey: float = 40.0,
    predators: float = 9.0,
    *,
    r: float = 1.1,
    a: float = 0.4,
    b: float = 0.1,
    d: float = 0.4,
    K: float | None = None,
) -> System:
    """Return prey and predators whose populations oscillate through predation.

    Args:
        prey:
            The initial number of prey.

        predators:
            The initial number of predators.

        r:
            The growth rate of the prey.

        a:
            The rate at which predation removes prey.

        b:
            The rate at which predation feeds predators.

        d:
            The mortality rate of the predators.

        K:
            The carrying capacity of the prey, or `None` for exponential growth of the prey.

    Returns:
        system:
            The model, with the populations `prey` and `predators`.

    Examples:
        >>> from skcomplex.models.biology.ecology import lotka_volterra
        >>> [element.label for element in lotka_volterra().elements]
        ['prey', 'predators']
    """
    growth: BaseInteraction = (
        ExponentialGrowthInteraction('growth', target='prey', r=r)
        if K is None
        else LogisticGrowthInteraction('growth', target='prey', r=r, K=K)
    )
    return System(
        elements=[PopulationElement('prey', n=prey), PopulationElement('predators', n=predators)],
        interactions=[
            growth,
            PredationInteraction('predation', predator='predators', prey='prey', a=a, b=b),
            MortalityInteraction('death', target='predators', d=d),
        ],
    )

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 ExponentialGrowthInteraction(GrowthInteraction):
    """Growth at a constant per capita rate, `dn/dt = r n`, the exponential form of a growth without harvest.

    Args:
        label:
            The label of the interaction.

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

        r:
            The per capita growth rate, negative for decay.

        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 ExponentialGrowthInteraction
        >>> system = System(
        ...     elements=[CompartmentElement('rabbits', n=1.0)],
        ...     interactions=[ExponentialGrowthInteraction('growth', target='rabbits', r=0.5)],
        ... )
        >>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
        >>> round(float(trajectory.values('rabbits.n')[-1]), 3)
        2.718
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        r: Annotated[float, ParameterSpec()] = 1.0,
        state: str = 'n',
    ) -> None:
        super().__init__(label, target=target, r=r, form='exponential', state=state)

    def _get_rate(self, namespace: SimpleNamespace) -> Expression:
        """Return the growth `r n`, 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 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)]

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

Equations

prey.n: diff(n_prey(t), t) == -a_predation*n_predators(t)*n_prey(t) + r_growth*n_prey(t)
predators.n: diff(n_predators(t), t) == b_predation*n_predators(t)*n_prey(t) - d_death*n_predators(t)