Systems biology

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

AtomMapSpec(reaction, reactants, products) dataclass

Bases: BaseSpec

Atom transitions of a reaction, a string of atom letters for every occurrence of a metabolite.

Parameters:

Name	Type	Description	Default
reaction	str	The name of the reaction.	required
reactants	Sequence[tuple[str, str]]	Every occurrence of a reactant with its letters, as pairs of a metabolite and a string.	required
products	Sequence[tuple[str, str]]	Every occurrence of a product with its letters, as pairs of a metabolite and a string.	required

Raises:

Type	Description
ValidationError	A product names a letter that no reactant carries, or a reactant repeats a letter, raised by prepare.

Examples:

>>> from skcomplex.systems.biology.systems_biology import AtomMapSpec
>>> atoms = AtomMapSpec('v1', [('AcCoA', 'ab'), ('OAC', 'cdef')], [('Cit', 'fedbac')])
>>> atoms.origins[('Cit', 1)], atoms.origins[('Cit', 5)]
(('OAC', 4), ('AcCoA', 1))

origins property

Reactant atom of every product atom, as pairs of a metabolite and an atom counted from 1.

BackboneGeometry(structure, system=0)

Bases: BaseConfig[None, Array]

Chain of a structure rebuilt from the angles phi and psi of its backbone, with all other geometry kept.

The angles of a chain of n residues are one array of 2 (n - 1) values in degrees: phi of the residues from the second to the last, then psi of the residues from the first to the one before the last. The bond lengths, the bond angles, omega, and the side chains stay as they are in the system. Biopython warns on every rebuild that NumPy reads where without out, and the class filters this one warning.

Parameters:

Name	Type	Description	Default
structure	StructureSpec	The structure, of one chain.	required
system	int	The index of the system that gives the geometry.	0

Attributes:

Name	Type	Description
angles_		The angles of the system, in degrees.
positions_		The positions of the atoms after the last run, atoms by 3.

Raises:

Type	Description
ValidationError	The structure holds more than one chain, or the angles have another length.

Examples:

>>> import numpy as np
>>> from skcomplex.systems.biology.systems_biology import BackboneGeometry, read_structure
>>> geometry = BackboneGeometry(read_structure('1L2Y')).prepare()
>>> geometry.angles_.shape
(38,)
>>> geometry.run(geometry.angles_).positions_.shape
(304, 3)

Source code in src/skcomplex/systems/biology/systems_biology/_backbone_geometry.py

def __init__(self, structure: StructureSpec, system: int = 0) -> None:
    self.structure = structure
    self.system = system

read_effective_atoms()

Return the mean atom of every residue of the last run.

Returns:

Name	Type	Description
positions	Array	Residues by 3.

Examples:

>>> from skcomplex.systems.biology.systems_biology import BackboneGeometry, read_structure
>>> geometry = BackboneGeometry(read_structure('1L2Y')).prepare()
>>> geometry.run(geometry.angles_).read_effective_atoms().shape
(20, 3)

Source code in src/skcomplex/systems/biology/systems_biology/_backbone_geometry.py

def read_effective_atoms(self) -> Array:
    """Return the mean atom of every residue of the last run.

    Returns:
        positions:
            Residues by 3.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import BackboneGeometry, read_structure
        >>> geometry = BackboneGeometry(read_structure('1L2Y')).prepare()
        >>> geometry.run(geometry.angles_).read_effective_atoms().shape
        (20, 3)
    """
    names = self.structure
    return StructureSpec(
        self.positions_[None],
        names.atom_names,
        names.residue_names,
        names.residue_numbers,
        names.chain_names,
    ).effective_atoms[0]

BackboneVisual(position, structure=None, radius=None, states=(), *, trail=0, strip=True, path=None, recording=None)

Bases: PointsVisual

Visual that logs the effective atoms of a chain as points in three dimensions, joined by a line strip in order.

A state of coordinates holds one row of three values per residue. With a structure, the state holds the angles of a BackboneGeometry instead, and the visual rebuilds the chain from them at every frame.

Parameters:

Name	Type	Description	Default
position	str	The name of the lattice state of coordinates or angles.	required
structure	StructureSpec | None	The structure whose chain the angles rebuild, or None for a state of coordinates.	None
radius	float | None	The radius of the points, or None for a share of the extent of the positions.	None
states	Sequence[str] | None	The names of the other states to log as time series.	()
trail	int	The number of recent positions the trail of each point holds, or 0 for no trail.	0
strip	bool	Whether to join the points in member order by one line strip.	True
path	str | None	The entity path of the group to draw, which also names it in the problem.	None
recording	RecordingStream | None	The Rerun recording to log to, or None for the active recording of the thread.	None

Attributes:

Name	Type	Description
geometry_		The chain the angles rebuild, or None for a state of coordinates.

Examples:

>>> from skcomplex.systems.biology.systems_biology import BackboneVisual
>>> BackboneVisual(position='x').strip
True

Source code in src/skcomplex/systems/biology/systems_biology/_fuzzy_oil_drop.py

def __init__(
    self,
    position: str,
    structure: StructureSpec | None = None,
    radius: float | None = None,
    states: Sequence[str] | None = (),
    *,
    trail: int = 0,
    strip: bool = True,
    path: str | None = None,
    recording: rr.RecordingStream | None = None,
) -> None:
    super().__init__(position, radius, states=states, trail=trail, strip=strip, path=path, recording=recording)
    self.structure = structure

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.systems.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/systems/biology/common/networks/_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.systems.biology.systems_biology import BooleanNodeElement
>>> BooleanNodeElement('Sic1', x=1).x
1

Source code in src/skcomplex/systems/biology/common/networks/_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

ChainCoordinateElement(label, x=0.0)

Bases: Element

Coordinate of an effective atom of a protein chain.

Parameters:

Name	Type	Description	Default
label	str	The label of the element.	required
x	Annotated[float, StateSpec()]	The coordinate in ångström.	0.0

Examples:

>>> from skcomplex.systems.biology.systems_biology import ChainCoordinateElement
>>> ChainCoordinateElement('site').x
0.0

Source code in src/skcomplex/systems/biology/systems_biology/_fuzzy_oil_drop.py

def __init__(self, label: str, x: Annotated[float, StateSpec()] = 0.0) -> None:
    super().__init__(label)
    self.x = x

DihedralAngleElement(label, angle=0.0)

Bases: Element

Backbone dihedral angle of a protein chain.

Parameters:

Name	Type	Description	Default
label	str	The label of the element.	required
angle	Annotated[float, StateSpec()]	The angle in degrees.	0.0

Examples:

>>> from skcomplex.systems.biology.systems_biology import DihedralAngleElement
>>> DihedralAngleElement('site').angle
0.0

Source code in src/skcomplex/systems/biology/systems_biology/_fuzzy_oil_drop.py

def __init__(self, label: str, angle: Annotated[float, StateSpec()] = 0.0) -> None:
    super().__init__(label)
    self.angle = angle

EllipsePath(long_diagonal=LONG_DIAGONAL, short_diagonal=SHORT_DIAGONAL, cells=72)

Bases: BaseConfig[None, Array]

Ellipse path on the Ramachandran map, with the projection of dihedral angles by the shortest periodic distance.

The path is phi = A cos(t) - B sin(t), psi = -A cos(t) - B sin(t), a centred ellipse whose long diagonal runs from the upper left to the lower right of the map, with t = 0 in the lower right and t growing clockwise. Konieczny et al. (2023 §3.7.1, Fig. 3.43d and e) draw this path and print no equation or number for it. The defaults, 83° and 58°, are readings of a figure and not published numbers: a fit to panel e of that figure gives 84.5° and 56.3° with the centre at the origin, which agrees within the 2.7° of one of its pixels. The origin and the direction of t are a convention of this class, since the book states neither. The equation and the parameters of Brylinski et al. (2004, Bioinformatics 20:199) have not been read.

Parameters:

Name	Type	Description	Default
long_diagonal	float	The value A in radians.	LONG_DIAGONAL
short_diagonal	float	The value B in radians.	SHORT_DIAGONAL
cells	int	The number of equal cells of the parameter t.	72

Attributes:

Name	Type	Description
parameters_		The parameter t of the projection of every pair of the last run, in [0, 2 pi).
angles_		The projected phi and psi of the last run in radians, pairs by 2.
cells_		The cell of the path of every pair of the last run.

Raises:

Type	Description
ValidationError	A diagonal or the number of cells is not positive, or the angles are not pairs.

Examples:

>>> import numpy as np
>>> from skcomplex.systems.biology.systems_biology import EllipsePath
>>> path = EllipsePath().prepare()
>>> path.read_points(np.array([0.0])).shape
(1, 2)
>>> path.run(np.array([[1.0, -1.0]])).cells_.tolist()
[0]

Source code in src/skcomplex/systems/biology/systems_biology/_early_stage.py

def __init__(
    self,
    long_diagonal: float = LONG_DIAGONAL,
    short_diagonal: float = SHORT_DIAGONAL,
    cells: int = 72,
) -> None:
    self.long_diagonal = long_diagonal
    self.short_diagonal = short_diagonal
    self.cells = cells

read_points(parameters)

Return the points of the path at some values of its parameter.

Parameters:

Name	Type	Description	Default
parameters	Array	The values of t.	required

Returns:

Name	Type	Description
points	Array	The phi and psi in radians, values by 2.

Examples:

>>> import numpy as np
>>> from skcomplex.systems.biology.systems_biology import EllipsePath
>>> phi, psi = EllipsePath(1.0, 0.5).read_points(np.array([0.0]))[0]
>>> float(phi), float(psi)
(1.0, -1.0)

Source code in src/skcomplex/systems/biology/systems_biology/_early_stage.py

def read_points(self, parameters: Array) -> Array:
    """Return the points of the path at some values of its parameter.

    Args:
        parameters:
            The values of `t`.

    Returns:
        points:
            The `phi` and `psi` in radians, values by 2.

    Examples:
        >>> import numpy as np
        >>> from skcomplex.systems.biology.systems_biology import EllipsePath
        >>> phi, psi = EllipsePath(1.0, 0.5).read_points(np.array([0.0]))[0]
        >>> float(phi), float(psi)
        (1.0, -1.0)
    """
    along, across = self.long_diagonal * np.cos(parameters), self.short_diagonal * np.sin(parameters)
    return np.column_stack([along - across, -along - across])

EnsembleReplayInteraction(label, *, target, state, ensemble)

Bases: Interaction

Drive that sets a lattice state to one member of an ensemble per unit of time, holding the last after the end.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the lattice group.	required
state	str	The name of the state.	required
ensemble	Array	The members, members by the lattice shape.	required

Examples:

>>> import numpy as np
>>> from skcomplex.systems.biology.systems_biology import EnsembleReplayInteraction
>>> EnsembleReplayInteraction('replay', target='chain', state='x', ensemble=np.zeros((2, 3))).state
'x'

Source code in src/skcomplex/systems/biology/systems_biology/_fuzzy_oil_drop.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    state: str,
    ensemble: Array,
) -> None:
    super().__init__(label)
    self.target = target
    self.state = state
    self.ensemble = ensemble

EnzymeNodeElement(label, visual=None, *, X=0.1)

Bases: Element

Enzyme of a signalling circuit, the active fraction X of a pool that other enzymes activate and deactivate.

Parameters:

Name	Type	Description	Default
label	str	The label of the enzyme.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
X	Annotated[float, StateSpec(domain=closed(0, 1))]	The initial active fraction.	0.1

Examples:

>>> from skcomplex.systems.biology.systems_biology import EnzymeNodeElement
>>> EnzymeNodeElement('A', X=0.2).X
0.2

Source code in src/skcomplex/systems/biology/systems_biology/_adaptation.py

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

EnzymeRegulationInteraction(label, *, target, source, sign=1, k=1.0, K=0.5)

Bases: Interaction

Michaelis-Menten activation or deactivation of an enzyme by another enzyme of a circuit.

An activating source converts the inactive fraction of the target at the rate k X_s (1 - X) / (1 - X + K), and a deactivating source the active fraction at the rate k X_s X / (X + K), the kinetics of Ma et al. (2009).

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the regulated enzyme.	required
source	Annotated[str, RoleSpec()]	The path of the regulating enzyme.	required
sign	int	1 for activation or -1 for deactivation.	1
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The catalytic rate.	1.0
K	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant.	0.5

Raises:

Type	Description
ValidationError	The sign is neither 1 nor -1.

Examples:

>>> from skcomplex.systems.biology.systems_biology import EnzymeRegulationInteraction
>>> EnzymeRegulationInteraction('ab', target='B', source='A', sign=-1).sign
-1

Source code in src/skcomplex/systems/biology/systems_biology/_adaptation.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    source: Annotated[str, RoleSpec()],
    sign: int = 1,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
) -> None:
    super().__init__(label)
    self.target = target
    self.source = source
    self.sign = sign
    self.k = k
    self.K = K

EnzymeStimulusInteraction(label, *, target, sign=1, k=1.0, K=0.5, level=0.5, change=0.0, at=0.0)

Bases: Interaction

Activation or deactivation of an enzyme by an outside stimulus or a basal enzyme of fixed concentration.

The regulator has the concentration level, which rises by change at the time at, as the input I of a circuit, and stays fixed for a basal enzyme F with no change. It activates the target at the rate k I (1 - X) / (1 - X + K) or deactivates it at the rate k I X / (X + K).

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the regulated enzyme.	required
sign	int	1 for activation or -1 for deactivation.	1
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The catalytic rate.	1.0
K	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Michaelis constant.	0.5
level	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The concentration of the regulator before the change.	0.5
change	float	The step of the concentration at the time at.	0.0
at	float	The time of the step.	0.0

Examples:

>>> from skcomplex.systems.biology.systems_biology import EnzymeStimulusInteraction
>>> EnzymeStimulusInteraction('input', target='A', level=0.5, change=0.1, at=300.0).change
0.1

Source code in src/skcomplex/systems/biology/systems_biology/_adaptation.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    sign: int = 1,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    level: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    change: float = 0.0,
    at: float = 0.0,
) -> None:
    super().__init__(label)
    self.target = target
    self.sign = sign
    self.k = k
    self.K = K
    self.level = level
    self.change = change
    self.at = at

FateLandscapeElement(label, visual=None, *, x=0.0, y=0.0, a=-0.2, theta1=0.0, theta2=0.0, v1=0.0, v2=0.0)

Bases: Element

Cell state in the compactified elliptic umbilic, the minimal landscape of a decision among three fates.

The state (x, y) descends the potential f = x^4 + y^4 + x^3 - 2 x y^2 + a (x^2 + y^2) + theta1 x + theta2 y of Rand et al. (2021, figure 5D), dx/dt = -df/dx and dy/dt = -df/dy. The signals theta1 and theta2 change in time at the rates v1 and v2, and the sign of a makes the central rest point an attractor or a repellor.

