Evolution

Evolution: evolutionary games, quasispecies, hypercycles, finite populations, landscapes, and coevolution.

AdaptiveClimbingInteraction(label, *, target, mutation=0.1, length=20, seed=0)

Bases: EvolutionInteraction

One generation of a population that climbs a random fitness landscape as a whole, under strong selection.

All individuals share one genome. Each mutates one random locus with probability u, to a genome of random fitness in [0, 1], with probability mutation. If a mutant is fitter, the whole population takes the fittest mutant. If every individual mutated to a less fit genome, the population escapes to the fittest of them. Otherwise it stays. It is the best-mutant replacement of EvolutionInteraction on genomes stored as the integers their bits spell.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of individuals.	required
mutation	Annotated[float, ParameterSpec(domain=closed(0, 1))]	The probability that an individual mutates in a generation.	0.1
length	int	The number of loci of a genome.	20
seed	int	The seed of the fitness landscape.	0

Examples:

>>> from skcomplex.systems.biology.evolution import AdaptiveClimbingInteraction
>>> AdaptiveClimbingInteraction('climbing', target='population', mutation=0.1).replacement
'best-mutant'

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    mutation: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 0.1,
    length: int = 20,
    seed: int = 0,
) -> None:
    super().__init__(label, target=target, strategy='trait', replacement='best-mutant', mutation=mutation)
    self.length = length
    self.seed = seed

AlleleFrequencyElement(label, visual=None, *, x=0.5, Ns=0.0)

Bases: Element

Frequency of an allele in a large population under selection and genetic drift, Kimura's diffusion.

The frequency follows the Itô equation dx = N s x (1 - x) dt + sqrt(2 x (1 - x)) dW, the continuum limit of the Moran process with the scaled selection N s, in units of N generations. The frequencies 0 and 1, loss and fixation, are absorbing, and a boundary keeps the frequency between them.

Parameters:

Name	Type	Description	Default
label	str	The label of the allele or of the group of replicate populations.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=closed(0, 1))]	The initial frequency.	0.5
Ns	Annotated[float, ParameterSpec()]	The scaled selection, the population size times the excess fitness.	0.0

Examples:

>>> from skcomplex.systems.biology.evolution import AlleleFrequencyElement
>>> AlleleFrequencyElement('allele', Ns=2.0).Ns
2.0

Source code in src/skcomplex/systems/biology/evolution/_drift.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.5,
    Ns: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.Ns = Ns

AssemblyLevelElement(label, visual=None, *, N=0.0, k=1.0, alpha=1.0, N0=1.0)

Bases: Element

Objects of one assembly index, discovered from the objects one assembly step below them.

In a group ordered by assembly index, the objects of index a + 1 grow as dN_{a+1}/dt = k (N_a)^alpha, Eq. (3) of Sharma et al. (2022), where the first member reads the constant pool N_0 of building blocks. The selectivity alpha = 1 reuses every object built before, and alpha < 1 only a subset, the signature of selection.

Parameters:

Name	Type	Description	Default
label	str	The label of the level or of the group of levels.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
N	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial number of objects.	0.0
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of discovery.	1.0
alpha	Annotated[float, ParameterSpec(domain=closed(0, 1))]	The selectivity of reuse.	1.0
N0	float	The constant number of building blocks below the first level.	1.0

Examples:

>>> from skcomplex.systems.biology.evolution import AssemblyLevelElement
>>> AssemblyLevelElement('levels', alpha=0.5).alpha
0.5

Source code in src/skcomplex/systems/biology/evolution/_assembly.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    N: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 1.0,
    N0: float = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.N = N
    self.k = k
    self.alpha = alpha
    self.N0 = N0

CatalysisInteraction(label, *, target, omega=1.0, topology=((0.0, 1.0), (1.0, 0.0)))

Bases: ReplicatorInteraction

Catalysed replication of the molecules of a group, with a common outflow that keeps the total at 1.

Molecule i replicates at the rate f_i = lambda_i + omega sum_j T_ij x_j, with its own growth lambda_i and the 0 or 1 topology T of the catalytic links read as the weights of a total over the molecules. The outflow is the total replication, so the concentrations follow the replicator equation dx_i/dt = x_i (f_i - sum_j x_j f_j).

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of molecules.	required
omega	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The catalytic rate of a link.	1.0
topology	Sequence[Sequence[float]]	The catalytic links, 1 where molecule j catalyses molecule i.	((0.0, 1.0), (1.0, 0.0))

Examples:

>>> from skcomplex.systems.biology.evolution import CatalysisInteraction
>>> CatalysisInteraction('catalysis', target='soup', topology=((0, 1), (1, 0))).topology
((0, 1), (1, 0))

Source code in src/skcomplex/systems/biology/evolution/_prebiotic.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    omega: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    topology: Sequence[Sequence[float]] = ((0.0, 1.0), (1.0, 0.0)),
) -> None:
    super().__init__(label, target=target)
    self.omega = omega
    self.topology = topology

ChemicalElement(label, visual=None, *, x=0.0, tie=0.5, fresh=1.0, h=0.5, p=0.02, period=10.0, k=1)

Bases: Element

Molecular species of a pool whose relative population grows with the populations of its catalysts.

The populations follow the catalytic network: every species gains h times the populations of the species that catalyse it, and the populations are normalised, so they approach the Perron eigenvector of the network. Every period steps the least populated species, with ties broken at random, is replaced by a new species with k random catalysts, and every species catalyses the new one with probability p.

Parameters:

Name	Type	Description	Default
label	str	The label of the species or of the group of species.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The relative population.	0.0
tie	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The random key that breaks ties between the least populated species.	0.5
fresh	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	Whether the species appeared at the previous graph update.	1.0
h	Annotated[float, ParameterSpec(domain=open(0, oo))]	The growth of a step from the catalysts.	0.5
p	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability of every catalytic link of a new species.	0.02
period	Annotated[float, ParameterSpec(domain=open(0, oo))]	The number of steps between two graph updates.	10.0
k	int	The number of catalysts of a new species.	1

Examples:

>>> from skcomplex.systems.biology.evolution import ChemicalElement
>>> ChemicalElement('species').period
10.0

Source code in src/skcomplex/systems/biology/evolution/_networks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    tie: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    fresh: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 1.0,
    h: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
    period: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    k: int = 1,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.tie = tie
    self.fresh = fresh
    self.h = h
    self.p = p
    self.period = period
    self.k = k

ClimberElement(label, visual=VisualSpec('agents', {'color_by': 'fitness', 'color_range': (0.0, 1.0)}), *, genome=0.0, fitness=0.0, n=10.0, k=2.0, seed=0.0)

Bases: Element

Agent that walks an NK fitness landscape, flipping one random bit of its genome when that raises its fitness.

The fitness of a genome of n bits is the mean of the contributions of its bits, each a random function of the bit and the next k bits, fixed by seed.

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('agents', {'color_by': 'fitness', 'color_range': (0.0, 1.0)})
genome	Annotated[float, StateSpec(domain=NN)]	The bits of the genome.	0.0
fitness	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The fitness of the genome.	0.0
n	Annotated[float, ParameterSpec(domain=open(0, oo))]	The number of bits.	10.0
k	Annotated[float, ParameterSpec(domain=NN)]	The number of other bits each bit interacts with.	2.0
seed	Annotated[float, ParameterSpec(domain=NN)]	The seed of the landscape.	0.0

Examples:

>>> from skcomplex.core.components.elements import GroupElement
>>> from skcomplex.systems.biology.evolution import ClimberElement
>>> GroupElement('climbers', element=ClimberElement, size=10).size
10

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('agents', {'color_by': 'fitness', 'color_range': (0.0, 1.0)}),
    *,
    genome: Annotated[float, StateSpec(domain=NN)] = 0.0,
    fitness: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    k: Annotated[float, ParameterSpec(domain=NN)] = 2.0,
    seed: Annotated[float, ParameterSpec(domain=NN)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.genome = genome
    self.fitness = fitness
    self.n = n
    self.k = k
    self.seed = seed

climb(agent, system)

Flip a random bit, and keep the flip when it raises the fitness.

Parameters:

Name	Type	Description	Default
agent	ClimberView	The climber.	required
system	SystemView	The system.	required

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def climb(self, agent: ClimberView, system: SystemView) -> None:
    """Flip a random bit, and keep the flip when it raises the fitness.

    Args:
        agent:
            The climber.

        system:
            The system.
    """
    n, k, seed = int(agent.n), int(agent.k), int(agent.seed)
    mutant = int(agent.genome) ^ (1 << int(system.rng.integers(n)))
    fitness = ClimberElement._nk_fitness(mutant, n, k, seed)
    if fitness > agent.fitness:
        agent.genome = float(mutant)
        agent.fitness = fitness

CoevolutionaryClimbingInteraction(label, *, species, n=12, k=2, c=1, s=2, seed=0)

Bases: Interaction

Hill climbing of species on coupled NK landscapes, one point mutation per species and step, the NKCS model.

Every gene of a species contributes to its fitness through a random table of its own state, the states of the next k genes of the species, and the states of c random genes in each of s other species. The fitness of a species is the mean contribution of its n genes. At every step the species are picked in turn, and each flips one random bit and keeps the flip when it raises its fitness against the current genomes of the others. The tables, the partner species, and the coupled genes are fixed by seed.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
species	Annotated[str, RoleSpec()]	The path of the group of species.	required
n	int	The number of genes of every species.	12
k	int	The number of genes of its own species each gene depends on.	2
c	int	The number of genes of every partner species each gene depends on.	1
s	int	The number of partner species of every species.	2
seed	int	The seed of the landscapes and the couplings.	0

Examples:

>>> from skcomplex.systems.biology.evolution import CoevolutionaryClimbingInteraction
>>> CoevolutionaryClimbingInteraction('climbing', species='species', k=4).k
4

Source code in src/skcomplex/systems/biology/evolution/_nkcs.py

def __init__(
    self,
    label: str,
    *,
    species: Annotated[str, RoleSpec()],
    n: int = 12,
    k: int = 2,
    c: int = 1,
    s: int = 2,
    seed: int = 0,
) -> None:
    super().__init__(label)
    self.species = species
    self.n = n
    self.k = k
    self.c = c
    self.s = s
    self.seed = seed

climb(system)

Let every species in turn try one point mutation, and update the fitness of all.

Parameters:

Name	Type	Description	Default
system	SystemView	The system.	required

Source code in src/skcomplex/systems/biology/evolution/_nkcs.py

def climb(self, system: SystemView) -> None:
    """Let every species in turn try one point mutation, and update the fitness of all.

    Args:
        system:
            The system.
    """
    species = system.agents(self.species, CoevolvingSpeciesView)
    genomes = [int(item.genome) for item in species]
    for index, item in enumerate(species):
        mutant = genomes[index] ^ (1 << int(system.rng.integers(self.n)))
        trial = [*genomes[:index], mutant, *genomes[index + 1 :]]
        accepted = self._fitness(index, trial) > self._fitness(index, genomes)
        if accepted:
            genomes[index] = mutant
        item.moves = float(accepted)
    for index, item in enumerate(species):
        item.genome = float(genomes[index])
        item.fitness = self._fitness(index, genomes)

CoevolvingSpeciesElement(label, visual=None, *, genome=0.0, fitness=0.0, moves=0.0)

Bases: Element

Species of an NKCS ecosystem, a genome whose fitness depends on its own genes and on genes of other species.

A CoevolutionaryClimbingInteraction mutates the genome and sets the fitness and the number of accepted moves of the last step.

Parameters:

Name	Type	Description	Default
label	str	The label of the species or of the group of species.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
genome	Annotated[float, StateSpec(domain=NN)]	The bits of the genome, as an integer.	0.0
fitness	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The fitness of the genome against the current genomes of the other species.	0.0
moves	Annotated[float, StateSpec(domain=NN)]	Whether the species accepted a mutation in the last step.	0.0

Examples:

>>> from skcomplex.systems.biology.evolution import CoevolvingSpeciesElement
>>> CoevolvingSpeciesElement('species', genome=5.0).genome
5.0

Source code in src/skcomplex/systems/biology/evolution/_nkcs.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    genome: Annotated[float, StateSpec(domain=NN)] = 0.0,
    fitness: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    moves: Annotated[float, StateSpec(domain=NN)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.genome = genome
    self.fitness = fitness
    self.moves = moves

CompetingSpeciesElement(label, visual=None, *, x=0.0, links=0.0)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the species or of the group of species.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The relative population.	0.0
links	Annotated[object, StateSpec(vector=True)]	The effects of every species on this one.	0.0

Examples:

>>> from skcomplex.systems.biology.evolution import CompetingSpeciesElement
>>> CompetingSpeciesElement('species', x=0.1).x
0.1

Source code in src/skcomplex/systems/biology/evolution/_networks.py

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

CompetitiveSelectionInteraction(label, *, species, p=0.05, duration=20.0, step=0.05)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
species	Annotated[str, RoleSpec()]	The path of the group of species.	required
p	float	The probability of every link of a new species.	0.05
duration	float	The time of the population dynamics between replacements.	20.0
step	float	The Euler step of the population dynamics.	0.05

Examples:

>>> from skcomplex.systems.biology.evolution import CompetitiveSelectionInteraction
>>> CompetitiveSelectionInteraction('selection', species='species', p=0.1).p
0.1

Source code in src/skcomplex/systems/biology/evolution/_networks.py

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

select(system)

Run the population dynamics, then replace the least populated species.

Parameters:

Name	Type	Description	Default
system	SystemView	The system.	required

Source code in src/skcomplex/systems/biology/evolution/_networks.py

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]

CyclicSpeciesElement(label, visual=None, *, strategy=0)

Bases: Element

Site of a lattice of rock, paper, and scissors, 0, 1, and 2, where each strategy beats the one before it.

The sites update one at a time in random order by a LocalRuleInteraction, and a site adopts the strategy of a random neighbour that beats it.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
strategy	Annotated[float, StateSpec(domain=NN)]	The initial strategy, 0, 1, or 2.	0

Examples:

>>> from skcomplex.systems.biology.evolution import CyclicSpeciesElement
>>> CyclicSpeciesElement('site', strategy=2).strategy
2

Source code in src/skcomplex/systems/biology/evolution/_games.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    strategy: Annotated[float, StateSpec(domain=NN)] = 0,
) -> None:
    super().__init__(label, visual=visual)
    self.strategy = strategy

