Systems biology

Systems biology: gene circuits, Boolean networks, circadian clocks, signalling cascades, and Min-protein waves.

BooleanFunctionInteraction(label, *, target, inputs, table=(0,))

Bases: MapInteraction

Synchronous update of a node by a truth table of any number of input nodes.

A table of 2^k entries reads the k inputs as the bits of a binary number, from the most significant, and gives the next value of the target. The update is the sum, over the patterns whose entry is 1, of the product that is 1 exactly on that pattern.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the updated node.	required
inputs	Annotated[list[str], RoleSpec(variadic=True)]	The paths of the input nodes, in the order of the bits of the table.	required
table	tuple[int, ...]	The next value of the target for each of the 2^k input patterns.	(0,)

Raises:

Type	Description
ValidationError	The table does not have one entry per pattern of the inputs.

Examples:

>>> from skcomplex.models.biology.systems_biology import BooleanFunctionInteraction
>>> BooleanFunctionInteraction('and', target='c', inputs=['a', 'b'], table=(0, 0, 0, 1)).table
(0, 0, 0, 1)

Source code in src/skcomplex/models/biology/systems_biology/_boolean_function.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    inputs: Annotated[list[str], RoleSpec(variadic=True)],
    table: tuple[int, ...] = (0,),
) -> None:
    super().__init__(label, target=target, state='x', function=self._get_next)
    self.inputs = inputs
    self.table = table

BooleanNodeElement(label, visual=None, *, x=0)

Bases: Element

Node of a Boolean network, a gene or protein that is on or off.

Parameters:

Name	Type	Description	Default
label	str	The label of the node.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial value, 0 or 1.	0

Examples:

>>> from skcomplex.models.biology.systems_biology import BooleanNodeElement
>>> BooleanNodeElement('Sic1', x=1).x
1

Source code in src/skcomplex/models/biology/systems_biology/_boolean_node.py

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

GProteinElement(label, visual=None, *, g=0.0, k_on=1.0, k_gtp=0.5, hormone=1.0, start=1.0, width=1.0)

Bases: Element

G protein activated by a hormone pulse and inactivated by its own GTPase, the timer that ends the signal.

The active fraction follows dg/dt = k_on H(t) (1 - g) - k_gtp g, where the hormone H(t) is a rectangular pulse of height hormone from start for width. Once the hormone is gone, the signal ends after the time 1 / k_gtp.

The element holds the states and parameters and has no laws of its own: the interactions of hormone_cascade_amplifier read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the G protein.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
g	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial active fraction.	0.0
k_on	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The activation rate per unit of hormone.	1.0
k_gtp	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of GTP hydrolysis, which inactivates the protein.	0.5
hormone	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The height of the hormone pulse.	1.0
start	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The time the hormone pulse starts.	1.0
width	Annotated[float, ParameterSpec(domain=open(0, oo))]	The duration of the hormone pulse.	1.0

Examples:

>>> from skcomplex.models.biology.systems_biology import GProteinElement
>>> GProteinElement('g_protein').k_gtp
0.5

Source code in src/skcomplex/models/biology/systems_biology/_cascades.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    g: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    k_on: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    k_gtp: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    hormone: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    start: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    width: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.g = g
    self.k_on = k_on
    self.k_gtp = k_gtp
    self.hormone = hormone
    self.start = start
    self.width = width

GeneElement(label, visual=None, *, m=0.0, p=0.0)

Bases: Element

Gene with a count m of messenger RNA and a count p of protein.

A gene has no laws of its own. Reactions transcribe, translate, and degrade its molecules.

Parameters:

Name	Type	Description	Default
label	str	The label of the gene.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
m	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial number of mRNA molecules.	0.0
p	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial number of protein molecules.	0.0

Examples:

>>> from skcomplex.models.biology.systems_biology import GeneElement
>>> GeneElement('lacI', p=20.0).p
20.0

Source code in src/skcomplex/models/biology/systems_biology/_genes.py

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

GeneralizedMassActionInteraction(label, *, species, stoichiometry=((0.0,),), gamma=(0.0,), orders=((0.0,),))

Bases: Interaction

Rates of a biochemical system in the generalized mass action form of biochemical systems theory.

Flux j is the power law v_j = gamma_j prod_k x_k^f_jk, and species i changes as dx_i/dt = sum_j s_ij v_j. The kinetic orders f_jk are real numbers, negative for inhibition.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
species	Annotated[list[str], RoleSpec(variadic=True)]	The paths of the species, in the order of the rows of the stoichiometry.	required
stoichiometry	Sequence[Sequence[float]]	The stoichiometric coefficient s_ij of every species in every flux, one row per species.	((0.0,),)
gamma	Sequence[float]	The rate constant of every flux.	(0.0,)
orders	Sequence[Sequence[float]]	The kinetic order f_jk of every species in every flux, one row per flux.	((0.0,),)

Raises:

Type	Description
ValidationError	The matrices do not have one row per flux or species and one entry per species or flux.

Examples:

>>> from skcomplex.models.biology.systems_biology import GeneralizedMassActionInteraction
>>> v = GeneralizedMassActionInteraction(
...     'v', species=['X'], stoichiometry=((1, -1),), gamma=(1, 2), orders=((0,), (1,))
... )
>>> v.gamma
(1, 2)

Source code in src/skcomplex/models/biology/systems_biology/_generalized_mass_action.py

def __init__(
    self,
    label: str,
    *,
    species: Annotated[list[str], RoleSpec(variadic=True)],
    stoichiometry: Sequence[Sequence[float]] = ((0.0,),),
    gamma: Sequence[float] = (0.0,),
    orders: Sequence[Sequence[float]] = ((0.0,),),
) -> None:
    super().__init__(label)
    self.species = species
    self.stoichiometry = stoichiometry
    self.gamma = gamma
    self.orders = orders

MinSkeletonElement(label, visual=None, *, md=0.0, mde=0.0, cdt=0.0, cdd=0.0, ce=0.0, Dm=0.013, Dc=60.0, k_d=0.1, k_dd=0.1, k_de=0.5, k_off=0.4, lam=6.0)

Bases: Element

Site of the Min skeleton model, with membrane MinD and MinDE and cytosolic MinD-ATP, MinD-ADP, and MinE.

The reactions are the attachment and self-recruitment of MinD-ATP to the membrane, (k_d + k_dd md) cdt, the recruitment of MinE by membrane MinD, k_de md ce, the dissociation of MinDE complexes, k_off mde, which releases MinD-ADP and MinE to the cytosol, and the nucleotide exchange lam cdd from MinD-ADP to MinD-ATP. They conserve the totals of MinD, md + mde + cdt + cdd, and of MinE, mde + ce.

The element holds the states and parameters and has no laws of its own: the interactions of min_protein_skeleton read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
md	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial membrane MinD.	0.0
mde	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial membrane MinDE complexes.	0.0
cdt	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial cytosolic MinD-ATP.	0.0
cdd	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial cytosolic MinD-ADP.	0.0
ce	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial cytosolic MinE.	0.0
Dm	Annotated[float, ParameterSpec(domain=open(0, oo))]	The diffusivity on the membrane.	0.013
Dc	Annotated[float, ParameterSpec(domain=open(0, oo))]	The diffusivity in the cytosol.	60.0
k_d	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The attachment rate of MinD-ATP.	0.1
k_dd	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The recruitment rate of MinD-ATP by membrane MinD.	0.1
k_de	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The recruitment rate of MinE by membrane MinD.	0.5
k_off	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The dissociation rate of MinDE complexes.	0.4
lam	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The nucleotide exchange rate.	6.0

Examples:

>>> from skcomplex.models.biology.systems_biology import MinSkeletonElement
>>> MinSkeletonElement('site').lam
6.0

Source code in src/skcomplex/models/biology/systems_biology/_min_proteins.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    md: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    mde: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    cdt: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    cdd: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    ce: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    Dm: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.013,
    Dc: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 60.0,
    k_d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    k_dd: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    k_de: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    k_off: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.4,
    lam: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 6.0,
) -> None:
    super().__init__(label, visual=visual)
    self.md = md
    self.mde = mde
    self.cdt = cdt
    self.cdd = cdd
    self.ce = ce
    self.Dm = Dm
    self.Dc = Dc
    self.k_d = k_d
    self.k_dd = k_dd
    self.k_de = k_de
    self.k_off = k_off
    self.lam = lam

PerClockElement(label, visual=None, *, M=0.5, P0=0.5, P1=0.5, P2=0.5, PN=0.5, v_s=0.76, v_m=0.65, K_m=0.5, k_s=0.38, v_d=0.95, K_d=0.2, k1=1.9, k2=1.3, K_I=1.0, n=4.0, V1=3.2, V2=1.58, V3=5.0, V4=2.5, K=2.0)

Bases: Element

PER protein clock of Drosophila, after Goldbeter: per mRNA, three phosphorylation forms, and nuclear PER.

The mRNA M is transcribed under repression by nuclear PER and degraded. Cytosolic PER is translated as P0, phosphorylated reversibly to P1 and P2, and P2 is degraded or enters the nucleus as PN. Concentrations are in µM and time in hours.

The element holds the states and parameters and has no laws of its own: the interactions of circadian_clock read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the clock.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
M	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial per mRNA concentration.	0.5
P0	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial unphosphorylated PER.	0.5
P1	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial monophosphorylated PER.	0.5
P2	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial bisphosphorylated PER.	0.5
PN	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial nuclear PER.	0.5
v_s	Annotated[float, ParameterSpec(domain=open(0, oo))]	The maximal transcription rate.	0.76
v_m	Annotated[float, ParameterSpec(domain=open(0, oo))]	The maximal mRNA degradation rate.	0.65
K_m	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant of mRNA degradation.	0.5
k_s	Annotated[float, ParameterSpec(domain=open(0, oo))]	The translation rate.	0.38
v_d	Annotated[float, ParameterSpec(domain=open(0, oo))]	The maximal degradation rate of P2.	0.95
K_d	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant of the degradation of P2.	0.2
k1	Annotated[float, ParameterSpec(domain=open(0, oo))]	The rate of nuclear entry.	1.9
k2	Annotated[float, ParameterSpec(domain=open(0, oo))]	The rate of nuclear exit.	1.3
K_I	Annotated[float, ParameterSpec(domain=open(0, oo))]	The threshold of repression by nuclear PER.	1.0
n	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Hill coefficient of repression.	4.0
V1	Annotated[float, ParameterSpec(domain=open(0, oo))]	The maximal rate of the first phosphorylation.	3.2
V2	Annotated[float, ParameterSpec(domain=open(0, oo))]	The maximal rate of the first dephosphorylation.	1.58
V3	Annotated[float, ParameterSpec(domain=open(0, oo))]	The maximal rate of the second phosphorylation.	5.0
V4	Annotated[float, ParameterSpec(domain=open(0, oo))]	The maximal rate of the second dephosphorylation.	2.5
K	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant of the phosphorylations and dephosphorylations.	2.0