Parameters:

Name	Type	Description	Default
label	str	The label of the cell or of the group of cells.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial first coordinate.	0.0
y	Annotated[float, StateSpec()]	The initial second coordinate.	0.0
a	Annotated[float, ParameterSpec()]	The curvature of the potential at the centre.	-0.2
theta1	Annotated[float, ParameterSpec()]	The first signal at the start.	0.0
theta2	Annotated[float, ParameterSpec()]	The second signal at the start.	0.0
v1	Annotated[float, ParameterSpec()]	The rate of change of the first signal.	0.0
v2	Annotated[float, ParameterSpec()]	The rate of change of the second signal.	0.0

Examples:

>>> from skcomplex.systems.biology.systems_biology import FateLandscapeElement
>>> FateLandscapeElement('cell', a=-0.2).a
-0.2

Source code in src/skcomplex/systems/biology/systems_biology/_fates.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    y: Annotated[float, StateSpec()] = 0.0,
    a: Annotated[float, ParameterSpec()] = -0.2,
    theta1: Annotated[float, ParameterSpec()] = 0.0,
    theta2: Annotated[float, ParameterSpec()] = 0.0,
    v1: Annotated[float, ParameterSpec()] = 0.0,
    v2: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.a = a
    self.theta1 = theta1
    self.theta2 = theta2
    self.v1 = v1
    self.v2 = v2

FragmentSpec(metabolite, atoms, name='') dataclass

Bases: BaseSpec

Fragment of a metabolite whose mass isotopomer distribution a simulation gives or an experiment measures.

Parameters:

Name	Type	Description	Default
metabolite	str	The metabolite.	required
atoms	Sequence[int]	The atoms of the fragment, counted from 1.	required
name	str	The name of the fragment, or the empty string for the name of the metabolite.	''

Examples:

>>> from skcomplex.systems.biology.systems_biology import FragmentSpec
>>> FragmentSpec('Glu', (1, 2, 3, 4, 5)).label
'Glu'

label property

Name of the fragment, the name of its metabolite without one of its own.

FuzzyOilDropAnalysis(structure, kind='relative-distance', scale=None, cutoff=CUTOFF, *, observed=None, state=None)

Bases: BaseAnalysis

Distributions of hydrophobicity of a structure and their relative distance, after the fuzzy oil drop system.

The theoretical distribution T is a three-dimensional Gaussian over the effective atoms along the principal axes of the chain (Eq. 1.1). The observed distribution O is Levitt's sum over the pairs of residues within a cutoff (Eq. 1.2). The uniform distribution R gives every residue the same value. Then RD = D_KL(O|T) / (D_KL(O|T) + D_KL(O|R)) in bits (Eqs. 1.3 to 1.6), and K minimises D_KL(O|M) for the profile M = [T + K [T_max - T]_n]_n (Eq. 1.9), where [...]_n normalises to a sum of one, by SciPy's bounded scalar minimiser. The equations follow Konieczny et al. (2023 §1.3.1), and the curves of Fig. 1.41 give back its K = 0.5. Each sigma is a third of the largest absolute coordinate along its axis, one reading of the book's Gaussian that reaches near zero at 3 sigma (p. 33 and Fig. 1.41), which gives no formula for sigma. The book names no scale, no orientation of the axes, no fit of K, and no upper bound of it, for which it reports values above 3 (p. 41), so the Kyte-Doolittle scale, the principal axes, the minimisation, and the bound 5 are choices of this class, as is leaving a residue out of its own sum in Eq. 1.2, where the book speaks of interactions between residues.

The bundled entries do not reach the RD < 0.5 and K below 0.5 that the book expects of fast-folding and single-chain proteins (p. 40): the mean RD and K are 0.77 and 2.2 for 1L2Y, 0.70 and 0.98 for 1UBQ, and 0.77 and 2.2 for 5PTI. The gap comes from sigma, since a wider Gaussian closes it and no other choice does.

Parameters:

Name	Type	Description	Default
structure	StructureSpec	The structure that gives the residues, and the systems when the source is None.	required
kind	str	'profile', 'relative-distance', or 'environment-factor'.	'relative-distance'
scale	Mapping[str, float] | None	The hydrophobicity by the three-letter name of a residue, or None for the Kyte-Doolittle scale rescaled to [0, 1].	None
cutoff	float	The cutoff distance of the observed distribution in ångström.	CUTOFF
observed	ArrayLike | None	The observed distribution to use in place of Levitt's sum, one value per residue, or None.	None
state	str | None	The dotted name of the lattice state of a trajectory that holds the coordinates, or None for the first state.	None

Attributes:

Name	Type	Description
positions_		The effective atoms, frames by residues by 3.
theoretical_		The distribution T, frames by residues.
observed_		The distribution O, frames by residues.
uniform_		The distribution R, one value per residue.
divergences_		D_KL(O|T) and D_KL(O|R) in bits, frames by 2.
values_		For the kind 'relative-distance' the RD, and for the kind 'environment-factor' the K, per frame.

Raises:

Type	Description
AnalysisError	The kind is unknown, the cutoff is not positive, the structure holds a residue the scale does not cover, or the positions do not hold one row per residue.

Examples:

>>> from skcomplex.systems.biology.systems_biology import FuzzyOilDropAnalysis, read_structure
>>> values = FuzzyOilDropAnalysis(read_structure('1L2Y')).prepare().run().values_
>>> values.shape, bool(values.min() >= 0 and values.max() <= 1)
((38,), True)

Source code in src/skcomplex/systems/biology/systems_biology/_fuzzy_oil_drop.py

def __init__(
    self,
    structure: StructureSpec,
    kind: str = 'relative-distance',
    scale: Mapping[str, float] | None = None,
    cutoff: float = CUTOFF,
    *,
    observed: npt.ArrayLike | None = None,
    state: str | None = None,
) -> None:
    self.structure = structure
    self.kind = kind
    self.scale = scale
    self.cutoff = cutoff
    self.observed = observed
    self.state = state

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.systems.biology.systems_biology import GProteinElement
>>> GProteinElement('g_protein').k_gtp
0.5

Source code in src/skcomplex/systems/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.systems.biology.systems_biology import GeneElement
>>> GeneElement('lacI', p=20.0).p
20.0

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

GlucoseSensorInteraction(label, *, node='node', source='source', threshold=None)

Bases: Interaction

Boolean node that reads whether a substrate lies above a threshold, or the previous node of a chain.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
node	Annotated[str, RoleSpec()]	The path of the Boolean node.	'node'
source	Annotated[str, RoleSpec()]	The path of what the node reads: a substrate with a threshold, or a Boolean node without one.	'source'
threshold	float | None	The concentration above which the node turns on, or None for a node that copies another.	None

Examples:

>>> from skcomplex.systems.biology.systems_biology import GlucoseSensorInteraction
>>> GlucoseSensorInteraction('sensor', node='a', source='glucose', threshold=0.1).threshold
0.1

Source code in src/skcomplex/systems/biology/systems_biology/_regulatory_fba.py

def __init__(
    self,
    label: str,
    *,
    node: Annotated[str, RoleSpec()] = 'node',
    source: Annotated[str, RoleSpec()] = 'source',
    threshold: float | None = None,
) -> None:
    super().__init__(label)
    self.node = node
    self.source = source
    self.threshold = threshold

GlycolysisElement(label, visual=None, *, variables=None)

Bases: DeclaredElement

Yeast glycolysis in a flow reactor, the fifteen metabolites of the model of Nielsen and coworkers.

The model of Nielsen, Sørensen, Hynne, and Busse (1998), BioModels BIOMD0000000042, follows glucose GLC through hexokinase, phosphofructokinase with its inhibition by the ratio of ATP to AMP, aldolase, glyceraldehyde phosphate dehydrogenase, phosphoglycerate kinase, pyruvate kinase, pyruvate decarboxylase, alcohol dehydrogenase, adenylate kinase, and a glycerol branch to P. The reactor is fed with ATP, ADP, NADH, NAD, and glucose at fixed concentrations and washed out at the rate flow. The concentrations, in mM, and the rate constants, in mM and per minute, are declared variables named as in the source, and variables replaces any of them.

Parameters:

Name	Type	Description	Default
label	str	The label of the cell extract.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
variables	Mapping[str, DeclaredVariableSpec] | None	The declared states and parameters that replace those of the source, by name, or None for none.	None

Examples:

>>> from skcomplex.systems.biology.systems_biology import GlycolysisElement
>>> GlycolysisElement('extract').get_value('V4')
10.0

Source code in src/skcomplex/core/components/elements/_declared.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    variables: Mapping[str, DeclaredVariableSpec] | None = None,
) -> None:
    super().__init__(label, visual=visual)
    self.variables = variables

LabelMeasurementSpec(experiment, fragment, values, deviation=0.002) dataclass

Bases: BaseSpec

Measured mass isotopomer distribution of a fragment in one tracer experiment.

Parameters:

Name	Type	Description	Default
experiment	str	The name of the tracer experiment.	required
fragment	str	The name of the fragment.	required
values	ArrayLike	The measured distribution.	required
deviation	float	The standard deviation of every value.	0.002

Examples:

>>> from skcomplex.systems.biology.systems_biology import LabelMeasurementSpec
>>> LabelMeasurementSpec('a', 'Glu', (0.5, 0.5), 0.002).deviation
0.002

distribution property

Measured distribution as an array.

LabellingAnalysis(kind='simulation', *, atoms=(), tracers=None, fragments=(), measurements=(), fluxes=None, measured_fluxes=None, symmetric=(), starts=5, intervals='linear', random_state=None)

Bases: BaseAnalysis

Isotope labelling of a metabolic network at metabolic and isotopic steady state, of one kind.

The analysis takes a system that holds a metabolic network. It walks the atom maps backward from the fragments to collect the elementary metabolite units (Antoniewicz et al. 2007), and solves per unit size one sparse linear system with scipy.sparse.linalg.splu.

• 'simulation': the mass isotopomer distributions of the fragments for given fluxes.
• 'estimation': the fluxes by weighted least squares over measured distributions and measured fluxes, from seeded feasible starts that the core solve draws. The net fluxes run over a null space basis of the balances, with one exchange flux per reversible reaction, so every balance holds exactly. A 95% interval of every flux and a chi-square test of the fit follow.

Parameters:

Name	Type	Description	Default
kind	str	'simulation' or 'estimation'.	'simulation'
atoms	Sequence[AtomMapSpec]	The atom map of every reaction that moves traced atoms.	()
tracers	Mapping[str, Sequence[TracerSpec]] | None	The labelled substrates of every parallel experiment, by the name of the experiment.	None
fragments	Sequence[FragmentSpec]	The fragments to simulate or that measurements name.	()
measurements	Sequence[LabelMeasurementSpec]	The measured distributions, for 'estimation'.	()
fluxes	Mapping[str, float | tuple[float, float]] | None	The flux of every reaction for 'simulation', a number, or a pair of a forward and a backward flux for a reversible reaction.	None
measured_fluxes	Mapping[str, tuple[float, float]] | None	The value and the standard deviation of every measured flux for 'estimation', by reaction. A deviation of zero fixes the flux.	None
symmetric	Sequence[str]	The metabolites whose atoms a producing reaction scrambles, half in each direction.	()
starts	int	The number of feasible starts of the fit.	5
intervals	str	'linear' for intervals of the linearised fit, or 'profile' for profiled ones (Antoniewicz et al. 2006).	'linear'
random_state	int | Generator | None	The seed or generator of the starts.	None

Attributes:

Name	Type	Description
distributions_		The distribution of every fragment of every experiment.
fluxes_		The estimated net flux of every reaction, and the exchange <reaction>.exchange of a reversible one.
intervals_		The 95% interval of every estimated flux.
identifiable_		Whether the data identify every estimated flux, false for an interval that reaches a bound.
residual_		The sum of the squares of the weighted residuals at the optimum.
chi2_		The range of the sum of squares that the test accepts at the 5% level.
accepted_		Whether the residual lies in that range.

Raises:

Type	Description
AnalysisError	A setting is unknown or a metabolite has no source of labels, raised by prepare.
ValidationError	An atom map is not valid, or two maps give one metabolite different numbers of atoms, raised by prepare.

Examples:

>>> from skcomplex.core.components.spaces import StaticTime
>>> from skcomplex.core.components.system import System
>>> from skcomplex.systems.biology.common.metabolism import ReactionSpec, metabolic_network_type
>>> from skcomplex.systems.biology.systems_biology import (
...     AtomMapSpec,
...     FragmentSpec,
...     LabellingAnalysis,
...     TracerSpec,
... )
>>> reactions = {'in': ReactionSpec({'B': 1.0}), 'out': ReactionSpec({'B': -1.0})}
>>> network = metabolic_network_type('ChainElement', reactions)('net')
>>> analysis = LabellingAnalysis(
...     atoms=[AtomMapSpec('in', [('A', 'ab')], [('B', 'ba')])],
...     tracers={'first': [TracerSpec('A', {'00': 0.5, '01': 0.5})]},
...     fragments=[FragmentSpec('B', (1, 2))],
...     fluxes={'in': 1.0, 'out': 1.0},
... )
>>> system = System(elements=[network], time=StaticTime())
>>> analysis.prepare(system).run().distributions_['first']['B'].tolist()
[0.5, 0.5, 0.0]

Source code in src/skcomplex/systems/biology/systems_biology/_carbon_13_mfa.py

def __init__(
    self,
    kind: str = 'simulation',
    *,
    atoms: Sequence[AtomMapSpec] = (),
    tracers: Mapping[str, Sequence[TracerSpec]] | None = None,
    fragments: Sequence[FragmentSpec] = (),
    measurements: Sequence[LabelMeasurementSpec] = (),
    fluxes: Mapping[str, float | tuple[float, float]] | None = None,
    measured_fluxes: Mapping[str, tuple[float, float]] | None = None,
    symmetric: Sequence[str] = (),
    starts: int = 5,
    intervals: str = 'linear',
    random_state: int | np.random.Generator | None = None,
) -> None:
    self.kind = kind
    self.atoms = atoms
    self.tracers = tracers
    self.fragments = fragments
    self.measurements = measurements
    self.fluxes = fluxes
    self.measured_fluxes = measured_fluxes
    self.symmetric = symmetric
    self.starts = starts
    self.intervals = intervals
    self.random_state = random_state

LateralInhibitionInteraction(label, *, cells, kernel=((0.0,),))

Bases: Interaction

Lateral inhibition of the neural fate of cells by their neighbours, through a kernel of the distance.

Following Eq. (2) of Rand et al. (2021), da_i/dt = sigma(a_i - h_i) - a_i with the inhibition h_i = sum_{j != i} K_ij d(a_j), where sigma(u) = 1 / (1 + exp(-u / width)) and d(a) = sigma(a - 1 / 2). A cell under an inhibition below 0.37 can only turn neural, above 0.63 only epidermal, and in between it is bistable.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
cells	Annotated[str, RoleSpec()]	The path of the group of cells.	required
kernel	Sequence[Sequence[float]]	The inhibition K_ij of every cell by every other, with a zero diagonal.	((0.0,),)