DependencyFieldInteraction(label, *, species, couplings)

Bases: Interaction

Field of every species, the sum of the couplings w of its edges to the species it depends on.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
species	Annotated[str, RoleSpec()]	The path of the group of species on the nodes.	required
couplings	Annotated[str, RoleSpec()]	The path of the variable group on the edges.	required

Examples:

>>> from skcomplex.systems.biology.evolution import DependencyFieldInteraction
>>> DependencyFieldInteraction('field', species='species', couplings='couplings').couplings
'couplings'

Source code in src/skcomplex/systems/biology/evolution/_networks.py

def __init__(
    self,
    label: str,
    *,
    species: Annotated[str, RoleSpec()],
    couplings: Annotated[str, RoleSpec()],
) -> None:
    super().__init__(label)
    self.species = species
    self.couplings = couplings

GenotypeClassElement(label, visual=None, *, x=0.0)

Bases: Element

Share of a population whose genomes lie a given number of point mutations away from the master sequence.

Parameters:

Name	Type	Description	Default
label	str	The label of the class.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial share.	0.0

Examples:

>>> from skcomplex.systems.biology.evolution import GenotypeClassElement
>>> GenotypeClassElement('k1', x=0.1).x
0.1

Source code in src/skcomplex/systems/biology/evolution/_prebiotic.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x

IndividualElement(label, visual=VisualSpec('agents', {'color_by': 'trait'}), *, trait=0.0, fitness=1.0)

Bases: Element

Individual of a finite population, with an integer trait or genome.

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('agents', {'color_by': 'trait'})
trait	Annotated[float, StateSpec(domain=ZZ)]	The trait, a type index or the bits of a genome.	0.0
fitness	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The fitness of the trait.	1.0

Examples:

>>> from skcomplex.core.components.elements import GroupElement
>>> from skcomplex.systems.biology.evolution import IndividualElement
>>> GroupElement('population', element=IndividualElement, size=10, values={'trait': 0.0}).size
10

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('agents', {'color_by': 'trait'}),
    *,
    trait: Annotated[float, StateSpec(domain=ZZ)] = 0.0,
    fitness: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.trait = trait
    self.fitness = fitness

LocusElement(label, visual=VisualSpec('image', {'low': -1.0, 'high': 1.0, 'history': 64}), *, s=-1.0, h=0.0)

Bases: Element

Locus of the genome shared by a monomorphic population, with an allele of -1 or 1.

Parameters:

Name	Type	Description	Default
label	str	The label of the locus or the group of loci.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('image', {'low': -1.0, 'high': 1.0, 'history': 64})
s	Annotated[float, StateSpec(domain=ZZ)]	The allele, -1 or 1.	-1.0
h	Annotated[float, ParameterSpec()]	The random field of the locus, the fitness gained by the allele 1 over the allele -1, halved.	0.0

Examples:

>>> from skcomplex.core.components.elements import GroupElement
>>> from skcomplex.systems.biology.evolution import LocusElement
>>> GroupElement('genome', element=LocusElement, size=10, values={'h': 0.05}).size
10

Source code in src/skcomplex/systems/biology/evolution/_substitutions.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('image', {'low': -1.0, 'high': 1.0, 'history': 64}),
    *,
    s: Annotated[float, StateSpec(domain=ZZ)] = -1.0,
    h: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.s = s
    self.h = h

MasterSequenceElement(label, visual=None, *, x=0.0, u=0.3, sigma=0.5)

Bases: GenotypeClassElement

Share of the population with the fittest genome, in a sharp-peak landscape where every other genome is less fit.

The master sequence has fitness 1 and every other genome 1 - sigma. A genome mutates at one locus with probability u per generation. The master sequence holds the parameters of the landscape that the maps of every class read.

Parameters:

Name	Type	Description	Default
label	str	The label of the class.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	float	The initial share.	0.0
u	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The mutation rate per genome and generation.	0.3
sigma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The selective disadvantage of the other genomes.	0.5

Examples:

>>> from skcomplex.systems.biology.evolution import MasterSequenceElement
>>> MasterSequenceElement('master', u=0.2).u
0.2

Source code in src/skcomplex/systems/biology/evolution/_prebiotic.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: float = 0.0,
    u: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
) -> None:
    super().__init__(label, visual=visual, x=x)
    self.u = u
    self.sigma = sigma

MoleculeElement(label, visual=None, *, x=0.5, growth=0.0)

Bases: Element

Self-replicating molecule of a prebiotic soup, with its concentration x and its own growth rate.

Parameters:

Name	Type	Description	Default
label	str	The label of the molecule or of the group of molecules.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[object, StateSpec(domain=closed_open(0, oo))]	The initial concentration.	0.5
growth	Annotated[object, ParameterSpec()]	The rate of replication without catalysis.	0.0

Examples:

>>> from skcomplex.systems.biology.evolution import MoleculeElement
>>> MoleculeElement('molecule', growth=0.1).growth
0.1

Source code in src/skcomplex/systems/biology/evolution/_prebiotic.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    growth: Annotated[object, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.growth = growth

MoranStepInteraction(label, *, target, advantage=0.0, mutation=0.0, phase=None)

Bases: EvolutionInteraction

One generation of the Moran process: as many birth and death events as there are individuals.

At every event an individual reproduces with probability proportional to its fitness (1 + advantage)^trait, its offspring mutates to the trait one above or one below with probability mutation / 2 each, and it replaces an individual chosen uniformly at random, so the size of the population stays constant. It is the Moran replacement of EvolutionInteraction with proportional selection and steps of the trait.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of individuals.	required
advantage	Annotated[float, ParameterSpec()]	The selective advantage of one unit of the trait.	0.0
mutation	Annotated[float, ParameterSpec(domain=closed(0, 1))]	The probability that an offspring mutates.	0.0
phase	str | None	The phase of the step the generation runs in, or None for the default phase.	None

Examples:

>>> from skcomplex.systems.biology.evolution import MoranStepInteraction
>>> MoranStepInteraction('moran', target='population', advantage=0.1).replacement
'moran'

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    advantage: Annotated[float, ParameterSpec()] = 0.0,
    mutation: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 0.0,
    phase: str | None = None,
) -> None:
    super().__init__(
        label,
        target=target,
        strategy='trait',
        selection='proportional',
        replacement='moran',
        mutation=mutation,
        alphabet=None,
        change='step',
        phase=phase,
    )
    self.advantage = advantage

NeutralNetworkAnalysis(phenotype_of, kind='network', alphabet='ACGU', state=None)

Bases: BaseAnalysis

Neutral networks of a map from sequences to phenotypes, with robustness and evolvability as Wagner defines them.

The neighbours of a genotype are the sequences one point mutation away. Its robustness is the number of its neighbours of the same phenotype, and its evolvability the number of different phenotypes among its neighbours. The robustness of a phenotype is the mean robustness of its genotypes, and its evolvability the number of different phenotypes among the neighbours of all its genotypes. The definitions follow Wagner (2008), as Raman (2021, §13.3) states them, so the phenotype of a genotype itself counts when a neighbour has it.

Parameters:

Name	Type	Description	Default
phenotype_of	Callable[[str], Hashable]	The function from a sequence to a phenotype key.	required
kind	str	'network', 'genotype-robustness', 'genotype-evolvability', 'phenotype-robustness', or 'phenotype-evolvability'.	'network'
alphabet	Sequence[str]	The letters of the sequences.	'ACGU'
state	str | None	The dotted name of the lattice state of a trajectory that holds one sequence per row, or None for the first state.	None

Attributes:

Name	Type	Description
genotypes_		The sequences of the source, in an array of their shape without the letters.
values_		For a genotype kind, one value per sequence of the source. For a phenotype kind, the value of every phenotype of the source, by key.
networks_		For the kind 'network', the graph per phenotype of the sequences of the source that one point mutation joins.

Raises:

Type	Description
AnalysisError	The kind is unknown, the mapping is not callable, or a sequence holds a letter outside the alphabet.

Examples:

>>> import itertools
>>> from skcomplex.systems.biology.evolution import NeutralNetworkAnalysis
>>> genotypes = [''.join(letters) for letters in itertools.product('01', repeat=4)]
>>> analysis = NeutralNetworkAnalysis(lambda sequence: sequence[0], 'genotype-robustness', alphabet='01')
>>> set(analysis.prepare(genotypes).run().values_.tolist())
{3.0}

Source code in src/skcomplex/systems/biology/evolution/_neutral_network_analysis.py

def __init__(
    self,
    phenotype_of: Callable[[str], Hashable],
    kind: str = 'network',
    alphabet: Sequence[str] = 'ACGU',
    state: str | None = None,
) -> None:
    self.phenotype_of = phenotype_of
    self.kind = kind
    self.alphabet = alphabet
    self.state = state

NucleotideElement(label, base='A')

Bases: Element

Site of an RNA sequence that holds one base.

Parameters:

Name	Type	Description	Default
label	str	The label of the element.	required
base	Annotated[str, StateSpec(domain=FiniteEnumeratedSet(list(BASES)))]	The base of the site.	'A'

Examples:

>>> from skcomplex.systems.biology.evolution import NucleotideElement
>>> NucleotideElement('site').base
'A'

Source code in src/skcomplex/systems/biology/evolution/_rna_walks.py

def __init__(
    self,
    label: str,
    base: Annotated[str, StateSpec(domain=FiniteEnumeratedSet(list(BASES)))] = 'A',
) -> None:
    super().__init__(label)
    self.base = base

PolymerElement(label, visual=None, *, length=1.0, born=0.0)

Bases: Element

Linear polymer of an assembly pool, with its length and the step it was made.

Parameters:

