Competitive Jain-Krishna model with signed couplings
The Jain-Krishna model with signed interactions, whose surviving species come to cooperate.
- Name:
competitive_jain_krishna - Field:
biology/evolution - Formalism: agent-based system
- References: Jain and Krishna 2002; Thurner et al. 2018
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('competitive_jain_krishna')
entry = describe('competitive_jain_krishna')
trajectory = create_simulation('competitive_jain_krishna').prepare().run().trajectory_
Factory
@register_system(
field='biology/evolution',
tags=('autocatalytic set', 'competition', 'cooperation', 'punctuated equilibrium', 'agents'),
references=('Jain and Krishna 2002', 'Thurner et al. 2018'),
defaults={'until': 400},
)
def competitive_jain_krishna(
n: int = 50,
*,
p: float = 0.05,
duration: float = 20.0,
seed: int | None = 1,
) -> System:
"""Return the Jain-Krishna model with signed interactions, whose surviving species come to cooperate.
Every link of the interaction matrix is present with the probability `p` and drawn from `[-1, 1]`, and every
species degrades at a rate drawn from `[-1, 0]`. The populations follow the linear catalytic dynamics with a floor
at zero, and at every step the least populated species is replaced by a new random one. Selection keeps the species
that catalyse each other, so once an autocatalytic set forms the links among the surviving species turn mostly
positive, the emergence of cooperation, and diversity grows until a keystone species is lost.
Args:
n:
The number of species.
p:
The probability of every link.
duration:
The time of the population dynamics between replacements.
seed:
The seed of the initial matrix.
Returns:
system:
The system, with the populations `x` and the rows `links` of the interaction matrix of the group `species`.
Examples:
>>> from skcomplex.systems.biology.evolution import competitive_jain_krishna
>>> competitive_jain_krishna(10).elements[0].prepare().size_
10
"""
rng = np.random.default_rng(seed)
matrix = rng.uniform(-1, 1, (n, n)) * (rng.random((n, n)) < p)
np.fill_diagonal(matrix, -rng.random(n))
species = GroupElement(
'species',
element=CompetingSpeciesElement,
size=n,
values={'x': np.full(n, 1 / n), 'links': matrix},
)
selection = CompetitiveSelectionInteraction('selection', species='species', p=p, duration=duration)
return System(elements=[species], interactions=[selection], time=DiscreteTime())
Components
Elements:
CompetingSpeciesElement, fromElement
Interactions:
CompetitiveSelectionInteraction, fromInteraction
class CompetingSpeciesElement(Element):
"""Species of the competitive Jain-Krishna model, a relative population and the signed effects of all species on it.
The links are the row of the interaction matrix of the species, `M_ij` for every species `j`, between -1 and 1, with
the degradation `M_ii` between -1 and 0 on the diagonal.
Args:
label:
The label of the species or of the group of species.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The relative population.
links:
The effects of every species on this one.
Examples:
>>> from skcomplex.systems.biology.evolution import CompetingSpeciesElement
>>> CompetingSpeciesElement('species', x=0.1).x
0.1
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
links: Annotated[object, StateSpec(vector=True)] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.links = links
@classmethod
def _get_view(cls) -> type[CompetingSpeciesView]:
"""Return the view of a species, `CompetingSpeciesView`."""
return CompetingSpeciesView
class CompetitiveSelectionInteraction(Interaction):
"""Population dynamics of signed catalytic interactions, with the least populated species replaced at every step.
Between replacements the relative populations follow Eq. (5.37) of Thurner et al. (2018), `dx_i/dt = f_i` with
`f_i = sum_j M_ij x_j - x_i phi`, where `phi = sum_jk M_jk x_k` keeps the total at one, except that a species at
zero stays there while `f_i < 0`. The rule integrates them for the time `duration` by Euler steps of `step`, then
replaces the least populated species, ties broken at random, by a new one whose links to and from every other
species are drawn from `[-1, 1]` with the probability `p`, whose degradation is drawn from `[-1, 0]`, and whose
population starts at `1 / N`.
Args:
label:
The label of the interaction.
species:
The path of the group of species.
p:
The probability of every link of a new species.
duration:
The time of the population dynamics between replacements.
step:
The Euler step of the population dynamics.
Examples:
>>> from skcomplex.systems.biology.evolution import CompetitiveSelectionInteraction
>>> CompetitiveSelectionInteraction('selection', species='species', p=0.1).p
0.1
"""
def __init__(
self,
label: str,
*,
species: Annotated[str, RoleSpec()],
p: float = 0.05,
duration: float = 20.0,
step: float = 0.05,
) -> None:
super().__init__(label)
self.species = species
self.p = p
self.duration = duration
self.step = step
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the selection rule.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return [RuleLaw(self.select)]
def select(self, system: SystemView) -> None:
"""Run the population dynamics, then replace the least populated species.
Args:
system:
The system.
"""
species = system.agents(self.species, CompetingSpeciesView)
matrix = np.array([np.asarray(item.links, dtype=float) for item in species])
populations = self._relax(matrix, np.array([item.x for item in species]))
size = len(species)
weakest = np.flatnonzero(populations <= populations.min() + 1e-12)
k = int(system.rng.choice(weakest))
drawn = system.rng.uniform(-1, 1, (2, size)) * (system.rng.random((2, size)) < self.p)
matrix[k, :], matrix[:, k] = drawn[0], drawn[1]
matrix[k, k] = -system.rng.random()
populations[k] = 1 / size
populations /= populations.sum()
for index, item in enumerate(species):
item.x = float(populations[index])
item.links = matrix[index]
def _relax(self, matrix: Array, populations: Array) -> Array:
"""Return the populations after the dynamics of Eq.
(5.37) for the time `duration`.
"""
x = populations.copy()
for _ in range(round(self.duration / self.step)):
drive = matrix @ x
f = drive - x * drive.sum()
f[(x <= 0) & (f < 0)] = 0.0
x = np.clip(x + self.step * f, 0.0, None)
total = x.sum()
x = x / total if total > 0 else np.full_like(x, 1 / x.size)
return x