Examples:

>>> from skcomplex.systems.biology.systems_biology import LateralInhibitionInteraction
>>> LateralInhibitionInteraction('inhibition', cells='cells', kernel=((0.0, 1.0), (1.0, 0.0))).kernel
((0.0, 1.0), (1.0, 0.0))

Source code in src/skcomplex/systems/biology/systems_biology/_fates.py

def __init__(
    self,
    label: str,
    *,
    cells: Annotated[str, RoleSpec()],
    kernel: Sequence[Sequence[float]] = ((0.0,),),
) -> None:
    super().__init__(label)
    self.cells = cells
    self.kernel = kernel

LogicalRegulationInteraction(label, *, target, regulators=(), thresholds=(), parameters=None)

Bases: Interaction

Regulation of a logical variable, which steps one level at a time towards the focal level its regulators set.

In the formalism of Thomas, every regulator acts when its level reaches its threshold. The logical parameters map every combination of acting regulators, a tuple of truth values in the order of the regulators, to the focal level of the target. The target rises or falls by one level towards its focal level at the unit rate, so the transitions of all variables are asynchronous random events.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the regulated variable.	required
regulators	Annotated[Sequence[str], RoleSpec(variadic=True)]	The paths of the regulating variables.	()
thresholds	Sequence[int]	The level at which every regulator acts, in the order of the regulators.	()
parameters	Mapping[tuple[bool, ...], int] | None	The focal level of the target for every combination of acting regulators.	None

Examples:

>>> from skcomplex.systems.biology.systems_biology import LogicalRegulationInteraction
>>> regulation = LogicalRegulationInteraction(
...     'x_regulation', target='x', regulators=['y'], thresholds=(1,), parameters={(False,): 2, (True,): 0}
... )
>>> regulation.parameters[(True,)]
0

Source code in src/skcomplex/systems/biology/systems_biology/_logical_networks.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    regulators: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
    thresholds: Sequence[int] = (),
    parameters: Mapping[tuple[bool, ...], int] | None = None,
) -> None:
    super().__init__(label)
    self.target = target
    self.regulators = regulators
    self.thresholds = thresholds
    self.parameters = parameters

LogicalVariableElement(label, visual=None, *, level=0.0)

Bases: Element

Variable of a generalised logical network, a discrete level from zero to its largest level.

The levels stand for qualitative concentrations, such as absent, basal, medium, or high, cut at the thresholds at which the variable acts on its targets.

Parameters:

Name	Type	Description	Default
label	str	The label of the variable.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
level	Annotated[float, StateSpec(domain=NN)]	The initial level.	0.0

Examples:

>>> from skcomplex.systems.biology.systems_biology import LogicalVariableElement
>>> LogicalVariableElement('x', level=2.0).level
2.0

Source code in src/skcomplex/systems/biology/systems_biology/_logical_networks.py

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

MajorityNodeElement(label, visual=None, *, sigma=1.0, eta=0.0)

Bases: Element

Node of a noisy Boolean network, a state of plus or minus one that follows the majority of its inputs.

Parameters:

Name	Type	Description	Default
label	str	The label of the node or of the group of nodes.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
sigma	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([-1, 1]))]	The state, -1 or 1.	1.0
eta	Annotated[float, ParameterSpec(domain=closed(0, 1))]	The probability that the node takes the state opposite to the majority.	0.0

Examples:

>>> from skcomplex.systems.biology.systems_biology import MajorityNodeElement
>>> MajorityNodeElement('node', eta=0.1).eta
0.1

Source code in src/skcomplex/systems/biology/systems_biology/_boolean_systems.py

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

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 system, 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.systems.biology.systems_biology import MinSkeletonElement
>>> MinSkeletonElement('site').lam
6.0

Source code in src/skcomplex/systems/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

NoisyMajorityInteraction(label, *, nodes, inputs=((0.0,),))

Bases: Interaction

Update of every node of a Boolean network to the sign of the sum of its inputs, flipped with a probability.

Node i takes sign(sum_j sigma_ij) over its inputs with the probability 1 - eta and the opposite sign with the probability eta, all nodes at once, the network of Huepe and Aldana-González (2002).

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
nodes	Annotated[str, RoleSpec()]	The path of the group of nodes.	required
inputs	Sequence[Sequence[float]] | sparray	The input matrix, dense or sparse, one row per node with 1 for every input, an odd number per row.	((0.0,),)

Examples:

>>> from skcomplex.systems.biology.systems_biology import NoisyMajorityInteraction
>>> NoisyMajorityInteraction('majority', nodes='nodes', inputs=((0.0, 1.0), (1.0, 0.0))).nodes
'nodes'

Source code in src/skcomplex/systems/biology/systems_biology/_boolean_systems.py

def __init__(
    self,
    label: str,
    *,
    nodes: Annotated[str, RoleSpec()],
    inputs: Sequence[Sequence[float]] | scipy.sparse.sparray = ((0.0,),),
) -> None:
    super().__init__(label)
    self.nodes = nodes
    self.inputs = inputs

PathCellElement(label, cell=0.0, amino=0.0)

Bases: Element

Residue of a protein chain whose dihedral angles lie in one cell of an ellipse path.

Parameters:

Name	Type	Description	Default
label	str	The label of the element.	required
cell	Annotated[float, StateSpec()]	The cell of the path.	0.0
amino	Annotated[float, ParameterSpec()]	The index of the amino acid of the residue.	0.0

Examples:

>>> from skcomplex.systems.biology.systems_biology import PathCellElement
>>> PathCellElement('residue').cell
0.0

Source code in src/skcomplex/systems/biology/systems_biology/_early_stage.py

def __init__(
    self,
    label: str,
    cell: Annotated[float, StateSpec()] = 0.0,
    amino: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label)
    self.cell = cell
    self.amino = amino

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.systems.biology.systems_biology import PerClockElement
>>> PerClockElement('clock').v_s
0.76

Source code in src/skcomplex/systems/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.systems.biology.systems_biology import PhosphorylationLoopElement
>>> PhosphorylationLoopElement('loop', A=0.5).A
0.5

Source code in src/skcomplex/systems/biology/systems_biology/_networks.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.systems.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/systems/biology/systems_biology/_probabilistic.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.systems.biology.systems_biology import ProbabilisticNodeElement
>>> ProbabilisticNodeElement('x1', x=1).x
1

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

ProneuralCellElement(label, visual=None, *, a=0.0, width=0.127)

Bases: Element

Cell of an epithelium that adopts the neural fate when the inhibition from its neighbours is weak.

The state a runs from 0, the epidermal fate, to 1, the neural fate of a sensory bristle. The sigmoids of the lateral inhibition share the width width.

Parameters:

Name	Type	Description	Default
label	str	The label of the cell or of the group of cells.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
a	Annotated[float, StateSpec(domain=closed(0, 1))]	The initial neural state.	0.0
width	Annotated[float, ParameterSpec(domain=open(0, oo))]	The width of the sigmoids.	0.127

Examples:

>>> from skcomplex.systems.biology.systems_biology import ProneuralCellElement
>>> ProneuralCellElement('cell', a=0.1).a
0.1

Source code in src/skcomplex/systems/biology/systems_biology/_fates.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    a: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.0,
    width: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.127,
) -> None:
    super().__init__(label, visual=visual)
    self.a = a
    self.width = width

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.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.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).prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('mRNA.x')[-1]), 3)
1.0

Source code in src/skcomplex/systems/biology/common/kinetics/_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.systems.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/systems/biology/systems_biology/_networks.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.systems.biology.systems_biology import SpeciesElement
>>> SpeciesElement('camp', x=2.0).x
2.0

Source code in src/skcomplex/systems/biology/common/kinetics/_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

StructureSpec(coordinates, atom_names, residue_names, residue_numbers, chain_names) dataclass

Bases: BaseSpec

Spec of a molecular structure, every system of an ensemble with the same atoms.

Parameters:

Name	Type	Description	Default
coordinates	Array	The positions in ångström, systems by atoms by 3.	required
atom_names	tuple[str, ...]	The name of every atom.	required
residue_names	tuple[str, ...]	The name of the residue of every atom.	required
residue_numbers	tuple[int, ...]	The number of the residue of every atom in its chain.	required
chain_names	tuple[str, ...]	The name of the chain of every atom.	required

Raises:

Type	Description
ValidationError	The coordinates are not systems by atoms by 3, or the names do not hold one entry per atom.

Examples:

>>> import numpy as np
>>> from skcomplex.systems.biology.systems_biology import StructureSpec
>>> atoms = StructureSpec(np.zeros((1, 2, 3)), ('N', 'CA'), ('GLY', 'GLY'), (1, 1), ('A', 'A'))
>>> atoms.residues
(('A', 1, 'GLY'),)

effective_atoms property

Mean atom of every residue, systems by residues by 3.

residues property

Chain, number, and name of every residue, in the order of the atoms.

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.systems.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/systems/biology/systems_biology/_boolean.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

ToyNetworkElement(label, visual=None, *, objective=None, parsimonious=False, infeasible='raise', name='ToyNetwork', reactions=TOY_REACTIONS, default_objective=TOY_OBJECTIVE, variables=None)

Bases: MetabolicNetworkElement

The toy network of Table 8.2 of Raman (2021), which grows on A alone.

Parameters:

Name	Type	Description	Default
label	str	The label of the element.	required
visual	VisualSpec | None	How the element is drawn, or None for the bars of its program.	None
objective	Mapping[str, float] | None	The coefficient of every reaction in the objective, or None for default_objective.	None
parsimonious	bool	Whether to minimise the total absolute flux after the objective.	False
infeasible	str	The infeasibility policy of the program.	'raise'
name	str | None	The name of the network.	'ToyNetwork'
reactions	Mapping[str, ReactionSpec]	The reactions by name.	TOY_REACTIONS
default_objective	Mapping[str, float]	The objective of the network when objective is None.	TOY_OBJECTIVE
variables	Mapping[str, DeclaredVariableSpec] | None	The declared states and parameters that replace the ones of the reactions, by name, or None for none.	None

Examples:

>>> from skcomplex.systems.biology.common.metabolism import ToyNetworkElement
>>> ToyNetworkElement('cell').get_value('upper_v1')
1000.0

Source code in src/skcomplex/systems/biology/common/metabolism/_toy_network.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    objective: Mapping[str, float] | None = None,
    parsimonious: bool = False,
    infeasible: str = 'raise',
    name: str | None = 'ToyNetwork',
    reactions: Mapping[str, ReactionSpec] = TOY_REACTIONS,
    default_objective: Mapping[str, float] = TOY_OBJECTIVE,
    variables: Mapping[str, DeclaredVariableSpec] | None = None,
) -> None:
    super().__init__(
        label,
        visual,
        objective=objective,
        parsimonious=parsimonious,
        infeasible=infeasible,
        name=name,
        reactions=reactions,
        default_objective=default_objective,
        variables=variables,
    )

TracerSpec(metabolite, fractions) dataclass

Bases: BaseSpec

Labelling of a substrate, the fraction of every isotopomer as a bit string with atom 1 first.

Parameters:

Name	Type	Description	Default
metabolite	str	The metabolite.	required
fractions	Mapping[str, float]	The fraction of every isotopomer, such as {'00': 0.5, '01': 0.25, '11': 0.25}, where 1 is a labelled atom. The fractions sum to 1.	required

Raises:

Type	Description
ValidationError	The fractions do not sum to 1, or their bit strings differ in length, raised by prepare.

Examples:

>>> from skcomplex.systems.biology.systems_biology import TracerSpec
>>> TracerSpec('AcCoA', {'00': 0.5, '11': 0.5}).size
2

size property

Number of atoms of the metabolite.

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 system, with the concentrations X1.x, X2.x, and so on.

Raises:

Type	Description
ValidationError	The form is unknown.

Examples:

>>> from skcomplex.systems.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/systems/biology/systems_biology/_networks.py