Name	Type	Description	Default
label	str	The label of the polymer or of the group of polymers.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
length	Annotated[float, StateSpec(domain=NN)]	The number of monomers.	1.0
born	Annotated[float, StateSpec(domain=NN)]	The step at which the polymer was made.	0.0

Examples:

>>> from skcomplex.systems.biology.evolution import PolymerElement
>>> PolymerElement('polymer', length=3.0).length
3.0

Source code in src/skcomplex/systems/biology/evolution/_assembly.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    length: Annotated[float, StateSpec(domain=NN)] = 1.0,
    born: Annotated[float, StateSpec(domain=NN)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.length = length
    self.born = born

PolymerJoiningInteraction(label, *, polymers, mode='undirected')

Bases: Interaction

Joining of two polymers of an assembly pool into a new one, at random or by reuse of the newest polymer.

At every step, undirected exploration joins two polymers drawn at random from the pool, and directed exploration joins the polymer made most recently with one drawn at random. The joined polymer, of the summed length, is added to the pool, and its parts stay in it.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
polymers	Annotated[str, RoleSpec()]	The path of the variable group of polymers.	required
mode	str	'undirected' or 'directed'.	'undirected'

Raises:

Type	Description
ValidationError	The mode is unknown.

Examples:

>>> from skcomplex.systems.biology.evolution import PolymerJoiningInteraction
>>> PolymerJoiningInteraction('joining', polymers='pool', mode='directed').mode
'directed'

Source code in src/skcomplex/systems/biology/evolution/_assembly.py

def __init__(self, label: str, *, polymers: Annotated[str, RoleSpec()], mode: str = 'undirected') -> None:
    super().__init__(label)
    self.polymers = polymers
    self.mode = mode

join(system)

Join two polymers and add the product to the pool.

Parameters:

Name	Type	Description	Default
system	SystemView	The system.	required

Source code in src/skcomplex/systems/biology/evolution/_assembly.py

def join(self, system: SystemView) -> None:
    """Join two polymers and add the product to the pool.

    Args:
        system:
            The system.
    """
    pool = system.agents(self.polymers, PolymerView)
    partner = pool[int(system.rng.integers(len(pool)))]
    if self.mode == 'directed':
        first = max(pool, key=lambda polymer: polymer.born)
    else:
        first = pool[int(system.rng.integers(len(pool)))]
    born = max(polymer.born for polymer in pool) + 1
    system.create(self.polymers, length=first.length + partner.length, born=born)

ReplicatorInteraction(label, *, target, fitness=None, state='x')

Bases: RateInteraction

Replicator equation of the shares of a group, dx_i/dt = x_i (f_i - sum_j x_j f_j), with a fitness f.

Every member of the group holds the share of a type, such as a strategy or a molecule. Its share grows with its fitness above the mean fitness of the group, so the total of the shares stays constant.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of types.	required
fitness	Callable[[SimpleNamespace], object] | None	The function of the namespace that returns the fitness of a member, such as lambda s: total(s.target.x, weights=payoff) for the payoff of a strategy against the others.	None
state	str	The name of the state of the members that holds the share.	'x'

Raises:

Type	Description
ValidationError	The fitness is not callable, raised by prepare.

Examples:

>>> from skcomplex.core.components.elements import GroupElement
>>> from skcomplex.core.components.laws.utils import total
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.biology.evolution import ReplicatorInteraction, StrategyElement
>>> population = GroupElement('population', element=StrategyElement, size=2, values={'x': [0.5, 0.5]})
>>> payoff = ((1.0, 1.0), (0.0, 0.0))
>>> replication = ReplicatorInteraction(
...     'replication', target='population', fitness=lambda s: total(s.target.x, weights=payoff)
... )
>>> system = System(elements=[population], interactions=[replication])
>>> trajectory = Simulation(system, dt=0.1).prepare().run(until=20.0).trajectory_
>>> [round(float(share), 3) for share in trajectory.values('population.x')[-1]]
[1.0, 0.0]

Source code in src/skcomplex/systems/biology/common/populations/_replicator.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    fitness: Callable[[SimpleNamespace], object] | None = None,
    state: str = 'x',
) -> None:
    super().__init__(label, target=target, state=state)
    self.fitness = fitness

SpeciesElement(label, visual=None, *, fitness=0.5, updates=200, neighbors='ring')

Bases: Element

Species of the Bak-Sneppen system, with a fitness in the unit interval that depends on its neighbours.

Every update replaces the fitness of the least fit species and of its neighbours by new uniform random numbers, and is recorded as an extinction event with the replaced species and its fitness. One step holds updates updates.

Parameters:

Name	Type	Description	Default
label	str	The label of the species or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
fitness	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial fitness.	0.5
updates	int	The number of updates per step.	200
neighbors	str	'ring' for the two neighbours on a ring or 'graph' for the neighbours of a node.	'ring'

Examples:

>>> from skcomplex.systems.biology.evolution import SpeciesElement
>>> SpeciesElement('species').updates
200

Source code in src/skcomplex/systems/biology/evolution/_coevolution.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    fitness: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    updates: int = 200,
    neighbors: str = 'ring',
) -> None:
    super().__init__(label, visual=visual)
    self.fitness = fitness
    self.updates = updates
    self.neighbors = neighbors

StrategyElement(label, visual=None, *, x=0.5)

Bases: Element

Strategy of an evolutionary game, with the share x of the population that plays it.

Parameters:

Name	Type	Description	Default
label	str	The label of the strategy or of the group of strategies.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[object, StateSpec(domain=closed_open(0, oo))]	The initial share.	0.5

Examples:

>>> from skcomplex.systems.biology.evolution import StrategyElement
>>> StrategyElement('hawks', x=0.4).x
0.4

Source code in src/skcomplex/systems/biology/evolution/_games.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x

SubstitutionInteraction(label, *, target, N=100.0, mu=0.01, F0=1.0, couplings=None, state='s')

Bases: TransitionInteraction

Origin and fixation of mutants that flip one locus of the genome of a monomorphic population.

Every locus of each of the N individuals mutates at the rate mu, so mutants of a given locus arise at the rate N mu. A mutant from the fitness f to f' fixes with the Moran probability (1 - f / f') / (1 - (f / f')^N), so the allele s of the locus moves between -1 and 1 by transitions at the product of the two. The fitness is F0 + sum_i h_i s_i + sum_ij J_ij s_i s_j / 2, and a flip of locus i changes it by -2 s_i (h_i + sum_j J_ij s_j).

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of loci.	required
N	Annotated[float, ParameterSpec(domain=open(0, oo))]	The size of the population.	100.0
mu	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate at which one locus of an individual mutates.	0.01
F0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The fitness of a genome without fields and couplings, above the sum of the |h_i| and |J_ij| so that every fitness is positive.	1.0
couplings	ArrayLike | None	The symmetric matrix J of couplings between loci, with a zero diagonal, or None for no couplings.	None
state	str	The name of the state of the loci that holds the allele.	's'

Examples:

>>> from skcomplex.systems.biology.evolution import SubstitutionInteraction
>>> SubstitutionInteraction('substitution', target='genome', N=50.0).N
50.0

Source code in src/skcomplex/systems/biology/evolution/_substitutions.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    N: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
    mu: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    F0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    couplings: npt.ArrayLike | None = None,
    state: str = 's',
) -> None:
    super().__init__(label, target=target, state=state, source=-1.0, destination=1.0)
    self.N = N
    self.mu = mu
    self.F0 = F0
    self.couplings = couplings

TangledNatureLifeCycleInteraction(label, *, target, k=5.0, mu=0.01, p_mut=0.01, p_kill=0.2, length=8, seed=0, phase='generation')

Bases: LifeCycleInteraction

Life cycle of the organisms of the Tangled Nature system, run by the visits of a phase drawn with replacement.

A visited organism of type S reproduces with probability 1 / (1 + exp(-H(S))), where H(S) = k / N sum_S' J(S, S') n(S') - mu N reads the composition n of the population of N organisms, and is replaced by two offspring, whose genome bits flip with probability p_mut each. Then a random organism dies with probability p_kill. The couplings J are uniform in [-1, 1] between different types, 0 within a type, and fixed by seed.

Parameters:

Name	Type	Description	Default
label	str	The label of the life cycle.	required
target	Annotated[str, RoleSpec()]	The path of the variable group of organisms.	required
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the couplings.	5.0
mu	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The limit of resources, which lowers reproduction in large populations.	0.01
p_mut	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability that a bit of an offspring's genome flips.	0.01
p_kill	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability that a visit kills a random organism.	0.2
length	int	The number of bits of a genome.	8
seed	int	The seed of the couplings.	0
phase	str | None	The phase whose visits run the life cycle.	'generation'

Examples:

>>> from skcomplex.systems.biology.evolution import TangledNatureLifeCycleInteraction
>>> TangledNatureLifeCycleInteraction('evolution', target='organisms').offspring
2

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0,
    mu: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    p_mut: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    p_kill: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    length: int = 8,
    seed: int = 0,
    phase: str | None = 'generation',
) -> None:
    super().__init__(
        label,
        target=target,
        birth=self.compute_birth,
        offspring=2,
        inherit=('trait',),
        mutate=self.mutate_genome,
        schedule='member',
        phase=phase,
        replace=True,
        cull=self.compute_kill,
        pick='random',
    )
    self.k = k
    self.mu = mu
    self.p_mut = p_mut
    self.p_kill = p_kill
    self.length = length
    self.seed = seed

compute_birth(agent, system)

Return the probability 1 / (1 + exp(-H(S))) that an organism of type S reproduces.

Parameters:

Name	Type	Description	Default
agent	IndividualView	The organism.	required
system	SystemView	The system view of the life cycle.	required

Returns:

Name	Type	Description
probability	float	The probability of reproduction.

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def compute_birth(self, agent: IndividualView, system: SystemView) -> float:
    """Return the probability `1 / (1 + exp(-H(S)))` that an organism of type `S` reproduces.

    Args:
        agent:
            The organism.

        system:
            The system view of the life cycle.

    Returns:
        probability:
            The probability of reproduction.
    """
    composition = system.composition('target', 'trait')
    size = sum(composition.values())
    row = self._get_couplings(self.seed, self.length)[int(agent.trait)]
    coupling = sum(row[int(np.asarray(genome))] * count for genome, count in composition.items())
    return float(expit(system.parameter('k') / size * coupling - system.parameter('mu') * size))

compute_kill(agent, system)

Return the probability p_kill that a visit kills a random organism.

Parameters:

Name	Type	Description	Default
agent	IndividualView	The visited organism.	required
system	SystemView	The system view of the life cycle.	required

Returns:

Name	Type	Description
probability	float	The probability of a death.

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def compute_kill(self, agent: IndividualView, system: SystemView) -> float:
    """Return the probability `p_kill` that a visit kills a random organism.

    Args:
        agent:
            The visited organism.

        system:
            The system view of the life cycle.

    Returns:
        probability:
            The probability of a death.
    """
    return system.parameter('p_kill')

mutate_genome(values, system)

Return the values of an offspring whose every genome bit flips with probability p_mut.

Parameters:

Name	Type	Description	Default
values	dict[str, object]	The values of the offspring, with the genome of its parent as trait.	required
system	SystemViewProtocol	The system view of the life cycle.	required

Returns:

Name	Type	Description
values	dict[str, object]	The values with the mutated genome.

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

def mutate_genome(self, values: dict[str, object], system: SystemViewProtocol) -> dict[str, object]:
    """Return the values of an offspring whose every genome bit flips with probability `p_mut`.

    Args:
        values:
            The values of the offspring, with the genome of its parent as `trait`.

        system:
            The system view of the life cycle.

    Returns:
        values:
            The values with the mutated genome.
    """
    genome = int(np.asarray(values['trait'], dtype=float))
    flips = np.flatnonzero(system.rng.random(self.length) < system.parameter('p_mut'))
    return {**values, 'trait': float(genome ^ sum(1 << int(locus) for locus in flips))}

TaxonElement(label, visual=None, *, h=1.0, age=0.0)