Examples:

>>> from skcomplex.models.biology.systems_biology import PerClockElement
>>> PerClockElement('clock').v_s
0.76

Source code in src/skcomplex/models/biology/systems_biology/_clocks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    M: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    P0: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    P1: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    P2: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    PN: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    v_s: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.76,
    v_m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.65,
    K_m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    k_s: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.38,
    v_d: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.95,
    K_d: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.2,
    k1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.9,
    k2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.3,
    K_I: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 4.0,
    V1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 3.2,
    V2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.58,
    V3: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 5.0,
    V4: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.5,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
) -> None:
    super().__init__(label, visual=visual)
    self.M = M
    self.P0 = P0
    self.P1 = P1
    self.P2 = P2
    self.PN = PN
    self.v_s = v_s
    self.v_m = v_m
    self.K_m = K_m
    self.k_s = k_s
    self.v_d = v_d
    self.K_d = K_d
    self.k1 = k1
    self.k2 = k2
    self.K_I = K_I
    self.n = n
    self.V1 = V1
    self.V2 = V2
    self.V3 = V3
    self.V4 = V4
    self.K = K

PhosphorylationLoopElement(label, visual=None, *, A=1.0, B=0.0, C=0.0, kp1=1.0, kd1=2.0, kp2=2.0, kd2=1.0, kp3=2.0, kd3=1.0, K1=0.01, K2=0.01, K3=0.01, A_T=1.0, B_T=1.0, C_T=1.0)

Bases: Element

Three proteins, each phosphorylated and dephosphorylated by Michaelis-Menten kinetics, that regulate each other.

The states are the unphosphorylated amounts, and each total is conserved. With u_T = A_T - A, and so on:

• dA/dt = kp1 u_T / (K1 + u_T) - kd1 C A / (K1 + A),
• dB/dt = kp2 A (B_T - B) / (K2 + B_T - B) - kd2 B / (K2 + B),
• dC/dt = kp3 B (C_T - C) / (K3 + C_T - C) - kd3 C / (K3 + C).

The element holds the states and parameters and has no laws of its own: the interactions of phosphorylation_oscillator read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the loop.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
A	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial amount of unphosphorylated A.	1.0
B	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial amount of unphosphorylated B.	0.0
C	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial amount of unphosphorylated C.	0.0
kp1	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The maximal rate that restores A.	1.0
kd1	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The catalytic rate of the conversion of A by C.	2.0
kp2	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The catalytic rate of the conversion of B by A.	2.0
kd2	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The maximal rate that removes B.	1.0
kp3	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The catalytic rate of the conversion of C by B.	2.0
kd3	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The maximal rate that removes C.	1.0
K1	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant of both reactions of A.	0.01
K2	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant of both reactions of B.	0.01
K3	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant of both reactions of C.	0.01
A_T	Annotated[float, ParameterSpec(domain=open(0, oo))]	The total amount of A.	1.0
B_T	Annotated[float, ParameterSpec(domain=open(0, oo))]	The total amount of B.	1.0
C_T	Annotated[float, ParameterSpec(domain=open(0, oo))]	The total amount of C.	1.0

Examples:

>>> from skcomplex.models.biology.systems_biology import PhosphorylationLoopElement
>>> PhosphorylationLoopElement('loop', A=0.5).A
0.5

Source code in src/skcomplex/models/biology/systems_biology/_phosphorylation_loop.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    A: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    B: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    C: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    kp1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    kd1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    kp2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    kd2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    kp3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    kd3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    K1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
    K2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
    K3: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
    A_T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    B_T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    C_T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.A = A
    self.B = B
    self.C = C
    self.kp1 = kp1
    self.kd1 = kd1
    self.kp2 = kp2
    self.kd2 = kd2
    self.kp3 = kp3
    self.kd3 = kd3
    self.K1 = K1
    self.K2 = K2
    self.K3 = K3
    self.A_T = A_T
    self.B_T = B_T
    self.C_T = C_T

ProbabilisticFunctionInteraction(label, *, target, inputs, tables=((0,),), probabilities=(1.0,))

Bases: MapInteraction

Synchronous map of a node by one of several truth tables, drawn with given probabilities at every step.

Each table of 2^k entries reads the k inputs as the bits of a binary number, from the most significant, and gives the next value x of the target. A second update of the same step records the index of the drawn table in the state function of the target, from the same draw.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the updated node, a ProbabilisticNodeElement.	required
inputs	Annotated[list[str], RoleSpec(variadic=True)]	The paths of the input nodes, in the order of the bits of the tables.	required
tables	Sequence[Sequence[int]]	The candidate truth tables.	((0,),)
probabilities	Sequence[float]	The probability of every table.	(1.0,)

Raises:

Type	Description
ValidationError	The tables and the probabilities differ in number, or a table does not have one entry per pattern.

Examples:

>>> from skcomplex.models.biology.systems_biology import ProbabilisticFunctionInteraction
>>> function = ProbabilisticFunctionInteraction(
...     'f', target='c', inputs=['a', 'b'], tables=((0, 0, 0, 1), (0, 1, 1, 1)), probabilities=(0.7, 0.3)
... )
>>> function.probabilities
(0.7, 0.3)

Source code in src/skcomplex/models/biology/systems_biology/_probabilistic_function.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    inputs: Annotated[list[str], RoleSpec(variadic=True)],
    tables: Sequence[Sequence[int]] = ((0,),),
    probabilities: Sequence[float] = (1.0,),
) -> None:
    super().__init__(label, target=target, state='x', function=self._get_drawn_value)
    self.inputs = inputs
    self.tables = tables
    self.probabilities = probabilities

ProbabilisticNodeElement(label, visual=None, *, x=0, function=0)

Bases: Element

Node of a probabilistic Boolean network, on or off, with the index of the function that set it last.

Parameters:

Name	Type	Description	Default
label	str	The label of the node.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial value, 0 or 1.	0
function	Annotated[float, StateSpec(domain=NN)]	The index of the candidate function of the last update, 0 at the start.	0

Examples:

>>> from skcomplex.models.biology.systems_biology import ProbabilisticNodeElement
>>> ProbabilisticNodeElement('x1', x=1).x
1

Source code in src/skcomplex/models/biology/systems_biology/_probabilistic_node.py

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

RepressionInteraction(label, *, reactants=None, products=None, modifiers=None, alpha=1.0, alpha0=0.0, K=1.0, n=1.0, mode='rate')

Bases: ReactionInteraction

Synthesis of a product at the rate alpha0 + alpha (1 - hill(y, K, n)), repressed by a modifier y.

The product, such as the mRNA of a gene, gains units at the leaky rate alpha0 plus the rate alpha lowered by the decreasing Hill function 1 / (1 + (y / K)^n) of the one modifier, the repressor. The repressor enters the rate and does not change. In the mode 'jump' every synthesis is a random event that adds one unit.

Parameters:

Name	Type	Description	Default
label	str	The label of the repression.	required
reactants	Annotated[Stoichiometry | None, RoleSpec(variadic=True)]	The coefficients of the reactants, or None for a synthesis from a source.	None
products	Annotated[Stoichiometry | None, RoleSpec(variadic=True)]	The coefficient of the synthesized state, such as {('gene', 'm'): 1}.	None
modifiers	Annotated[Stoichiometry | None, RoleSpec(variadic=True)]	The coefficient of the repressor state, such as {('repressor', 'p'): 1}.	None
alpha	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The synthesis rate without repressor.	1.0
alpha0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The leaky synthesis rate at saturating repressor.	0.0
K	Annotated[float, ParameterSpec(domain=open(0, oo))]	The repressor amount at half repression.	1.0
n	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Hill coefficient of the repression.	1.0
mode	str	The form of the laws, 'rate', 'jump', or 'langevin'.	'rate'

Raises:

Type	Description
ValidationError	The mode is unknown, the reaction has neither reactants nor products, or an entry is not a pair of a path and a state name with a positive coefficient.

Examples:

>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.systems_biology import RepressionInteraction, SpeciesElement
>>> repression = RepressionInteraction(
...     'repression', products={('mRNA', 'x'): 1}, modifiers={('repressor', 'x'): 1}, alpha=2.0, K=1.0
... )
>>> species = [SpeciesElement('mRNA'), SpeciesElement('repressor', x=1.0)]
>>> trajectory = Simulation(System(elements=species, interactions=[repression]), dt=0.1).check().prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('mRNA.x')[-1]), 3)
1.0

Source code in src/skcomplex/models/biology/systems_biology/_repression.py

def __init__(
    self,
    label: str,
    *,
    reactants: Annotated[Stoichiometry | None, RoleSpec(variadic=True)] = None,
    products: Annotated[Stoichiometry | None, RoleSpec(variadic=True)] = None,
    modifiers: Annotated[Stoichiometry | None, RoleSpec(variadic=True)] = None,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    alpha0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    mode: str = 'rate',
) -> None:
    super().__init__(label, reactants=reactants, products=products, modifiers=modifiers, mode=mode)
    self.alpha = alpha
    self.alpha0 = alpha0
    self.K = K
    self.n = n

SSystemInteraction(label, *, species, alpha=(), beta=(), g=((0.0,),), h=((0.0,),))

Bases: GeneralizedMassActionInteraction

Rates of a biochemical system in the S-system form of biochemical systems theory.

Every species has one aggregated production and one aggregated degradation, both power laws: dx_i/dt = alpha_i prod_j x_j^g_ij - beta_i prod_j x_j^h_ij. It is the generalized mass action form with two fluxes per species.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
species	Annotated[list[str], RoleSpec(variadic=True)]	The paths of the species.	required
alpha	Annotated[Sequence[float], ParameterSpec(per=species, domain=closed_open(0, oo))]	The rate constant of the production of every species.	()
beta	Annotated[Sequence[float], ParameterSpec(per=species, domain=closed_open(0, oo))]	The rate constant of the degradation of every species.	()
g	Sequence[Sequence[float]]	The kinetic orders of the productions, one row per species.	((0.0,),)
h	Sequence[Sequence[float]]	The kinetic orders of the degradations, one row per species.	((0.0,),)

Raises:

Type	Description
ValidationError	The kinetic orders do not have one row and one column per species.

Examples:

>>> from skcomplex.models.biology.systems_biology import SSystemInteraction
>>> SSystemInteraction('rates', species=['X'], alpha=[1.0], beta=[1.0], g=((0.0,),), h=((1.0,),)).h
((1.0,),)

