Random Lotka-Volterra with higher-order interactions
A random Lotka-Volterra community with third-order interactions.
- Name:
higher_order_random_lotka_volterra - Field:
biology/ecology - Formalism: ordinary differential equations
- References: Sidhom and Galla 2024
Factory
@register_model(
field='biology/ecology',
tags=('random community', 'Lotka-Volterra', 'higher-order interactions', 'coupling tensor', 'ode'),
references=('Sidhom and Galla 2024',),
defaults={'until': 1000.0, 'dt': 10.0},
)
def higher_order_random_lotka_volterra(
size: int = 200,
*,
mu: float = -4.0,
sigma: float = 1.0,
gamma: float = 0.0,
k: float = 1.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return a random Lotka-Volterra community with third-order interactions.
The abundances follow `dx_i/dt = x_i (k - x_i + sum_(j<k) a_ijk x_j x_k)` (Sidhom and Galla 2024, Eq. 2), where the
sum runs over the pairs of other species. The coefficients are Gaussian with mean `mu / binom(N - 1, 2)` and
variance `3 sigma^2 / (2 binom(N - 1, 2))`, and the three coefficients of a set of three species have the
correlation `gamma` (Eqs. 4 and 5). The defaults are the competitive community of the paper's Fig. 2, 200 species
with uncorrelated coefficients at `mu = -4` and `sigma = 1`, in its phase with a unique stable fixed point, where
the theory gives a fraction `0.928` of surviving species and a mean abundance `0.392`. The paper does not state the
distribution of the initial abundances, which are drawn uniformly from the interval from 0 to 1.
Args:
size:
The number of species.
mu:
The mean of the interactions, negative for competition.
sigma:
The heterogeneity of the interactions.
gamma:
The correlation of the three coefficients of a set of three species, from `-1/2` to 1, where 1 gives every
species of a set the same coefficient.
k:
The growth rate at low abundance, which is also the carrying capacity without interactions.
random_state:
The seed or the generator that draws the tensor and the initial abundances.
Returns:
system:
The model, with the abundances `species.n` of the members of the group `species`, the growth `growth`, and
the interactions `triplets`, whose tensor is the parameter `triplets.a`.
Examples:
>>> from skcomplex.models.biology.ecology import higher_order_random_lotka_volterra
>>> higher_order_random_lotka_volterra(5).interactions[1].a.shape
(5, 5, 5)
"""
rng = np.random.default_rng(random_state)
tensor = _draw_tensor(size, mu, sigma, gamma, rng)
species = GroupElement('species', element=PopulationElement, size=size, values={'n': rng.random(size)})
return System(
elements=[species],
interactions=[
LogisticGrowthInteraction('growth', target='species', r=k, K=k),
TripletInteraction('triplets', target='species', a=tensor, state='n'),
],
)
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 TripletInteraction(RateInteraction):
"""Interactions of every species with the pairs of other species, `dx_i/dt += x_i sum_(j<k) a_ijk x_j x_k`.
The coefficient `a_ijk` is the effect of the combination of the species `j` and `k` on the growth of the species
`i` (Sidhom and Galla 2024, Eq. 2). The tensor holds it for both orders of the pair, `a_ijk = a_ikj`, so the sum
over the sets `j < k` is half the contraction of the tensor with the abundances, and it holds 0 where indices
coincide.
Args:
label:
The label of the interaction.
target:
The path of the group of species.
a:
The coupling tensor of order 3, with one index per species, or a number for every entry.
state:
The name of the state that holds the abundance of a species.
Examples:
>>> from skcomplex.models.biology.ecology import TripletInteraction
>>> TripletInteraction('triplets', target='species', a=0.0, state='n').check().state
'n'
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
a: Annotated[object, ParameterSpec(shape=('target', 'target', 'target'))] = 0.0,
state: str = 'x',
) -> None:
super().__init__(label, target=target, state=state, rate=self._triplets)
self.a = a
def _triplets(self, namespace: SimpleNamespace) -> object:
"""Return the growth of a species from the pairs of other species."""
abundance = getattr(namespace.target, self.state)
return abundance * contraction(namespace.a, abundance, abundance) / 2
Equations
species.n[i]: diff(n_species(t), t) == -r_growth*(n_species(t)/K_growth - 1)*n_species(t) + 1/2*_reduced_0*n_species(t)