Bases: Element

Species of an ecosystem with its field, the sum of the couplings of the species it depends on, and its age.

A species whose field is negative goes extinct, and a new species takes its niche.

Parameters:

Name	Type	Description	Default
label	str	The label of the species or of the group of species.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
h	Annotated[float, StateSpec()]	The field, the sum of the couplings of its dependencies.	1.0
age	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The number of steps since the species appeared.	0.0

Examples:

>>> from skcomplex.systems.biology.evolution import TaxonElement
>>> TaxonElement('species').h
1.0

Source code in src/skcomplex/systems/biology/evolution/_networks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    h: Annotated[float, StateSpec()] = 1.0,
    age: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.h = h
    self.age = age

adaptive_climbing(size=10, *, u=0.1, length=20, seed=0)

Return a finite population that climbs a random fitness landscape by mutation and strong selection.

The population climbs to a local optimum in an adaptive walk, ever slower near the top. When every individual mutates in the same generation, which is likely only in small populations with a large u, it escapes to a lower genome and starts a new walk, the wandering regime.

Parameters:

Name	Type	Description	Default
size	int	The number of individuals.	10
u	float	The probability that an individual mutates in a generation.	0.1
length	int	The number of loci of a genome.	20
seed	int	The seed of the fitness landscape.	0

Returns:

Name	Type	Description
system	System	The system, with the group population.

Examples:

>>> from skcomplex.systems.biology.evolution import adaptive_climbing
>>> adaptive_climbing(5).elements[0].prepare().size_
5

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

@register_system(
    field='biology/evolution',
    tags=('fitness landscape', 'agents', 'random energy system', 'stochastic escape'),
    references=('Gros 2015',),
    defaults={'until': 200},
)
def adaptive_climbing(
    size: int = 10,
    *,
    u: float = 0.1,
    length: int = 20,
    seed: int = 0,
) -> System:
    """Return a finite population that climbs a random fitness landscape by mutation and strong selection.

    The population climbs to a local optimum in an adaptive walk, ever slower near the top. When every individual
    mutates in the same generation, which is likely only in small populations with a large `u`, it escapes to a lower
    genome and starts a new walk, the wandering regime.

    Args:
        size:
            The number of individuals.

        u:
            The probability that an individual mutates in a generation.

        length:
            The number of loci of a genome.

        seed:
            The seed of the fitness landscape.

    Returns:
        system:
            The system, with the group `population`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import adaptive_climbing
        >>> adaptive_climbing(5).elements[0].prepare().size_
        5
    """
    fitness = AdaptiveClimbingInteraction._random_fitness(0, seed)
    return System(
        elements=[
            GroupElement(
                'population',
                element=IndividualElement,
                size=size,
                values={'trait': 0.0, 'fitness': fitness},
                variable=True,
                visual=VisualSpec('agents', {'color_by': 'fitness'}),
            ),
        ],
        interactions=[
            AdaptiveClimbingInteraction('climbing', target='population', mutation=u, length=length, seed=seed),
        ],
    )

assembly_growth(levels=10, *, k=1.0, alpha=1.0)

Return the numbers of objects at successive assembly indices, each discovered from the level below.

The objects of assembly index a + 1 grow as dN_{a+1}/dt = k (N_a)^alpha from a constant pool of building blocks. Without selection, alpha = 1, every object built before is reused and N_a(t) = (k t)^a / a!, the exponential expansion of the assembly possible. Selection, alpha < 1, reuses only a subset, so the discovery of new objects slows down and the deeper levels fill far later.

Parameters:

Name	Type	Description	Default
levels	int	The number of assembly indices above the building blocks.	10
k	float	The rate of discovery.	1.0
alpha	float	The selectivity of reuse, 1 for no selection.	1.0

Returns:

Name	Type	Description
system	System	The system, with the numbers levels.N of the group levels, ordered by assembly index from 1.

Examples:

>>> from skcomplex.systems.biology.evolution import assembly_growth
>>> assembly_growth(5).elements[0].prepare().size_
5

Source code in src/skcomplex/systems/biology/evolution/_assembly.py

@register_system(
    field='biology/evolution',
    tags=('assembly theory', 'selection', 'open-ended evolution', 'ode'),
    references=('Sharma et al. 2023',),
    defaults={'until': 10.0, 'dt': 0.1},
)
def assembly_growth(levels: int = 10, *, k: float = 1.0, alpha: float = 1.0) -> System:
    """Return the numbers of objects at successive assembly indices, each discovered from the level below.

    The objects of assembly index `a + 1` grow as `dN_{a+1}/dt = k (N_a)^alpha` from a constant pool of building
    blocks. Without selection, `alpha = 1`, every object built before is reused and `N_a(t) = (k t)^a / a!`, the
    exponential expansion of the assembly possible. Selection, `alpha < 1`, reuses only a subset, so the discovery of
    new objects slows down and the deeper levels fill far later.

    Args:
        levels:
            The number of assembly indices above the building blocks.

        k:
            The rate of discovery.

        alpha:
            The selectivity of reuse, 1 for no selection.

    Returns:
        system:
            The system, with the numbers `levels.N` of the group `levels`, ordered by assembly index from 1.

    Examples:
        >>> from skcomplex.systems.biology.evolution import assembly_growth
        >>> assembly_growth(5).elements[0].prepare().size_
        5
    """
    group = GroupElement('levels', element=AssemblyLevelElement, size=levels, values={'k': k, 'alpha': alpha})
    return System(elements=[group])

bak_sneppen(size=200, *, graph=None, updates=None, random_state=0)

Return the Bak-Sneppen system, whose fitnesses self-organise above a critical value in avalanches of extinction.

On a ring the fitnesses gather above the critical value 0.667, and the minimum fitness rises towards it through avalanches of every size.

Parameters:

Name	Type	Description	Default
size	int	The number of species on the ring, unused on a graph.	200
graph	Graph | None	The graph whose nodes are the species, or None for a ring.	None
updates	int | None	The number of updates per step, or None for one update per species.	None
random_state	int | Generator | None	The seed or generator of the initial fitnesses.	0

Returns:

Name	Type	Description
system	System	The system, with the fitnesses species.fitness, and extinction events with the site and the minimum fitness of every update.

Examples:

>>> from skcomplex.systems.biology.evolution import bak_sneppen
>>> bak_sneppen(50).elements[0].values['updates']
50

Source code in src/skcomplex/systems/biology/evolution/_coevolution.py