Source code in src/skcomplex/models/biology/systems_biology/_s_system.py

def __init__(
    self,
    label: str,
    *,
    species: Annotated[list[str], RoleSpec(variadic=True)],
    alpha: Annotated[Sequence[float], ParameterSpec(per='species', domain=RealSet.closed_open(0, oo))] = (),
    beta: Annotated[Sequence[float], ParameterSpec(per='species', domain=RealSet.closed_open(0, oo))] = (),
    g: Sequence[Sequence[float]] = ((0.0,),),
    h: Sequence[Sequence[float]] = ((0.0,),),
) -> None:
    super().__init__(label, species=species)
    self.alpha = alpha
    self.beta = beta
    self.g = g
    self.h = h

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

Bases: Element

Molecular species, an mRNA, an enzyme, a messenger, or a metabolite, with a nonnegative concentration x.

A species has no laws of its own, and reactions synthesize, convert, and degrade it.

Parameters:

Name	Type	Description	Default
label	str	The label of the 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 initial concentration.	0.0

Examples:

>>> from skcomplex.models.biology.systems_biology import SpeciesElement
>>> SpeciesElement('camp', x=2.0).x
2.0

Source code in src/skcomplex/models/biology/systems_biology/_species.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

ThresholdFunctionInteraction(label, *, target, inputs, weights=(), decays=False)

Bases: MapInteraction

Synchronous update of a node by the sign of the weighted sum of its active inputs, after Li et al.

The node turns on when the sum w_0 x_0 + w_1 x_1 + ... over its inputs is positive and off when it is negative. When the sum is zero, a node that decays turns off and any other node keeps its value.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the updated node.	required
inputs	Annotated[list[str], RoleSpec(variadic=True)]	The paths of the input nodes.	required
weights	Annotated[Sequence[float], ParameterSpec(per=inputs)]	The weight of every input, positive for activation and negative for inhibition.	()
decays	bool	Whether the node turns off without net input.	False

Examples:

>>> from skcomplex.models.biology.systems_biology import ThresholdFunctionInteraction
>>> ThresholdFunctionInteraction('sic1', target='Sic1', inputs=['Swi5', 'Cln1_2'], weights=[1.0, -1.0]).weights
[1.0, -1.0]

Source code in src/skcomplex/models/biology/systems_biology/_threshold_function.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    inputs: Annotated[list[str], RoleSpec(variadic=True)],
    weights: Annotated[Sequence[float], ParameterSpec(per='inputs')] = (),
    decays: bool = False,
) -> None:
    super().__init__(label, target=target, state='x', function=self._get_next)
    self.inputs = inputs
    self.weights = weights
    self.decays = decays

biochemical_systems_theory(x=(0.5, 1.0, 4.0), *, form='s-system', alpha=PATHWAY_ALPHA, beta=PATHWAY_BETA, g=PATHWAY_G, h=PATHWAY_H, stoichiometry=PATHWAY_STOICHIOMETRY, gamma=PATHWAY_GAMMA, orders=PATHWAY_ORDERS)

Return a biochemical system of power-law rates, in the S-system or the generalized mass action form.

The S-system form gives every species one production and one degradation, dx_i/dt = alpha_i prod_j x_j^g_ij - beta_i prod_j x_j^h_ij. The generalized mass action form sums power-law fluxes, dx_i/dt = sum_j s_ij gamma_j prod_k x_k^f_jk. The book states only the two forms. The default is a choice, an unbranched pathway X1 -> X2 -> X3 whose end product inhibits the supply of X1, written in both forms with the same rates, dX1/dt = 2 X3^-0.5 - X1^0.5, dX2/dt = X1^0.5 - X2^0.5, and dX3/dt = X2^0.5 - X3^0.5. The steady state of an S-system solves the linear system (g - h) log x = log(beta / alpha), here x = (2, 2, 2).

Parameters:

Name	Type	Description	Default
x	Sequence[float]	The initial concentrations of the species.	(0.5, 1.0, 4.0)
form	str	's-system' or 'gma'.	's-system'
alpha	Sequence[float]	The rate constants of the productions of the S-system.	PATHWAY_ALPHA
beta	Sequence[float]	The rate constants of the degradations of the S-system.	PATHWAY_BETA
g	Sequence[Sequence[float]]	The kinetic orders of the productions of the S-system, one row per species.	PATHWAY_G
h	Sequence[Sequence[float]]	The kinetic orders of the degradations of the S-system, one row per species.	PATHWAY_H
stoichiometry	Sequence[Sequence[float]]	The stoichiometric coefficients of the generalized mass action form, one row per species.	PATHWAY_STOICHIOMETRY
gamma	Sequence[float]	The rate constants of the fluxes of the generalized mass action form.	PATHWAY_GAMMA
orders	Sequence[Sequence[float]]	The kinetic orders of the fluxes of the generalized mass action form, one row per flux.	PATHWAY_ORDERS

Returns:

Name	Type	Description
system	System	The model, with the concentrations X1.x, X2.x, and so on.

Raises:

Type	Description
ValidationError	The form is unknown.

Examples:

>>> from skcomplex.models.biology.systems_biology import biochemical_systems_theory
>>> biochemical_systems_theory(form='gma').interactions[0].gamma
(2.0, 1.0, 1.0, 1.0)

Source code in src/skcomplex/models/biology/systems_biology/_networks.py

@register_model(
    field='biology/systems_biology',
    tags=('biochemical systems theory', 'power law', 'S-system', 'generalized mass action', 'ode'),
    references=(RAMAN, 'Savageau 1976', 'Voit 2000'),
    defaults={'until': 20.0, 'dt': 0.05},
)
def biochemical_systems_theory(
    x: Sequence[float] = (0.5, 1.0, 4.0),
    *,
    form: str = 's-system',
    alpha: Sequence[float] = PATHWAY_ALPHA,
    beta: Sequence[float] = PATHWAY_BETA,
    g: Sequence[Sequence[float]] = PATHWAY_G,
    h: Sequence[Sequence[float]] = PATHWAY_H,
    stoichiometry: Sequence[Sequence[float]] = PATHWAY_STOICHIOMETRY,
    gamma: Sequence[float] = PATHWAY_GAMMA,
    orders: Sequence[Sequence[float]] = PATHWAY_ORDERS,
) -> System:
    """Return a biochemical system of power-law rates, in the S-system or the generalized mass action form.

    The S-system form gives every species one production and one degradation,
    `dx_i/dt = alpha_i prod_j x_j^g_ij - beta_i prod_j x_j^h_ij`. The generalized mass action form sums power-law
    fluxes, `dx_i/dt = sum_j s_ij gamma_j prod_k x_k^f_jk`. The book states only the two forms. The default is a
    choice, an unbranched pathway `X1 -> X2 -> X3` whose end product inhibits the supply of `X1`, written in both forms
    with the same rates, `dX1/dt = 2 X3^-0.5 - X1^0.5`, `dX2/dt = X1^0.5 - X2^0.5`, and `dX3/dt = X2^0.5 - X3^0.5`.
    The steady state of an S-system solves the linear system `(g - h) log x = log(beta / alpha)`, here `x = (2, 2, 2)`.

    Args:
        x:
            The initial concentrations of the species.

        form:
            `'s-system'` or `'gma'`.

        alpha:
            The rate constants of the productions of the S-system.

        beta:
            The rate constants of the degradations of the S-system.

        g:
            The kinetic orders of the productions of the S-system, one row per species.

        h:
            The kinetic orders of the degradations of the S-system, one row per species.

        stoichiometry:
            The stoichiometric coefficients of the generalized mass action form, one row per species.

        gamma:
            The rate constants of the fluxes of the generalized mass action form.

        orders:
            The kinetic orders of the fluxes of the generalized mass action form, one row per flux.

    Returns:
        system:
            The model, with the concentrations `X1.x`, `X2.x`, and so on.

    Raises:
        ValidationError:
            The form is unknown.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import biochemical_systems_theory
        >>> biochemical_systems_theory(form='gma').interactions[0].gamma
        (2.0, 1.0, 1.0, 1.0)
    """
    species = [SpeciesElement(f'X{index}', x=float(value)) for index, value in enumerate(x, start=1)]
    paths = [item.label for item in species]
    if form == 's-system':
        rates: BaseInteraction = SSystemInteraction(
            'rates',
            species=paths,
            alpha=[float(value) for value in alpha],
            beta=[float(value) for value in beta],
            g=tuple(tuple(row) for row in g),
            h=tuple(tuple(row) for row in h),
        )
    elif form == 'gma':
        rates = GeneralizedMassActionInteraction(
            'rates',
            species=paths,
            stoichiometry=tuple(tuple(row) for row in stoichiometry),
            gamma=tuple(gamma),
            orders=tuple(tuple(row) for row in orders),
        )
    else:
        error_msg = f'Unknown form {form!r}. Use s-system or gma.'
        raise ValidationError(error_msg)
    return System(elements=species, interactions=[rates])

boolean_repressilator(x=(1, 0, 0), *, delay=0)

Return the Boolean repressilator, where LacI represses TetR, TetR represses cI, and cI represses LacI.

The transfer functions are TetR* = NOT LacI, cI* = NOT TetR, and LacI* = NOT cI. From a state with one or two genes on, the genes switch on and off in turn with period 6, and the states 000 and 111 alternate with period 2. A repression delayed by delay steps passes through that many relay nodes, <repressor>_1 and so on, which copy the repressor one step later and start at its initial value. The oscillation remains, with period 6 (delay + 1).

Parameters:

Name	Type	Description	Default
x	Sequence[int]	The initial values of LacI, TetR, and cI.	(1, 0, 0)
delay	int	The number of steps by which every repression is delayed.	0

Returns:

Name	Type	Description
system	System	The model, with the values LacI.x, TetR.x, cI.x, and the relays <repressor>_<k>.x.

Raises:

Type	Description
ValidationError	The delay is negative.

Examples:

>>> from skcomplex.models.biology.systems_biology import boolean_repressilator
>>> [node.label for node in boolean_repressilator(delay=1).elements]
['LacI', 'TetR', 'cI', 'LacI_1', 'TetR_1', 'cI_1']

Source code in src/skcomplex/models/biology/systems_biology/_boolean_models.py