@register_system(
    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 system, with the concentrations `X1.x`, `X2.x`, and so on.

    Raises:
        ValidationError:
            The form is unknown.

    Examples:
        >>> from skcomplex.systems.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 system, 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.systems.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/systems/biology/systems_biology/_boolean_systems.py

@register_system(
    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 system, with the values `LacI.x`, `TetR.x`, `cI.x`, and the relays `<repressor>_<k>.x`.

    Raises:
        ValidationError:
            The delay is negative.

    Examples:
        >>> from skcomplex.systems.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)

bristle_lateral_inhibition(cells=40, *, strength=1.0, length=2.5, width=0.127, noise=0.05, seed=0)

Return a ring of cells that inhibit the neural fate of their neighbours, leaving sparse, evenly spaced bristles.

Following Eq. (2) of Rand et al. (2021), every cell relaxes to sigma(a_i - h_i) under the inhibition h_i = sum_j K_ij d(a_j), with the Gaussian kernel K_ij = strength exp(-d_ij^2 / (2 length^2)) of the distance d_ij along the ring. All cells start near the epidermal state and rise together, until small differences break the symmetry at a saddle on the diagonal: the leading cells turn neural and silence their neighbours, so the neural cells end several cell diameters apart.

Parameters:

Name	Type	Description	Default
cells	int	The number of cells on the ring.	40
strength	float	The inhibition between neighbouring cells at zero distance.	1.0
length	float	The range of the inhibition, in cell diameters.	2.5
width	float	The width of the sigmoids.	0.127
noise	float	The largest initial neural state, drawn uniformly from zero.	0.05
seed	int | None	The seed of the initial states.	0

Returns:

Name	Type	Description
system	System	The system, with the neural states cells.a.

Examples:

>>> from skcomplex.systems.biology.systems_biology import bristle_lateral_inhibition
>>> bristle_lateral_inhibition(10).elements[0].prepare().size_
10

Source code in src/skcomplex/systems/biology/systems_biology/_fates.py

@register_system(
    field='biology/systems_biology',
    tags=('lateral inhibition', 'pattern formation', 'bristles', 'notch-delta', 'ode'),
    references=('Corson et al. 2017', 'Rand et al. 2021'),
    defaults={'until': 60.0, 'dt': 0.5},
)
def bristle_lateral_inhibition(
    cells: int = 40,
    *,
    strength: float = 1.0,
    length: float = 2.5,
    width: float = 0.127,
    noise: float = 0.05,
    seed: int | None = 0,
) -> System:
    """Return a ring of cells that inhibit the neural fate of their neighbours, leaving sparse, evenly spaced bristles.

    Following Eq. (2) of Rand et al. (2021), every cell relaxes to `sigma(a_i - h_i)` under the inhibition `h_i = sum_j
    K_ij d(a_j)`, with the Gaussian kernel `K_ij = strength exp(-d_ij^2 / (2 length^2))` of the distance `d_ij` along
    the ring. All cells start near the epidermal state and rise together, until small differences break the symmetry
    at a saddle on the diagonal: the leading cells turn neural and silence their neighbours, so the neural cells end
    several cell diameters apart.

    Args:
        cells:
            The number of cells on the ring.

        strength:
            The inhibition between neighbouring cells at zero distance.

        length:
            The range of the inhibition, in cell diameters.

        width:
            The width of the sigmoids.

        noise:
            The largest initial neural state, drawn uniformly from zero.

        seed:
            The seed of the initial states.

    Returns:
        system:
            The system, with the neural states `cells.a`.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import bristle_lateral_inhibition
        >>> bristle_lateral_inhibition(10).elements[0].prepare().size_
        10
    """
    index = np.arange(cells)
    distance = np.abs(index[:, None] - index[None, :])
    distance = np.minimum(distance, cells - distance)
    kernel = strength * np.exp(-(distance**2) / (2 * length**2))
    np.fill_diagonal(kernel, 0.0)
    states = np.random.default_rng(seed).uniform(0, noise, cells)
    group = GroupElement('cells', element=ProneuralCellElement, size=cells, values={'a': states, 'width': width})
    inhibition = LateralInhibitionInteraction('inhibition', cells='cells', kernel=tuple(map(tuple, kernel.tolist())))
    return System(elements=[group], interactions=[inhibition])

carbon_13_mfa(network=None)

Return a metabolic network that LabellingAnalysis fits to labelling data, the TCA example by default.

Labelled atoms of a tracer spread through the network along its atom transitions, so the measured mass isotopomer distributions of its products carry the fluxes (Raman 2021, Section 9.7). tca_labelling_analysis gives the analysis of the TCA example of Antoniewicz et al. (2007) with data simulated from its known fluxes.

Parameters:

Name	Type	Description	Default
network	MetabolicNetworkElement | None	The metabolic network, or None for the TCA example.	None

Returns:

Name	Type	Description
system	System	The system.

Examples:

>>> from skcomplex.systems.biology.systems_biology import carbon_13_mfa
>>> carbon_13_mfa().elements[0].label
'tca'

Source code in src/skcomplex/systems/biology/systems_biology/_carbon_13_mfa.py

@register_system(
    field='biology/systems_biology',
    tags=('carbon-13 flux analysis', 'isotope labelling', 'elementary metabolite units', 'metabolism', 'static'),
    references=('Raman 2021', 'Antoniewicz et al. 2007', 'Antoniewicz et al. 2006'),
    defaults={'until': 0.0},
)
def carbon_13_mfa(network: MetabolicNetworkElement | None = None) -> System:
    """Return a metabolic network that `LabellingAnalysis` fits to labelling data, the TCA example by default.

    Labelled atoms of a tracer spread through the network along its atom transitions, so the measured mass isotopomer
    distributions of its products carry the fluxes (Raman 2021, Section 9.7). `tca_labelling_analysis` gives the
    analysis of the TCA example of Antoniewicz et al. (2007) with data simulated from its known fluxes.

    Args:
        network:
            The metabolic network, or `None` for the TCA example.

    Returns:
        system:
            The system.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import carbon_13_mfa
        >>> carbon_13_mfa().elements[0].label
        'tca'
    """
    return System(elements=[TcaNetworkElement('tca') if network is None else network], time=StaticTime())

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

Return Goldbeter's system 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 system, with the concentrations clock.M, clock.P0, clock.P1, clock.P2, and clock.PN.

Examples:

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

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

@register_system(
    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 system 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 system, with the concentrations `clock.M`, `clock.P0`, `clock.P1`, `clock.P2`, and `clock.PN`.

    Examples:
        >>> from skcomplex.systems.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)

dynamic_fba(approach='static', *, biomass=0.01, glucose=10.0, vmax=10.0, km=0.015, interval=SYNCHRONISATION)

Return a batch culture of the E. coli core system, which needs the extra metabolism.

The program solves for the fluxes with Michaelis-Menten glucose uptake and then the smallest total flux. By the static optimisation approach it solves in a discrete partition at every synchronisation point, and the biomass and the glucose integrate with the fluxes held (Raman 2021, Eqs. 9.8 and 9.9). By the direct approach it solves at every evaluation of the drift. Without glucose the program is infeasible, and the policy 'zero' stops growth, so biomass grows in one phase and stays.

Parameters:

Name	Type	Description	Default
approach	str	'static' for the static optimisation approach or 'direct' for the direct approach.	'static'
biomass	float	The initial biomass in gDW per L.	0.01
glucose	float	The initial glucose in mM.	10.0
vmax	float	The largest glucose uptake in mmol per gDW per h.	10.0
km	float	The half-saturation concentration of glucose in mM.	0.015
interval	float	The synchronisation interval of the static optimisation approach in h.	SYNCHRONISATION

Returns:

Name	Type	Description
system	System	The system.

Raises:

Type	Description
ValidationError	The approach is unknown.
MissingExtraError	The extra is not installed.

Examples:

>>> from skcomplex.systems.biology.systems_biology import dynamic_fba
>>> [system.label for system in dynamic_fba().systems]
['fba', 'culture']

Source code in src/skcomplex/systems/biology/systems_biology/_dynamic_fba.py

@register_system(
    field='biology/systems_biology',
    tags=('dynamic flux balance', 'linear program', 'metabolism', 'batch culture', 'co-simulation'),
    references=('Raman 2021', 'Mahadevan et al. 2002', 'Varma and Palsson 1994'),
    defaults={'until': 10.0, 'dt': SYNCHRONISATION},
    settings=_settings,
)
def dynamic_fba(
    approach: str = 'static',
    *,
    biomass: float = 0.01,
    glucose: float = 10.0,
    vmax: float = 10.0,
    km: float = 0.015,
    interval: float = SYNCHRONISATION,
) -> System:
    """Return a batch culture of the E. coli core system, which needs the extra `metabolism`.

    The program solves for the fluxes with Michaelis-Menten glucose uptake and then the smallest total flux. By the
    static optimisation approach it solves in a discrete partition at every synchronisation point, and the biomass and
    the glucose integrate with the fluxes held (Raman 2021, Eqs. 9.8 and 9.9). By the direct approach it solves at every
    evaluation of the drift. Without glucose the program is infeasible, and the policy `'zero'` stops growth, so biomass
    grows in one phase and stays.

    Args:
        approach:
            `'static'` for the static optimisation approach or `'direct'` for the direct approach.

        biomass:
            The initial biomass in gDW per L.

        glucose:
            The initial glucose in mM.

        vmax:
            The largest glucose uptake in mmol per gDW per h.

        km:
            The half-saturation concentration of glucose in mM.

        interval:
            The synchronisation interval of the static optimisation approach in h.

    Returns:
        system:
            The system.

    Raises:
        ValidationError:
            The approach is unknown.

        MissingExtraError:
            The extra is not installed.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import dynamic_fba
        >>> [system.label for system in dynamic_fba().systems]
        ['fba', 'culture']
    """
    if approach not in APPROACHES:
        error_msg = f'Unknown approach {approach!r}. Use one of {", ".join(APPROACHES)}.'
        raise ValidationError(error_msg)
    cell = load_network('textbook', 'cell', needed_by='dynamic_fba')
    cell.parsimonious, cell.infeasible = True, 'zero'
    amounts = [BiomassElement('biomass', x=biomass), SubstrateElement('glucose', c=glucose)]
    static = approach == 'static'
    cell_path = 'fba/cell' if static else 'cell'
    biomass_path = 'culture/biomass' if static else 'biomass'
    interactions = [
        ExchangeInteraction(
            'uptake',
            cell=cell_path,
            biomass=biomass_path,
            substrate='culture/glucose' if static else 'glucose',
            reaction=GLUCOSE_EXCHANGE,
            vmax=vmax,
            km=km,
            interval=interval if static else None,
        ),
        GrowthInteraction('growth', cell=cell_path, biomass=biomass_path, reaction=GROWTH),
    ]
    if static:
        return System(
            systems=[System('fba', elements=[cell], time=DiscreteTime()), System('culture', elements=amounts)],
            interactions=interactions,
        )
    return System(elements=[cell, *amounts], interactions=interactions)

early_stage_folding_intermediate(entry=ENTRY, *, cells=CELLS, temperature=1.0)

Return the residues of a protein sampling the cells of an ellipse path of the Ramachandran map by heat bath.

Each residue takes a cell of EllipsePath, the limited conformational space of the early stage after Konieczny et al. (2023 §3.7.1), with the local energy -log p(cell | amino acid) of the probabilities of read_cell_probabilities. The heat bath and the 72 cells of equal t are choices of this system: the 72 of the book counts 5° steps of one angle of the whole map (p. 147). Each unit of time makes as many draws as the chain has residues, each drawing the cell of one random residue from p(cell | amino acid) raised to one over the temperature, renormalised, so at temperature 1 the cells follow the probabilities. The local energies enter through a row score, their sum over the chain, since the change of the score with the cell of one residue is the change of its local energy.

Parameters:

Name	Type	Description	Default
entry	str	The bundled entry whose sequence the system samples.	ENTRY
cells	int	The number of cells of the path.	CELLS
temperature	float	The temperature, in the units of the local energy.	1.0

Returns:

Name	Type	Description
system	System	The system, with the lattice residues of the state cell.

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.
ValidationError	The temperature is not positive.

Examples:

>>> from skcomplex.systems.biology.systems_biology import early_stage_folding_intermediate
>>> early_stage_folding_intermediate().elements[0].space.shape
(1, 58)

Source code in src/skcomplex/systems/biology/systems_biology/_early_stage.py

@register_system(
    field='biology/systems_biology',
    tags=('protein', 'ramachandran', 'early stage', 'heat bath', 'monte carlo'),
    references=('Konieczny et al. 2023', 'Brylinski et al. 2004', 'Brylinski et al. 2005'),
    defaults={'until': 300},
)
def early_stage_folding_intermediate(
    entry: str = ENTRY,
    *,
    cells: int = CELLS,
    temperature: float = 1.0,
) -> System:
    """Return the residues of a protein sampling the cells of an ellipse path of the Ramachandran map by heat bath.

    Each residue takes a cell of `EllipsePath`, the limited conformational space of the early stage after Konieczny et
    al. (2023 §3.7.1), with the local energy `-log p(cell | amino acid)` of the probabilities of
    `read_cell_probabilities`. The heat bath and the 72 cells of equal `t` are choices of this system: the 72 of the
    book counts 5° steps of one angle of the whole map (p. 147).
    Each unit of time makes as many draws as the chain has residues, each drawing the cell of one random residue from
    `p(cell | amino acid)` raised to one over the temperature, renormalised, so at temperature 1 the cells follow the
    probabilities. The local energies enter through a row score, their sum over the chain, since the change of the score
    with the cell of one residue is the change of its local energy.

    Args:
        entry:
            The bundled entry whose sequence the system samples.

        cells:
            The number of cells of the path.

        temperature:
            The temperature, in the units of the local energy.

    Returns:
        system:
            The system, with the lattice `residues` of the state `cell`.

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

        ValidationError:
            The temperature is not positive.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import early_stage_folding_intermediate
        >>> early_stage_folding_intermediate().elements[0].space.shape
        (1, 58)
    """
    import_extra('Bio.PDB', 'early_stage_folding_intermediate', 'structures')
    if temperature <= 0:
        error_msg = f'The temperature of the system must be positive, but it is {temperature}.'
        raise ValidationError(error_msg)
    names = [residue for _, _, residue in read_structure(entry).residues]
    amino = np.array([AMINO_ACIDS.index(name) for name in names], dtype=float)
    energies = -np.log(read_cell_probabilities(cells))
    return System(
        elements=[
            LatticeElement(
                'residues',
                element=PathCellElement,
                space=LatticeSpace((1, len(names)), boundary='fixed'),
                values={'amino': amino[None]},
                visual=VisualSpec('image', {'states': ('cell',), 'low': 0.0, 'high': float(cells - 1)}),
            ),
        ],
        interactions=[
            SamplingInteraction(
                'sampling',
                target='residues',
                state='cell',
                objective=RowScoreSpec(partial(_score_cells, energies), proposals=len(names)),
                move='random-value',
                intensity=1 / temperature,
                rule='heat-bath',
                values=tuple(float(cell) for cell in range(cells)),
            ),
        ],
        time=DiscreteTime(),
    )

elementary_flux_modes(knockouts=())

Return the toy network of Raman whose steady-state cone ElementaryModeAnalysis enumerates.

Every steady-state flux of the network is a non-negative combination of the elementary modes of its cone (Raman 2021, Section 9.8). The toy network has two modes, through v1 and through v2 and v3, both with growth yield 1/12 on A, so the bound 1000 on v6 limits the growth to 50 whichever route carries the flux.

Parameters:

Name	Type	Description	Default
knockouts	Sequence[str]	The reactions to fix at zero.	()

Returns:

Name	Type	Description
system	System	The system.

Examples:

>>> from skcomplex.core.analysis import ElementaryModeAnalysis
>>> from skcomplex.systems.biology.systems_biology import elementary_flux_modes
>>> analysis = ElementaryModeAnalysis(normalize='cell.v_b1').prepare(elementary_flux_modes()).run()
>>> analysis.modes_.shape[0]
2

Source code in src/skcomplex/systems/biology/systems_biology/_elementary_flux_modes.py

@register_system(
    field='biology/systems_biology',
    tags=('elementary flux modes', 'extreme pathways', 'linear program', 'metabolism', 'static'),
    references=('Raman 2021', 'Schuster and Hilgetag 1994', 'Schilling et al. 2000', 'Gagneur and Klamt 2004'),
    defaults={'until': 0.0},
)
def elementary_flux_modes(knockouts: Sequence[str] = ()) -> System:
    """Return the toy network of Raman whose steady-state cone `ElementaryModeAnalysis` enumerates.

    Every steady-state flux of the network is a non-negative combination of the elementary modes of its cone (Raman
    2021, Section 9.8). The toy network has two modes, through `v1` and through `v2` and `v3`, both with growth yield
    `1/12` on `A`, so the bound 1000 on `v6` limits the growth to 50 whichever route carries the flux.

    Args:
        knockouts:
            The reactions to fix at zero.

    Returns:
        system:
            The system.

    Examples:
        >>> from skcomplex.core.analysis import ElementaryModeAnalysis
        >>> from skcomplex.systems.biology.systems_biology import elementary_flux_modes
        >>> analysis = ElementaryModeAnalysis(normalize='cell.v_b1').prepare(elementary_flux_modes()).run()
        >>> analysis.modes_.shape[0]
        2
    """
    closed = {f'{bound}_{name}': 0.0 for name in knockouts for bound in ('lower', 'upper')}
    return System(elements=[ToyNetworkElement('cell').set_params(**closed)], time=StaticTime())

enzyme_adaptation_circuit(links=NEGATIVE_FEEDBACK_BUFFER, *, k=RATES, K=CONSTANTS, k_basal=BASAL_RATES, K_basal=BASAL_CONSTANTS, level=0.5, change=0.1, at=300.0)

Return a circuit of three enzymes A, B, and C of Ma et al., whose output C may adapt to a step of the input.

The input I activates A, and links[i][j], 1, -1, or 0, says whether enzyme j activates, deactivates, or leaves enzyme i, with Michaelis-Menten kinetics of the rate k[j][i] and the constant K[j][i], the last row of k and K for the input. A node with no activator gets a basal activating enzyme and one with no deactivator a basal deactivating enzyme, both of concentration 0.5 with the constants k_basal and K_basal. The input steps from level by change at the time at. Ma et al. enumerated the 16,038 topologies of three nodes whose input reaches C and found that adaptation, a sensitive transient followed by a precise return, needs a negative feedback loop with a buffer node or an incoherent feedforward loop with a proportioner node. The default is the negative feedback loop A to C to B, B deactivating A, with a parameter set that adapts: the buffer B works in saturation, so its steady state fixes the output.

Parameters:

Name	Type	Description	Default
links	Sequence[Sequence[int]]	The signs of the regulations, one row per regulated node and one column per regulating node.	NEGATIVE_FEEDBACK_BUFFER
k	Sequence[Sequence[float]]	The catalytic rates, one row per regulating node and the input last, one column per regulated node.	RATES
K	Sequence[Sequence[float]]	The Michaelis constants, in the layout of k.	CONSTANTS
k_basal	Sequence[float]	The catalytic rates of the basal enzymes of every node.	BASAL_RATES
K_basal	Sequence[float]	The Michaelis constants of the basal enzymes of every node.	BASAL_CONSTANTS
level	float	The input before the step.	0.5
change	float	The step of the input.	0.1
at	float	The time of the step.	300.0

Returns:

Name	Type	Description
system	System	The system, with the active fractions A.X, B.X, and C.X.

Raises:

Type	Description
ValidationError	The links are not a three by three matrix of 1, -1, and 0.

Examples:

>>> from skcomplex.systems.biology.systems_biology import enzyme_adaptation_circuit
>>> [interaction.label for interaction in enzyme_adaptation_circuit().interactions][:3]
['input', 'B_A', 'C_B']

Source code in src/skcomplex/systems/biology/systems_biology/_adaptation.py

@register_system(
    field='biology/systems_biology',
    tags=('adaptation', 'network motifs', 'negative feedback', 'incoherent feedforward', 'ode'),
    references=('Ma et al. 2009', 'Raman 2021'),
    defaults={'until': 800.0, 'dt': 1.0},
)
def enzyme_adaptation_circuit(
    links: Sequence[Sequence[int]] = NEGATIVE_FEEDBACK_BUFFER,
    *,
    k: Sequence[Sequence[float]] = RATES,
    K: Sequence[Sequence[float]] = CONSTANTS,
    k_basal: Sequence[float] = BASAL_RATES,
    K_basal: Sequence[float] = BASAL_CONSTANTS,
    level: float = 0.5,
    change: float = 0.1,
    at: float = 300.0,
) -> System:
    """Return a circuit of three enzymes A, B, and C of Ma et al., whose output C may adapt to a step of the input.

    The input `I` activates A, and `links[i][j]`, 1, -1, or 0, says whether enzyme `j` activates, deactivates, or
    leaves enzyme `i`, with Michaelis-Menten kinetics of the rate `k[j][i]` and the constant `K[j][i]`, the last row of
    `k` and `K` for the input. A node with no activator gets a basal activating enzyme and one with no deactivator a
    basal deactivating enzyme, both of concentration 0.5 with the constants `k_basal` and `K_basal`. The input steps
    from `level` by `change` at the time `at`. Ma et al. enumerated the 16,038 topologies of three nodes whose input
    reaches C and found that adaptation, a sensitive transient followed by a precise return, needs a negative feedback
    loop with a buffer node or an incoherent feedforward loop with a proportioner node. The default is the negative
    feedback loop A to C to B, B deactivating A, with a parameter set that adapts: the buffer B works in saturation, so
    its steady state fixes the output.

    Args:
        links:
            The signs of the regulations, one row per regulated node and one column per regulating node.

        k:
            The catalytic rates, one row per regulating node and the input last, one column per regulated node.

        K:
            The Michaelis constants, in the layout of `k`.

        k_basal:
            The catalytic rates of the basal enzymes of every node.

        K_basal:
            The Michaelis constants of the basal enzymes of every node.

        level:
            The input before the step.

        change:
            The step of the input.

        at:
            The time of the step.

    Returns:
        system:
            The system, with the active fractions `A.X`, `B.X`, and `C.X`.

    Raises:
        ValidationError:
            The links are not a three by three matrix of 1, -1, and 0.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import enzyme_adaptation_circuit
        >>> [interaction.label for interaction in enzyme_adaptation_circuit().interactions][:3]
        ['input', 'B_A', 'C_B']
    """
    signs = np.asarray(links)
    if signs.shape != (3, 3) or not np.isin(signs, (-1, 0, 1)).all():
        error_msg = f'The links must be a 3 by 3 matrix of 1, -1, and 0, not {links!r}.'
        raise ValidationError(error_msg)
    rates, constants = np.asarray(k, dtype=float), np.asarray(K, dtype=float)
    interactions: list[BaseInteraction] = [
        EnzymeStimulusInteraction(
            'input',
            target='A',
            k=float(rates[3, 0]),
            K=float(constants[3, 0]),
            level=level,
            change=change,
            at=at,
        ),
    ]
    for target, source in zip(*np.nonzero(signs), strict=True):
        interactions.append(
            EnzymeRegulationInteraction(
                f'{NODES[source]}_{NODES[target]}',
                target=NODES[target],
                source=NODES[source],
                sign=int(signs[target, source]),
                k=float(rates[source, target]),
                K=float(constants[source, target]),
            ),
        )
    for index, node in enumerate(NODES):
        activated = index == 0 or (signs[index] > 0).any()
        for sign, present in ((1, activated), (-1, (signs[index] < 0).any())):
            if not present:
                interactions.append(
                    EnzymeStimulusInteraction(
                        f'basal_{"activation" if sign > 0 else "deactivation"}_{node}',
                        target=node,
                        sign=sign,
                        k=float(k_basal[index]),
                        K=float(K_basal[index]),
                    ),
                )
    return System(elements=[EnzymeNodeElement(node) for node in NODES], interactions=interactions)

flux_balance(network='toy', *, oxygen=1000.0, parsimonious=False)

Return flux balance analysis of a metabolic network, which sets its fluxes to maximise growth.

The steady-state balances of the metabolites and the bounds of the reactions confine the fluxes, and the optimum of the growth flux sets them (Raman 2021, Section 8.6). The toy network of Table 8.2 reaches growth 50, set by the bound 1000 on v6. The E. coli core system needs the extra metabolism, and grows at 0.8739 h^-1 on glucose with oxygen and at 0.2117 h^-1 without.

Parameters:

Name	Type	Description	Default
network	str	'toy' for Raman's toy network or 'ecoli' for the E. coli core system.	'toy'
oxygen	float	The largest oxygen uptake of the E. coli core system.	1000.0
parsimonious	bool	Whether to minimise the total absolute flux after growth.	False

Returns:

Name	Type	Description
system	System	The system.

Raises:

Type	Description
ValidationError	The network is unknown.
MissingExtraError	The E. coli core system needs the extra, which is not installed.

Examples:

>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.biology.systems_biology import flux_balance
>>> trajectory = Simulation(flux_balance()).prepare().run(0.0).trajectory_
>>> round(float(trajectory.values('cell.v_vbio')[0]), 6)
50.0

Source code in src/skcomplex/systems/biology/systems_biology/_flux_balance.py

@register_system(
    field='biology/systems_biology',
    tags=('flux balance', 'linear program', 'metabolism', 'static'),
    references=('Raman 2021', 'Orth et al. 2010'),
    defaults={'until': 0.0},
)
def flux_balance(
    network: str = 'toy',
    *,
    oxygen: float = 1000.0,
    parsimonious: bool = False,
) -> System:
    """Return flux balance analysis of a metabolic network, which sets its fluxes to maximise growth.

    The steady-state balances of the metabolites and the bounds of the reactions confine the fluxes, and the optimum of
    the growth flux sets them (Raman 2021, Section 8.6). The toy network of Table 8.2 reaches growth 50, set by the
    bound 1000 on `v6`. The E. coli core system needs the extra `metabolism`, and grows at 0.8739 h^-1 on glucose with
    oxygen and at 0.2117 h^-1 without.

    Args:
        network:
            `'toy'` for Raman's toy network or `'ecoli'` for the E. coli core system.

        oxygen:
            The largest oxygen uptake of the E. coli core system.

        parsimonious:
            Whether to minimise the total absolute flux after growth.

    Returns:
        system:
            The system.

    Raises:
        ValidationError:
            The network is unknown.

        MissingExtraError:
            The E. coli core system needs the extra, which is not installed.

    Examples:
        >>> from skcomplex.core.engine.runtime import Simulation
        >>> from skcomplex.systems.biology.systems_biology import flux_balance
        >>> trajectory = Simulation(flux_balance()).prepare().run(0.0).trajectory_
        >>> round(float(trajectory.values('cell.v_vbio')[0]), 6)
        50.0
    """
    if network not in NETWORKS:
        error_msg = f'Unknown network {network!r}. Use one of {", ".join(NETWORKS)}.'
        raise ValidationError(error_msg)
    if network == 'toy':
        return System(elements=[ToyNetworkElement('cell', parsimonious=parsimonious)], time=StaticTime())
    element = load_network('textbook', 'cell', needed_by='flux_balance')
    element.parsimonious = parsimonious
    identifier = f'lower_{element._identifier(OXYGEN_EXCHANGE)}'
    element.set_params(**{identifier: -float(oxygen)})
    return System(elements=[element], time=StaticTime())

fuzzy_oil_drop(entry=ENTRY, *, dynamics='replay', proposals=1, width=60.0, temperature=0.0)

Return the fuzzy oil drop of a protein, an ensemble replayed one member per frame or a chain folded by moves.

The default replays the systems of an entry, which sets the effective atoms of the chain, the mean atom of each residue, to one system per unit of time. With dynamics='folding', the chain starts extended, with all phi = -120° and psi = 130°, and each unit of time proposes moves of its dihedral angles, accepted by the row score RD plus one for every pair of effective atoms closer than 3 Å and at least two residues apart. The observables RD and K follow Konieczny et al. (2023 §1.3.1) through FuzzyOilDropAnalysis. The book gives no folding dynamics, so the extended start, the clash rule, the rotation width of 60°, and the temperature are choices of this system.

Parameters:

Name	Type	Description	Default
entry	str	The bundled entry, '1L2Y', '1UBQ', or '5PTI'.	ENTRY
dynamics	str	'replay' or 'folding'.	'replay'
proposals	int	The number of dihedral moves per unit of time of the folding.	1
width	float	The largest rotation of a dihedral move, in degrees.	60.0
temperature	float	The temperature of the folding, one over the intensity of its sampling, where 0 accepts only moves that do not raise the score.	0.0

Returns:

Name	Type	Description
system	System	The system, with the lattice chain of the state x or angle.

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.
ValidationError	The dynamics are unknown.

Examples:

>>> from skcomplex.systems.biology.systems_biology import fuzzy_oil_drop
>>> fuzzy_oil_drop().elements[0].space.shape
(20, 3)

Source code in src/skcomplex/systems/biology/systems_biology/_fuzzy_oil_drop.py

@register_system(
    field='biology/systems_biology',
    tags=('protein', 'hydrophobicity', 'divergence entropy', 'structures', 'monte carlo'),
    references=('Konieczny et al. 2023', 'Levitt 1976'),
    defaults={'until': 37},
    observables=_read_observables,
)
def fuzzy_oil_drop(
    entry: str = ENTRY,
    *,
    dynamics: str = 'replay',
    proposals: int = 1,
    width: float = 60.0,
    temperature: float = 0.0,
) -> System:
    """Return the fuzzy oil drop of a protein, an ensemble replayed one member per frame or a chain folded by moves.

    The default replays the systems of an entry, which sets the effective atoms of the chain, the mean atom of each
    residue, to one system per unit of time. With `dynamics='folding'`, the chain starts extended, with all `phi =
    -120°` and `psi = 130°`, and each unit of time proposes moves of its dihedral angles, accepted by the row score `RD`
    plus one for every pair of effective atoms closer than 3 Å and at least two residues apart. The observables `RD` and
    `K` follow Konieczny et al. (2023 §1.3.1) through `FuzzyOilDropAnalysis`. The book gives no folding dynamics, so the
    extended start, the clash rule, the rotation width of 60°, and the temperature are choices of this system.

    Args:
        entry:
            The bundled entry, `'1L2Y'`, `'1UBQ'`, or `'5PTI'`.

        dynamics:
            `'replay'` or `'folding'`.

        proposals:
            The number of dihedral moves per unit of time of the folding.

        width:
            The largest rotation of a dihedral move, in degrees.

        temperature:
            The temperature of the folding, one over the intensity of its sampling, where 0 accepts only moves that do
            not raise the score.

    Returns:
        system:
            The system, with the lattice `chain` of the state `x` or `angle`.

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

        ValidationError:
            The dynamics are unknown.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import fuzzy_oil_drop
        >>> fuzzy_oil_drop().elements[0].space.shape
        (20, 3)
    """
    import_extra('Bio.PDB', 'fuzzy_oil_drop', 'structures')
    if dynamics not in DYNAMICS:
        error_msg = f'Unknown dynamics {dynamics!r}. Use one of {", ".join(DYNAMICS)}.'
        raise ValidationError(error_msg)
    structure = read_structure(entry)
    if dynamics == 'replay':
        ensemble = structure.effective_atoms
        return System(
            label=entry,
            elements=[
                LatticeElement(
                    'chain',
                    element=ChainCoordinateElement,
                    space=LatticeSpace((ensemble.shape[1], COORDINATES), boundary='fixed'),
                    values={'x': ensemble[0]},
                    visual=VisualSpec(BackboneVisual, {'position': 'x'}),
                ),
            ],
            interactions=[EnsembleReplayInteraction('replay', target='chain', state='x', ensemble=ensemble)],
            time=DiscreteTime(),
        )
    angles = np.full(BackboneGeometry(structure).prepare().angles_.size, EXTENDED_PHI)
    angles[angles.size // 2 :] = EXTENDED_PSI
    return System(
        label=entry,
        elements=[
            LatticeElement(
                'chain',
                element=DihedralAngleElement,
                space=LatticeSpace((1, angles.size), boundary='fixed'),
                values={'angle': angles[None]},
                visual=VisualSpec(BackboneVisual, {'position': 'angle', 'structure': structure}),
            ),
        ],
        interactions=[
            SamplingInteraction(
                'folding',
                target='chain',
                state='angle',
                objective=RowScoreSpec(partial(_score_angles, entry), proposals=proposals),
                move=partial(_rotate_dihedral, width),
                intensity=float('inf') if temperature == 0 else 1 / temperature,
            ),
        ],
        time=DiscreteTime(),
    )

generalized_logical_network(regulations=None, levels=None)

Return a generalised logical network of Thomas, whose variables step asynchronously towards focal levels.

Every variable takes discrete levels, and each of its regulators acts once its level reaches a threshold. The logical parameters give the focal level of the variable for every combination of acting regulators, and the variable rises or falls one level at a time towards it, each step a random event, so the dynamics explores the asynchronous state transition graph. The default is the negative circuit of a three-level x that activates y at its second threshold and a Boolean y that represses x, whose state graph is a cycle, a sustained oscillation. A positive circuit, two mutual repressors, has two stable states instead.

Parameters:

Name	Type	Description	Default
regulations	Regulations | None	For every variable, its regulators, their thresholds, and the logical parameters of its focal level, or None for the negative circuit.	None
levels	Mapping[str, int] | None	The initial level of every variable, or None for 0.	None

Returns:

Name	Type	Description
system	System	The system, with the levels <variable>.level.

Examples:

>>> from skcomplex.systems.biology.systems_biology import generalized_logical_network
>>> [element.label for element in generalized_logical_network().elements]
['x', 'y']

Source code in src/skcomplex/systems/biology/systems_biology/_logical_networks.py

@register_system(
    field='biology/systems_biology',
    tags=('logical network', 'multi-valued', 'asynchronous', 'thomas', 'jumps'),
    references=('Thomas 1991', 'Raman 2021'),
    defaults={'until': 50.0, 'dt': 0.25},
)
def generalized_logical_network(
    regulations: Regulations | None = None,
    levels: Mapping[str, int] | None = None,
) -> System:
    """Return a generalised logical network of Thomas, whose variables step asynchronously towards focal levels.

    Every variable takes discrete levels, and each of its regulators acts once its level reaches a threshold. The
    logical parameters give the focal level of the variable for every combination of acting regulators, and the
    variable rises or falls one level at a time towards it, each step a random event, so the dynamics explores the
    asynchronous state transition graph. The default is the negative circuit of a three-level `x` that activates `y` at
    its second threshold and a Boolean `y` that represses `x`, whose state graph is a cycle, a sustained oscillation.
    A positive circuit, two mutual repressors, has two stable states instead.

    Args:
        regulations:
            For every variable, its regulators, their thresholds, and the logical parameters of its focal level, or
            `None` for the negative circuit.

        levels:
            The initial level of every variable, or `None` for 0.

    Returns:
        system:
            The system, with the levels `<variable>.level`.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import generalized_logical_network
        >>> [element.label for element in generalized_logical_network().elements]
        ['x', 'y']
    """
    network = NEGATIVE_CIRCUIT if regulations is None else regulations
    start = levels or {}
    return System(
        elements=[LogicalVariableElement(name, level=float(start.get(name, 0))) for name in network],
        interactions=[
            LogicalRegulationInteraction(
                f'{name}_regulation',
                target=name,
                regulators=list(regulators),
                thresholds=tuple(thresholds),
                parameters=dict(parameters),
            )
            for name, (regulators, thresholds, parameters) in network.items()
        ],
    )

glycolysis_oscillations(*, flow=0.011, V4=10.0)

Return the glycolysis of a yeast extract in a flow reactor, whose NADH oscillates with a period of minutes.

The model of Nielsen et al. (1998), BioModels BIOMD0000000042, has fifteen metabolites, the inhibition of phosphofructokinase by ATP and its activation by AMP, and the exchange with the feed of the reactor at the rate flow. With the published parameters the reactor settles on sustained oscillations of NADH and of all the other metabolites. Raman (2021, exercise 5.7) varies V4, K4GAP, K4NAD, k8f, and k8b by half to study their effect on NAD and NADH.

Parameters:

Name	Type	Description	Default
flow	float	The flow rate of the reactor, per minute.	0.011
V4	float	The largest rate of glyceraldehyde phosphate dehydrogenase, in mM per minute.	10.0

Returns:

Name	Type	Description
system	System	The system, with the concentrations extract.<metabolite>, such as extract.NADH, in mM.

Examples:

>>> from skcomplex.systems.biology.systems_biology import glycolysis_oscillations
>>> glycolysis_oscillations(flow=0.02).elements[0].get_value('flow')
0.02

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

@register_system(
    field='biology/systems_biology',
    tags=('glycolysis', 'metabolic oscillations', 'nadh', 'biomodels', 'ode'),
    references=('Nielsen et al. 1998', 'Raman 2021'),
    defaults={'until': 300.0, 'dt': 0.5},
)
def glycolysis_oscillations(*, flow: float = 0.011, V4: float = 10.0) -> System:
    """Return the glycolysis of a yeast extract in a flow reactor, whose NADH oscillates with a period of minutes.

    The model of Nielsen et al. (1998), BioModels BIOMD0000000042, has fifteen metabolites, the inhibition of
    phosphofructokinase by `ATP` and its activation by `AMP`, and the exchange with the feed of the reactor at the rate
    `flow`. With the published parameters the reactor settles on sustained oscillations of `NADH` and of all the
    other metabolites. Raman (2021, exercise 5.7) varies `V4`, `K4GAP`, `K4NAD`, `k8f`, and `k8b` by half to study
    their effect on `NAD` and `NADH`.

    Args:
        flow:
            The flow rate of the reactor, per minute.

        V4:
            The largest rate of glyceraldehyde phosphate dehydrogenase, in mM per minute.

    Returns:
        system:
            The system, with the concentrations `extract.<metabolite>`, such as `extract.NADH`, in mM.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import glycolysis_oscillations
        >>> glycolysis_oscillations(flow=0.02).elements[0].get_value('flow')
        0.02
    """
    variables = {
        'flow': DeclaredVariableSpec(flow, ParameterSpec(domain=RealSet.closed_open(0, oo))),
        'V4': DeclaredVariableSpec(V4, ParameterSpec(domain=RealSet.closed_open(0, oo))),
    }
    return System(elements=[GlycolysisElement('extract', variables=variables)])

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 system 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 system, 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.systems.biology.systems_biology import goodwin_oscillator
>>> [stage.label for stage in goodwin_oscillator(4).elements]
['x1', 'x2', 'x3', 'x4']

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

@register_system(
    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 system 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 system, 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.systems.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 system, 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.systems.biology.systems_biology import hormone_cascade_amplifier
>>> [element.label for element in hormone_cascade_amplifier().elements][-1]
'phosphorylase'

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

@register_system(
    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 system, 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.systems.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 system 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 system, with the values <protein>.x.

Raises:

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

Examples:

>>> from skcomplex.systems.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/systems/biology/systems_biology/_boolean_systems.py

@register_system(
    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 system 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 system, with the values `<protein>.x`.

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

    Examples:
        >>> from skcomplex.systems.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)

metabolic_perturbations(genes=())

Return the E. coli core system with some genes knocked out, which needs the extra metabolism.

A gene that is off closes the reactions whose rules need it (Raman 2021, Section 9.1). CobraAnalysis scans the deletions and adjusts the fluxes by MoMA and ROOM, StrainDesignAnalysis finds OptKnock designs, and ProductionScanAnalysis scans fluxes with an enforced product.

Parameters:

Name	Type	Description	Default
genes	Sequence[str]	The names of the genes to knock out.	()

Returns:

Name	Type	Description
system	System	The system.

Raises:

Type	Description
MissingExtraError	The extra is not installed.

Examples:

>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.biology.systems_biology import metabolic_perturbations
>>> system = metabolic_perturbations(genes=('b4025',))
>>> trajectory = Simulation(system).prepare().run(0.0).trajectory_
>>> round(float(trajectory.values('cell.objective_0')[0]), 3)
0.863

Source code in src/skcomplex/systems/biology/systems_biology/_metabolic_perturbations.py

@register_system(
    field='biology/systems_biology',
    tags=('flux balance', 'gene knock-out', 'linear program', 'metabolism', 'static'),
    references=('Raman 2021', 'Segre et al. 2002', 'Shlomi et al. 2005', 'Burgard et al. 2003'),
    defaults={'until': 0.0},
)
def metabolic_perturbations(genes: Sequence[str] = ()) -> System:
    """Return the E. coli core system with some genes knocked out, which needs the extra `metabolism`.

    A gene that is off closes the reactions whose rules need it (Raman 2021, Section 9.1). `CobraAnalysis` scans the
    deletions and adjusts the fluxes by MoMA and ROOM, `StrainDesignAnalysis` finds OptKnock designs, and
    `ProductionScanAnalysis` scans fluxes with an enforced product.

    Args:
        genes:
            The names of the genes to knock out.

    Returns:
        system:
            The system.

    Raises:
        MissingExtraError:
            The extra is not installed.

    Examples:
        >>> from skcomplex.core.engine.runtime import Simulation
        >>> from skcomplex.systems.biology.systems_biology import metabolic_perturbations
        >>> system = metabolic_perturbations(genes=('b4025',))
        >>> trajectory = Simulation(system).prepare().run(0.0).trajectory_
        >>> round(float(trajectory.values('cell.objective_0')[0]), 3)
        0.863
    """
    element = load_network('textbook', 'cell', needed_by='metabolic_perturbations')
    element.set_params(**{element.gene_name(gene): 0.0 for gene in genes})
    return System(elements=[element], time=StaticTime())

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

Return the skeleton system 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 system, with the membrane fields cell.md and cell.mde and the cytosolic fields cell.cdt, cell.cdd, and cell.ce.

Examples:

>>> from skcomplex.systems.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/systems/biology/systems_biology/_min_proteins.py

@register_system(
    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 system 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 system, with the membrane fields `cell.md` and `cell.mde` and the cytosolic fields `cell.cdt`,
            `cell.cdd`, and `cell.ce`.

    Examples:
        >>> from skcomplex.systems.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 system is a jump process.	True

Returns:

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

Examples:

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

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

@register_system(
    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 system is a jump process.

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

    Examples:
        >>> from skcomplex.systems.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])

noisy_boolean_network(n=1000, *, k=3, eta=0.1, seed=0)

Return the Boolean network of Huepe and Aldana-González, nodes that follow the noisy majority of k inputs.

Every node takes the sign of the sum of its k random inputs, flipped with the probability eta. The order parameter Psi, the time average of the absolute mean state, undergoes a continuous phase transition at the noise eta_c = (1 - 1 / Phi'(0)) / 2 of mean-field theory, with Phi'(0) = k C(k - 1, (k - 1) / 2) / 2^(k - 1), which is 1/6 for k = 3. Below it the network orders, above it the mean state falls to fluctuations of order 1 / sqrt(n).

Parameters:

Name	Type	Description	Default
n	int	The number of nodes.	1000
k	int	The number of inputs of every node, odd.	3
eta	float	The noise, the probability of a flip.	0.1
seed	int | None	The seed of the inputs.	0

Returns:

Name	Type	Description
system	System	The system, with the states nodes.sigma.

Raises:

Type	Description
ValidationError	The number of inputs is even.

Examples:

>>> from skcomplex.systems.biology.systems_biology import noisy_boolean_network
>>> noisy_boolean_network(10).elements[0].prepare().size_
10

Source code in src/skcomplex/systems/biology/systems_biology/_boolean_systems.py

@register_system(
    field='biology/systems_biology',
    tags=('boolean network', 'noise', 'phase transition', 'mean field', 'stochastic map'),
    references=('Huepe and Aldana-González 2002', 'Gros 2015'),
    defaults={'until': 200, 'dt': 1},
)
def noisy_boolean_network(n: int = 1000, *, k: int = 3, eta: float = 0.1, seed: int | None = 0) -> System:
    """Return the Boolean network of Huepe and Aldana-González, nodes that follow the noisy majority of `k` inputs.

    Every node takes the sign of the sum of its `k` random inputs, flipped with the probability `eta`. The order
    parameter `Psi`, the time average of the absolute mean state, undergoes a continuous phase transition at the noise
    `eta_c = (1 - 1 / Phi'(0)) / 2` of mean-field theory, with `Phi'(0) = k C(k - 1, (k - 1) / 2) / 2^(k - 1)`, which is
    1/6 for `k = 3`. Below it the network orders, above it the mean state falls to fluctuations of order `1 /
    sqrt(n)`.

    Args:
        n:
            The number of nodes.

        k:
            The number of inputs of every node, odd.

        eta:
            The noise, the probability of a flip.

        seed:
            The seed of the inputs.

    Returns:
        system:
            The system, with the states `nodes.sigma`.

    Raises:
        ValidationError:
            The number of inputs is even.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import noisy_boolean_network
        >>> noisy_boolean_network(10).elements[0].prepare().size_
        10
    """
    if k % 2 == 0:
        error_msg = f'The number of inputs must be odd so that the majority is defined, not {k}.'
        raise ValidationError(error_msg)
    rng = np.random.default_rng(seed)
    inputs = np.zeros((n, n))
    for node in range(n):
        inputs[node, rng.choice(np.delete(np.arange(n), node), size=k, replace=False)] = 1.0
    nodes = GroupElement(
        'nodes',
        element=MajorityNodeElement,
        size=n,
        values={'sigma': np.ones(n), 'eta': eta},
    )
    majority = NoisyMajorityInteraction('majority', nodes='nodes', inputs=scipy.sparse.csr_array(inputs))
    return System(elements=[nodes], interactions=[majority], time=DiscreteTime())

petri_net(places=None, transitions=None, *, stochastic=True)

Return a Petri net, places that hold tokens and transitions that consume and produce them when they fire.

Every transition is a pair of input and output arcs with weights and a rate. In a stochastic Petri net it fires at random times with the mass-action propensity of its input places, consuming the tokens of its input arcs and producing those of its output arcs, and otherwise it is a continuous Petri net of rate equations. Every place invariant, a weighted sum of tokens that no transition changes, is conserved. The default net is the expression of a gene that switches on, is transcribed and translated, and is switched off by its own protein, with the invariant gene_on + gene_off = 1.

Parameters:

Name	Type	Description	Default
places	Mapping[str, int] | None	The initial tokens of every place, or None for the gene expression net.	None
transitions	Transitions | None	For every transition, its label, input arcs, output arcs, and rate, or None for the gene expression net.	None
stochastic	bool	Whether the transitions fire at random times.	True

Returns:

Name	Type	Description
system	System	The system, with the tokens <place>.n of every place.

Examples:

>>> from skcomplex.systems.biology.systems_biology import petri_net
>>> [element.label for element in petri_net().elements]
['gene_off', 'gene_on', 'mRNA', 'protein']

Source code in src/skcomplex/systems/biology/systems_biology/_logical_networks.py

@register_system(
    field='biology/systems_biology',
    tags=('petri net', 'stochastic simulation', 'gene expression', 'invariants', 'jumps'),
    references=('Petri 1962', 'Reddy et al. 1993', 'Raman 2021'),
    defaults={'until': 100.0, 'dt': 0.5},
)
def petri_net(
    places: Mapping[str, int] | None = None,
    transitions: Transitions | None = None,
    *,
    stochastic: bool = True,
) -> System:
    """Return a Petri net, places that hold tokens and transitions that consume and produce them when they fire.

    Every transition is a pair of input and output arcs with weights and a rate. In a stochastic Petri net it fires at
    random times with the mass-action propensity of its input places, consuming the tokens of its input arcs and
    producing those of its output arcs, and otherwise it is a continuous Petri net of rate equations. Every place
    invariant, a weighted sum of tokens that no transition changes, is conserved. The default net is the expression
    of a gene that switches on, is transcribed and translated, and is switched off by its own protein, with the
    invariant `gene_on + gene_off = 1`.

    Args:
        places:
            The initial tokens of every place, or `None` for the gene expression net.

        transitions:
            For every transition, its label, input arcs, output arcs, and rate, or `None` for the gene expression net.

        stochastic:
            Whether the transitions fire at random times.

    Returns:
        system:
            The system, with the tokens `<place>.n` of every place.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import petri_net
        >>> [element.label for element in petri_net().elements]
        ['gene_off', 'gene_on', 'mRNA', 'protein']
    """
    marking = GENE_EXPRESSION_PLACES if places is None else places
    arcs = GENE_EXPRESSION_TRANSITIONS if transitions is None else transitions
    mode = 'jump' if stochastic else 'rate'
    return System(
        elements=[CompartmentElement(name, n=float(tokens)) for name, tokens in marking.items()],
        interactions=[
            ReactionInteraction(
                label,
                reactants={(place, 'n'): weight for place, weight in inputs.items()} or None,
                products={(place, 'n'): weight for place, weight in outputs.items()} or None,
                k=rate,
                mode=mode,
            )
            for label, inputs, outputs, rate in arcs
        ],
    )

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 system, with the amounts loop.A, loop.B, and loop.C.

Examples:

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

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

@register_system(
    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 system, with the amounts `loop.A`, `loop.B`, and `loop.C`.

    Examples:
        >>> from skcomplex.systems.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 system, with the values <node>.x and the indices of the drawn functions <node>.function.

Examples:

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

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

@register_system(
    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 system, with the values `<node>.x` and the indices of the drawn functions `<node>.function`.

    Examples:
        >>> from skcomplex.systems.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 system, with the values node_<i>.x.

Raises:

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

Examples:

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

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

@register_system(
    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 system, with the values `node_<i>.x`.

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

    Examples:
        >>> from skcomplex.systems.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)

read_biopython_structure(structure, systems=None)

Return a structure as an object of Biopython.

Parameters:

Name	Type	Description	Default
structure	StructureSpec	The structure.	required
systems	tuple[int, ...] | None	The indices of the systems to keep, or None for every system.	None

Returns:

Name	Type	Description
object	BiopythonStructureProtocol	The Bio.PDB.Structure of the atoms in their order.

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.

Examples:

>>> from skcomplex.systems.biology.systems_biology import read_biopython_structure, read_structure
>>> len(list(read_biopython_structure(read_structure('1L2Y'), (0, 1)).get_chains()))
2

Source code in src/skcomplex/systems/biology/systems_biology/_structures.py

def read_biopython_structure(
    structure: StructureSpec,
    systems: tuple[int, ...] | None = None,
) -> BiopythonStructureProtocol:
    """Return a structure as an object of Biopython.

    Args:
        structure:
            The structure.

        systems:
            The indices of the systems to keep, or `None` for every system.

    Returns:
        object:
            The `Bio.PDB.Structure` of the atoms in their order.

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

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import read_biopython_structure, read_structure
        >>> len(list(read_biopython_structure(read_structure('1L2Y'), (0, 1)).get_chains()))
        2
    """
    pdb = import_extra('Bio.PDB', STRUCTURE_SYSTEMS, 'structures')
    builder = pdb.StructureBuilder.StructureBuilder()
    builder.init_structure('structure')
    chosen = range(len(structure.coordinates)) if systems is None else systems
    for serial, index in enumerate(chosen):
        builder.init_model(serial, serial + 1)
        chain, number = None, None
        builder.init_seg(' ')
        for atom, (name, residue, position, owner) in enumerate(
            zip(
                structure.atom_names,
                structure.residue_names,
                structure.residue_numbers,
                structure.chain_names,
                strict=True,
            ),
        ):
            if owner != chain:
                chain, number = owner, None
                builder.init_chain(owner)
            if position != number:
                number = position
                builder.init_residue(residue, POLYMER, position, POLYMER)
            builder.init_atom(
                name,
                structure.coordinates[index, atom].astype(np.float32),
                0.0,
                1.0,
                POLYMER,
                name if len(name) == FULL_NAME else f' {name:<3}',
                atom + 1,
                name[0],
            )
    return builder.get_structure()

read_cell_probabilities(cells=CELLS) cached

Return the probability of every cell of the ellipse path per amino acid, from the native bundled structures.

The projected native angles of the bundled entries are counted per amino acid and cell, with one added to every count. The bundled entries and the added one are choices of this function, since Konieczny et al. (2023 §3.7.1) print entropies of such probabilities (Table 3.3) and not the probabilities or how they are estimated.

Parameters:

Name	Type	Description	Default
cells	int	The number of cells of the path.	CELLS

Returns:

Name	Type	Description
probabilities	Array	One row per amino acid, in the order of AMINO_ACIDS, and one column per cell.

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.

Examples:

>>> from skcomplex.systems.biology.systems_biology import read_cell_probabilities
>>> read_cell_probabilities().sum(axis=1).round(9).tolist() == [1.0] * 20
True

Source code in src/skcomplex/systems/biology/systems_biology/_early_stage.py

@lru_cache(maxsize=4)
def read_cell_probabilities(cells: int = CELLS) -> Array:
    """Return the probability of every cell of the ellipse path per amino acid, from the native bundled structures.

    The projected native angles of the bundled entries are counted per amino acid and cell, with one added to every
    count. The bundled entries and the added one are choices of this function, since Konieczny et al. (2023 §3.7.1)
    print entropies of such probabilities (Table 3.3) and not the probabilities or how they are estimated.

    Args:
        cells:
            The number of cells of the path.

    Returns:
        probabilities:
            One row per amino acid, in the order of `AMINO_ACIDS`, and one column per cell.

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

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import read_cell_probabilities
        >>> read_cell_probabilities().sum(axis=1).round(9).tolist() == [1.0] * 20
        True
    """
    counts = np.ones((len(AMINO_ACIDS), cells))
    path = EllipsePath(cells=cells).prepare()
    for entry in BUNDLED:
        names, pairs = read_native_pairs(read_structure(entry))
        for name, cell in zip(names, path.run(pairs).cells_, strict=True):
            counts[AMINO_ACIDS.index(name), cell] += 1
    return counts / counts.sum(axis=1, keepdims=True)

read_kyte_doolittle_scale()

Return the Kyte-Doolittle hydrophobicity of the twenty amino acids of Biopython, rescaled to [0, 1].

Returns:

Name	Type	Description
scale	dict[str, float]	The hydrophobicity by the three-letter name of a residue.

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.

Examples:

>>> from skcomplex.systems.biology.systems_biology import read_kyte_doolittle_scale
>>> scale = read_kyte_doolittle_scale()
>>> scale['ILE'], scale['ARG']
(1.0, 0.0)

Source code in src/skcomplex/systems/biology/systems_biology/_fuzzy_oil_drop.py

def read_kyte_doolittle_scale() -> dict[str, float]:
    """Return the Kyte-Doolittle hydrophobicity of the twenty amino acids of Biopython, rescaled to `[0, 1]`.

    Returns:
        scale:
            The hydrophobicity by the three-letter name of a residue.

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

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import read_kyte_doolittle_scale
        >>> scale = read_kyte_doolittle_scale()
        >>> scale['ILE'], scale['ARG']
        (1.0, 0.0)
    """
    system = 'fuzzy_oil_drop'
    letters = import_extra('Bio.SeqUtils.ProtParamData', system, 'structures').kd
    names = import_extra('Bio.Data.IUPACData', system, 'structures').protein_letters_1to3
    low, high = min(letters.values()), max(letters.values())
    return {names[letter].upper(): (value - low) / (high - low) for letter, value in letters.items()}

read_native_pairs(structure)

Return the residues with both a phi and a psi in the first system of a structure, and these angles.

Parameters:

Name	Type	Description	Default
structure	StructureSpec	The structure, of one chain.	required

Returns:

Name	Type	Description
pairs	tuple[list[str], Array]	The names of the residues, and their phi and psi in radians, residues by 2.

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.

Examples:

>>> from skcomplex.systems.biology.systems_biology import read_native_pairs, read_structure
>>> names, angles = read_native_pairs(read_structure('1L2Y'))
>>> len(names), angles.shape
(18, (18, 2))

Source code in src/skcomplex/systems/biology/systems_biology/_early_stage.py

def read_native_pairs(structure: StructureSpec) -> tuple[list[str], Array]:
    """Return the residues with both a `phi` and a `psi` in the first system of a structure, and these angles.

    Args:
        structure:
            The structure, of one chain.

    Returns:
        pairs:
            The names of the residues, and their `phi` and `psi` in radians, residues by 2.

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

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import read_native_pairs, read_structure
        >>> names, angles = read_native_pairs(read_structure('1L2Y'))
        >>> len(names), angles.shape
        (18, (18, 2))
    """
    angles = BackboneGeometry(structure).prepare().angles_
    split = angles.size // 2
    names = [residue for _, _, residue in structure.residues][1:-1]
    pairs = np.column_stack([angles[: split - 1], angles[split + 1 :]])
    return names, np.radians(pairs)

read_structure(source)

Return the polymer atoms of every system of a bundled entry, or of a PDB or an mmCIF file.

Parameters:

Name	Type	Description	Default
source	str | Path	The name of a bundled entry, '1L2Y', '1UBQ', or '5PTI', or the path of a .pdb, .cif, or .mmcif file.	required

Returns:

Name	Type	Description
structure	StructureSpec	The structure.

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.
ValidationError	The source is no bundled entry and no file of a known format.

Examples:

>>> from skcomplex.systems.biology.systems_biology import read_structure
>>> read_structure('1L2Y').coordinates.shape
(38, 304, 3)

Source code in src/skcomplex/systems/biology/systems_biology/_structures.py

def read_structure(source: str | Path) -> StructureSpec:
    """Return the polymer atoms of every system of a bundled entry, or of a PDB or an mmCIF file.

    Args:
        source:
            The name of a bundled entry, `'1L2Y'`, `'1UBQ'`, or `'5PTI'`, or the path of a `.pdb`, `.cif`, or `.mmcif`
            file.

    Returns:
        structure:
            The structure.

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

        ValidationError:
            The source is no bundled entry and no file of a known format.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import read_structure
        >>> read_structure('1L2Y').coordinates.shape
        (38, 304, 3)
    """
    if source in BUNDLED:
        with gzip.open(str(files(__package__) / f'{source}.cif.gz'), 'rt') as handle:
            return _read_systems(handle, is_pdb=False)
    path = Path(source)
    if path.suffix == SUFFIX_PDB or path.suffix in SUFFIXES_MMCIF:
        return _read_systems(str(path), is_pdb=path.suffix == SUFFIX_PDB)
    error_msg = f'{source!r} is no bundled entry, {", ".join(BUNDLED)}, and no {SUFFIX_PDB} or .cif file.'
    raise ValidationError(error_msg)

read_trajectory(path, topology)

Return the frames of a molecular dynamics trajectory as the systems of a structure, through MDTraj.

Parameters:

Name	Type	Description	Default
path	str | Path	The path of a DCD, XTC, TRR, NetCDF, or HDF5 trajectory.	required
topology	str | Path	The path of a file that gives the atoms, such as a PDB file.	required

Returns:

Name	Type	Description
structure	StructureSpec	The structure, with one system per frame, in ångström.

Raises:

Type	Description
MissingExtraError	The extra trajectories is not installed.

Source code in src/skcomplex/systems/biology/systems_biology/_structures.py

def read_trajectory(path: str | Path, topology: str | Path) -> StructureSpec:
    """Return the frames of a molecular dynamics trajectory as the systems of a structure, through MDTraj.

    Args:
        path:
            The path of a DCD, XTC, TRR, NetCDF, or HDF5 trajectory.

        topology:
            The path of a file that gives the atoms, such as a PDB file.

    Returns:
        structure:
            The structure, with one system per frame, in ångström.

    Raises:
        MissingExtraError:
            The extra `trajectories` is not installed.
    """
    trajectory = import_extra('mdtraj', STRUCTURE_SYSTEMS, 'trajectories').load(str(path), top=str(topology))
    atoms = list(trajectory.topology.atoms)
    return StructureSpec(
        np.asarray(trajectory.xyz, dtype=float) * NM_TO_ANGSTROM,
        tuple(atom.name for atom in atoms),
        tuple(atom.residue.name for atom in atoms),
        tuple(int(atom.residue.resSeq) for atom in atoms),
        tuple(str(atom.residue.chain.chain_id or atom.residue.chain.index) for atom in atoms),
    )._validated()

regulatory_fba(*, biomass=0.01, glucose=10.0, acetate=10.0, threshold=0.1, lag=10, acetate_vmax=10.0, km=0.015)

Return a batch culture on glucose and acetate with Boolean catabolite repression, which needs the extra.

Boolean rules alternate with the program (Raman 2021, Eq. 9.7). A node reads whether glucose lies above a threshold, a chain of nodes delays it, and while the last node is on it holds the genes of acetate kinase off through update laws. So the culture leaves acetate alone while glucose lasts, stops at depletion, and grows again on acetate after the lag of the chain.

The alternation is that of Raman, but the rules are an illustration and not those of Covert et al. (2001): the genes of acetate kinase, the threshold, and the chain of nodes come from no source, and the diauxie of Raman is on glucose and lactose, which the E. coli core system cannot take up.

Parameters:

Name	Type	Description	Default
biomass	float	The initial biomass in gDW per L.	0.01
glucose	float	The initial glucose in mM.	10.0
acetate	float	The initial acetate in mM.	10.0
threshold	float	The glucose above which the sensor is on, in mM.	0.1
lag	int	The number of nodes after the sensor, each a synchronisation interval of delay.	10
acetate_vmax	float	The largest acetate uptake in mmol per gDW per h.	10.0
km	float	The half-saturation concentration of both substrates in mM.	0.015

Returns:

Name	Type	Description
system	System	The system.

Raises:

Type	Description
MissingExtraError	The extra is not installed.

Examples:

>>> from skcomplex.systems.biology.systems_biology import regulatory_fba
>>> [system.label for system in regulatory_fba(lag=1).systems]
['control', 'culture']

Source code in src/skcomplex/systems/biology/systems_biology/_regulatory_fba.py

@register_system(
    field='biology/systems_biology',
    tags=('regulatory flux balance', 'boolean network', 'linear program', 'metabolism', 'diauxie', 'co-simulation'),
    references=('Raman 2021', 'Covert et al. 2001'),
    defaults={
        'until': 8.0,
        'dt': SYNCHRONISATION,
        'solver': CoSimulationSolver(sync_dt=SYNCHRONISATION, solvers={'culture': RK4Solver(dt=ODE_STEP)}),
    },
)
def regulatory_fba(
    *,
    biomass: float = 0.01,
    glucose: float = 10.0,
    acetate: float = 10.0,
    threshold: float = 0.1,
    lag: int = 10,
    acetate_vmax: float = 10.0,
    km: float = 0.015,
) -> System:
    """Return a batch culture on glucose and acetate with Boolean catabolite repression, which needs the extra.

    Boolean rules alternate with the program (Raman 2021, Eq. 9.7). A node reads whether glucose lies above a threshold,
    a chain of nodes delays it, and while the last node is on it holds the genes of acetate kinase off through update
    laws. So the culture leaves acetate alone while glucose lasts, stops at depletion, and grows again on acetate after
    the lag of the chain.

    The alternation is that of Raman, but the rules are an illustration and not those of Covert et al. (2001): the genes
    of acetate kinase, the threshold, and the chain of nodes come from no source, and the diauxie of Raman is on glucose
    and lactose, which the E. coli core system cannot take up.

    Args:
        biomass:
            The initial biomass in gDW per L.

        glucose:
            The initial glucose in mM.

        acetate:
            The initial acetate in mM.

        threshold:
            The glucose above which the sensor is on, in mM.

        lag:
            The number of nodes after the sensor, each a synchronisation interval of delay.

        acetate_vmax:
            The largest acetate uptake in mmol per gDW per h.

        km:
            The half-saturation concentration of both substrates in mM.

    Returns:
        system:
            The system.

    Raises:
        MissingExtraError:
            The extra is not installed.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import regulatory_fba
        >>> [system.label for system in regulatory_fba(lag=1).systems]
        ['control', 'culture']
    """
    cell = load_network('textbook', 'cell', needed_by='regulatory_fba')
    cell.parsimonious, cell.infeasible = True, 'zero'
    identifier = cell._identifier(ACETATE_EXCHANGE)
    cell.set_params(
        **{f'lower_{identifier}': -float(acetate_vmax)},
        **{cell.gene_name(gene): 0.0 for gene in ACETATE_GENES},
    )
    nodes = [BooleanNodeElement(f'node_{position}', x=1) for position in range(lag + 1)]
    chain = [
        GlucoseSensorInteraction(
            f'delay_{position}',
            node=f'control/node_{position}',
            source=f'control/node_{position - 1}',
        )
        for position in range(1, lag + 1)
    ]
    return System(
        systems=[
            System('control', elements=[cell, *nodes], time=DiscreteTime()),
            System(
                'culture',
                elements=[
                    BiomassElement('biomass', x=biomass),
                    SubstrateElement('glucose', c=glucose),
                    SubstrateElement('acetate', c=acetate),
                ],
            ),
        ],
        interactions=[
            GlucoseSensorInteraction(
                'sensor',
                node='control/node_0',
                source='culture/glucose',
                threshold=threshold,
            ),
            *chain,
            SwitchInteraction('repression', node=f'control/node_{lag}', cell='control/cell', genes=ACETATE_GENES),
            ExchangeInteraction(
                'glucose_uptake',
                cell='control/cell',
                biomass='culture/biomass',
                substrate='culture/glucose',
                reaction=GLUCOSE_EXCHANGE,
                km=km,
                interval=SYNCHRONISATION,
            ),
            ExchangeInteraction(
                'acetate_uptake',
                cell='control/cell',
                biomass='culture/biomass',
                substrate='culture/acetate',
                reaction=ACETATE_EXCHANGE,
                vmax=acetate_vmax,
                km=km,
                interval=SYNCHRONISATION,
            ),
            GrowthInteraction('growth', cell='control/cell', biomass='culture/biomass', reaction=GROWTH),
        ],
    )

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 system is a jump process.	False

Returns:

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

Examples:

>>> from skcomplex.systems.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/systems/biology/systems_biology/_genes.py

@register_system(
    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 system is a jump process.

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

    Examples:
        >>> from skcomplex.systems.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)

tca_labelling_analysis(kind='estimation', *, noise=DEVIATION, random_state=0, starts=5, intervals='linear')

Return the labelling analysis of the TCA example, with data simulated from its known fluxes plus seeded noise.

Parameters:

Name	Type	Description	Default
kind	str	'simulation' for the published tracer at the true fluxes, or 'estimation' for a fit to two parallel tracer experiments, 50% [1-13C] and 50% [1,2-13C] acetyl-CoA, with v1 fixed at 100.	'estimation'
noise	float	The standard deviation of the noise of the data.	DEVIATION
random_state	int	The seed of the noise and of the starts.	0
starts	int	The number of feasible starts of the fit.	5
intervals	str	'linear' or 'profile'.	'linear'

Returns:

Name	Type	Description
analysis	LabellingAnalysis	The analysis, not checked yet.

Examples:

>>> from skcomplex.systems.biology.systems_biology import carbon_13_mfa, tca_labelling_analysis
>>> analysis = tca_labelling_analysis('simulation').prepare(carbon_13_mfa()).run()
>>> [round(float(item), 4) for item in analysis.distributions_['published']['Glu']]
[0.3464, 0.2695, 0.2708, 0.0807, 0.0286, 0.0039]

Source code in src/skcomplex/systems/biology/systems_biology/_carbon_13_mfa.py

def tca_labelling_analysis(
    kind: str = 'estimation',
    *,
    noise: float = DEVIATION,
    random_state: int = 0,
    starts: int = 5,
    intervals: str = 'linear',
) -> LabellingAnalysis:
    """Return the labelling analysis of the TCA example, with data simulated from its known fluxes plus seeded noise.

    Args:
        kind:
            `'simulation'` for the published tracer at the true fluxes, or `'estimation'` for a fit to two parallel
            tracer experiments, 50% [1-13C] and 50% [1,2-13C] acetyl-CoA, with `v1` fixed at 100.

        noise:
            The standard deviation of the noise of the data.

        random_state:
            The seed of the noise and of the starts.

        starts:
            The number of feasible starts of the fit.

        intervals:
            `'linear'` or `'profile'`.

    Returns:
        analysis:
            The analysis, not checked yet.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import carbon_13_mfa, tca_labelling_analysis
        >>> analysis = tca_labelling_analysis('simulation').prepare(carbon_13_mfa()).run()
        >>> [round(float(item), 4) for item in analysis.distributions_['published']['Glu']]
        [0.3464, 0.2695, 0.2708, 0.0807, 0.0286, 0.0039]
    """
    if kind == 'simulation':
        return LabellingAnalysis(
            'simulation',
            atoms=TCA_ATOMS,
            tracers={'published': PUBLISHED_TRACERS},
            fragments=(GLUTAMATE,),
            fluxes=TCA_TRUE_FLUXES,
            symmetric=TCA_SYMMETRIC,
        )
    truth = LabellingAnalysis(
        'simulation',
        atoms=TCA_ATOMS,
        tracers=PARALLEL_TRACERS,
        fragments=(GLUTAMATE,),
        fluxes=TCA_TRUE_FLUXES,
        symmetric=TCA_SYMMETRIC,
    )
    simulated = truth.prepare(carbon_13_mfa()).run().distributions_
    rng = np.random.default_rng(random_state)
    measurements = []
    for experiment, values in simulated.items():
        clean = values[GLUTAMATE.label]
        noisy = clean + rng.normal(0.0, noise, clean.size) if noise else clean
        measurements.append(LabelMeasurementSpec(experiment, GLUTAMATE.label, noisy, DEVIATION))
    return LabellingAnalysis(
        'estimation',
        atoms=TCA_ATOMS,
        tracers=PARALLEL_TRACERS,
        fragments=(GLUTAMATE,),
        measurements=measurements,
        measured_fluxes={'v1': (100.0, 0.0)},
        symmetric=TCA_SYMMETRIC,
        starts=starts,
        intervals=intervals,
        random_state=random_state,
    )

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 system, with the values g1.x, g2.x, and g3.x.

Examples:

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

Source code in src/skcomplex/systems/biology/systems_biology/_boolean_systems.py

@register_system(
    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 system, with the values `g1.x`, `g2.x`, and `g3.x`.

    Examples:
        >>> from skcomplex.systems.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)

three_way_fate_decision(cells=30, *, a=-0.2, theta1=0.0, theta2=0.0, v1=0.0, v2=0.0, spread=0.6, seed=0)

Return cells that choose among three fates on the compactified elliptic umbilic landscape of Rand et al.

Every cell descends the potential x^4 + y^4 + x^3 - 2 x y^2 + a (x^2 + y^2) + theta1 x + theta2 y. With a < 0 and no signal the centre is a repellor and three attractors, the fates, lie around it, separated by three saddles, so cells that start near the centre split among the three fates. Signals that drift at the rates v1 and v2 move the system across the fold curves of its parameter space, where a fate disappears and its cells switch to another.

Parameters:

Name	Type	Description	Default
cells	int	The number of cells.	30
a	float	The curvature of the potential at the centre.	-0.2
theta1	float	The first signal at the start.	0.0
theta2	float	The second signal at the start.	0.0
v1	float	The rate of change of the first signal.	0.0
v2	float	The rate of change of the second signal.	0.0
spread	float	The half-width of the uniform initial states around the centre.	0.6
seed	int | None	The seed of the initial states.	0

Returns:

Name	Type	Description
system	System	The system, with the states cells.x and cells.y.

Examples:

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

Source code in src/skcomplex/systems/biology/systems_biology/_fates.py

@register_system(
    field='biology/systems_biology',
    tags=('cell fate', 'morse-smale', 'elliptic umbilic', 'landscape', 'ode'),
    references=('Rand et al. 2021',),
    defaults={'until': 40.0, 'dt': 0.2},
)
def three_way_fate_decision(
    cells: int = 30,
    *,
    a: float = -0.2,
    theta1: float = 0.0,
    theta2: float = 0.0,
    v1: float = 0.0,
    v2: float = 0.0,
    spread: float = 0.6,
    seed: int | None = 0,
) -> System:
    """Return cells that choose among three fates on the compactified elliptic umbilic landscape of Rand et al.

    Every cell descends the potential `x^4 + y^4 + x^3 - 2 x y^2 + a (x^2 + y^2) + theta1 x + theta2 y`. With `a < 0`
    and no signal the centre is a repellor and three attractors, the fates, lie around it, separated by three saddles,
    so cells that start near the centre split among the three fates. Signals that drift at the rates `v1` and `v2`
    move the system across the fold curves of its parameter space, where a fate disappears and its cells switch to
    another.

    Args:
        cells:
            The number of cells.

        a:
            The curvature of the potential at the centre.

        theta1:
            The first signal at the start.

        theta2:
            The second signal at the start.

        v1:
            The rate of change of the first signal.

        v2:
            The rate of change of the second signal.

        spread:
            The half-width of the uniform initial states around the centre.

        seed:
            The seed of the initial states.

    Returns:
        system:
            The system, with the states `cells.x` and `cells.y`.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import three_way_fate_decision
        >>> three_way_fate_decision(5).elements[0].prepare().size_
        5
    """
    rng = np.random.default_rng(seed)
    start = rng.uniform(-spread, spread, (2, cells))
    values = {'x': start[0], 'y': start[1], 'a': a, 'theta1': theta1, 'theta2': theta2, 'v1': v1, 'v2': v2}
    return System(elements=[GroupElement('cells', element=FateLandscapeElement, size=cells, values=values)])

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 system is a jump process.	False

Returns:

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

Examples:

>>> from skcomplex.systems.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/systems/biology/systems_biology/_genes.py

@register_system(
    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 system is a jump process.

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

    Examples:
        >>> from skcomplex.systems.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)

write_structure(structure, path)

Write every system of a structure to a .pdb, .cif, or .mmcif file.

Parameters:

Name	Type	Description	Default
structure	StructureSpec	The structure.	required
path	str | Path	The path of the file, whose suffix gives the format.	required

Raises:

Type	Description
MissingExtraError	The extra structures is not installed.
ValidationError	The suffix of the path names no known format.

Examples:

>>> import tempfile
>>> from pathlib import Path
>>> from skcomplex.systems.biology.systems_biology import read_structure, write_structure
>>> with tempfile.TemporaryDirectory() as folder:
...     write_structure(read_structure('1UBQ'), Path(folder) / 'copy.cif')
...     read_structure(Path(folder) / 'copy.cif').coordinates.shape
(1, 602, 3)

Source code in src/skcomplex/systems/biology/systems_biology/_structures.py

def write_structure(structure: StructureSpec, path: str | Path) -> None:
    """Write every system of a structure to a `.pdb`, `.cif`, or `.mmcif` file.

    Args:
        structure:
            The structure.

        path:
            The path of the file, whose suffix gives the format.

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

        ValidationError:
            The suffix of the path names no known format.

    Examples:
        >>> import tempfile
        >>> from pathlib import Path
        >>> from skcomplex.systems.biology.systems_biology import read_structure, write_structure
        >>> with tempfile.TemporaryDirectory() as folder:
        ...     write_structure(read_structure('1UBQ'), Path(folder) / 'copy.cif')
        ...     read_structure(Path(folder) / 'copy.cif').coordinates.shape
        (1, 602, 3)
    """
    target = Path(path)
    pdb = import_extra('Bio.PDB', STRUCTURE_SYSTEMS, 'structures')
    if target.suffix != SUFFIX_PDB and target.suffix not in SUFFIXES_MMCIF:
        error_msg = f'The structure file {target} has no suffix of a known format, {SUFFIX_PDB} or .cif.'
        raise ValidationError(error_msg)
    writer = pdb.PDBIO() if target.suffix == SUFFIX_PDB else pdb.MMCIFIO()
    writer.set_structure(read_biopython_structure(structure))
    writer.save(str(target))

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 system, with the values <protein>.x.

Examples:

>>> from skcomplex.systems.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/systems/biology/systems_biology/_boolean.py

@register_system(
    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 system, with the values `<protein>.x`.

    Examples:
        >>> from skcomplex.systems.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)