@register_system(
    field='biology/evolution',
    tags=('lattice', 'network', 'self-organized criticality', 'extremal dynamics', 'stochastic'),
    references=('Bak and Sneppen 1993', 'Gros 2015', 'Jensen 2023'),
    defaults={'until': 2000, 'dt': 10},
)
def bak_sneppen(
    size: int = 200,
    *,
    graph: nx.Graph | None = None,
    updates: int | None = None,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Bak-Sneppen system, whose fitnesses self-organise above a critical value in avalanches of extinction.

    On a ring the fitnesses gather above the critical value 0.667, and the minimum fitness rises towards it through
    avalanches of every size.

    Args:
        size:
            The number of species on the ring, unused on a graph.

        graph:
            The graph whose nodes are the species, or `None` for a ring.

        updates:
            The number of updates per step, or `None` for one update per species.

        random_state:
            The seed or generator of the initial fitnesses.

    Returns:
        system:
            The system, with the fitnesses `species.fitness`, and `extinction` events with the site and the minimum
            fitness of every update.

    Examples:
        >>> from skcomplex.systems.biology.evolution import bak_sneppen
        >>> bak_sneppen(50).elements[0].values['updates']
        50
    """
    space = (LatticeSpace(size) if graph is None else NetworkSpace(graph)).prepare()
    count = int(space.size_ or 0)
    fitness = np.random.default_rng(random_state).random(count)
    species = LatticeElement(
        'species',
        element=SpeciesElement,
        space=space,
        values={
            'fitness': fitness,
            'updates': count if updates is None else updates,
            'neighbors': 'ring' if graph is None else 'graph',
        },
        visual=FITNESS_VISUAL if graph is None else FITNESS_SITES,
    )
    return System(elements=[species])

competitive_jain_krishna(n=50, *, p=0.05, duration=20.0, seed=1)

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.

Parameters:

Name	Type	Description	Default
n	int	The number of species.	50
p	float	The probability of every link.	0.05
duration	float	The time of the population dynamics between replacements.	20.0
seed	int | None	The seed of the initial matrix.	1

Returns:

Name	Type	Description
system	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

Source code in src/skcomplex/systems/biology/evolution/_networks.py

@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())

eigen_quasispecies(classes=12, *, u=0.3, sigma=0.5)

Return the shares of a population at every number of mutations from the master sequence, the linear chain system.

Below the error threshold, u < sigma, a quasispecies forms around the master sequence, whose share settles at x0 = 1 - u / sigma. Above it the master sequence is lost, the error catastrophe. Back mutations and multiple mutations are neglected, so the chain holds only forward mutations, and the last class collects the far mutants.

Parameters:

Name	Type	Description	Default
classes	int	The number of classes of mutants after the master sequence.	12
u	float	The mutation rate per genome and generation.	0.3
sigma	float	The selective disadvantage of the mutants.	0.5

Returns:

Name	Type	Description
system	System	The system, with the shares master.x and k1.x to k<classes>.x.

Examples:

>>> from skcomplex.systems.biology.evolution import eigen_quasispecies
>>> [element.label for element in eigen_quasispecies(2).elements]
['master', 'k1', 'k2']

Source code in src/skcomplex/systems/biology/evolution/_prebiotic.py

@register_system(
    field='biology/evolution',
    tags=('quasispecies', 'error threshold', 'discrete map', 'mutation', 'selection'),
    references=('Eigen 1971', 'Gros 2015'),
    defaults={'until': 100, 'dt': 1},
)
def eigen_quasispecies(classes: int = 12, *, u: float = 0.3, sigma: float = 0.5) -> System:
    """Return the shares of a population at every number of mutations from the master sequence, the linear chain system.

    Below the error threshold, `u < sigma`, a quasispecies forms around the master sequence, whose share settles at
    `x0 = 1 - u / sigma`. Above it the master sequence is lost, the error catastrophe. Back mutations and multiple
    mutations are neglected, so the chain holds only forward mutations, and the last class collects the far mutants.

    Args:
        classes:
            The number of classes of mutants after the master sequence.

        u:
            The mutation rate per genome and generation.

        sigma:
            The selective disadvantage of the mutants.

    Returns:
        system:
            The system, with the shares `master.x` and `k1.x` to `k<classes>.x`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import eigen_quasispecies
        >>> [element.label for element in eigen_quasispecies(2).elements]
        ['master', 'k1', 'k2']
    """
    labels = ['master', *(f'k{k}' for k in range(1, classes + 1))]
    elements = [
        MasterSequenceElement('master', x=1.0, u=u, sigma=sigma),
        *(GenotypeClassElement(label) for label in labels[1:]),
    ]
    selection = MapInteraction('selection', target='master', state='x', function=_replicate_master)
    mutations = [
        MapInteraction(
            f'mutation_{target}',
            target=target,
            state='x',
            function=partial(_replicate_mutants, first=source == 'master', last=target == labels[-1]),
            sources=(source, 'master'),
        )
        for source, target in itertools.pairwise(labels)
    ]
    interactions = [selection, *mutations]
    return System(elements=elements, interactions=interactions)

fold_sequence(sequence) cached

Return the minimum free energy secondary structure of an RNA sequence at 37 °C and its energy, memoised.

Parameters:

Name	Type	Description	Default
sequence	str	The sequence over the letters ACGU.	required

Returns:

Name	Type	Description
fold	tuple[str, float]	The structure in dot-bracket notation and its free energy in kcal/mol.

Raises:

Type	Description
MissingExtraError	The extra rna is not installed.

Examples:

>>> from skcomplex.systems.biology.evolution import fold_sequence
>>> structure, energy = fold_sequence('GGGGAAACCCC')
>>> structure, energy < 0
('((((...))))', True)

Source code in src/skcomplex/systems/biology/evolution/_rna_walks.py

@lru_cache(maxsize=FOLD_CACHE)
def fold_sequence(sequence: str) -> tuple[str, float]:
    """Return the minimum free energy secondary structure of an RNA sequence at 37 °C and its energy, memoised.

    Args:
        sequence:
            The sequence over the letters `ACGU`.

    Returns:
        fold:
            The structure in dot-bracket notation and its free energy in kcal/mol.

    Raises:
        MissingExtraError:
            The extra `rna` is not installed.

    Examples:
        >>> from skcomplex.systems.biology.evolution import fold_sequence
        >>> structure, energy = fold_sequence('GGGGAAACCCC')
        >>> structure, energy < 0
        ('((((...))))', True)
    """
    structure, energy = import_extra('RNA', SYSTEM, EXTRA).fold_compound(sequence).mfe()
    return str(structure), float(energy)

hypercycles(molecules=4, *, alpha=0.1, omega=1.0, topology='cycle')

Return a network of self-replicating molecules whose catalysis keeps a set of them alive, a hypercycle.

The growth rates rise as lambda_i = alpha i. In a 'cycle', every molecule catalyses the next one. Up to four molecules the hypercycle keeps all of them at a fixed point, and from five some members die out and the others settle at a fixed point. In a 'complete' network every molecule catalyses every other, and about sqrt(2 omega / alpha) of the fastest replicators survive while the others die out.

Parameters:

Name	Type	Description	Default
molecules	int	The number of molecules.	4
alpha	float	The step of the autocatalytic growth rates.	0.1
omega	float	The catalytic rate of a link.	1.0
topology	str	'cycle' or 'complete'.	'cycle'

Returns:

Name	Type	Description
system	System	The system, with the concentrations soup.x of the molecules, members of the group soup.

Raises:

Type	Description
ValidationError	The topology is unknown.

Examples:

>>> from skcomplex.systems.biology.evolution import hypercycles
>>> hypercycles(3).interactions[0].topology
((0.0, 0.0, 1.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0))

Source code in src/skcomplex/systems/biology/evolution/_prebiotic.py

@register_system(
    field='biology/evolution',
    tags=('hypercycle', 'autocatalysis', 'prebiotic', 'ode'),
    references=('Eigen and Schuster 1977', 'Gros 2015'),
    defaults={'until': 50.0, 'dt': 0.5},
)
def hypercycles(
    molecules: int = 4,
    *,
    alpha: float = 0.1,
    omega: float = 1.0,
    topology: str = 'cycle',
) -> System:
    """Return a network of self-replicating molecules whose catalysis keeps a set of them alive, a hypercycle.

    The growth rates rise as `lambda_i = alpha i`. In a `'cycle'`, every molecule catalyses the next one. Up to four
    molecules the hypercycle keeps all of them at a fixed point, and from five some members die out and the others
    settle at a fixed point. In a `'complete'` network every molecule catalyses every other, and about
    `sqrt(2 omega / alpha)` of the fastest replicators survive while the others die out.

    Args:
        molecules:
            The number of molecules.

        alpha:
            The step of the autocatalytic growth rates.

        omega:
            The catalytic rate of a link.

        topology:
            `'cycle'` or `'complete'`.

    Returns:
        system:
            The system, with the concentrations `soup.x` of the molecules, members of the group `soup`.

    Raises:
        ValidationError:
            The topology is unknown.

    Examples:
        >>> from skcomplex.systems.biology.evolution import hypercycles
        >>> hypercycles(3).interactions[0].topology
        ((0.0, 0.0, 1.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0))
    """
    if topology == 'cycle':
        links = np.roll(np.eye(molecules), 1, axis=0)
    elif topology == 'complete':
        links = 1 - np.eye(molecules)
    else:
        error_msg = f'Unknown topology {topology!r}. Use cycle or complete.'
        raise ValidationError(error_msg)
    soup = GroupElement(
        'soup',
        element=MoleculeElement,
        size=molecules,
        values={
            'x': [1 / molecules] * molecules,
            'growth': [float(alpha * index) for index in range(1, molecules + 1)],
        },
    )
    catalysis = CatalysisInteraction(
        'catalysis',
        target='soup',
        omega=omega,
        topology=tuple(tuple(float(value) for value in row) for row in links),
    )
    return System(elements=[soup], interactions=[catalysis])

jain_krishna(n=40, *, p=0.02, h=0.5, period=10, seed=1)

Return the system of Jain and Krishna, where selection of the fittest builds autocatalytic sets of molecules.

The populations of n molecular species approach the Perron eigenvector of their catalytic network between graph updates, every period steps. At every graph update the least populated species is replaced by a new one, linked to the others with probability p. Once a cycle of catalysts forms by chance, the species it supports win, and the autocatalytic set grows until it spans the network, so the number of links and the populated species jump.

Parameters:

Name	Type	Description	Default
n	int	The number of species.	40
p	float	The probability of every catalytic link.	0.02
h	float	The growth of a step from the catalysts.	0.5
period	int	The number of steps between two graph updates.	10
seed	int | None	The seed of the initial network.	1

Returns:

Name	Type	Description
system	System	The system, with the populations species.x of the variable group species on the nodes and the variable group links on the edges, where an edge from i to j means that j catalyses i.

Examples:

>>> from skcomplex.systems.biology.evolution import jain_krishna
>>> round(float(jain_krishna(n=10).elements[0].values['x'].sum()), 6)
1.0

Source code in src/skcomplex/systems/biology/evolution/_networks.py

@register_system(
    field='biology/evolution',
    tags=('adaptive network', 'autocatalytic set', 'Perron eigenvalue', 'origin of life', 'stochastic map'),
    references=('Jain and Krishna 1998', 'Jain and Krishna 2001', 'Thurner et al. 2018'),
    defaults={'until': 1500, 'dt': 5.0},
)
def jain_krishna(
    n: int = 40,
    *,
    p: float = 0.02,
    h: float = 0.5,
    period: int = 10,
    seed: int | None = 1,
) -> System:
    """Return the system of Jain and Krishna, where selection of the fittest builds autocatalytic sets of molecules.

    The populations of `n` molecular species approach the Perron eigenvector of their catalytic network between graph
    updates, every `period` steps. At every graph update the least populated species is replaced by a new one, linked
    to the others with probability `p`. Once a cycle of catalysts forms by chance, the species it supports win, and the
    autocatalytic set grows until it spans the network, so the number of links and the populated species jump.

    Args:
        n:
            The number of species.

        p:
            The probability of every catalytic link.

        h:
            The growth of a step from the catalysts.

        period:
            The number of steps between two graph updates.

        seed:
            The seed of the initial network.

    Returns:
        system:
            The system, with the populations `species.x` of the variable group `species` on the nodes and the variable
            group `links` on the edges, where an edge from `i` to `j` means that `j` catalyses `i`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import jain_krishna
        >>> round(float(jain_krishna(n=10).elements[0].values['x'].sum()), 6)
        1.0
    """
    rng = np.random.default_rng(seed)
    graph = nx.gnp_random_graph(n, p, seed=int(rng.integers(2**31)), directed=True)
    space = NetworkSpace(graph)
    species = GroupElement(
        'species',
        element=ChemicalElement,
        space=space,
        on='sites',
        variable=True,
        capacity=2 * n,
        values={
            'x': np.full(n, 1 / n),
            'tie': rng.random(n),
            'fresh': 0.0,
            'h': h,
            'p': p,
            'period': float(period),
            'k': max(1, round(p * n)),
        },
        visual=VisualSpec('sites', {'states': ['x']}),
    )
    links = GroupElement(
        'links',
        element=Element,
        space=space,
        on='edges',
        variable=True,
        capacity=n * n,
        visual=HIDDEN,
    )
    return System(elements=[species, links], time=DiscreteTime())

kimura_diffusion(x=0.1, *, Ns=2.0, replicates=500)

Return replicate populations whose allele frequencies diffuse by drift until the allele is lost or fixed.

Each replicate follows dx = N s x (1 - x) dt + sqrt(2 x (1 - x)) dW, whose density obeys Kimura's diffusion equation dp/dt = -N s d/dx [x (1 - x) p] + d^2/dx^2 [x (1 - x) p]. Every replicate ends at loss or fixation, and an allele that starts at the frequency x fixes with the probability (1 - exp(-N s x)) / (1 - exp(-N s)), or x for a neutral allele.

Parameters:

Name	Type	Description	Default
x	float	The initial frequency of the allele.	0.1
Ns	float	The scaled selection, the population size times the excess fitness.	2.0
replicates	int	The number of independent replicate populations.	500

Returns:

Name	Type	Description
system	System	The system, with the frequencies alleles.x of the replicates.

Examples:

>>> from skcomplex.systems.biology.evolution import kimura_diffusion
>>> kimura_diffusion(replicates=10).elements[0].prepare().size_
10

Source code in src/skcomplex/systems/biology/evolution/_drift.py

@register_system(
    field='biology/evolution',
    tags=('population genetics', 'genetic drift', 'fixation', 'sde', 'diffusion'),
    references=('Kimura 1955', 'Bertin 2021'),
    defaults={'until': 10.0, 'dt': 0.05, 'solver': EulerMaruyamaSolver(dt=0.001)},
)
def kimura_diffusion(x: float = 0.1, *, Ns: float = 2.0, replicates: int = 500) -> System:
    """Return replicate populations whose allele frequencies diffuse by drift until the allele is lost or fixed.

    Each replicate follows `dx = N s x (1 - x) dt + sqrt(2 x (1 - x)) dW`, whose density obeys Kimura's diffusion
    equation `dp/dt = -N s d/dx [x (1 - x) p] + d^2/dx^2 [x (1 - x) p]`. Every replicate ends at loss or fixation, and
    an allele that starts at the frequency `x` fixes with the probability `(1 - exp(-N s x)) / (1 - exp(-N s))`, or `x`
    for a neutral allele.

    Args:
        x:
            The initial frequency of the allele.

        Ns:
            The scaled selection, the population size times the excess fitness.

        replicates:
            The number of independent replicate populations.

    Returns:
        system:
            The system, with the frequencies `alleles.x` of the replicates.

    Examples:
        >>> from skcomplex.systems.biology.evolution import kimura_diffusion
        >>> kimura_diffusion(replicates=10).elements[0].prepare().size_
        10
    """
    alleles = GroupElement('alleles', element=AlleleFrequencyElement, size=replicates, values={'x': x, 'Ns': Ns})
    return System(elements=[alleles])

moran_process(size=50, *, mutants=1, selection=0.0, mutation=0.0)

Return the Moran process, a population of constant size where fitter traits fix by selection and drift.

The population starts with mutants individuals of trait 1 among individuals of trait 0. Without mutation one trait eventually fixes, the mutant with probability (1 - 1 / r^m) / (1 - 1 / r^N) for the relative fitness r = 1 + selection. With mutation the traits diffuse along the trait axis, in clusters rather than as a spreading Gaussian, which distinguishes evolutionary diffusion from a random walk.

Parameters:

Name	Type	Description	Default
size	int	The number of individuals.	50
mutants	int	The initial number of individuals of trait 1.	1
selection	float	The selective advantage of one unit of the trait.	0.0
mutation	float	The probability that an offspring's trait mutates by one unit.	0.0

Returns:

Name	Type	Description
system	System	The system, with the group population.

Examples:

>>> from skcomplex.systems.biology.evolution import moran_process
>>> float(moran_process(10, mutants=3).elements[0].values['trait'].sum())
3.0

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

@register_system(
    field='biology/evolution',
    tags=('population genetics', 'agents', 'fixation', 'genetic drift', 'mutation'),
    references=('Moran 1958', 'Nowak 2006', 'Jensen 2023'),
    defaults={'until': 100},
)
def moran_process(
    size: int = 50,
    *,
    mutants: int = 1,
    selection: float = 0.0,
    mutation: float = 0.0,
) -> System:
    """Return the Moran process, a population of constant size where fitter traits fix by selection and drift.

    The population starts with `mutants` individuals of trait 1 among individuals of trait 0. Without mutation one trait
    eventually fixes, the mutant with probability `(1 - 1 / r^m) / (1 - 1 / r^N)` for the relative fitness
    `r = 1 + selection`. With mutation the traits diffuse along the trait axis, in clusters rather than as a spreading
    Gaussian, which distinguishes evolutionary diffusion from a random walk.

    Args:
        size:
            The number of individuals.

        mutants:
            The initial number of individuals of trait 1.

        selection:
            The selective advantage of one unit of the trait.

        mutation:
            The probability that an offspring's trait mutates by one unit.

    Returns:
        system:
            The system, with the group `population`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import moran_process
        >>> float(moran_process(10, mutants=3).elements[0].values['trait'].sum())
        3.0
    """
    traits = np.zeros(size)
    traits[:mutants] = 1.0
    return System(
        elements=[
            GroupElement(
                'population',
                element=IndividualElement,
                size=size,
                values={'trait': traits, 'fitness': (1 + selection) ** traits},
                variable=True,
            ),
        ],
        interactions=[MoranStepInteraction('moran', target='population', advantage=selection, mutation=mutation)],
    )

mutate_point(row, random)

Return a row of indices of bases with the base at one random site replaced by one of the other three.

Parameters:

Name	Type	Description	Default
row	ndarray	The indices of the bases, which the function changes.	required
random	Generator	The random generator.	required

Returns:

Name	Type	Description
row	ndarray	The row after the mutation.

Examples:

>>> import numpy as np
>>> from skcomplex.systems.biology.evolution import mutate_point
>>> row = mutate_point(np.zeros(4), np.random.default_rng(0))
>>> int((row != 0).sum())
1

Source code in src/skcomplex/systems/biology/evolution/_rna_walks.py

def mutate_point(row: np.ndarray, random: np.random.Generator) -> np.ndarray:
    """Return a row of indices of bases with the base at one random site replaced by one of the other three.

    Args:
        row:
            The indices of the bases, which the function changes.

        random:
            The random generator.

    Returns:
        row:
            The row after the mutation.

    Examples:
        >>> import numpy as np
        >>> from skcomplex.systems.biology.evolution import mutate_point
        >>> row = mutate_point(np.zeros(4), np.random.default_rng(0))
        >>> int((row != 0).sum())
        1
    """
    site = random.integers(row.size)
    row[site] = (row[site] + random.integers(1, len(BASES))) % len(BASES)
    return row

nk_landscape(climbers=50, *, n=10, k=2, seed=0, random_state=0)

Return agents that climb Kauffman's NK fitness landscape from random genomes until they reach local peaks.

The landscape grows more rugged with k, from a single peak at k = 0 to a random landscape at k = n - 1, so the climbers end on more peaks of lower fitness.

Parameters:

Name	Type	Description	Default
climbers	int	The number of agents.	50
n	int	The number of bits of a genome.	10
k	int	The number of other bits each bit interacts with.	2
seed	int	The seed of the landscape.	0
random_state	int | Generator | None	The seed or generator of the initial genomes.	0

Returns:

Name	Type	Description
system	System	The system, with the group climbers.

Examples:

>>> from skcomplex.systems.biology.evolution import nk_landscape
>>> nk_landscape(5).elements[0].prepare().size_
5

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

@register_system(
    field='biology/evolution',
    tags=('fitness landscape', 'agents', 'NK system', 'hill climbing', 'ruggedness'),
    references=('Kauffman and Levin 1987', 'Downey 2018', 'Thurner et al. 2018'),
    defaults={'until': 50},
)
def nk_landscape(
    climbers: int = 50,
    *,
    n: int = 10,
    k: int = 2,
    seed: int = 0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return agents that climb Kauffman's NK fitness landscape from random genomes until they reach local peaks.

    The landscape grows more rugged with `k`, from a single peak at `k = 0` to a random landscape at `k = n - 1`, so
    the climbers end on more peaks of lower fitness.

    Args:
        climbers:
            The number of agents.

        n:
            The number of bits of a genome.

        k:
            The number of other bits each bit interacts with.

        seed:
            The seed of the landscape.

        random_state:
            The seed or generator of the initial genomes.

    Returns:
        system:
            The system, with the group `climbers`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import nk_landscape
        >>> nk_landscape(5).elements[0].prepare().size_
        5
    """
    genomes = np.random.default_rng(random_state).integers(0, 2**n, size=climbers)
    fitness = [ClimberElement._nk_fitness(int(genome), n, k, seed) for genome in genomes]
    group = GroupElement(
        'climbers',
        element=ClimberElement,
        size=climbers,
        values={
            'genome': genomes.astype(float),
            'fitness': np.asarray(fitness),
            'n': float(n),
            'k': float(k),
            'seed': float(seed),
        },
        variable=True,
    )
    return System(elements=[group])

nkcs(n_species=10, *, n=12, k=2, c=1, s=2, seed=0, random_state=0)

Return the NKCS model of Kauffman and Johnsen, species climbing NK landscapes coupled to each other's genes.

Every gene of a species depends on k genes of its own species and on c genes in each of s partner species. A species keeps a point mutation when it raises its fitness, which changes the landscapes of the species that depend on it, so adaptations trigger coevolutionary avalanches. Rugged landscapes, k above about c s, freeze the ecosystem once every species sits on a peak, and smoother ones keep it in a chaotic regime of endless adaptation.

Parameters:

Name	Type	Description	Default
n_species	int	The number of species.	10
n	int	The number of genes of every species.	12
k	int	The number of genes of its own species each gene depends on.	2
c	int	The number of genes of every partner species each gene depends on.	1
s	int	The number of partner species of every species.	2
seed	int	The seed of the landscapes and the couplings.	0
random_state	int | Generator | None	The seed or generator of the initial genomes.	0

Returns:

Name	Type	Description
system	System	The system, with the genomes species.genome, fitnesses species.fitness, and accepted mutations species.moves of the group species.

Examples:

>>> from skcomplex.systems.biology.evolution import nkcs
>>> nkcs(5).elements[0].prepare().size_
5

Source code in src/skcomplex/systems/biology/evolution/_nkcs.py

@register_system(
    field='biology/evolution',
    tags=('coevolution', 'fitness landscape', 'avalanches', 'criticality', 'agents'),
    references=('Kauffman and Johnsen 1991', 'Thurner et al. 2018'),
    defaults={'until': 200},
)
def nkcs(
    n_species: int = 10,
    *,
    n: int = 12,
    k: int = 2,
    c: int = 1,
    s: int = 2,
    seed: int = 0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the NKCS model of Kauffman and Johnsen, species climbing NK landscapes coupled to each other's genes.

    Every gene of a species depends on `k` genes of its own species and on `c` genes in each of `s` partner species.
    A species keeps a point mutation when it raises its fitness, which changes the landscapes of the species that
    depend on it, so adaptations trigger coevolutionary avalanches. Rugged landscapes, `k` above about `c s`, freeze
    the ecosystem once every species sits on a peak, and smoother ones keep it in a chaotic regime of endless
    adaptation.

    Args:
        n_species:
            The number of species.

        n:
            The number of genes of every species.

        k:
            The number of genes of its own species each gene depends on.

        c:
            The number of genes of every partner species each gene depends on.

        s:
            The number of partner species of every species.

        seed:
            The seed of the landscapes and the couplings.

        random_state:
            The seed or generator of the initial genomes.

    Returns:
        system:
            The system, with the genomes `species.genome`, fitnesses `species.fitness`, and accepted mutations
            `species.moves` of the group `species`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import nkcs
        >>> nkcs(5).elements[0].prepare().size_
        5
    """
    genomes = np.random.default_rng(random_state).integers(0, 2**n, size=n_species).astype(float)
    species = GroupElement('species', element=CoevolvingSpeciesElement, size=n_species, values={'genome': genomes})
    climbing = CoevolutionaryClimbingInteraction('climbing', species='species', n=n, k=k, c=c, s=s, seed=seed)
    return System(elements=[species], interactions=[climbing], time=DiscreteTime())

origin_fixation_walk(loci=96, *, N=100.0, mu=0.01, field=0.01, coupling=0.0, F0=1.0, random_state=0)

Return a monomorphic population that climbs a fitness landscape as mutants arise and fix one at a time.

Under strong selection and weak mutation the population carries one genome of loci alleles -1 or 1, and every substitution of a mutant moves it to a neighbouring genome. The walk climbs the landscape and settles on the stationary distribution Q ~ f^(N - 1), a Boltzmann distribution of the energy -ln f at the temperature 1 / (N - 1). Random fields alone, coupling = 0, give a landscape with one peak. Couplings J_ij of variance coupling^2 / loci give the rugged landscape of a spin glass, where a small population escapes local peaks.

Parameters:

Name	Type	Description	Default
loci	int	The number of loci of the genome.	96
N	float	The size of the population.	100.0
mu	float	The rate at which one locus of an individual mutates.	0.01
field	float	The standard deviation of the random fields h_i.	0.01
coupling	float	The scale of the random couplings J_ij, or 0 for none.	0.0
F0	float	The fitness of a genome without fields and couplings, above the sum of the |h_i| and |J_ij| so that every fitness is positive.	1.0
random_state	int | Generator | None	The seed or generator of the fields and couplings.	0

Returns:

Name	Type	Description
system	System	The system, with the alleles genome.s on a one-dimensional lattice of loci, starting from the genome of all -1, and the fitness fitness of the genome.

Examples:

>>> from skcomplex.systems.biology.evolution import origin_fixation_walk
>>> origin_fixation_walk(8).elements[0].prepare().size_
8

Source code in src/skcomplex/systems/biology/evolution/_substitutions.py

@register_system(
    field='biology/evolution',
    tags=('population genetics', 'fitness landscape', 'jump process', 'fixation', 'strong selection weak mutation'),
    references=('Bertin 2021', 'Sella and Hirsh 2005', 'Drossel 2001'),
    defaults={'until': 640.0, 'dt': 10.0},
    observables=_walk_observables,
)
def origin_fixation_walk(
    loci: int = 96,
    *,
    N: float = 100.0,
    mu: float = 0.01,
    field: float = 0.01,
    coupling: float = 0.0,
    F0: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a monomorphic population that climbs a fitness landscape as mutants arise and fix one at a time.

    Under strong selection and weak mutation the population carries one genome of `loci` alleles `-1` or `1`, and
    every substitution of a mutant moves it to a neighbouring genome. The walk climbs the landscape and settles on the
    stationary distribution `Q ~ f^(N - 1)`, a Boltzmann distribution of the energy `-ln f` at the temperature
    `1 / (N - 1)`. Random fields alone, `coupling = 0`, give a landscape with one peak. Couplings `J_ij` of variance
    `coupling^2 / loci` give the rugged landscape of a spin glass, where a small population escapes local peaks.

    Args:
        loci:
            The number of loci of the genome.

        N:
            The size of the population.

        mu:
            The rate at which one locus of an individual mutates.

        field:
            The standard deviation of the random fields `h_i`.

        coupling:
            The scale of the random couplings `J_ij`, or 0 for none.

        F0:
            The fitness of a genome without fields and couplings, above the sum of the `|h_i|` and `|J_ij|` so that
            every fitness is positive.

        random_state:
            The seed or generator of the fields and couplings.

    Returns:
        system:
            The system, with the alleles `genome.s` on a one-dimensional lattice of loci, starting from the genome of
            all `-1`, and the fitness `fitness` of the genome.

    Examples:
        >>> from skcomplex.systems.biology.evolution import origin_fixation_walk
        >>> origin_fixation_walk(8).elements[0].prepare().size_
        8
    """
    rng = np.random.default_rng(random_state)
    fields = rng.normal(0.0, field, loci)
    upper = np.triu(rng.normal(0.0, coupling / np.sqrt(loci), (loci, loci)), 1)
    matrix = upper + upper.T
    couplings = matrix if coupling else None
    genome = LatticeElement(
        'genome',
        element=LocusElement,
        space=LatticeSpace(loci),
        values={'s': -np.ones(loci), 'h': fields},
    )
    return System(
        elements=[genome],
        interactions=[SubstitutionInteraction('substitution', target='genome', N=N, mu=mu, F0=F0, couplings=couplings)],
    )

polymer_assembly_pool(mode='undirected', *, capacity=2000)

Return an assembly pool of linear polymers grown from one monomer by joining pairs of its polymers.

Undirected exploration joins two random polymers of the pool, and directed exploration, a model of selection, joins the newest polymer with a random one. Undirected joining fills the pool with many short polymers, while directed joining builds ever longer chains, of a higher assembly index, about the base 2 logarithm of the length.

Parameters:

Name	Type	Description	Default
mode	str	'undirected' or 'directed'.	'undirected'
capacity	int	The largest number of polymers the pool holds.	2000

Returns:

Name	Type	Description
system	System	The system, with the lengths pool.length of the variable group pool.

Examples:

>>> from skcomplex.systems.biology.evolution import polymer_assembly_pool
>>> polymer_assembly_pool('directed').interactions[0].mode
'directed'

Source code in src/skcomplex/systems/biology/evolution/_assembly.py

@register_system(
    field='biology/evolution',
    tags=('assembly theory', 'selection', 'polymers', 'agents'),
    references=('Sharma et al. 2023',),
    defaults={'until': 200},
)
def polymer_assembly_pool(mode: str = 'undirected', *, capacity: int = 2000) -> System:
    """Return an assembly pool of linear polymers grown from one monomer by joining pairs of its polymers.

    Undirected exploration joins two random polymers of the pool, and directed exploration, a model of selection,
    joins the newest polymer with a random one. Undirected joining fills the pool with many short polymers, while
    directed joining builds ever longer chains, of a higher assembly index, about the base 2 logarithm of the length.

    Args:
        mode:
            `'undirected'` or `'directed'`.

        capacity:
            The largest number of polymers the pool holds.

    Returns:
        system:
            The system, with the lengths `pool.length` of the variable group `pool`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import polymer_assembly_pool
        >>> polymer_assembly_pool('directed').interactions[0].mode
        'directed'
    """
    pool = GroupElement('pool', element=PolymerElement, size=1, variable=True, capacity=capacity)
    return System(
        elements=[pool],
        interactions=[PolymerJoiningInteraction('joining', polymers='pool', mode=mode)],
        time=DiscreteTime(),
    )

read_base_pair_distance(first, second)

Return the number of base pairs that one of two secondary structures of one length holds and the other does not.

Parameters:

Name	Type	Description	Default
first	str	The first structure in dot-bracket notation.	required
second	str	The second structure in dot-bracket notation.	required

Returns:

Name	Type	Description
distance	int	The base-pair distance.

Raises:

Type	Description
MissingExtraError	The extra rna is not installed.

Examples:

>>> from skcomplex.systems.biology.evolution import read_base_pair_distance
>>> read_base_pair_distance('((((...))))', '.(((...))).')
1

Source code in src/skcomplex/systems/biology/evolution/_rna_walks.py

def read_base_pair_distance(first: str, second: str) -> int:
    """Return the number of base pairs that one of two secondary structures of one length holds and the other does not.

    Args:
        first:
            The first structure in dot-bracket notation.

        second:
            The second structure in dot-bracket notation.

    Returns:
        distance:
            The base-pair distance.

    Raises:
        MissingExtraError:
            The extra `rna` is not installed.

    Examples:
        >>> from skcomplex.systems.biology.evolution import read_base_pair_distance
        >>> read_base_pair_distance('((((...))))', '.(((...))).')
        1
    """
    return int(import_extra('RNA', SYSTEM, EXTRA).bp_distance(first, second))

read_plane_layout(structure)

Return the positions in a plane of the bases of a secondary structure, in the simple layout of ViennaRNA.

Parameters:

Name	Type	Description	Default
structure	str	The structure in dot-bracket notation.	required

Returns:

Name	Type	Description
positions	NDArray[float64]	One row of two coordinates per base.

Raises:

Type	Description
MissingExtraError	The extra rna is not installed.

Examples:

>>> from skcomplex.systems.biology.evolution import read_plane_layout
>>> read_plane_layout('((((...))))').shape
(11, 2)

Source code in src/skcomplex/systems/biology/evolution/_rna_walks.py

def read_plane_layout(structure: str) -> npt.NDArray[np.float64]:
    """Return the positions in a plane of the bases of a secondary structure, in the simple layout of ViennaRNA.

    Args:
        structure:
            The structure in dot-bracket notation.

    Returns:
        positions:
            One row of two coordinates per base.

    Raises:
        MissingExtraError:
            The extra `rna` is not installed.

    Examples:
        >>> from skcomplex.systems.biology.evolution import read_plane_layout
        >>> read_plane_layout('((((...))))').shape
        (11, 2)
    """
    layout = import_extra('RNA', SYSTEM, EXTRA).simple_xy_coordinates(structure)
    return np.array([[point.X, point.Y] for point in layout[: len(structure)]], dtype=float)

read_sequence(row)

Return the sequence of a row of indices of bases.

Parameters:

Name	Type	Description	Default
row	ndarray	The indices of the bases.	required

Returns:

Name	Type	Description
sequence	str	The sequence over the letters ACGU.

Examples:

>>> import numpy as np
>>> from skcomplex.systems.biology.evolution import read_sequence
>>> read_sequence(np.array([0.0, 3.0, 2.0]))
'AUG'

Source code in src/skcomplex/systems/biology/evolution/_rna_walks.py

def read_sequence(row: np.ndarray) -> str:
    """Return the sequence of a row of indices of bases.

    Args:
        row:
            The indices of the bases.

    Returns:
        sequence:
            The sequence over the letters `ACGU`.

    Examples:
        >>> import numpy as np
        >>> from skcomplex.systems.biology.evolution import read_sequence
        >>> read_sequence(np.array([0.0, 3.0, 2.0]))
        'AUG'
    """
    return ''.join(BASES[int(value)] for value in row)

replicator_dynamics(x=None, *, payoff=None, V=2.0, C=3.0)

Return the replicator dynamics of a game, by default hawks and doves, which settle at a share V / C of hawks.

Every strategy is a member of the group population, and it grows with its payoff above the mean payoff. The hawks and doves game has the payoff matrix ((V - C) / 2, V), (0, V / 2) with hawks first. When the cost of a fight C exceeds the value V, hawks and doves coexist with V / C hawks, and otherwise hawks take over. A payoff matrix gives any other game, such as the prisoner's dilemma or rock, paper, scissors.

Parameters:

Name	Type	Description	Default
x	Sequence[float] | None	The initial shares of the strategies, or None for equal shares.	None
payoff	Sequence[Sequence[float]] | None	The payoff matrix of another game, one row per strategy, or None for hawks and doves.	None
V	float	The value of the resource contested by hawks and doves.	2.0
C	float	The cost of a fight between hawks.	3.0

Returns:

Name	Type	Description
system	System	The system, with the shares population.x of the strategies, hawks first.

Examples:

>>> from skcomplex.systems.biology.evolution import replicator_dynamics
>>> replicator_dynamics(payoff=((0, -1, 1), (1, 0, -1), (-1, 1, 0))).elements[0].prepare().size_
3

Source code in src/skcomplex/systems/biology/evolution/_games.py

@register_system(
    field='biology/evolution',
    tags=('evolutionary game', 'replicator', 'ode', 'game theory'),
    references=('Taylor and Jonker 1978', 'Maynard Smith 1982', 'Gros 2015', 'Szabó and Fáth 2007'),
    defaults={'until': 30.0, 'dt': 0.1},
)
def replicator_dynamics(
    x: Sequence[float] | None = None,
    *,
    payoff: Sequence[Sequence[float]] | None = None,
    V: float = 2.0,
    C: float = 3.0,
) -> System:
    """Return the replicator dynamics of a game, by default hawks and doves, which settle at a share `V / C` of hawks.

    Every strategy is a member of the group `population`, and it grows with its payoff above the mean payoff. The
    hawks and doves game has the payoff matrix `((V - C) / 2, V), (0, V / 2)` with hawks first. When the cost of a fight
    `C` exceeds the value `V`, hawks and doves coexist with `V / C` hawks, and otherwise hawks take over. A payoff
    matrix gives any other game, such as the prisoner's dilemma or rock, paper, scissors.

    Args:
        x:
            The initial shares of the strategies, or `None` for equal shares.

        payoff:
            The payoff matrix of another game, one row per strategy, or `None` for hawks and doves.

        V:
            The value of the resource contested by hawks and doves.

        C:
            The cost of a fight between hawks.

    Returns:
        system:
            The system, with the shares `population.x` of the strategies, hawks first.

    Examples:
        >>> from skcomplex.systems.biology.evolution import replicator_dynamics
        >>> replicator_dynamics(payoff=((0, -1, 1), (1, 0, -1), (-1, 1, 0))).elements[0].prepare().size_
        3
    """
    matrix = (((V - C) / 2, V), (0.0, V / 2)) if payoff is None else tuple(tuple(row) for row in payoff)
    strategies = len(matrix)
    shares = [1 / strategies] * strategies if x is None else list(x)
    population = GroupElement('population', element=StrategyElement, size=strategies, values={'x': shares})
    weights = np.asarray(matrix, dtype=float)
    replication = ReplicatorInteraction(
        'replication',
        target='population',
        fitness=lambda s: total(s.target.x, weights=weights),
    )
    return System(elements=[population], interactions=[replication])

rna_neutral_networks(walkers=8, length=30, *, start=START, target=TARGET, regime='neutral', seed=0)

Return walkers of RNA sequences that take one proposed point mutation per unit of time at temperature zero.

In the 'neutral' regime a mutation is accepted when the minimum free energy fold of the sequence stays on the target. In the 'adaptive' regime a mutation is accepted when the base-pair distance between the fold and the target does not grow. The neighbours of the visited sequences give the robustness and evolvability of Raman (2021 §13.3) through NeutralNetworkAnalysis. The folds come from ViennaRNA at 37 °C with its default energy parameters. Every walker records its score as walk.score.

Parameters:

Name	Type	Description	Default
walkers	int	The number of walkers, the rows of the lattice.	8
length	int	The length of the sequences, the columns of the lattice.	30
start	str | None	The sequence every walker starts on, or None for random sequences.	START
target	str	The target structure in dot-bracket notation, of the length of the sequences.	TARGET
regime	str	'neutral' or 'adaptive'.	'neutral'
seed	int	The seed of the random starting sequences.	0

Returns:

Name	Type	Description
system	System	The system, with the lattice walkers of the state base.

Raises:

Type	Description
MissingExtraError	The extra rna is not installed.
ValidationError	The regime is unknown, or the start or the target does not have the length of the sequences.

Examples:

>>> from skcomplex.systems.biology.evolution import rna_neutral_networks
>>> rna_neutral_networks(walkers=2).elements[0].space.shape
(2, 30)

Source code in src/skcomplex/systems/biology/evolution/_rna_walks.py

@register_system(
    field='biology/evolution',
    tags=('rna', 'genotype networks', 'neutral networks', 'robustness', 'evolvability', 'monte carlo'),
    references=('Raman 2021', 'Wagner 2008', 'Lorenz et al. 2011'),
    defaults={'until': 300},
)
def rna_neutral_networks(
    walkers: int = 8,
    length: int = 30,
    *,
    start: str | None = START,
    target: str = TARGET,
    regime: str = 'neutral',
    seed: int = 0,
) -> System:
    """Return walkers of RNA sequences that take one proposed point mutation per unit of time at temperature zero.

    In the `'neutral'` regime a mutation is accepted when the minimum free energy fold of the sequence stays on the
    target. In the `'adaptive'` regime a mutation is accepted when the base-pair distance between the fold and the
    target does not grow. The neighbours of the visited sequences give the robustness and evolvability of
    Raman (2021 §13.3) through `NeutralNetworkAnalysis`. The folds come from ViennaRNA at 37 °C with its default energy
    parameters. Every walker records its score as `walk.score`.

    Args:
        walkers:
            The number of walkers, the rows of the lattice.

        length:
            The length of the sequences, the columns of the lattice.

        start:
            The sequence every walker starts on, or `None` for random sequences.

        target:
            The target structure in dot-bracket notation, of the length of the sequences.

        regime:
            `'neutral'` or `'adaptive'`.

        seed:
            The seed of the random starting sequences.

    Returns:
        system:
            The system, with the lattice `walkers` of the state `base`.

    Raises:
        MissingExtraError:
            The extra `rna` is not installed.

        ValidationError:
            The regime is unknown, or the start or the target does not have the length of the sequences.

    Examples:
        >>> from skcomplex.systems.biology.evolution import rna_neutral_networks
        >>> rna_neutral_networks(walkers=2).elements[0].space.shape
        (2, 30)
    """
    import_extra('RNA', 'rna_neutral_networks', 'rna')
    if regime not in REGIMES:
        error_msg = f'Unknown regime {regime!r}. Use one of {", ".join(REGIMES)}.'
        raise ValidationError(error_msg)
    if len(target) != length or (start is not None and len(start) != length):
        error_msg = f'The start {start!r} and the target {target!r} must have the length {length} of the sequences.'
        raise ValidationError(error_msg)
    if start is None:
        bases = np.random.default_rng(seed).integers(0, len(BASES), size=(walkers, length))
    else:
        bases = np.tile([BASES.index(letter) for letter in start], (walkers, 1))
    return System(
        elements=[
            LatticeElement(
                'walkers',
                element=NucleotideElement,
                space=LatticeSpace((walkers, length), boundary='fixed'),
                values={'base': bases},
                visual=VisualSpec('image', {'low': 0.0, 'high': float(LAST_BASE)}),
            ),
        ],
        interactions=[
            SamplingInteraction(
                'walk',
                target='walkers',
                state='base',
                objective=RowScoreSpec(partial(_score_row, target, regime)),
                move=mutate_point,
                intensity=float('inf'),
            ),
        ],
        time=DiscreteTime(),
    )

rock_paper_scissors(shape=(60, 60), *, random_state=0)

Return rock, paper, and scissors on a lattice, whose cyclic invasions keep the three strategies coexisting.

Starting from random strategies, the strategies chase each other in domains and fronts, and each keeps about a third of the lattice.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns of the periodic lattice.	(60, 60)
random_state	int | Generator | None	The seed or generator of the initial strategies.	0

Returns:

Name	Type	Description
system	System	The system, with the strategies lattice.strategy.

Examples:

>>> from skcomplex.systems.biology.evolution import rock_paper_scissors
>>> rock_paper_scissors((10, 10)).elements[0].prepare().size_
100

Source code in src/skcomplex/systems/biology/evolution/_games.py

@register_system(
    field='biology/evolution',
    tags=('evolutionary game', 'lattice', 'cyclic dominance', 'spirals', 'stochastic'),
    references=('Szabó and Fáth 2007', 'Reichenbach, Mobilia, and Frey 2007'),
    defaults={'until': 100, 'dt': 1},
)
def rock_paper_scissors(
    shape: Sequence[int] = (60, 60),
    *,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return rock, paper, and scissors on a lattice, whose cyclic invasions keep the three strategies coexisting.

    Starting from random strategies, the strategies chase each other in domains and fronts, and each keeps about a
    third of the lattice.

    Args:
        shape:
            The number of rows and columns of the periodic lattice.

        random_state:
            The seed or generator of the initial strategies.

    Returns:
        system:
            The system, with the strategies `lattice.strategy`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import rock_paper_scissors
        >>> rock_paper_scissors((10, 10)).elements[0].prepare().size_
        100
    """
    strategies = np.random.default_rng(random_state).integers(0, STRATEGIES, size=tuple(shape)).astype(float)
    return System(
        elements=[
            LatticeElement(
                'lattice',
                element=CyclicSpeciesElement,
                space=LatticeSpace(tuple(shape)),
                values={'strategy': strategies},
            ),
        ],
        interactions=[
            LocalRuleInteraction(
                'invasion',
                target='lattice',
                state='strategy',
                rule=_invade,
                neighborhood=NEAREST,
                schedule='random-sequential',
                stochastic=True,
            ),
        ],
    )

sole_manrubia(n=100, k=5, *, noise=0.05, seed=1)

Return the ecosystem of Solé and Manrubia, whose species go extinct when their dependencies turn against them.

Every species depends on k species through couplings drawn from [-1, 1], and its field is the sum of these couplings. The couplings drift at random. A species whose field turns negative goes extinct, the species that depend on it lose a coupling, and their fields may turn negative in turn, so extinctions come in avalanches with a broad distribution of sizes. A new species with k random dependencies replaces every extinct one.

Parameters:

Name	Type	Description	Default
n	int	The number of species.	100
k	int	The number of dependencies of every species.	5
noise	float	The standard deviation of the drift of every coupling at a step.	0.05
seed	int | None	The seed of the network and of the initial couplings.	1

Returns:

Name	Type	Description
system	System	The system, with the fields species.h and ages species.age of the variable group species and the couplings couplings.w of the variable group couplings on the edges.

Examples:

>>> from skcomplex.systems.biology.evolution import sole_manrubia
>>> [component.label for component in sole_manrubia(n=20).elements]
['species', 'couplings']

Source code in src/skcomplex/systems/biology/evolution/_networks.py

@register_system(
    field='biology/evolution',
    tags=('adaptive network', 'extinction', 'avalanches', 'coevolution', 'stochastic map'),
    references=('Solé and Manrubia 1996', 'Thurner et al. 2018'),
    defaults={'until': 300, 'dt': 1.0},
)
def sole_manrubia(
    n: int = 100,
    k: int = 5,
    *,
    noise: float = 0.05,
    seed: int | None = 1,
) -> System:
    """Return the ecosystem of Solé and Manrubia, whose species go extinct when their dependencies turn against them.

    Every species depends on `k` species through couplings drawn from `[-1, 1]`, and its field is the sum of these
    couplings. The couplings drift at random. A species whose field turns negative goes extinct, the species that
    depend on it lose a coupling, and their fields may turn negative in turn, so extinctions come in avalanches with a
    broad distribution of sizes. A new species with `k` random dependencies replaces every extinct one.

    Args:
        n:
            The number of species.

        k:
            The number of dependencies of every species.

        noise:
            The standard deviation of the drift of every coupling at a step.

        seed:
            The seed of the network and of the initial couplings.

    Returns:
        system:
            The system, with the fields `species.h` and ages `species.age` of the variable group `species` and the
            couplings `couplings.w` of the variable group `couplings` on the edges.

    Examples:
        >>> from skcomplex.systems.biology.evolution import sole_manrubia
        >>> [component.label for component in sole_manrubia(n=20).elements]
        ['species', 'couplings']
    """
    rng = np.random.default_rng(seed)
    graph = nx.DiGraph()
    graph.add_nodes_from(range(n))
    for node in range(n):
        others = rng.choice(np.delete(np.arange(n), node), size=k, replace=False)
        graph.add_edges_from((node, int(other)) for other in others)
    space = NetworkSpace(graph).prepare()
    edges = graph.number_of_edges()
    couplings = rng.uniform(-1, 1, edges)
    sources, _ = space.edges_[None]
    fields = np.bincount(sources, weights=couplings, minlength=n)
    species = GroupElement(
        'species',
        element=TaxonElement,
        space=space,
        on='sites',
        variable=True,
        capacity=2 * n,
        values={'h': fields},
        visual=VisualSpec('sites', {'states': ['age']}),
    )
    edge_group = GroupElement(
        'couplings',
        element=WeightedEdgeElement,
        space=space,
        on='edges',
        variable=True,
        capacity=4 * edges,
        values={'w': couplings, 'fresh': 0.0},
        visual=HIDDEN,
    )
    return System(
        elements=[species, edge_group],
        interactions=[
            DependencyFieldInteraction('field', species='species', couplings='couplings'),
            EdgeDrawInteraction('drift', target='species', edges='couplings', noise=noise),
            MapInteraction('ageing', target='species', state='age', function=lambda s: 1, mode='increment'),
            TopologyInteraction(
                'extinction',
                target='species',
                action='replace-node',
                condition=lambda s: s.source.h < 0,
                edges=k,
            ),
        ],
        time=DiscreteTime(),
    )

tangled_nature(organisms=100, *, k=5.0, mu=0.01, p_mut=0.01, p_kill=0.2, length=8, seed=0, random_state=0)

Return the Tangled Nature system, where organisms whose types interact through random couplings coevolve.

Organisms of one random type start the run. Mutants spread through the space of genomes, and the population settles into quasi-stable configurations of a few dominant types, the q-ESS, that end in abrupt transitions, a punctuated equilibrium. The trait of an organism holds the bits of its genome. A generation is N / p_kill visits to organisms drawn with replacement, where N is the population at its start.

Parameters:

Name	Type	Description	Default
organisms	int	The initial number of organisms.	100
k	float	The strength of the couplings.	5.0
mu	float	The limit of resources, which lowers reproduction in large populations.	0.01
p_mut	float	The probability that a bit of an offspring's genome flips.	0.01
p_kill	float	The probability that a chosen organism dies.	0.2
length	int	The number of bits of a genome.	8
seed	int	The seed of the couplings.	0
random_state	int | Generator | None	The seed or generator of the initial type.	0

Returns:

Name	Type	Description
system	System	The system, with the group organisms.

Examples:

>>> from skcomplex.systems.biology.evolution import tangled_nature
>>> tangled_nature(20).elements[0].prepare().size_
20

Source code in src/skcomplex/systems/biology/evolution/_finite_populations.py

@register_system(
    field='biology/evolution',
    tags=('coevolution', 'agents', 'punctuated equilibrium', 'ecosystems'),
    references=('Christensen et al. 2002', 'Jensen 2023'),
    defaults={'until': 50},
)
def tangled_nature(
    organisms: int = 100,
    *,
    k: float = 5.0,
    mu: float = 0.01,
    p_mut: float = 0.01,
    p_kill: float = 0.2,
    length: int = 8,
    seed: int = 0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Tangled Nature system, where organisms whose types interact through random couplings coevolve.

    Organisms of one random type start the run. Mutants spread through the space of genomes, and the population
    settles into quasi-stable configurations of a few dominant types, the q-ESS, that end in abrupt transitions, a
    punctuated equilibrium. The trait of an organism holds the bits of its genome. A generation is `N / p_kill`
    visits to organisms drawn with replacement, where `N` is the population at its start.

    Args:
        organisms:
            The initial number of organisms.

        k:
            The strength of the couplings.

        mu:
            The limit of resources, which lowers reproduction in large populations.

        p_mut:
            The probability that a bit of an offspring's genome flips.

        p_kill:
            The probability that a chosen organism dies.

        length:
            The number of bits of a genome.

        seed:
            The seed of the couplings.

        random_state:
            The seed or generator of the initial type.

    Returns:
        system:
            The system, with the group `organisms`.

    Examples:
        >>> from skcomplex.systems.biology.evolution import tangled_nature
        >>> tangled_nature(20).elements[0].prepare().size_
        20
    """
    genome = float(np.random.default_rng(random_state).integers(2**length))
    life_cycle = TangledNatureLifeCycleInteraction(
        'evolution',
        target=ORGANISMS,
        k=k,
        mu=mu,
        p_mut=p_mut,
        p_kill=p_kill,
        length=length,
        seed=seed,
    )
    visual = VisualSpec('histograms', {'color_by': 'trait', 'shown': ('distinct',)})
    return System(
        elements=[
            GroupElement(
                ORGANISMS,
                element=IndividualElement,
                size=organisms,
                values={'trait': genome},
                variable=True,
                visual=visual,
            ),
        ],
        interactions=[life_cycle],
        time=DiscreteTime(
            schedule=PhasesScheduleSpec([PhaseSpec('generation', activation='visits', visits=1 / p_kill)]),
        ),
    )