@register_model(
    field='biology/systems_biology',
    tags=('boolean network', 'repressilator', 'oscillator', 'discrete map', 'delay'),
    references=(RAMAN, 'Elowitz and Leibler 2000'),
    defaults={'until': 24, 'dt': 1},
)
def boolean_repressilator(x: Sequence[int] = (1, 0, 0), *, delay: int = 0) -> System:
    """Return the Boolean repressilator, where LacI represses TetR, TetR represses cI, and cI represses LacI.

    The transfer functions are `TetR* = NOT LacI`, `cI* = NOT TetR`, and `LacI* = NOT cI`. From a state with one or two
    genes on, the genes switch on and off in turn with period 6, and the states `000` and `111` alternate with period
    2. A repression delayed by `delay` steps passes through that many relay nodes, `<repressor>_1` and so on, which
    copy the repressor one step later and start at its initial value. The oscillation remains, with period
    `6 (delay + 1)`.

    Args:
        x:
            The initial values of LacI, TetR, and cI.

        delay:
            The number of steps by which every repression is delayed.

    Returns:
        system:
            The model, with the values `LacI.x`, `TetR.x`, `cI.x`, and the relays `<repressor>_<k>.x`.

    Raises:
        ValidationError:
            The delay is negative.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import boolean_repressilator
        >>> [node.label for node in boolean_repressilator(delay=1).elements]
        ['LacI', 'TetR', 'cI', 'LacI_1', 'TetR_1', 'cI_1']
    """
    if delay < 0:
        error_msg = f'The delay must be a nonnegative number of steps, not {delay}.'
        raise ValidationError(error_msg)
    genes = [BooleanNodeElement(name, x=int(value)) for name, value in zip(REPRESSILATOR_GENES, x, strict=True)]
    relays = [BooleanNodeElement(f'{gene.label}_{step}', x=gene.x) for step in range(1, delay + 1) for gene in genes]
    functions: list[BooleanFunctionInteraction] = []
    for repressor, target in zip(REPRESSILATOR_GENES, (*REPRESSILATOR_GENES[1:], REPRESSILATOR_GENES[0]), strict=True):
        chain = [repressor, *(f'{repressor}_{step}' for step in range(1, delay + 1))]
        functions.extend(
            BooleanFunctionInteraction(f'relay_{relay}', target=relay, inputs=[source], table=(0, 1))
            for source, relay in itertools.pairwise(chain)
        )
        functions.append(
            BooleanFunctionInteraction(f'repression_{target}', target=target, inputs=[chain[-1]], table=(1, 0)),
        )
    return System(elements=[*genes, *relays], interactions=functions)

circadian_clock(*, v_s=0.76, v_d=0.95, k_s=0.38, n=4.0)

Return Goldbeter's model of the Drosophila circadian clock, PER that represses its own transcription.

Multiple phosphorylation and nuclear entry delay the negative feedback, which yields a limit cycle with a period near 24 hours for the published parameters. Larger degradation v_d lengthens the period, from 22.6 hours at v_d = 0.8 to 29.3 hours at v_d = 1.6.

Parameters:

Name	Type	Description	Default
v_s	float	The maximal transcription rate in µM/h.	0.76
v_d	float	The maximal degradation rate of bisphosphorylated PER in µM/h.	0.95
k_s	float	The translation rate per hour.	0.38
n	float	The Hill coefficient of repression.	4.0

Returns:

Name	Type	Description
system	System	The model, with the concentrations clock.M, clock.P0, clock.P1, clock.P2, and clock.PN.

Examples:

>>> from skcomplex.models.biology.systems_biology import circadian_clock
>>> circadian_clock(v_d=1.2).elements[0].v_d
1.2

Source code in src/skcomplex/models/biology/systems_biology/_clocks.py

@register_model(
    field='biology/systems_biology',
    tags=('circadian rhythm', 'oscillations', 'gene regulation'),
    references=('Goldbeter 1995', 'Trzaska 2023'),
    defaults={'until': 240.0, 'dt': 0.5},
)
def circadian_clock(
    *,
    v_s: float = 0.76,
    v_d: float = 0.95,
    k_s: float = 0.38,
    n: float = 4.0,
) -> System:
    """Return Goldbeter's model of the Drosophila circadian clock, PER that represses its own transcription.

    Multiple phosphorylation and nuclear entry delay the negative feedback, which yields a limit cycle with a period
    near 24 hours for the published parameters. Larger degradation `v_d` lengthens the period, from 22.6 hours at
    `v_d = 0.8` to 29.3 hours at `v_d = 1.6`.

    Args:
        v_s:
            The maximal transcription rate in µM/h.

        v_d:
            The maximal degradation rate of bisphosphorylated PER in µM/h.

        k_s:
            The translation rate per hour.

        n:
            The Hill coefficient of repression.

    Returns:
        system:
            The model, with the concentrations `clock.M`, `clock.P0`, `clock.P1`, `clock.P2`, and `clock.PN`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import circadian_clock
        >>> circadian_clock(v_d=1.2).elements[0].v_d
        1.2
    """
    clock = PerClockElement('clock', v_s=v_s, v_d=v_d, k_s=k_s, n=n)
    reactions = [
        ReactionInteraction(
            'transcription',
            target='clock',
            products={'M': 1},
            propensity=lambda s: s.target.v_s * (1 - hill(s.target.PN, s.target.K_I, s.target.n)),
        ),
        MichaelisMentenInteraction(
            'mrna_degradation',
            target='clock',
            reactants={'M': 1},
            substrate='M',
            vmax=lambda s: s.target.v_m,
            km=lambda s: s.target.K_m,
        ),
        ReactionInteraction(
            'translation',
            target='clock',
            products={'P0': 1},
            propensity=lambda s: s.target.k_s * s.target.M,
        ),
        MichaelisMentenInteraction(
            'first_phosphorylation',
            target='clock',
            reactants={'P0': 1},
            products={'P1': 1},
            substrate='P0',
            vmax=lambda s: s.target.V1,
            km=lambda s: s.target.K,
        ),
        MichaelisMentenInteraction(
            'first_dephosphorylation',
            target='clock',
            reactants={'P1': 1},
            products={'P0': 1},
            substrate='P1',
            vmax=lambda s: s.target.V2,
            km=lambda s: s.target.K,
        ),
        MichaelisMentenInteraction(
            'second_phosphorylation',
            target='clock',
            reactants={'P1': 1},
            products={'P2': 1},
            substrate='P1',
            vmax=lambda s: s.target.V3,
            km=lambda s: s.target.K,
        ),
        MichaelisMentenInteraction(
            'second_dephosphorylation',
            target='clock',
            reactants={'P2': 1},
            products={'P1': 1},
            substrate='P2',
            vmax=lambda s: s.target.V4,
            km=lambda s: s.target.K,
        ),
        ReactionInteraction(
            'nuclear_entry',
            target='clock',
            reactants={'P2': 1},
            products={'PN': 1},
            propensity=lambda s: s.target.k1 * s.target.P2,
        ),
        ReactionInteraction(
            'nuclear_exit',
            target='clock',
            reactants={'PN': 1},
            products={'P2': 1},
            propensity=lambda s: s.target.k2 * s.target.PN,
        ),
        MichaelisMentenInteraction(
            'per_degradation',
            target='clock',
            reactants={'P2': 1},
            substrate='P2',
            vmax=lambda s: s.target.v_d,
            km=lambda s: s.target.K_d,
        ),
    ]
    return System(elements=[clock], interactions=reactions)

goodwin_oscillator(variables=3, x=None, *, a=100.0, b=1.0, K=1.0, n=12.0, k=1.0)

Return the Goodwin oscillator, a chain of stages whose last one represses the synthesis of the first.

The first stage, an mRNA, is synthesized at a / (1 + (x_N / K)^n) and every later stage at k times the stage before it, and every stage decays at rate b: dx_1/dt = a / (1 + (x_N / K)^n) - b x_1 and dx_i/dt = k x_(i-1) - b x_i. The book points to BioModels and prints no equations, so the model takes the standard form of Goodwin with the degradation of Griffith. With equal degradation rates the fixed point loses stability when the loop gain passes sec(pi / N)^N, 8 for three stages. When a is large the repressor saturates and the gain is close to n, so three stages oscillate for Hill coefficients above about 8 and settle for lower ones. The concentrations and the time are dimensionless.

Parameters:

Name	Type	Description	Default
variables	int	The number of stages, at least 2.	3
x	Sequence[float] | None	The initial concentrations of the stages, or None for 0.1, 0.2, and so on.	None
a	float	The maximal synthesis rate of the first stage.	100.0
b	float	The degradation rate of every stage.	1.0
K	float	The concentration of the last stage at half repression.	1.0
n	float	The Hill coefficient of the repression.	12.0
k	float	The synthesis rate of every later stage per unit of the stage before it.	1.0

Returns:

Name	Type	Description
system	System	The model, with the concentrations x1.x to x<variables>.x.

Raises:

Type	Description
ValidationError	There are fewer than two stages, or the initial concentrations do not match them.

Examples:

>>> from skcomplex.models.biology.systems_biology import goodwin_oscillator
>>> [stage.label for stage in goodwin_oscillator(4).elements]
['x1', 'x2', 'x3', 'x4']

Source code in src/skcomplex/models/biology/systems_biology/_networks.py

@register_model(
    field='biology/systems_biology',
    tags=('oscillator', 'negative feedback', 'Hill function', 'gene regulation', 'ode'),
    references=(RAMAN, 'Goodwin 1965', 'Griffith 1968'),
    defaults={'until': 60.0, 'dt': 0.1},
)
def goodwin_oscillator(
    variables: int = 3,
    x: Sequence[float] | None = None,
    *,
    a: float = 100.0,
    b: float = 1.0,
    K: float = 1.0,
    n: float = 12.0,
    k: float = 1.0,
) -> System:
    """Return the Goodwin oscillator, a chain of stages whose last one represses the synthesis of the first.

    The first stage, an mRNA, is synthesized at `a / (1 + (x_N / K)^n)` and every later stage at `k` times the stage
    before it, and every stage decays at rate `b`: `dx_1/dt = a / (1 + (x_N / K)^n) - b x_1` and
    `dx_i/dt = k x_(i-1) - b x_i`. The book points to BioModels and prints no equations, so the model takes the standard
    form of Goodwin with the degradation of Griffith. With equal degradation rates the fixed point loses stability when
    the loop gain passes `sec(pi / N)^N`, 8 for three stages. When `a` is large the repressor saturates and the gain is
    close to `n`, so three stages oscillate for Hill coefficients above about 8 and settle for lower ones. The
    concentrations and the time are dimensionless.

    Args:
        variables:
            The number of stages, at least 2.

        x:
            The initial concentrations of the stages, or `None` for 0.1, 0.2, and so on.

        a:
            The maximal synthesis rate of the first stage.

        b:
            The degradation rate of every stage.

        K:
            The concentration of the last stage at half repression.

        n:
            The Hill coefficient of the repression.

        k:
            The synthesis rate of every later stage per unit of the stage before it.

    Returns:
        system:
            The model, with the concentrations `x1.x` to `x<variables>.x`.

    Raises:
        ValidationError:
            There are fewer than two stages, or the initial concentrations do not match them.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import goodwin_oscillator
        >>> [stage.label for stage in goodwin_oscillator(4).elements]
        ['x1', 'x2', 'x3', 'x4']
    """
    initial = [0.1 * index for index in range(1, variables + 1)] if x is None else list(x)
    if variables < MIN_GOODWIN_STAGES or len(initial) != variables:
        error_msg = f'The Goodwin loop needs at least 2 stages and one initial value per stage, not {variables} stages.'
        raise ValidationError(error_msg)
    labels = [f'x{index}' for index in range(1, variables + 1)]
    stages = [SpeciesElement(label, x=float(value)) for label, value in zip(labels, initial, strict=True)]
    repression = RepressionInteraction(
        'repression',
        products={(labels[0], 'x'): 1},
        modifiers={(labels[-1], 'x'): 1},
        alpha=a,
        K=K,
        n=n,
    )
    syntheses = [
        ReactionInteraction(f'synthesis_{target}', products={(target, 'x'): 1}, modifiers={(source, 'x'): 1}, k=k)
        for source, target in itertools.pairwise(labels)
    ]
    degradations = [ReactionInteraction(f'degradation_{label}', reactants={(label, 'x'): 1}, k=b) for label in labels]
    return System(elements=stages, interactions=[repression, *syntheses, *degradations])

