Skip to content

Mutualism

Two species that grow logistically and feed each other, above their carrying capacities.

  • Name: mutualism
  • Field: biology/ecology
  • Formalism: ordinary differential equations
  • References: May 1981; Wright 1989; Argun et al. 2021

Factory

@register_model(
    field='biology/ecology',
    tags=('mutualism', 'ode', 'handling time'),
    references=('May 1981', 'Wright 1989', 'Argun et al. 2021'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def mutualism(
    first: float = 10.0,
    second: float = 10.0,
    *,
    r1: float = 1.0,
    r2: float = 1.0,
    K1: float = 50.0,
    K2: float = 50.0,
    a: float = 0.1,
    b: float = 0.5,
    handling: float = 0.5,
) -> System:
    """Return two species that grow logistically and feed each other, above their carrying capacities.

    With a positive handling time the benefit saturates and the populations settle above their capacities. Without it,
    strong mutualism grows without bound.

    Args:
        first:
            The initial size of the first population.

        second:
            The initial size of the second population.

        r1:
            The growth rate of the first population.

        r2:
            The growth rate of the second population.

        K1:
            The carrying capacity of the first population alone.

        K2:
            The carrying capacity of the second population alone.

        a:
            The effective search rate.

        b:
            The rate of encounters.

        handling:
            The handling time of the resource of the partner, 0 for type I mutualism.

    Returns:
        system:
            The model, with the populations `first` and `second`.

    Examples:
        >>> from skcomplex.models.biology.ecology import mutualism
        >>> [interaction.label for interaction in mutualism().interactions]
        ['growth_first', 'growth_second', 'mutualism']
    """
    return System(
        elements=[PopulationElement('first', n=first), PopulationElement('second', n=second)],
        interactions=[
            LogisticGrowthInteraction('growth_first', target='first', r=r1, K=K1),
            LogisticGrowthInteraction('growth_second', target='second', r=r2, K=K2),
            MutualismInteraction('mutualism', first='first', second='second', a=a, b=b, handling=handling),
        ],
    )

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 MutualismInteraction(Interaction):
    """Mutualism with a handling time, where encounters feed both species, `a b n1 n2 / (1 + a T n)` each.

    The benefit to each species saturates with the abundance of its partner, as the handling time `T` limits how much
    of the partner's resource it can use. A handling time of 0 gives the unbounded type I mutualism `a b n1 n2`.

    Args:
        label:
            The label of the interaction.

        first:
            The path of the first population.

        second:
            The path of the second population.

        a:
            The effective search rate.

        b:
            The rate of encounters.

        handling:
            The handling time `T` of the resource of the partner.

    Examples:
        >>> from skcomplex.models.biology.ecology import MutualismInteraction
        >>> MutualismInteraction('pollination', first='orchids', second='bees', handling=0.5).handling
        0.5
    """

    def __init__(
        self,
        label: str,
        *,
        first: Annotated[str, RoleSpec()],
        second: 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.1,
        handling: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    ) -> None:
        super().__init__(label)
        self.first = first
        self.second = second
        self.a = a
        self.b = b
        self.handling = handling

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the saturating benefits `a b n1 n2 / (1 + a T n2)` and `a b n1 n2 / (1 + a T n1)`.

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

        Returns:
            laws:
                The laws of the component.
        """
        first, second = s.first.n, s.second.n
        encounters = s.a * s.b * first * second
        return [
            RateLaw(first, encounters / (1 + s.a * s.handling * second)),
            RateLaw(second, encounters / (1 + s.a * s.handling * first)),
        ]

Equations

first.n: diff(n_first(t), t) == -r_growth_first*(n_first(t)/K_growth_first - 1)*n_first(t) + a_mutualism*b_mutualism*n_first(t)*n_second(t)/(a_mutualism*handling_mutualism*n_second(t) + 1)
second.n: diff(n_second(t), t) == -r_growth_second*(n_second(t)/K_growth_second - 1)*n_second(t) + a_mutualism*b_mutualism*n_first(t)*n_second(t)/(a_mutualism*handling_mutualism*n_first(t) + 1)