hormone_cascade_amplifier(alpha=(10.0, 100.0, 1.0, 10.0, 10.0), beta=(1.0, 2.0, 0.5, 1.0, 1.0), *, k_on=1.0, k_gtp=0.5, hormone=1.0, start=1.0, width=1.0)

Return the adenylate cyclase cascade, which amplifies a hormone pulse and ends it by the G-protein timer.

A hormone pulse activates the G protein, which activates adenylate cyclase, which makes cyclic AMP, which activates protein kinase A, phosphorylase kinase, and glycogen phosphorylase in turn (Konieczny et al. 4.4.3). Every stage is activated at the rate alpha_i times the active amount of the stage above it and inactivated at the rate beta_i, so the integral of every stage over the signal is the integral of the stage above it times the gain alpha_i / beta_i, and the cascade multiplies the integral of the G-protein signal by the product of the gains, 100,000 with the defaults. The G protein inactivates itself by GTP hydrolysis at the rate k_gtp, the timer of Konieczny et al. 4.4.6, so blocking it, as cholera toxin does, prolongs the signal.

Parameters:

Name	Type	Description	Default
alpha	Sequence[float]	The activation rates of adenylate cyclase, cyclic AMP, protein kinase A, phosphorylase kinase, and phosphorylase.	(10.0, 100.0, 1.0, 10.0, 10.0)
beta	Sequence[float]	The inactivation rates of the five stages.	(1.0, 2.0, 0.5, 1.0, 1.0)
k_on	float	The activation rate of the G protein per unit of hormone.	1.0
k_gtp	float	The rate of GTP hydrolysis of the G protein.	0.5
hormone	float	The height of the hormone pulse.	1.0
start	float	The time the hormone pulse starts.	1.0
width	float	The duration of the hormone pulse.	1.0

Returns:

Name	Type	Description
system	System	The model, with the active G protein g_protein.g and the active stages adenylate_cyclase.x, camp.x, protein_kinase_a.x, phosphorylase_kinase.x, and phosphorylase.x.

Raises:

Type	Description
ValidationError	The activation and inactivation rates do not give one value per stage.

Examples:

>>> from skcomplex.models.biology.systems_biology import hormone_cascade_amplifier
>>> [element.label for element in hormone_cascade_amplifier().elements][-1]
'phosphorylase'

Source code in src/skcomplex/models/biology/systems_biology/_cascades.py

@register_model(
    field='biology/systems_biology',
    tags=('signalling', 'cascade', 'amplification', 'hormone', 'g protein', 'ode'),
    references=('Konieczny et al. 2023', 'Heinrich, Neel, and Rapoport 2002'),
    defaults={'until': 30.0, 'dt': 0.2},
)
def hormone_cascade_amplifier(
    alpha: Sequence[float] = (10.0, 100.0, 1.0, 10.0, 10.0),
    beta: Sequence[float] = (1.0, 2.0, 0.5, 1.0, 1.0),
    *,
    k_on: float = 1.0,
    k_gtp: float = 0.5,
    hormone: float = 1.0,
    start: float = 1.0,
    width: float = 1.0,
) -> System:
    """Return the adenylate cyclase cascade, which amplifies a hormone pulse and ends it by the G-protein timer.

    A hormone pulse activates the G protein, which activates adenylate cyclase, which makes cyclic AMP, which activates
    protein kinase A, phosphorylase kinase, and glycogen phosphorylase in turn (Konieczny et al. 4.4.3). Every stage is
    activated at the rate `alpha_i` times the active amount of the stage above it and inactivated at the rate
    `beta_i`, so the integral of every stage over the signal is the integral of the stage above it times the gain
    `alpha_i / beta_i`, and the cascade multiplies the integral of the G-protein signal by the product of the gains,
    100,000 with the defaults. The G protein inactivates itself by GTP hydrolysis at the rate `k_gtp`, the timer of
    Konieczny et al. 4.4.6, so blocking it, as cholera toxin does, prolongs the signal.

    Args:
        alpha:
            The activation rates of adenylate cyclase, cyclic AMP, protein kinase A, phosphorylase kinase, and
            phosphorylase.

        beta:
            The inactivation rates of the five stages.

        k_on:
            The activation rate of the G protein per unit of hormone.

        k_gtp:
            The rate of GTP hydrolysis of the G protein.

        hormone:
            The height of the hormone pulse.

        start:
            The time the hormone pulse starts.

        width:
            The duration of the hormone pulse.

    Returns:
        system:
            The model, with the active G protein `g_protein.g` and the active stages `adenylate_cyclase.x`, `camp.x`,
            `protein_kinase_a.x`, `phosphorylase_kinase.x`, and `phosphorylase.x`.

    Raises:
        ValidationError:
            The activation and inactivation rates do not give one value per stage.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import hormone_cascade_amplifier
        >>> [element.label for element in hormone_cascade_amplifier().elements][-1]
        'phosphorylase'
    """
    if len(alpha) != len(STAGES) or len(beta) != len(STAGES):
        error_msg = f'Give one activation and one inactivation rate per stage, {len(STAGES)} of each.'
        raise ValidationError(error_msg)
    g_protein = GProteinElement('g_protein', k_on=k_on, k_gtp=k_gtp, hormone=hormone, start=start, width=width)
    stages = [SpeciesElement(label) for label in STAGES]
    upstreams = (('g_protein', 'g'), *((label, 'x') for label in STAGES[:-1]))
    activations = [
        ReactionInteraction(f'activation_of_{label}', products={(label, 'x'): 1}, modifiers={upstream: 1}, k=rate)
        for upstream, label, rate in zip(upstreams, STAGES, alpha, strict=True)
    ]
    inactivations = [
        ReactionInteraction(f'inactivation_of_{label}', reactants={(label, 'x'): 1}, k=rate)
        for label, rate in zip(STAGES, beta, strict=True)
    ]
    switching = [
        ReactionInteraction(
            'g_protein_activation',
            target='g_protein',
            products={'g': 1},
            propensity=lambda s: s.target.k_on
            * pulse(s.t, s.target.start, s.target.width, s.target.hormone)
            * (1 - s.target.g),
        ),
        ReactionInteraction(
            'gtp_hydrolysis',
            target='g_protein',
            reactants={'g': 1},
            propensity=lambda s: s.target.k_gtp * s.target.g,
        ),
    ]
    return System(elements=[g_protein, *stages], interactions=[*switching, *activations, *inactivations])

mammalian_cell_cycle(active=CELL_CYCLE_START, *, knockouts=())

Return the Boolean model of the mammalian cell cycle of Faure et al., ten proteins updated synchronously.

The rules are those printed by Raman for CycD, Rb, E2F, CycE, CycA, p27, Cdc20, Cdh1, UbcH10 (Ubc), and CycB. The growth signal CycD keeps its value. Without it the network rests in the quiescent state with Rb, p27, and Cdh1 on. With it the network leaves that state and runs through a cycle of period 7 in which the cyclins E, A, and B rise and fall in turn. A knocked-out protein stays off.

Parameters:

Name	Type	Description	Default
active	Sequence[str]	The proteins that are on at the start.	CELL_CYCLE_START
knockouts	Sequence[str]	The proteins that are knocked out.	()

Returns:

Name	Type	Description
system	System	The model, with the values <protein>.x.

Raises:

Type	Description
ValidationError	A name is not one of the ten proteins.

Examples:

>>> from skcomplex.models.biology.systems_biology import mammalian_cell_cycle
>>> [node.label for node in mammalian_cell_cycle().elements if node.x]
['CycD', 'Rb', 'p27', 'Cdh1']

Source code in src/skcomplex/models/biology/systems_biology/_boolean_models.py

@register_model(
    field='biology/systems_biology',
    tags=('boolean network', 'cell cycle', 'discrete map', 'attractors', 'knock-out'),
    references=(RAMAN, 'Faure et al. 2006'),
    defaults={'until': 21, 'dt': 1},
)
def mammalian_cell_cycle(
    active: Sequence[str] = CELL_CYCLE_START,
    *,
    knockouts: Sequence[str] = (),
) -> System:
    """Return the Boolean model of the mammalian cell cycle of Faure et al., ten proteins updated synchronously.

    The rules are those printed by Raman for CycD, Rb, E2F, CycE, CycA, p27, Cdc20, Cdh1, UbcH10 (`Ubc`), and CycB.
    The growth signal CycD keeps its value. Without it the network rests in the quiescent state with Rb, p27, and Cdh1
    on. With it the network leaves that state and runs through a cycle of period 7 in which the cyclins E, A, and B
    rise and fall in turn. A knocked-out protein stays off.

    Args:
        active:
            The proteins that are on at the start.

        knockouts:
            The proteins that are knocked out.

    Returns:
        system:
            The model, with the values `<protein>.x`.

    Raises:
        ValidationError:
            A name is not one of the ten proteins.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import mammalian_cell_cycle
        >>> [node.label for node in mammalian_cell_cycle().elements if node.x]
        ['CycD', 'Rb', 'p27', 'Cdh1']
    """
    return _rule_network(CELL_CYCLE_RULES, active, knockouts)

min_protein_skeleton(n_sites=128, *, length=100.0, total_D=200.0, total_E=120.0, random_state=0)

Return the skeleton model of the Min proteins on a periodic line, a mass-conserving reaction-diffusion system.

MinD and MinE cycle between the membrane, where they diffuse slowly, and the cytosol, where they diffuse fast, with the reactions of MinSkeletonElement. The reactions conserve the totals of MinD and MinE, so only diffusion moves mass, and the local reactive equilibria map the dynamics onto conserved fields with nonreciprocal cross-diffusion (Brauns and Marchetti, App. F). With the default totals, the uniform state is unstable to an oscillatory band of wavelengths near 17, and narrow bands of membrane MinD travel along the line and meet head-on. The fast cytosolic diffusion and the fast recruitment make the system stiff, so the catalog registers the implicit BDF method with the exact Jacobian.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites of the periodic line.	128
length	float	The length of the line.	100.0
total_D	float	The mean total density of MinD.	200.0
total_E	float	The mean total density of MinE.	120.0
random_state	int | Generator | None	The seed or generator of the 1% noise that moves MinD between the membrane and the cytosol.	0

Returns:

Name	Type	Description
system	System	The model, with the membrane fields cell.md and cell.mde and the cytosolic fields cell.cdt, cell.cdd, and cell.ce.

Examples:

>>> from skcomplex.models.biology.systems_biology import min_protein_skeleton
>>> values = min_protein_skeleton(32).elements[0].values
>>> round(float(values['ce'].mean() + values['mde'].mean()), 6)
120.0

Source code in src/skcomplex/models/biology/systems_biology/_min_proteins.py

@register_model(
    field='biology/systems_biology',
    tags=('field', 'min proteins', 'mass conserving', 'reaction-diffusion', 'traveling waves', 'pde'),
    references=('Brauns and Marchetti 2024', 'Halatek and Frey 2018', 'Huang, Meir, and Wingreen 2003'),
    defaults={'until': 400.0, 'dt': 4.0, 'solver': ScipyODESolver(method='BDF', jacobian='exact')},
)
def min_protein_skeleton(
    n_sites: int = 128,
    *,
    length: float = 100.0,
    total_D: float = 200.0,
    total_E: float = 120.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the skeleton model of the Min proteins on a periodic line, a mass-conserving reaction-diffusion system.

    MinD and MinE cycle between the membrane, where they diffuse slowly, and the cytosol, where they diffuse fast,
    with the reactions of `MinSkeletonElement`. The reactions conserve the totals of MinD and MinE, so only diffusion
    moves mass, and the local reactive equilibria map the dynamics onto conserved fields with nonreciprocal
    cross-diffusion (Brauns and Marchetti, App. F). With the default totals, the uniform state is unstable to an
    oscillatory band of wavelengths near 17, and narrow bands of membrane MinD travel along the line and meet head-on.
    The fast cytosolic diffusion and the fast recruitment make the system stiff, so the catalog registers the implicit
    `BDF` method with the exact Jacobian.

    Args:
        n_sites:
            The number of sites of the periodic line.

        length:
            The length of the line.

        total_D:
            The mean total density of MinD.

        total_E:
            The mean total density of MinE.

        random_state:
            The seed or generator of the 1% noise that moves MinD between the membrane and the cytosol.

    Returns:
        system:
            The model, with the membrane fields `cell.md` and `cell.mde` and the cytosolic fields `cell.cdt`,
            `cell.cdd`, and `cell.ce`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import min_protein_skeleton
        >>> values = min_protein_skeleton(32).elements[0].values
        >>> round(float(values['ce'].mean() + values['mde'].mean()), 6)
        120.0
    """
    rng = np.random.default_rng(random_state)
    uniform = _uniform_state(MinSkeletonElement('site'), total_D, total_E)
    values: dict[str, Array] = {name: np.full(n_sites, value) for name, value in uniform.items()}
    noise = 0.01 * uniform['md'] * rng.uniform(-1.0, 1.0, n_sites)
    values['md'] = values['md'] + noise
    values['cdt'] = values['cdt'] - noise
    space = LatticeSpace(n_sites, spacing=length / n_sites)
    cell = FieldElement('cell', element=MinSkeletonElement, space=space, values=values)
    reactions = [
        ReactionInteraction(
            'attachment',
            target='cell',
            reactants={'cdt': 1},
            products={'md': 1},
            propensity=lambda s: (s.target.k_d + s.target.k_dd * s.target.md) * s.target.cdt,
        ),
        ReactionInteraction(
            'recruitment',
            target='cell',
            reactants={'md': 1, 'ce': 1},
            products={'mde': 1},
            propensity=lambda s: s.target.k_de * s.target.md * s.target.ce,
        ),
        ReactionInteraction(
            'release',
            target='cell',
            reactants={'mde': 1},
            products={'cdd': 1, 'ce': 1},
            propensity=lambda s: s.target.k_off * s.target.mde,
        ),
        ReactionInteraction(
            'exchange',
            target='cell',
            reactants={'cdd': 1},
            products={'cdt': 1},
            propensity=lambda s: s.target.lam * s.target.cdd,
        ),
    ]
    membrane = [
        DiffusionInteraction(f'diffusion_{state}', target='cell', state=state, coefficient=lambda s: s.target.Dm)
        for state in ('md', 'mde')
    ]
    cytosol = [
        DiffusionInteraction(f'diffusion_{state}', target='cell', state=state, coefficient=lambda s: s.target.Dc)
        for state in ('cdt', 'cdd', 'ce')
    ]
    return System(elements=[cell], interactions=[*reactions, *membrane, *cytosol])

negative_autoregulation(p=0.0, *, alpha=2.0, K=50.0, n=1.0, b=10.0, beta=0.1, delta=1.0, stochastic=True)

Return a gene that represses its own transcription and makes its protein in bursts.

Every mRNA molecule yields on average b / delta proteins before it decays, so with fast mRNA decay the protein arrives in geometric bursts and its count fluctuates far more than a Poisson count. Negative feedback lowers the mean and the noise. A very large K removes the feedback.

Parameters:

Name	Type	Description	Default
p	float	The initial protein count, with no mRNA.	0.0
alpha	float	The transcription rate without repressor.	2.0
K	float	The protein count at half repression.	50.0
n	float	The Hill coefficient of repression.	1.0
b	float	The translation rate per mRNA molecule.	10.0
beta	float	The decay rate of the protein.	0.1
delta	float	The decay rate of mRNA.	1.0
stochastic	bool	Whether the model is a jump process.	True

Returns:

Name	Type	Description
system	System	The model, with the counts gene.m and gene.p.

Examples:

>>> from skcomplex.models.biology.systems_biology import negative_autoregulation
>>> negative_autoregulation().interactions[0].modifiers
{('gene', 'p'): 1}

Source code in src/skcomplex/models/biology/systems_biology/_genes.py

@register_model(
    field='biology/systems_biology',
    tags=('gene regulation', 'noise', 'bursts', 'jump process'),
    references=('Thattai and van Oudenaarden 2001', 'Becskei and Serrano 2000', 'arXiv:1608.06582'),
    defaults={'until': 200.0, 'dt': 1.0},
)
def negative_autoregulation(
    p: float = 0.0,
    *,
    alpha: float = 2.0,
    K: float = 50.0,
    n: float = 1.0,
    b: float = 10.0,
    beta: float = 0.1,
    delta: float = 1.0,
    stochastic: bool = True,
) -> System:
    """Return a gene that represses its own transcription and makes its protein in bursts.

    Every mRNA molecule yields on average `b / delta` proteins before it decays, so with fast mRNA decay the protein
    arrives in geometric bursts and its count fluctuates far more than a Poisson count. Negative feedback lowers the
    mean and the noise. A very large `K` removes the feedback.

    Args:
        p:
            The initial protein count, with no mRNA.

        alpha:
            The transcription rate without repressor.

        K:
            The protein count at half repression.

        n:
            The Hill coefficient of repression.

        b:
            The translation rate per mRNA molecule.

        beta:
            The decay rate of the protein.

        delta:
            The decay rate of mRNA.

        stochastic:
            Whether the model is a jump process.

    Returns:
        system:
            The model, with the counts `gene.m` and `gene.p`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import negative_autoregulation
        >>> negative_autoregulation().interactions[0].modifiers
        {('gene', 'p'): 1}
    """
    mode = 'jump' if stochastic else 'rate'
    feedback = _repress('autoregulation', 'gene', 'gene', {'alpha': alpha, 'alpha0': 0.0, 'K': K, 'n': n}, mode)
    expression = _express('gene', b, beta, delta, mode)
    return System(elements=[GeneElement('gene', p=p)], interactions=[feedback, *expression])

phosphorylation_oscillator(A=1.0, B=0.0, C=0.0, *, kp1=1.0, kd1=2.0, kp2=2.0, kd2=1.0, kp3=2.0, kd3=1.0, K=0.01)

Return three proteins whose phosphorylation cycles regulate each other in a loop that oscillates.

C converts A, A converts B, and B converts C, each reaction with Michaelis-Menten kinetics and each total conserved, so the loop is a negative feedback of three steps. The book prints the equations and no parameter values. The defaults are a choice: totals of 1, Michaelis constants K far below the totals, which make every cycle an ultrasensitive switch, and rates that put every switch at half the total of its enzyme. The proteins then switch on and off in turn on a limit cycle of period about 4.2. With large K the switches become graded and the loop settles on a fixed point. The amounts are fractions of the totals.

Parameters:

Name	Type	Description	Default
A	float	The initial amount of unphosphorylated A.	1.0
B	float	The initial amount of unphosphorylated B.	0.0
C	float	The initial amount of unphosphorylated C.	0.0
kp1	float	The maximal rate that restores A.	1.0
kd1	float	The catalytic rate of the conversion of A by C.	2.0
kp2	float	The catalytic rate of the conversion of B by A.	2.0
kd2	float	The maximal rate that removes B.	1.0
kp3	float	The catalytic rate of the conversion of C by B.	2.0
kd3	float	The maximal rate that removes C.	1.0
K	float	The Michaelis constant of every reaction.	0.01

Returns:

Name	Type	Description
system	System	The model, with the amounts loop.A, loop.B, and loop.C.

Examples:

>>> from skcomplex.models.biology.systems_biology import phosphorylation_oscillator
>>> phosphorylation_oscillator(K=0.1).elements[0].K2
0.1

Source code in src/skcomplex/models/biology/systems_biology/_networks.py

@register_model(
    field='biology/systems_biology',
    tags=('oscillator', 'phosphorylation', 'Michaelis-Menten', 'negative feedback', 'ode'),
    references=(RAMAN, 'Mogilner, Wollman, and Marshall 2006'),
    defaults={'until': 40.0, 'dt': 0.05},
)
def phosphorylation_oscillator(
    A: float = 1.0,
    B: float = 0.0,
    C: float = 0.0,
    *,
    kp1: float = 1.0,
    kd1: float = 2.0,
    kp2: float = 2.0,
    kd2: float = 1.0,
    kp3: float = 2.0,
    kd3: float = 1.0,
    K: float = 0.01,
) -> System:
    """Return three proteins whose phosphorylation cycles regulate each other in a loop that oscillates.

    C converts A, A converts B, and B converts C, each reaction with Michaelis-Menten kinetics and each total
    conserved, so the loop is a negative feedback of three steps. The book prints the equations and no parameter
    values. The defaults are a choice: totals of 1, Michaelis constants `K` far below the totals, which make every
    cycle an ultrasensitive switch, and rates that put every switch at half the total of its enzyme. The proteins then
    switch on and off in turn on a limit cycle of period about 4.2. With large `K` the switches become graded and the
    loop settles on a fixed point. The amounts are fractions of the totals.

    Args:
        A:
            The initial amount of unphosphorylated A.

        B:
            The initial amount of unphosphorylated B.

        C:
            The initial amount of unphosphorylated C.

        kp1:
            The maximal rate that restores A.

        kd1:
            The catalytic rate of the conversion of A by C.

        kp2:
            The catalytic rate of the conversion of B by A.

        kd2:
            The maximal rate that removes B.

        kp3:
            The catalytic rate of the conversion of C by B.

        kd3:
            The maximal rate that removes C.

        K:
            The Michaelis constant of every reaction.

    Returns:
        system:
            The model, with the amounts `loop.A`, `loop.B`, and `loop.C`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import phosphorylation_oscillator
        >>> phosphorylation_oscillator(K=0.1).elements[0].K2
        0.1
    """
    loop = PhosphorylationLoopElement(
        'loop',
        A=A,
        B=B,
        C=C,
        kp1=kp1,
        kd1=kd1,
        kp2=kp2,
        kd2=kd2,
        kp3=kp3,
        kd3=kd3,
        K1=K,
        K2=K,
        K3=K,
    )
    reactions = [
        MichaelisMentenInteraction(
            'a_restoration',
            target='loop',
            products={'A': 1},
            substrate=lambda s: s.target.A_T - s.target.A,
            vmax=lambda s: s.target.kp1,
            km=lambda s: s.target.K1,
        ),
        MichaelisMentenInteraction(
            'a_conversion',
            target='loop',
            reactants={'A': 1},
            substrate='A',
            vmax=lambda s: s.target.kd1,
            km=lambda s: s.target.K1,
            modifiers={'C': 1},
        ),
        MichaelisMentenInteraction(
            'b_conversion',
            target='loop',
            products={'B': 1},
            substrate=lambda s: s.target.B_T - s.target.B,
            vmax=lambda s: s.target.kp2,
            km=lambda s: s.target.K2,
            modifiers={'A': 1},
        ),
        MichaelisMentenInteraction(
            'b_removal',
            target='loop',
            reactants={'B': 1},
            substrate='B',
            vmax=lambda s: s.target.kd2,
            km=lambda s: s.target.K2,
        ),
        MichaelisMentenInteraction(
            'c_conversion',
            target='loop',
            products={'C': 1},
            substrate=lambda s: s.target.C_T - s.target.C,
            vmax=lambda s: s.target.kp3,
            km=lambda s: s.target.K3,
            modifiers={'B': 1},
        ),
        MichaelisMentenInteraction(
            'c_removal',
            target='loop',
            reactants={'C': 1},
            substrate='C',
            vmax=lambda s: s.target.kd3,
            km=lambda s: s.target.K3,
        ),
    ]
    return System(elements=[loop], interactions=reactions)

probabilistic_boolean_network(network=None, *, inputs=SHMULEVICH_INPUTS, active=('x1',))

Return a probabilistic Boolean network, whose nodes draw one of their candidate functions at every step.

Every node has candidate truth tables of the inputs, each with a probability, and at every step every node draws one of its tables independently and takes its value, all nodes together. The default network follows the three-gene example of Shmulevich et al., in which x1 and x3 have two candidate functions and x2 one. From one initial state the network reaches different attractors in different runs.

Parameters:

Name	Type	Description	Default
network	Mapping[str, Candidates] | None	A mapping from node name to its candidate tables with their probabilities, or None for the three-gene example.	None
inputs	Sequence[str]	The nodes every table reads, in the order of its bits, from the most significant.	SHMULEVICH_INPUTS
active	Sequence[str]	The nodes that are on at the start.	('x1',)

Returns:

Name	Type	Description
system	System	The model, with the values <node>.x and the indices of the drawn functions <node>.function.

Examples:

>>> from skcomplex.models.biology.systems_biology import probabilistic_boolean_network
>>> probabilistic_boolean_network().interactions[0].probabilities
(0.6, 0.4)

Source code in src/skcomplex/models/biology/systems_biology/_probabilistic.py

@register_model(
    field='biology/systems_biology',
    tags=('boolean network', 'gene regulation', 'discrete map', 'stochastic'),
    references=('Shmulevich et al. 2002', 'Raman 2021'),
    defaults={'until': 100, 'dt': 1},
)
def probabilistic_boolean_network(
    network: Mapping[str, Candidates] | None = None,
    *,
    inputs: Sequence[str] = SHMULEVICH_INPUTS,
    active: Sequence[str] = ('x1',),
) -> System:
    """Return a probabilistic Boolean network, whose nodes draw one of their candidate functions at every step.

    Every node has candidate truth tables of the inputs, each with a probability, and at every step every node draws
    one of its tables independently and takes its value, all nodes together. The default network follows the three-gene
    example of Shmulevich et al., in which `x1` and `x3` have two candidate functions and `x2` one. From one initial
    state the network reaches different attractors in different runs.

    Args:
        network:
            A mapping from node name to its candidate tables with their probabilities, or `None` for the three-gene
            example.

        inputs:
            The nodes every table reads, in the order of its bits, from the most significant.

        active:
            The nodes that are on at the start.

    Returns:
        system:
            The model, with the values `<node>.x` and the indices of the drawn functions `<node>.function`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import probabilistic_boolean_network
        >>> probabilistic_boolean_network().interactions[0].probabilities
        (0.6, 0.4)
    """
    chosen = SHMULEVICH_NETWORK if network is None else network
    nodes = [ProbabilisticNodeElement(name, x=int(name in active)) for name in chosen]
    functions = [
        ProbabilisticFunctionInteraction(
            f'functions_{name}',
            target=name,
            inputs=list(inputs),
            tables=tuple(tuple(table) for table, _ in candidates),
            probabilities=tuple(probability for _, probability in candidates),
        )
        for name, candidates in chosen.items()
    ]
    return System(elements=nodes, interactions=functions)

random_boolean_network(n_nodes=20, k=2, *, bias=0.5, random_state=0)

Return a Kauffman random Boolean network, nodes updated together by random functions of k random inputs.

Each node reads k distinct nodes, possibly itself, through a truth table whose entries are 1 with probability bias. The finite state space makes every run end on a cycle. With k = 2 and bias = 0.5 the network is critical, between frozen and chaotic.

Parameters:

Name	Type	Description	Default
n_nodes	int	The number of nodes.	20
k	int	The number of inputs per node, at most the number of nodes.	2
bias	float	The probability that an entry of a truth table is 1.	0.5
random_state	int | Generator | None	The seed or generator of the wiring, the tables, and the initial values.	0

Returns:

Name	Type	Description
system	System	The model, with the values node_<i>.x.

Raises:

Type	Description
ValidationError	k is not between 1 and the number of nodes.

Examples:

>>> from skcomplex.models.biology.systems_biology import random_boolean_network
>>> len(random_boolean_network(10, 3).interactions[0].table)
8

Source code in src/skcomplex/models/biology/systems_biology/_boolean.py

@register_model(
    field='biology/systems_biology',
    tags=('boolean network', 'gene regulation', 'discrete map', 'criticality'),
    references=('Kauffman 1969', 'Gros 2015'),
    defaults={'until': 50, 'dt': 1},
)
def random_boolean_network(
    n_nodes: int = 20,
    k: int = 2,
    *,
    bias: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a Kauffman random Boolean network, nodes updated together by random functions of `k` random inputs.

    Each node reads `k` distinct nodes, possibly itself, through a truth table whose entries are 1 with probability
    `bias`. The finite state space makes every run end on a cycle. With `k = 2` and `bias = 0.5` the network is
    critical, between frozen and chaotic.

    Args:
        n_nodes:
            The number of nodes.

        k:
            The number of inputs per node, at most the number of nodes.

        bias:
            The probability that an entry of a truth table is 1.

        random_state:
            The seed or generator of the wiring, the tables, and the initial values.

    Returns:
        system:
            The model, with the values `node_<i>.x`.

    Raises:
        ValidationError:
            `k` is not between 1 and the number of nodes.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import random_boolean_network
        >>> len(random_boolean_network(10, 3).interactions[0].table)
        8
    """
    if not 1 <= k <= n_nodes:
        error_msg = f'k = {k} inputs per node must be between 1 and the number of nodes, n_nodes = {n_nodes}.'
        raise ValidationError(error_msg)
    rng = np.random.default_rng(random_state)
    nodes = [BooleanNodeElement(f'node_{index}', x=int(rng.integers(2))) for index in range(n_nodes)]
    functions = [
        BooleanFunctionInteraction(
            f'function_{index}',
            target=f'node_{index}',
            inputs=[f'node_{source}' for source in rng.choice(n_nodes, size=k, replace=False)],
            table=tuple(int(entry) for entry in rng.random(2**k) < bias),
        )
        for index in range(n_nodes)
    ]
    return System(elements=nodes, interactions=functions)

repressilator(proteins=(20.0, 0.0, 0.0), *, alpha=40.0, alpha0=0.04, K=20.0, n=2.0, b=2.0, beta=0.2, delta=1.0, stochastic=False)

Return the repressilator of Elowitz and Leibler, three genes in a ring that each repress the next.

Gene i is repressed by gene i - 1, so the ring has an odd number of repressions and oscillates when transcription is strong and repression steep. Counts are molecule numbers, so the same rates run as differential equations or, with stochastic=True, as a jump process of single transcriptions, translations, and decays.

Parameters:

Name	Type	Description	Default
proteins	Sequence[float]	The initial protein counts of the three genes, with no mRNA.	(20.0, 0.0, 0.0)
alpha	float	The transcription rate without repressor.	40.0
alpha0	float	The leaky transcription rate.	0.04
K	float	The repressor count at half repression.	20.0
n	float	The Hill coefficient of repression.	2.0
b	float	The translation rate per mRNA molecule.	2.0
beta	float	The decay rate of the proteins.	0.2
delta	float	The decay rate of mRNA.	1.0
stochastic	bool	Whether the model is a jump process.	False

Returns:

Name	Type	Description
system	System	The model, with the counts gene_<i>.m and gene_<i>.p.

Examples:

>>> from skcomplex.models.biology.systems_biology import repressilator
>>> [interaction.modifiers for interaction in repressilator().interactions[:3]]
[{('gene_2', 'p'): 1}, {('gene_0', 'p'): 1}, {('gene_1', 'p'): 1}]

Source code in src/skcomplex/models/biology/systems_biology/_genes.py

@register_model(
    field='biology/systems_biology',
    tags=('gene regulation', 'oscillations', 'synthetic biology', 'jump process'),
    references=('Elowitz and Leibler 2000', 'arXiv:2111.04508'),
    defaults={'until': 200.0, 'dt': 0.5},
)
def repressilator(
    proteins: Sequence[float] = (20.0, 0.0, 0.0),
    *,
    alpha: float = 40.0,
    alpha0: float = 0.04,
    K: float = 20.0,
    n: float = 2.0,
    b: float = 2.0,
    beta: float = 0.2,
    delta: float = 1.0,
    stochastic: bool = False,
) -> System:
    """Return the repressilator of Elowitz and Leibler, three genes in a ring that each repress the next.

    Gene `i` is repressed by gene `i - 1`, so the ring has an odd number of repressions and oscillates when
    transcription is strong and repression steep. Counts are molecule numbers, so the same rates run as differential
    equations or, with `stochastic=True`, as a jump process of single transcriptions, translations, and decays.

    Args:
        proteins:
            The initial protein counts of the three genes, with no mRNA.

        alpha:
            The transcription rate without repressor.

        alpha0:
            The leaky transcription rate.

        K:
            The repressor count at half repression.

        n:
            The Hill coefficient of repression.

        b:
            The translation rate per mRNA molecule.

        beta:
            The decay rate of the proteins.

        delta:
            The decay rate of mRNA.

        stochastic:
            Whether the model is a jump process.

    Returns:
        system:
            The model, with the counts `gene_<i>.m` and `gene_<i>.p`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import repressilator
        >>> [interaction.modifiers for interaction in repressilator().interactions[:3]]
        [{('gene_2', 'p'): 1}, {('gene_0', 'p'): 1}, {('gene_1', 'p'): 1}]
    """
    mode = 'jump' if stochastic else 'rate'
    expression = {'b': b, 'beta': beta, 'delta': delta}
    repression = {'alpha': alpha, 'alpha0': alpha0, 'K': K, 'n': n}
    return _circuit(proteins, [(2, 0), (0, 1), (1, 2)], expression, repression, mode)

three_gene_boolean_network(x=(1, 0, 0))

Return the three-gene Boolean network of Raman, updated synchronously.

The genes follow g1* = g2 AND NOT g3, g2* = NOT g1, and g3* = NOT g3 OR g2. Its eight states fall into two basins: the point attractor 011 and the cycle 001 -> 010 -> 111 of period three, which the state 100 enters.

Parameters:

Name	Type	Description	Default
x	Sequence[int]	The initial values of g1, g2, and g3.	(1, 0, 0)

Returns:

Name	Type	Description
system	System	The model, with the values g1.x, g2.x, and g3.x.

Examples:

>>> from skcomplex.models.biology.systems_biology import three_gene_boolean_network
>>> three_gene_boolean_network().interactions[2].table
(1, 0, 1, 1)

Source code in src/skcomplex/models/biology/systems_biology/_boolean_models.py

@register_model(
    field='biology/systems_biology',
    tags=('boolean network', 'gene regulation', 'discrete map', 'attractors'),
    references=(RAMAN,),
    defaults={'until': 10, 'dt': 1},
)
def three_gene_boolean_network(x: Sequence[int] = (1, 0, 0)) -> System:
    """Return the three-gene Boolean network of Raman, updated synchronously.

    The genes follow `g1* = g2 AND NOT g3`, `g2* = NOT g1`, and `g3* = NOT g3 OR g2`. Its eight states fall into two
    basins: the point attractor `011` and the cycle `001 -> 010 -> 111` of period three, which the state `100` enters.

    Args:
        x:
            The initial values of `g1`, `g2`, and `g3`.

    Returns:
        system:
            The model, with the values `g1.x`, `g2.x`, and `g3.x`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import three_gene_boolean_network
        >>> three_gene_boolean_network().interactions[2].table
        (1, 0, 1, 1)
    """
    active = [name for name, value in zip(THREE_GENE_RULES, x, strict=True) if value]
    return _rule_network(THREE_GENE_RULES, active)

toggle_switch(proteins=(60.0, 0.0), *, alpha=20.0, alpha0=0.0, K=20.0, n=2.0, b=2.0, beta=0.2, delta=1.0, stochastic=False)

Return the genetic toggle switch of Gardner, Cantor, and Collins, two genes that repress each other.

With steep and strong repression the switch is bistable: one gene is expressed and silences the other, and the initial protein counts decide which. The stochastic mode lets noise flip the switch.

Parameters:

Name	Type	Description	Default
proteins	Sequence[float]	The initial protein counts of the two genes, with no mRNA.	(60.0, 0.0)
alpha	float	The transcription rate without repressor.	20.0
alpha0	float	The leaky transcription rate.	0.0
K	float	The repressor count at half repression.	20.0
n	float	The Hill coefficient of repression.	2.0
b	float	The translation rate per mRNA molecule.	2.0
beta	float	The decay rate of the proteins.	0.2
delta	float	The decay rate of mRNA.	1.0
stochastic	bool	Whether the model is a jump process.	False

Returns:

Name	Type	Description
system	System	The model, with the counts gene_0.m, gene_0.p, gene_1.m, and gene_1.p.

Examples:

>>> from skcomplex.models.biology.systems_biology import toggle_switch
>>> [element.p for element in toggle_switch((0.0, 50.0)).elements]
[0.0, 50.0]

Source code in src/skcomplex/models/biology/systems_biology/_genes.py

@register_model(
    field='biology/systems_biology',
    tags=('gene regulation', 'bistability', 'synthetic biology', 'jump process'),
    references=('Gardner, Cantor, and Collins 2000', 'arXiv:2111.04508'),
    defaults={'until': 100.0, 'dt': 0.5},
)
def toggle_switch(
    proteins: Sequence[float] = (60.0, 0.0),
    *,
    alpha: float = 20.0,
    alpha0: float = 0.0,
    K: float = 20.0,
    n: float = 2.0,
    b: float = 2.0,
    beta: float = 0.2,
    delta: float = 1.0,
    stochastic: bool = False,
) -> System:
    """Return the genetic toggle switch of Gardner, Cantor, and Collins, two genes that repress each other.

    With steep and strong repression the switch is bistable: one gene is expressed and silences the other, and the
    initial protein counts decide which. The stochastic mode lets noise flip the switch.

    Args:
        proteins:
            The initial protein counts of the two genes, with no mRNA.

        alpha:
            The transcription rate without repressor.

        alpha0:
            The leaky transcription rate.

        K:
            The repressor count at half repression.

        n:
            The Hill coefficient of repression.

        b:
            The translation rate per mRNA molecule.

        beta:
            The decay rate of the proteins.

        delta:
            The decay rate of mRNA.

        stochastic:
            Whether the model is a jump process.

    Returns:
        system:
            The model, with the counts `gene_0.m`, `gene_0.p`, `gene_1.m`, and `gene_1.p`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import toggle_switch
        >>> [element.p for element in toggle_switch((0.0, 50.0)).elements]
        [0.0, 50.0]
    """
    mode = 'jump' if stochastic else 'rate'
    expression = {'b': b, 'beta': beta, 'delta': delta}
    repression = {'alpha': alpha, 'alpha0': alpha0, 'K': K, 'n': n}
    return _circuit(proteins, [(1, 0), (0, 1)], expression, repression, mode)

yeast_cell_cycle(active=YEAST_START)

Return the Boolean network of the budding yeast cell cycle of Li et al., eleven proteins and their regulation.

Every protein turns on or off with the sign of the summed activations and inhibitions by the active proteins, and the proteins Cln3, Cln1,2, Swi5, Cdc20, and Mcm1 decay without input. Started from the G1 state with Cln3 on, the cell size signal, the network runs through the phases of the cycle into the stationary G1 state.

Parameters:

Name	Type	Description	Default
active	Sequence[str]	The proteins that are on at the start, among Cln3, MBF, SBF, Cln1_2, Cdh1, Swi5, Cdc20, Clb5_6, Sic1, Clb1_2, and Mcm1.	YEAST_START

Returns:

Name	Type	Description
system	System	The model, with the values <protein>.x.

Examples:

>>> from skcomplex.models.biology.systems_biology import yeast_cell_cycle
>>> [node.label for node in yeast_cell_cycle().elements if node.x]
['Cln3', 'Cdh1', 'Sic1']

Source code in src/skcomplex/models/biology/systems_biology/_boolean.py

@register_model(
    field='biology/systems_biology',
    tags=('boolean network', 'gene regulation', 'discrete map', 'cell cycle'),
    references=('Li et al. 2004', 'Gros 2015'),
    defaults={'until': 15, 'dt': 1},
)
def yeast_cell_cycle(active: Sequence[str] = YEAST_START) -> System:
    """Return the Boolean network of the budding yeast cell cycle of Li et al., eleven proteins and their regulation.

    Every protein turns on or off with the sign of the summed activations and inhibitions by the active proteins,
    and the proteins Cln3, Cln1,2, Swi5, Cdc20, and Mcm1 decay without input. Started from the G1 state with Cln3 on,
    the cell size signal, the network runs through the phases of the cycle into the stationary G1 state.

    Args:
        active:
            The proteins that are on at the start, among `Cln3`, `MBF`, `SBF`, `Cln1_2`, `Cdh1`, `Swi5`, `Cdc20`,
            `Clb5_6`, `Sic1`, `Clb1_2`, and `Mcm1`.

    Returns:
        system:
            The model, with the values `<protein>.x`.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import yeast_cell_cycle
        >>> [node.label for node in yeast_cell_cycle().elements if node.x]
        ['Cln3', 'Cdh1', 'Sic1']
    """
    nodes = [BooleanNodeElement(name, x=int(name in active)) for name in YEAST_NODES]
    functions = [
        ThresholdFunctionInteraction(
            f'regulation_{name}',
            target=name,
            inputs=[source for source, _ in YEAST_INPUTS[name]],
            weights=[float(weight) for _, weight in YEAST_INPUTS[name]],
            decays=name in YEAST_DECAYING,
        )
        for name in YEAST_NODES
    ]
    return System(elements=nodes, interactions=functions)