Skip to content

Physiology

Physiology: delayed feedback, excitable tissue, the circulation, cell exchange, regulation, adaptivity, tissues.

AdaptivePhaseCouplingInteraction(label, *, cells, links, N=1.0, alpha=0.0, beta=0.0, epsilon=0.01)

Bases: NeighborInteraction

Phase coupling of the cells of a layer through links that adapt to the phase differences they carry.

Along the link from cell j to cell i, dphi_i/dt += -(a_ij + kappa_ij) sin(phi_i - phi_j + alpha) / N and dkappa_ij/dt = -epsilon (kappa_ij + sin(phi_i - phi_j - beta)). The phase lag alpha sets the preferred phase difference of the coupling, and the age parameter beta the rule of the adaptation.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
cells Annotated[str, RoleSpec()]

The path of the group of cells on the nodes.

 required
links Annotated[str, RoleSpec()]

The path of the group of links on the edges of the layer.

 required
N Annotated[float, ParameterSpec(domain=open(0, oo))]

The number of cells of the layer, which normalizes the coupling.

 1.0
alpha Annotated[float, ParameterSpec()]

The phase lag of the coupling.

 0.0
beta Annotated[float, ParameterSpec()]

The age parameter of the adaptation.

 0.0
epsilon Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The adaptation rate.

 0.01

Examples:

>>> from skcomplex.systems.biology.physiology import AdaptivePhaseCouplingInteraction
>>> AdaptivePhaseCouplingInteraction('coupling', cells='cells', links='links', epsilon=0.3).epsilon
0.3
Source code in src/skcomplex/systems/biology/physiology/_disease.py 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
def __init__(
    self,
    label: str,
    *,
    cells: Annotated[str, RoleSpec()],
    links: Annotated[str, RoleSpec()],
    N: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    alpha: Annotated[float, ParameterSpec()] = 0.0,
    beta: Annotated[float, ParameterSpec()] = 0.0,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
) -> None:
    super().__init__(label)
    self.cells = cells
    self.links = links
    self.N = N
    self.alpha = alpha
    self.beta = beta
    self.epsilon = epsilon

AdaptiveResponderElement(label, visual=None, *, e=0.0, c=0.2, S=1.0, response_cost=0.0, restore_cost=0.0, repair_cost=0.0, prepare_cost=0.0, eta=0.5, anticipation=0.0, c0=0.2, tau=50.0, k0=0.1, kappa=0.2, price=1.0, theta=1.0, repair=5.0, upkeep=0.05, height=1.0, width=2.0, period=10.0, start=1.0)

Bases: Element

System that counters episodes of demand on a critical parameter and pays for its response in four costs.

The excursion of the critical parameter follows de/dt = (1 - a) s(t) - c e - k0 e, where the stimulus s is a train of episodes of height height and width width every period from start, c e is the response of a mechanism of capacity c, and a is the share of every episode that anticipation avoids. The capacity grows with the demand and decays to its baseline, dc/dt = eta e^2 - (c - c0) / tau, and the response draws the reserves, which are restored, dS/dt = -c e + kappa (1 - S). The four costs of Mobus accumulate as stocks: the response cost is the response effort c e, the restore cost the work price kappa (1 - S) of refilling the reserves, the repair cost repair (e - theta) while the excursion exceeds the damage threshold theta, and the preparedness cost the upkeep upkeep (c - c0 + a) of the extra capacity and of the anticipation.

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

Parameters:

Name Type Description Default
label str

The label of the system.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
e Annotated[float, StateSpec()]

The initial excursion of the critical parameter.

 0.0
c Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial capacity of the response.

 0.2
S Annotated[float, StateSpec()]

The initial reserves, 1 when full.

 1.0
response_cost Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial response cost.

 0.0
restore_cost Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial restore cost.

 0.0
repair_cost Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial repair cost.

 0.0
prepare_cost Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial preparedness cost.

 0.0
eta Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate at which demand builds capacity, 0 for a fixed stimulus-response.

 0.5
anticipation Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The share a of every episode that anticipation avoids.

 0.0
c0 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The baseline capacity.

 0.2
tau Annotated[float, ParameterSpec(domain=open(0, oo))]

The time over which unused capacity decays.

 50.0
k0 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The homeostatic relaxation rate of the excursion.

 0.1
kappa Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate of restoring the reserves.

 0.2
price Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The cost of restoring a unit of reserves.

 1.0
theta Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The excursion above which damage accrues.

 1.0
repair Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The repair cost per unit of excursion above the threshold and per unit time.

 5.0
upkeep Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The cost per unit time of a unit of extra capacity or of anticipation.

 0.05
height float

The height of the episodes.

 1.0
width float

The duration of the episodes.

 2.0
period float

The time between the starts of episodes.

 10.0
start float

The start of the first episode.

 1.0

Examples:

>>> from skcomplex.systems.biology.physiology import AdaptiveResponderElement
>>> AdaptiveResponderElement('organism').period
10.0
Source code in src/skcomplex/systems/biology/physiology/_adaptivity.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    e: Annotated[float, StateSpec()] = 0.0,
    c: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    S: Annotated[float, StateSpec()] = 1.0,
    response_cost: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    restore_cost: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    repair_cost: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    prepare_cost: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    anticipation: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    c0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 50.0,
    k0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    price: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    theta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    repair: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0,
    upkeep: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    height: float = 1.0,
    width: float = 2.0,
    period: float = 10.0,
    start: float = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.e = e
    self.c = c
    self.S = S
    self.response_cost = response_cost
    self.restore_cost = restore_cost
    self.repair_cost = repair_cost
    self.prepare_cost = prepare_cost
    self.eta = eta
    self.anticipation = anticipation
    self.c0 = c0
    self.tau = tau
    self.k0 = k0
    self.kappa = kappa
    self.price = price
    self.theta = theta
    self.repair = repair
    self.upkeep = upkeep
    self.height = height
    self.width = width
    self.period = period
    self.start = start

AlievPanfilovElement(label, visual=None, *, u=0.0, v=0.0, k=8.0, a=0.15, epsilon0=0.002, mu1=0.2, mu2=0.3)

Bases: Element

Cardiac tissue of Aliev and Panfilov, an excitable medium with the restitution of heart muscle.

The potential follows u' = -k u (u - a) (u - 1) - u v and the recovery v' = (epsilon0 + mu1 v / (u + mu2)) (-v - k u (u - a - 1)), and a diffusion of the potential couples the cells of the tissue. The recovery rate depends on the state, so the action potential lasts longer than in the FitzHugh-Nagumo medium and its duration adapts to the pacing rate.

Parameters:

Name Type Description Default
label str

The label of the tissue.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
u Annotated[float, StateSpec()]

The initial potential.

 0.0
v Annotated[float, StateSpec()]

The initial recovery.

 0.0
k Annotated[float, ParameterSpec(domain=open(0, oo))]

The strength of the excitation.

 8.0
a Annotated[float, ParameterSpec(domain=open(0, oo))]

The excitation threshold.

 0.15
epsilon0 Annotated[float, ParameterSpec(domain=open(0, oo))]

The recovery rate far from the excited state.

 0.002
mu1 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The growth of the recovery rate with the recovery.

 0.2
mu2 Annotated[float, ParameterSpec(domain=open(0, oo))]

The offset of the recovery rate.

 0.3

Examples:

>>> from skcomplex.systems.biology.physiology import AlievPanfilovElement
>>> AlievPanfilovElement('tissue').k
8.0
Source code in src/skcomplex/systems/biology/physiology/_cardiac.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec()] = 0.0,
    v: Annotated[float, StateSpec()] = 0.0,
    k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 8.0,
    a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.15,
    epsilon0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.002,
    mu1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    mu2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.3,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u
    self.v = v
    self.k = k
    self.a = a
    self.epsilon0 = epsilon0
    self.mu1 = mu1
    self.mu2 = mu2

BarkleyElement(label, visual=None, *, u=0.0, w=0.0, a=0.6, b=0.05, epsilon=0.02, gamma=1.0)

Bases: Element

Excitable medium of Barkley, a fast excitation u and a slow recovery w with a polynomial nullcline.

The excitation follows u' = u (1 - u) (u - (w + b) / a) and the recovery w' = epsilon (u^gamma - w). A diffusion of the excitation couples the sites of the medium. The threshold (w + b) / a rises with the recovery, so an excited site recovers before it can fire again.

Parameters:

Name Type Description Default
label str

The label of the medium.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
u Annotated[float, StateSpec()]

The initial excitation.

 0.0
w Annotated[float, StateSpec()]

The initial recovery.

 0.0
a Annotated[float, ParameterSpec(domain=open(0, oo))]

The inverse slope of the threshold.

 0.6
b Annotated[float, ParameterSpec()]

The offset of the threshold.

 0.05
epsilon Annotated[float, ParameterSpec(domain=open(0, oo))]

The ratio of the time scales of the recovery and the excitation.

 0.02
gamma Annotated[float, ParameterSpec(domain=open(0, oo))]

The exponent of the excitation in the recovery, 1 for Barkley's model and 3 for the variant of Bär and Eiswirth with spiral breakup.

 1.0

Examples:

>>> from skcomplex.systems.biology.physiology import BarkleyElement
>>> BarkleyElement('medium').a
0.6
Source code in src/skcomplex/systems/biology/physiology/_cardiac.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec()] = 0.0,
    w: Annotated[float, StateSpec()] = 0.0,
    a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.6,
    b: Annotated[float, ParameterSpec()] = 0.05,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.02,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u
    self.w = w
    self.a = a
    self.b = b
    self.epsilon = epsilon
    self.gamma = gamma

BranchingVesselElement(label, visual=None, *, r=1.0, l=1.0, volume=1.0, terminal=0.0, level=1.0, depth=1.0, n=2.0, r_c=1.0, l_c=1.0)

Bases: Element

Vessel of a space-filling branching network, sized by its level so that the network minimizes transport energy.

Every vessel branches into n daughters. Area-preserving branching gives the radii the ratio n^(-1/2) between levels and space filling the lengths the ratio n^(-1/3), and the capillaries at the last level depth have the same radius r_c and length l_c in every organism, so r = r_c n^((depth - level) / 2), l = l_c n^((depth - level) / 3), and the volume of the vessel is pi r^2 l.

Parameters:

Name Type Description Default
label str

The label of the vessel or of the group on the edges.

 required
visual VisualSpec | None

How the vessel is drawn, or None for default_visual_.

 None
r Annotated[float, StateSpec(domain=closed_open(0, oo))]

The radius.

 1.0
l Annotated[float, StateSpec(domain=closed_open(0, oo))]

The length.

 1.0
volume Annotated[float, StateSpec(domain=closed_open(0, oo))]

The volume.

 1.0
terminal Annotated[float, StateSpec(domain=closed(0, 1))]

Whether the vessel is a capillary.

 0.0
level Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The level of the vessel, 1 for the aorta's daughters and depth for the capillaries.

 1.0
depth Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The level of the capillaries.

 1.0
n Annotated[float, ParameterSpec(domain=open(1, oo))]

The branching ratio.

 2.0
r_c Annotated[float, ParameterSpec(domain=open(0, oo))]

The radius of a capillary.

 1.0
l_c Annotated[float, ParameterSpec(domain=open(0, oo))]

The length of a capillary.

 1.0

Examples:

>>> from skcomplex.systems.biology.physiology import BranchingVesselElement
>>> BranchingVesselElement('vessel', level=3.0).level
3.0
Source code in src/skcomplex/systems/biology/physiology/_distribution_networks.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    r: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    l: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    volume: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    terminal: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.0,
    level: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    depth: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(1, oo))] = 2.0,
    r_c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    l_c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.r = r
    self.l = l
    self.volume = volume
    self.terminal = terminal
    self.level = level
    self.depth = depth
    self.n = n
    self.r_c = r_c
    self.l_c = l_c

CellAdhesionInteraction(label, *, target, kinds, adhesion, temperature=1.0, area_strength=1.0, target_area=40.0, perimeter_strength=0.0, target_perimeter=0.0, state='cell', rule='metropolis', neighborhood=None, hits=1)

Bases: SamplingInteraction

Cellular Potts sampling of a tissue, whose cells move their boundaries under adhesion and area constraints.

The energy of the Graner-Glazier model is sum_<i,j> w_ij J(t(s_i), t(s_j)) [s_i != s_j] + lambda sum_s (a_s - A)^2 + mu sum_s (p_s - P)^2, where s_i is the cell of site i, t(s) the kind of cell s, J the adhesion energy per unit contact between kinds, w_ij the weight of the contact of two sites, a_s the area of a cell, and p_s its perimeter, the weighted number of its contacts with other cells. Every sweep proposes hits copies per site of the cell of a random neighbour into a random site, accepted by the Metropolis rule at the temperature. A heterotypic adhesion above the mean of the homotypic ones sorts the kinds, the differential adhesion of Steinberg.

Parameters:

Name Type Description Default
label str

The label of the interaction, unique among the parts of its system.

 required
target Annotated[str, RoleSpec()]

The path of the group of sites.

 required
kinds Sequence[int]

The kind of every cell, by label.

 required
adhesion ArrayLike

The adhesion energy per unit contact between cells of every pair of kinds, a symmetric matrix.

 required
temperature Annotated[float, ParameterSpec(domain=open(0, oo))]

The temperature of the copies, the amplitude of the fluctuations of the cell boundaries.

 1.0
area_strength Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The strength lambda of the area constraint.

 1.0
target_area Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The target area A of every cell.

 40.0
perimeter_strength Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The strength mu of the perimeter constraint.

 0.0
target_perimeter Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The target perimeter P of every cell.

 0.0
state str

The name of the state that holds the cell of a site.

 'cell'
rule str

The acceptance, 'metropolis' or 'glauber'.

 'metropolis'
neighborhood BaseNeighborhood | None

The neighbours of a site, or None for the default neighbourhood of the space.

 None
hits int

The number of copies proposed per site and sweep.

 1

Raises:

Type Description
ValidationError

The adhesion is not a symmetric square matrix, or a kind is not one of its rows.

Examples:

>>> from skcomplex.systems.biology.physiology import CellAdhesionInteraction
>>> adhesion = CellAdhesionInteraction('adhesion', target='tissue', kinds=(0, 1), adhesion=((1, 2), (2, 1)))
>>> adhesion.objective.contact(np.array([0, 0]), np.array([0, 1])).tolist()
[1.0, 2.0]
Source code in src/skcomplex/systems/biology/physiology/_cellular_potts.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    kinds: Sequence[int],
    adhesion: npt.ArrayLike,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    area_strength: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    target_area: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 40.0,
    perimeter_strength: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    target_perimeter: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    state: str = 'cell',
    rule: str = 'metropolis',
    neighborhood: BaseNeighborhood | None = None,
    hits: int = 1,
) -> None:
    super().__init__(
        label,
        target=target,
        state=state,
        objective=RegionScoreSpec(self._get_contact, self._get_area, self._get_perimeter, hits),
        move='copy',
        rule=rule,
        neighborhood=neighborhood,
    )
    self.kinds = kinds
    self.adhesion = adhesion
    self.temperature = temperature
    self.area_strength = area_strength
    self.target_area = target_area
    self.perimeter_strength = perimeter_strength
    self.target_perimeter = target_perimeter
    self.hits = hits

CellOscillatorElement(label, visual=None, *, phi=0.0, omega=0.0)

Bases: Element

Functional unit of a tissue layer reduced to a phase, a cell oscillator phi' = omega.

Healthy and pathological cells differ by their metabolism, which the natural frequency omega stands for. Couplings to other cells add to the rate of the phase.

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
phi Annotated[float, StateSpec(period=2 * pi)]

The initial phase.

 0.0
omega Annotated[float, ParameterSpec()]

The natural frequency.

 0.0

Examples:

>>> from skcomplex.systems.biology.physiology import CellOscillatorElement
>>> CellOscillatorElement('cell', omega=0.2).omega
0.2
Source code in src/skcomplex/systems/biology/physiology/_disease.py 
49
50
51
52
53
54
55
56
57
58
59
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    phi: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
    omega: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.phi = phi
    self.omega = omega

ChamberElement(label, visual=None, *, V=100.0, V0=0.0, E_min=1.0, E_max=1.0, period=0.8, systole=0.3)

Bases: Element

Elastic compartment of blood, a heart chamber or a vascular bed, whose pressure is elastance times volume.

The pressure is E(t) (V - V0). The elastance of a heart chamber rises from E_min to E_max and back as a half sine during the systole that starts every period, and the elastance of a vessel stays at E_min when E_max equals it. Volumes are in ml, pressures in mmHg, and times in s.

Parameters:

Name Type Description Default
label str

The label of the chamber.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
V Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial volume.

 100.0
V0 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The unstressed volume.

 0.0
E_min Annotated[float, ParameterSpec(domain=open(0, oo))]

The elastance in diastole, the inverse of the compliance.

 1.0
E_max Annotated[float, ParameterSpec(domain=open(0, oo))]

The largest elastance in systole.

 1.0
period Annotated[float, ParameterSpec(domain=open(0, oo))]

The duration of a heartbeat.

 0.8
systole Annotated[float, ParameterSpec(domain=open(0, oo))]

The duration of the contraction.

 0.3

Examples:

>>> from skcomplex.systems.biology.physiology import ChamberElement
>>> ChamberElement('arteries', V=150.0, E_min=0.8, E_max=0.8).E_max
0.8
Source code in src/skcomplex/systems/biology/physiology/_circulation.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    V: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 100.0,
    V0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    E_min: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    E_max: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    period: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.8,
    systole: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.3,
) -> None:
    super().__init__(label, visual=visual)
    self.V = V
    self.V0 = V0
    self.E_min = E_min
    self.E_max = E_max
    self.period = period
    self.systole = systole

ChemotacticTissueElement(label, visual=None, *, n=0.0, qx=0.0, qy=0.0, c=0.0, jx=0.0, jy=0.0, s=100.0, mu1=5000.0, mu2=5000.0, sigma=5000000.0, tau=1.0, alpha=0.33)

Bases: Element

Site of a tissue where cells and a chemical move with finite speeds, the hyperbolic limit of velocity jumps.

The cell density n and its flux q, and the chemical density c and its flux j, follow Eq. (6.50) of Bellomo et al. (2017) in two dimensions with the cells as the only source of the chemical, dn/dt = -div q, dq/dt = -(s^2 / 2) grad n - mu1 q + mu2 alpha n grad c, tau dc/dt = -div j + n, and tau dj/dt = -(s^2 / 2) grad c - sigma j. Cells and chemical move at the speed s and turn at the rates mu1 and sigma, and the cells turn towards the gradient of the chemical with the sensitivity alpha. In the limit of fast turning it is the Keller-Segel model with the diffusions s^2 / (2 mu1) and s^2 / (2 sigma).

Parameters:

Name Type Description Default
label str

The label of the tissue.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
n Annotated[float, StateSpec()]

The initial cell density.

 0.0
qx Annotated[float, StateSpec()]

The initial flux of cells along the first axis.

 0.0
qy Annotated[float, StateSpec()]

The initial flux of cells along the second axis.

 0.0
c Annotated[float, StateSpec()]

The initial chemical density.

 0.0
jx Annotated[float, StateSpec()]

The initial flux of the chemical along the first axis.

 0.0
jy Annotated[float, StateSpec()]

The initial flux of the chemical along the second axis.

 0.0
s Annotated[float, ParameterSpec(domain=open(0, oo))]

The speed of the cells and of the chemical.

 100.0
mu1 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The turning rate of the cells.

 5000.0
mu2 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate of turning of the cells towards the chemical.

 5000.0
sigma Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The turning rate of the chemical.

 5000000.0
tau Annotated[float, ParameterSpec(domain=open(0, oo))]

The time scale of the chemical.

 1.0
alpha Annotated[float, ParameterSpec()]

The chemotactic sensitivity of the cells.

 0.33

Examples:

>>> from skcomplex.systems.biology.physiology import ChemotacticTissueElement
>>> ChemotacticTissueElement('tissue').alpha
0.33
Source code in src/skcomplex/systems/biology/physiology/_chemotaxis.py 
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    n: Annotated[float, StateSpec()] = 0.0,
    qx: Annotated[float, StateSpec()] = 0.0,
    qy: Annotated[float, StateSpec()] = 0.0,
    c: Annotated[float, StateSpec()] = 0.0,
    jx: Annotated[float, StateSpec()] = 0.0,
    jy: Annotated[float, StateSpec()] = 0.0,
    s: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
    mu1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5000.0,
    mu2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5000.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0e6,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    alpha: Annotated[float, ParameterSpec()] = 0.33,
) -> None:
    super().__init__(label, visual=visual)
    self.n = n
    self.qx = qx
    self.qy = qy
    self.c = c
    self.jx = jx
    self.jy = jy
    self.s = s
    self.mu1 = mu1
    self.mu2 = mu2
    self.sigma = sigma
    self.tau = tau
    self.alpha = alpha

CoagulationElement(label, visual=None, *, T=0.0, C=0.0, F=0.0, P=0.0, ka=3.0, k_half=1.0, d_t=1.0, k_c=1.0, s_c=0.1, d_c=0.05, s_f=1.0, k_p=1.0, s_p=0.02, d_p=0.1, injury=1.0, injury_time=1.0, duration=1.0)

Bases: Element

Coagulation cascade: thrombin amplified by positive feedback, switched off by protein C, and fibrinolysis.

An injury of height injury from injury_time for duration activates thrombin, which activates more of itself through the feedback of the cascade, dT/dt = j(t) + ka T^2 / (k_half^2 + T^2) - (d_t + k_c C) T. Thrombin activates the inhibitor protein C, dC/dt = s_c T - d_c C, and converts fibrinogen to fibrin, dF/dt = s_f T - k_p P F. Fibrin activates plasmin, dP/dt = s_p F - d_p P, which dissolves the fibrin.

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

Parameters:

Name Type Description Default
label str

The label of the blood.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
T Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial thrombin.

 0.0
C Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial activated protein C.

 0.0
F Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial fibrin.

 0.0
P Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial plasmin.

 0.0
ka Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The maximal rate of the positive feedback.

 3.0
k_half Annotated[float, ParameterSpec(domain=open(0, oo))]

The thrombin at which the feedback is half its maximum.

 1.0
d_t Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The inactivation rate of thrombin by antithrombin.

 1.0
k_c Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The inactivation of thrombin per unit of protein C.

 1.0
s_c Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The activation of protein C per unit of thrombin.

 0.1
d_c Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The inactivation rate of protein C.

 0.05
s_f Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The production of fibrin per unit of thrombin.

 1.0
k_p Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The lysis of fibrin per unit of plasmin.

 1.0
s_p Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The activation of plasmin per unit of fibrin.

 0.02
d_p Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The inactivation rate of plasmin.

 0.1
injury Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate of thrombin activation by the injury.

 1.0
injury_time float

The start of the injury.

 1.0
duration float

The duration of the injury.

 1.0

Examples:

>>> from skcomplex.systems.biology.physiology import CoagulationElement
>>> CoagulationElement('blood').injury_time
1.0
Source code in src/skcomplex/systems/biology/physiology/_regulation.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    T: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    C: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    F: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    P: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    ka: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
    k_half: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    d_t: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    k_c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    s_c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    d_c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    s_f: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    k_p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    s_p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
    d_p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    injury: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    injury_time: float = 1.0,
    duration: float = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.T = T
    self.C = C
    self.F = F
    self.P = P
    self.ka = ka
    self.k_half = k_half
    self.d_t = d_t
    self.k_c = k_c
    self.s_c = s_c
    self.d_c = d_c
    self.s_f = s_f
    self.k_p = k_p
    self.s_p = s_p
    self.d_p = d_p
    self.injury = injury
    self.injury_time = injury_time
    self.duration = duration

ContourCellElement(label, visual=None, *, x=0.5, r=2.5, c=0.1, p=0.5)

Bases: Element

Logistic cell of a closed contour that takes up the signalling molecules of the next cell through a power map.

ContourUptakeInteraction iterates x_i <- (1 - c_i) r x_i (1 - x_i) + c_i x_{i+1}^{p_i}, where the logistic map is the intracellular dynamics, c_i the share of the exchange with the environment between the cells, and p_i the affinity of the cell for the uptake. The last cell of the contour takes up from the first.

Parameters:

Name Type Description Default
label str

The label of the group of cells.

 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 of signalling molecules in the cell, between 0 and 1.

 0.5
r Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The logistic parameter, the disposition of the environment to communication.

 2.5
c Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The coupling of the cell to the next one.

 0.1
p Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The affinity of the cell for the uptake of signalling molecules.

 0.5

Examples:

>>> from skcomplex.systems.biology.physiology import ContourCellElement
>>> ContourCellElement('cells').r
2.5
Source code in src/skcomplex/systems/biology/physiology/_exchange.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.5,
    c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.r = r
    self.c = c
    self.p = p

ContourUptakeInteraction(label, *, target)

Bases: CouplingInteraction

Uptake of every cell of a contour from the next cell, x_i <- (1 - c_i) r x_i (1 - x_i) + c_i x_{i+1}^p_i.

It is the map form of coupling to the next member of the group, with the power x_{i+1}^{p_i} of the next cell combined with the logistic map of the cell. The last cell of the contour takes up from the first.

Parameters:

Name Type Description Default
label str

The label of the interaction, unique among the parts of its system.

 required
target Annotated[str, RoleSpec()]

The path of the group of cells.

 required

Examples:

>>> from skcomplex.core.components.elements import GroupElement
>>> from skcomplex.core.components.spaces import DiscreteTime
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.biology.physiology import ContourCellElement, ContourUptakeInteraction
>>> cells = GroupElement('cells', element=ContourCellElement, size=2, values={'x': [0.25, 1.0], 'c': 0.5})
>>> uptake = ContourUptakeInteraction('uptake', target='cells')
>>> system = System(elements=[cells], interactions=[uptake], time=DiscreteTime())
>>> Simulation(system).prepare().run(until=1).trajectory_.values('cells.x')[-1].round(4).tolist()
[0.7344, 0.25]
Source code in src/skcomplex/systems/biology/physiology/_exchange.py 
169
170
171
172
173
174
175
176
177
178
179
def __init__(self, label: str, *, target: Annotated[str, RoleSpec()]) -> None:
    super().__init__(
        label,
        target=target,
        read='x',
        write='x',
        function=self._get_power,
        mode='next',
        form='map',
        combine=self._get_mix,
    )

CytokineLinkElement(label, visual=None, *, kappa=0.0, a=0.0)

Bases: Element

Link between two cells of a layer, a fixed connectivity a and an adaptive cytokine coupling kappa.

Parameters:

Name Type Description Default
label str

The label of the link or of the group on the edges.

 required
visual VisualSpec | None

How the link is drawn, or None for default_visual_.

 None
kappa Annotated[float, StateSpec(domain=closed(-1, 1))]

The initial adaptive coupling, between -1 and 1.

 0.0
a Annotated[float, ParameterSpec(domain=closed(0, 1))]

The fixed connectivity, between 0 and 1.

 0.0

Examples:

>>> from skcomplex.systems.biology.physiology import CytokineLinkElement
>>> CytokineLinkElement('link', a=1.0).a
1.0
Source code in src/skcomplex/systems/biology/physiology/_disease.py 
 97
 98
 99
100
101
102
103
104
105
106
107
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    kappa: Annotated[float, StateSpec(domain=RealSet.closed(-1, 1))] = 0.0,
    a: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.kappa = kappa
    self.a = a

DrugToxicityElement(label, visual=None, *, c=0.0, T=0.0, k=0.5, eta=0.05, dose=1.0, period=7.0, width=1.0, doses=6)

Bases: Element

Drug in the plasma, given in a train of infusions and eliminated, and the toxicity it accumulates.

The plasma concentration follows dc/dt = u(t) - k c, where the infusion u is dose for width time units every period, and the toxicity dT/dt = c - eta T, which the metabolism lowers at the rate eta.

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

Parameters:

Name Type Description Default
label str

The label of the patient.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
c Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial plasma concentration.

 0.0
T Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial toxicity.

 0.0
k Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The elimination rate of the drug.

 0.5
eta Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate at which the metabolism lowers the toxicity, 0 for a cumulative toxicity.

 0.05
dose Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The infusion rate during a dose.

 1.0
period float

The time between the starts of doses.

 7.0
width float

The duration of a dose.

 1.0
doses int

The number of doses.

 6

Examples:

>>> from skcomplex.systems.biology.physiology import DrugToxicityElement
>>> DrugToxicityElement('patient', doses=4).doses
4
Source code in src/skcomplex/systems/biology/physiology/_toxicity.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    c: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    T: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    dose: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    period: float = 7.0,
    width: float = 1.0,
    doses: int = 6,
) -> None:
    super().__init__(label, visual=visual)
    self.c = c
    self.T = T
    self.k = k
    self.eta = eta
    self.dose = dose
    self.period = period
    self.width = width
    self.doses = doses

ExchangingCellElement(label, visual=None, *, x=1.0, alpha=1.0, beta=10.0, g=0.1, xi=0.0)

Bases: Element

Cell that metabolizes a substance and exchanges it with its neighbours through gap junctions and nanotubes.

The concentration follows dx_i/dt = -alpha_i x_i / (beta_i + x_i) + Σ_j [g (x_j - x_i) + xi A_ij x_j], where the first term is Michaelis-Menten metabolism, the one law of the cell, g the exchange coefficient of the gap junctions, and xi the strength of the tunnelling nanotubes along the same links. The sum over the neighbours is a CouplingInteraction of a network of cells.

Parameters:

Name Type Description Default
label str

The label of the field of cells.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
x Annotated[float, StateSpec()]

The initial concentration of the substance.

 1.0
alpha Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The maximal rate of metabolism.

 1.0
beta Annotated[float, ParameterSpec(domain=open(0, oo))]

The half-saturation concentration of metabolism.

 10.0
g Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The exchange coefficient of the gap junctions.

 0.1
xi Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The strength of the nanotubes.

 0.0

Examples:

>>> from skcomplex.systems.biology.physiology import ExchangingCellElement
>>> ExchangingCellElement('cells').g
0.1
Source code in src/skcomplex/systems/biology/physiology/_exchange.py 
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 1.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    g: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    xi: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.alpha = alpha
    self.beta = beta
    self.g = g
    self.xi = xi

FeedbackLoopElement(label, visual=None, *, x=0.0, e=0.0, v=2.0, d=1.0, K=1.0, n=4.0, T=1.0, tau=0.5, steering=0.0, start=0.0, duration=0.0)

Bases: Element

Receptor-effector negative feedback loop that holds a product at a set point an outside signal steers.

The product x is made by the effector at the rate v e and consumed at the rate d x. The receptor reads the product after the transmission delay tau and inhibits the effector with the Hill function 1 / (1 + (x / K)^n), where the sensitivity K (1 + steering) rises while a steering signal acts, from start for duration. The effector follows the receptor with the response time T, de/dt = (1 / (1 + (x(t - tau) / K)^n) - e) / T.

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
x Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial product.

 0.0
e Annotated[float, StateSpec(domain=closed(0, 1))]

The initial effector activity.

 0.0
v Annotated[float, ParameterSpec(domain=open(0, oo))]

The largest rate of synthesis.

 2.0
d Annotated[float, ParameterSpec(domain=open(0, oo))]

The consumption rate of the product.

 1.0
K Annotated[float, ParameterSpec(domain=open(0, oo))]

The sensitivity of the receptor, the product at which it inhibits the effector by half.

 1.0
n Annotated[float, ParameterSpec(domain=open(0, oo))]

The steepness of the receptor.

 4.0
T Annotated[float, ParameterSpec(domain=open(0, oo))]

The response time of the effector.

 1.0
tau Annotated[float, ParameterSpec(domain=open(0, oo))]

The transmission delay between the receptor and the effector.

 0.5
steering Annotated[float, ParameterSpec()]

The relative shift of the sensitivity by the steering signal.

 0.0
start float

The time the steering signal starts.

 0.0
duration float

The duration of the steering signal.

 0.0

Examples:

>>> from skcomplex.systems.biology.physiology import FeedbackLoopElement
>>> FeedbackLoopElement('loop', tau=2.0).tau
2.0
Source code in src/skcomplex/systems/biology/physiology/_regulation.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    e: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.0,
    v: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    d: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 4.0,
    T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    steering: Annotated[float, ParameterSpec()] = 0.0,
    start: float = 0.0,
    duration: float = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.e = e
    self.v = v
    self.d = d
    self.K = K
    self.n = n
    self.T = T
    self.tau = tau
    self.steering = steering
    self.start = start
    self.duration = duration

FitzHughNagumoElement(label, visual=None, *, u=0.0, v=0.0)

Bases: Element

Cell of an excitable medium, with a potential u and a recovery v.

The cell has no laws of its own. FitzHugh-Nagumo kinetics and a diffusion of the potential change it.

Parameters:

Name Type Description Default
label str

The label of the cell.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
u Annotated[float, StateSpec()]

The initial potential.

 0.0
v Annotated[float, StateSpec()]

The initial recovery.

 0.0

Examples:

>>> from skcomplex.systems.biology.physiology import FitzHughNagumoElement
>>> FitzHughNagumoElement('cell', u=1.5).u
1.5
Source code in src/skcomplex/systems/biology/physiology/_physiology.py 
50
51
52
53
54
55
56
57
58
59
60
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec()] = 0.0,
    v: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u
    self.v = v

FlowAdaptationInteraction(label, *, vessels, edges=(), h=(), alpha1=1.0, alpha2=1.0, gamma=2.0)

Bases: Interaction

Adaptation of the conductances of a transport network to the pressure drops of its Kirchhoff flows.

The nodes accumulate no fluid, so the pressures p solve the mass balance h_i = sum_j C_ij (p_i - p_j) with the external flows h, which sum to zero, at every instant. Every conductance follows dC_ij/dt = C_ij (alpha1 |p_i - p_j|^gamma - alpha2): it grows with the power dissipated along its edge and decays at a constant rate. The balance is solved numerically with the pressure of the first node fixed at zero.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
vessels Annotated[str, RoleSpec()]

The path of the group of vessel segments on the edges.

 required
edges Sequence[tuple[int, int]]

The ends of every edge, in the order of the group.

 ()
h Sequence[float]

The external flow of every node, positive for sources and negative for sinks.

 ()
alpha1 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The growth of the conductance with the pressure drop.

 1.0
alpha2 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The decay rate of the conductance.

 1.0
gamma Annotated[float, ParameterSpec(domain=open(0, oo))]

The exponent of the pressure drop.

 2.0

Examples:

>>> from skcomplex.systems.biology.physiology import FlowAdaptationInteraction
>>> FlowAdaptationInteraction('adaptation', vessels='vessels', edges=((0, 1),), h=(1.0, -1.0)).gamma
2.0
Source code in src/skcomplex/systems/biology/physiology/_transport.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
def __init__(
    self,
    label: str,
    *,
    vessels: Annotated[str, RoleSpec()],
    edges: Sequence[tuple[int, int]] = (),
    h: Sequence[float] = (),
    alpha1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    alpha2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
) -> None:
    super().__init__(label)
    self.vessels = vessels
    self.edges = edges
    self.h = h
    self.alpha1 = alpha1
    self.alpha2 = alpha2
    self.gamma = gamma

GlucoseControlElement(label, visual=None, *, G=1.0, I=1.0, N=0.5, R=0.75, k0=0.5, k_i=0.5, k_n=1.0, s_i=2.0, s_n=1.0, d_i=1.0, d_n=1.0, K=1.0, n=4.0, meal=1.0, exercise=0.8, meal_time=2.0, exercise_time=12.0, duration=2.0)

Bases: Element

Blood glucose held near its norm by insulin, which lowers it, and glucagon, which raises it.

The glucose follows dG/dt = R - (k0 + k_i I) G + k_n N + u(t), where u is a meal of height meal from meal_time and an exercise of height -exercise from exercise_time, each lasting duration. High glucose secretes insulin, dI/dt = s_i h(G) - d_i I, and low glucose secretes glucagon, dN/dt = s_n (1 - h(G)) - d_n N, with the Hill switch h(G) = G^n / (G^n + K^n). The production R is set so that the norm G = K is a steady state.

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

Parameters:

Name Type Description Default
label str

The label of the organism.

 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 glucose, relative to the norm.

 1.0
I Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial insulin.

 1.0
N Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial glucagon.

 0.5
R Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The basal production of glucose.

 0.75
k0 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The uptake of glucose without insulin.

 0.5
k_i Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The uptake of glucose per unit of insulin.

 0.5
k_n Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The release of glucose per unit of glucagon, 0 without the glucagon loop.

 1.0
s_i Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The maximal secretion of insulin.

 2.0
s_n Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The maximal secretion of glucagon.

 1.0
d_i Annotated[float, ParameterSpec(domain=open(0, oo))]

The clearance rate of insulin.

 1.0
d_n Annotated[float, ParameterSpec(domain=open(0, oo))]

The clearance rate of glucagon.

 1.0
K Annotated[float, ParameterSpec(domain=open(0, oo))]

The norm of the glucose, where both hormones are half secreted.

 1.0
n Annotated[float, ParameterSpec(domain=open(0, oo))]

The steepness of the secretion switch.

 4.0
meal Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate of glucose intake during the meal.

 1.0
exercise Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate of extra glucose uptake during the exercise.

 0.8
meal_time float

The start of the meal.

 2.0
exercise_time float

The start of the exercise.

 12.0
duration float

The duration of the meal and of the exercise.

 2.0

Examples:

>>> from skcomplex.systems.biology.physiology import GlucoseControlElement
>>> GlucoseControlElement('organism').meal_time
2.0
Source code in src/skcomplex/systems/biology/physiology/_regulation.py 
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    G: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    I: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    N: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    R: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.75,
    k0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    k_i: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    k_n: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    s_i: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    s_n: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    d_i: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    d_n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 4.0,
    meal: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    exercise: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.8,
    meal_time: float = 2.0,
    exercise_time: float = 12.0,
    duration: float = 2.0,
) -> None:
    super().__init__(label, visual=visual)
    self.G = G
    self.I = I
    self.N = N
    self.R = R
    self.k0 = k0
    self.k_i = k_i
    self.k_n = k_n
    self.s_i = s_i
    self.s_n = s_n
    self.d_i = d_i
    self.d_n = d_n
    self.K = K
    self.n = n
    self.meal = meal
    self.exercise = exercise
    self.meal_time = meal_time
    self.exercise_time = exercise_time
    self.duration = duration

LayerCouplingInteraction(label, *, cells, sigma=1.0, layer=None)

Bases: NeighborInteraction

Coupling of every cell to its partner in another layer, dphi_i/dt += -sigma sin(phi_i - phi_j).

It runs on the edges of the layer of the network that joins the partners, both ways.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
cells Annotated[str, RoleSpec()]

The path of the group of cells on the nodes.

 required
sigma Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The strength of the coupling between the layers.

 1.0
layer str | None

The layer of the network whose edges join the partners.

 None

Examples:

>>> from skcomplex.systems.biology.physiology import LayerCouplingInteraction
>>> LayerCouplingInteraction('interlayer', cells='cells', sigma=0.5, layer='partners').layer
'partners'
Source code in src/skcomplex/systems/biology/physiology/_disease.py 
207
208
209
210
211
212
213
214
215
216
217
def __init__(
    self,
    label: str,
    *,
    cells: Annotated[str, RoleSpec()],
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    layer: str | None = None,
) -> None:
    super().__init__(label, layer=layer)
    self.cells = cells
    self.sigma = sigma

MackeyGlassElement(label, visual=None, *, x=0.5, beta=0.2, gamma=0.1, n=10.0, tau=17.0)

Bases: Element

Concentration x of circulating blood cells, made at a delayed Hill rate and lost at a constant rate.

The concentration follows x' = beta x(t - tau) / (1 + x(t - tau)^n) - gamma x.

Parameters:

Name Type Description Default
label str

The label of the element.

 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, also its history.

 0.5
beta Annotated[float, ParameterSpec()]

The production rate.

 0.2
gamma Annotated[float, ParameterSpec()]

The loss rate.

 0.1
n Annotated[float, ParameterSpec()]

The Hill coefficient of the production.

 10.0
tau Annotated[float, ParameterSpec(domain=open(0, oo))]

The delay of the production.

 17.0

Examples:

>>> from skcomplex.systems.biology.physiology import MackeyGlassElement
>>> MackeyGlassElement('blood', tau=6.0).tau
6.0
Source code in src/skcomplex/systems/biology/physiology/_physiology.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    beta: Annotated[float, ParameterSpec()] = 0.2,
    gamma: Annotated[float, ParameterSpec()] = 0.1,
    n: Annotated[float, ParameterSpec()] = 10.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 17.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.beta = beta
    self.gamma = gamma
    self.n = n
    self.tau = tau

ResistiveFlowInteraction(label, *, source, target, R=1.0, valve=False)

Bases: TransferInteraction

Flow of blood from one chamber to another through a resistance, Q = (P_source - P_target) / R.

A valve lets blood pass forward only, so the flow is zero when the target pressure is higher.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
source Annotated[str, RoleSpec()]

The path of the chamber the blood leaves.

 required
target Annotated[str, RoleSpec()]

The path of the chamber the blood enters.

 required
R Annotated[float, ParameterSpec(domain=open(0, oo))]

The resistance in mmHg s/ml.

 1.0
valve bool

Whether a valve prevents backward flow.

 False

Examples:

>>> from skcomplex.systems.biology.physiology import ResistiveFlowInteraction
>>> ResistiveFlowInteraction('aortic_valve', source='ventricle', target='arteries', valve=True).valve
True
Source code in src/skcomplex/systems/biology/physiology/_circulation.py 
108
109
110
111
112
113
114
115
116
117
118
119
def __init__(
    self,
    label: str,
    *,
    source: Annotated[str, RoleSpec()],
    target: Annotated[str, RoleSpec()],
    R: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    valve: bool = False,
) -> None:
    super().__init__(label, source=source, target=target, state='V', amount=self._get_flow)
    self.R = R
    self.valve = valve

TissueSiteElement(label, visual=None, *, cell=0.0)

Bases: Element

Site of a lattice tissue, which holds the label of the cell that covers it.

The sites of one label form a cell of the cellular Potts model, and CellAdhesionInteraction moves the boundaries between cells by copying the label of a neighbouring site.

Parameters:

Name Type Description Default
label str

The label of the site or the lattice group.

 required
visual VisualSpec | None

How the site is drawn, or None for default_visual_.

 None
cell Annotated[float, StateSpec(domain=NN)]

The initial label of the cell that covers the site, from 0.

 0.0

Attributes:

Name Type Description
default_visual_

The visual of the site when it sets none.

Examples:

>>> from skcomplex.systems.biology.physiology import TissueSiteElement
>>> TissueSiteElement('site', cell=3).cell
3
Source code in src/skcomplex/systems/biology/physiology/_cellular_potts.py 
69
70
71
72
73
74
75
76
77
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    cell: Annotated[float, StateSpec(domain=NN)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.cell = cell

TwoCompartmentElement(label, visual=None, *, Uc=1.0, Up=0.0)

Bases: Element

Drug concentrations of a patient in a central compartment Uc and a peripheral compartment Up.

The compartments have no laws of their own. Reactions exchange the drug between them and eliminate it.

Parameters:

Name Type Description Default
label str

The label of the patient.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
Uc Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial central concentration.

 1.0
Up Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial peripheral concentration.

 0.0

Examples:

>>> from skcomplex.systems.biology.physiology import TwoCompartmentElement
>>> TwoCompartmentElement('patient').Up
0.0
Source code in src/skcomplex/systems/biology/physiology/_pharmacokinetics.py 
42
43
44
45
46
47
48
49
50
51
52
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    Uc: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    Up: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.Uc = Uc
    self.Up = Up

VesselSegmentElement(label, visual=None, *, C=1.0)

Bases: Element

Vessel segment on an edge of a transport network, whose conductance adapts to the flow it carries.

Parameters:

Name Type Description Default
label str

The label of the segment or of the group on the edges.

 required
visual VisualSpec | None

How the segment is drawn, or None for default_visual_.

 None
C Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial conductance.

 1.0

Examples:

>>> from skcomplex.systems.biology.physiology import VesselSegmentElement
>>> VesselSegmentElement('vessel', C=2.0).C
2.0
Source code in src/skcomplex/systems/biology/physiology/_transport.py 
44
45
46
47
48
49
50
51
52
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    C: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.C = C

adaptive_flow_network(side=6, *, alpha1=1.0, alpha2=1.0, gamma=2.0, inflow=None, seed=0)

Return a square grid of vessels fed at one corner and drained at every other node, whose conductances adapt.

The pressures balance the flows at every node, h_i = sum_j C_ij (p_i - p_j), and every conductance grows with the dissipated power and decays, dC_ij/dt = C_ij (alpha1 |p_i - p_j|^gamma - alpha2). Under a steady supply the network settles where every surviving vessel carries the same pressure drop (alpha2 / alpha1)^(1 / gamma), so the conductances grow with the flow each vessel must carry: the vessels near the source thicken into trunks and the distal ones stay thin.

Parameters:

Name Type Description Default
side int

The number of nodes along each side of the grid.

 6
alpha1 float

The growth of the conductance with the pressure drop.

 1.0
alpha2 float

The decay rate of the conductance.

 1.0
gamma float

The exponent of the pressure drop.

 2.0
inflow float | None

The flow into the source at the first corner, which the other nodes drain equally, or None for the number of edges.

 None
seed int | None

The seed of the initial conductances.

 0

Returns:

Name Type Description
system System

The system, with the conductances vessels.C on the edges of the grid.

Examples:

>>> from skcomplex.systems.biology.physiology import adaptive_flow_network
>>> adaptive_flow_network(3).elements[0].prepare().size_
12
Source code in src/skcomplex/systems/biology/physiology/_transport.py 
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
@register_system(
    field='biology/physiology',
    tags=('adaptive network', 'vasculature', 'kirchhoff', 'remodelling', 'ode'),
    references=('Hu and Cai 2013', 'Martens and Klemm 2017', 'Berner et al. 2023'),
    defaults={'until': 20.0, 'dt': 0.1},
)
def adaptive_flow_network(
    side: int = 6,
    *,
    alpha1: float = 1.0,
    alpha2: float = 1.0,
    gamma: float = 2.0,
    inflow: float | None = None,
    seed: int | None = 0,
) -> System:
    """Return a square grid of vessels fed at one corner and drained at every other node, whose conductances adapt.

    The pressures balance the flows at every node, `h_i = sum_j C_ij (p_i - p_j)`, and every conductance grows with the
    dissipated power and decays, `dC_ij/dt = C_ij (alpha1 |p_i - p_j|^gamma - alpha2)`. Under a steady supply the
    network settles where every surviving vessel carries the same pressure drop `(alpha2 / alpha1)^(1 / gamma)`, so the
    conductances grow with the flow each vessel must carry: the vessels near the source thicken into trunks and the
    distal ones stay thin.

    Args:
        side:
            The number of nodes along each side of the grid.

        alpha1:
            The growth of the conductance with the pressure drop.

        alpha2:
            The decay rate of the conductance.

        gamma:
            The exponent of the pressure drop.

        inflow:
            The flow into the source at the first corner, which the other nodes drain equally, or `None` for the number
            of edges.

        seed:
            The seed of the initial conductances.

    Returns:
        system:
            The system, with the conductances `vessels.C` on the edges of the grid.

    Examples:
        >>> from skcomplex.systems.biology.physiology import adaptive_flow_network
        >>> adaptive_flow_network(3).elements[0].prepare().size_
        12
    """
    graph = nx.convert_node_labels_to_integers(nx.grid_2d_graph(side, side))
    space = NetworkSpace(graph).prepare()
    sources, targets = space.edges_[None]
    n_nodes, n_edges = graph.number_of_nodes(), len(sources)
    supply = float(n_edges) if inflow is None else inflow
    h = np.full(n_nodes, -supply / (n_nodes - 1))
    h[0] = supply
    rng = np.random.default_rng(seed)
    vessels = GroupElement(
        'vessels',
        element=VesselSegmentElement,
        space=space,
        on='edges',
        values={'C': rng.uniform(0.5, 1.5, n_edges)},
    )
    adaptation = FlowAdaptationInteraction(
        'adaptation',
        vessels='vessels',
        edges=tuple(zip(sources.tolist(), targets.tolist(), strict=True)),
        h=tuple(h.tolist()),
        alpha1=alpha1,
        alpha2=alpha2,
        gamma=gamma,
    )
    return System(elements=[vessels], interactions=[adaptation])

adaptive_response_costs(strategy='adaptive', *, height=1.0, width=2.0, period=10.0)

Return a system that meets recurring demands by stimulus-response, adaptive response, or anticipation.

Episodes of demand push a critical parameter away from its norm, and the system counters them, after Mobus. A stimulus-response keeps its baseline capacity, so large episodes overshoot the damage threshold and cost repairs. An adaptive response builds capacity under repeated demand and loses it when demand stops, at the price of its upkeep. An anticipatory response also predicts every episode and avoids half of it. Mobus's objective is the sum of the response, restore, and repair costs over the episodes, with the preparedness cost amortized over them: adaptive response lowers it below stimulus-response by avoiding repairs, and anticipation lowers it further by cutting the response it needs.

Parameters:

Name Type Description Default
strategy str

'stimulus-response', 'adaptive', or 'anticipatory'.

 'adaptive'
height float

The height of the episodes of demand.

 1.0
width float

The duration of the episodes.

 2.0
period float

The time between the starts of episodes.

 10.0

Returns:

Name Type Description
system System

The system, with the excursion organism.e, the capacity organism.c, the reserves organism.S, and the costs organism.response_cost, organism.restore_cost, organism.repair_cost, and organism.prepare_cost.

Raises:

Type Description
ValidationError

The strategy is unknown.

Examples:

>>> from skcomplex.systems.biology.physiology import adaptive_response_costs
>>> adaptive_response_costs('stimulus-response').elements[0].eta
0.0
Source code in src/skcomplex/systems/biology/physiology/_adaptivity.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
@register_system(
    field='biology/physiology',
    tags=('adaptation', 'anticipation', 'homeostasis', 'costs', 'stocks and flows', 'ode'),
    references=('Mobus 2022', 'Mobus 1999'),
    defaults={'until': 200.0, 'dt': 0.1},
)
def adaptive_response_costs(
    strategy: str = 'adaptive',
    *,
    height: float = 1.0,
    width: float = 2.0,
    period: float = 10.0,
) -> System:
    """Return a system that meets recurring demands by stimulus-response, adaptive response, or anticipation.

    Episodes of demand push a critical parameter away from its norm, and the system counters them, after Mobus. A
    stimulus-response keeps its baseline capacity, so large episodes overshoot the damage threshold and cost repairs.
    An adaptive response builds capacity under repeated demand and loses it when demand stops, at the price of its
    upkeep. An anticipatory response also predicts every episode and avoids half of it. Mobus's objective is the sum
    of the response, restore, and repair costs over the episodes, with the preparedness cost amortized over them:
    adaptive response lowers it below stimulus-response by avoiding repairs, and anticipation lowers it further by
    cutting the response it needs.

    Args:
        strategy:
            `'stimulus-response'`, `'adaptive'`, or `'anticipatory'`.

        height:
            The height of the episodes of demand.

        width:
            The duration of the episodes.

        period:
            The time between the starts of episodes.

    Returns:
        system:
            The system, with the excursion `organism.e`, the capacity `organism.c`, the reserves `organism.S`, and the
            costs `organism.response_cost`, `organism.restore_cost`, `organism.repair_cost`, and
            `organism.prepare_cost`.

    Raises:
        ValidationError:
            The strategy is unknown.

    Examples:
        >>> from skcomplex.systems.biology.physiology import adaptive_response_costs
        >>> adaptive_response_costs('stimulus-response').elements[0].eta
        0.0
    """
    if strategy not in STRATEGIES:
        error_msg = f'Unknown strategy {strategy!r}. Use one of {", ".join(STRATEGIES)}.'
        raise ValidationError(error_msg)
    eta, anticipation = STRATEGIES[strategy]
    organism = AdaptiveResponderElement(
        'organism',
        eta=eta,
        anticipation=anticipation,
        height=height,
        width=width,
        period=period,
    )
    interactions = [
        RateInteraction(
            'stimulus',
            target='organism',
            state='e',
            rate=lambda s: (1 - s.target.anticipation)
            * pulse_train(s.t, organism.period, organism.width, organism.height, start=organism.start),
        ),
        ReactionInteraction(
            'effort',
            target='organism',
            reactants={'e': 1, 'S': 1},
            products={'response_cost': 1},
            propensity=lambda s: s.target.c * s.target.e,
        ),
        ReactionInteraction(
            'effort_decay',
            target='organism',
            reactants={'e': 1},
            propensity=lambda s: s.target.k0 * s.target.e,
        ),
        RateInteraction('capacity_growth', target='organism', state='c', rate=lambda s: s.target.eta * s.target.e**2),
        RelaxationInteraction(
            'capacity_relaxation',
            target='organism',
            state='c',
            goal=lambda s: s.target.c0,
            time=lambda s: s.target.tau,
        ),
        RateInteraction('restoration', target='organism', state='S', rate=lambda s: s.target.kappa * (1 - s.target.S)),
        RateInteraction(
            'restoration_cost',
            target='organism',
            state='restore_cost',
            rate=lambda s: s.target.price * s.target.kappa * (1 - s.target.S),
        ),
        RateInteraction(
            'repair_cost',
            target='organism',
            state='repair_cost',
            rate=lambda s: s.target.repair * where(s.target.e > s.target.theta, s.target.e - s.target.theta, 0),
        ),
        RateInteraction(
            'preparation_cost',
            target='organism',
            state='prepare_cost',
            rate=lambda s: s.target.upkeep * (s.target.c - s.target.c0 + s.target.anticipation),
        ),
    ]
    return System(elements=[organism], interactions=interactions)

aliev_panfilov(shape=(100, 100), *, D=1.0, k=8.0, a=0.15, epsilon0=0.002, mu1=0.2, mu2=0.3)

Return the Aliev-Panfilov system of cardiac tissue, started from a broken wave that curls into a spiral.

A spiral wave in heart muscle is a re-entrant arrhythmia, a tachycardia, and its breakup into many waves systems fibrillation. The upper half of the tissue starts excited and its left half refractory, so the free end of the wave rotates.

Parameters:

Name Type Description Default
shape Sequence[int]

The number of rows and columns of the lattice with zero-flux edges.

 (100, 100)
D float

The diffusion of the potential.

 1.0
k float

The strength of the excitation.

 8.0
a float

The excitation threshold.

 0.15
epsilon0 float

The recovery rate far from the excited state.

 0.002
mu1 float

The growth of the recovery rate with the recovery.

 0.2
mu2 float

The offset of the recovery rate.

 0.3

Returns:

Name Type Description
system System

The system, with the potential tissue.u and the recovery tissue.v.

Examples:

>>> from skcomplex.systems.biology.physiology import aliev_panfilov
>>> aliev_panfilov((20, 20)).elements[0].values['k']
8.0
Source code in src/skcomplex/systems/biology/physiology/_cardiac.py 
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
@register_system(
    field='biology/physiology',
    tags=('field', 'excitable media', 'spirals', 'cardiac', 'pde'),
    references=('Aliev and Panfilov 1996', 'Datseris and Parlitz 2022'),
    defaults={'until': 200.0, 'dt': 2.0, 'solver': RK4Solver(dt=0.05)},
)
def aliev_panfilov(
    shape: Sequence[int] = (100, 100),
    *,
    D: float = 1.0,
    k: float = 8.0,
    a: float = 0.15,
    epsilon0: float = 0.002,
    mu1: float = 0.2,
    mu2: float = 0.3,
) -> System:
    """Return the Aliev-Panfilov system of cardiac tissue, started from a broken wave that curls into a spiral.

    A spiral wave in heart muscle is a re-entrant arrhythmia, a tachycardia, and its breakup into many waves systems
    fibrillation. The upper half of the tissue starts excited and its left half refractory, so the free end of the
    wave rotates.

    Args:
        shape:
            The number of rows and columns of the lattice with zero-flux edges.

        D:
            The diffusion of the potential.

        k:
            The strength of the excitation.

        a:
            The excitation threshold.

        epsilon0:
            The recovery rate far from the excited state.

        mu1:
            The growth of the recovery rate with the recovery.

        mu2:
            The offset of the recovery rate.

    Returns:
        system:
            The system, with the potential `tissue.u` and the recovery `tissue.v`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import aliev_panfilov
        >>> aliev_panfilov((20, 20)).elements[0].values['k']
        8.0
    """
    rows, columns = shape
    u, v = np.zeros(shape), np.zeros(shape)
    u[: rows // 2, :] = 1.0
    v[:, : columns // 2] = 2.0
    values = {'u': u, 'v': v, 'k': k, 'a': a, 'epsilon0': epsilon0, 'mu1': mu1, 'mu2': mu2}
    space = LatticeSpace((rows, columns), boundary='reflecting')
    tissue = FieldElement('tissue', element=AlievPanfilovElement, space=space, values=values)
    diffusion = DiffusionInteraction('diffusion', target='tissue', state='u', coefficient=D)
    return System(elements=[tissue], interactions=[diffusion])

barkley(shape=(64, 64), *, spacing=3.0, d=1.0, a=0.6, b=0.05, epsilon=0.02, gamma=1.0)

Return Barkley's excitable medium, started from orthogonal stripes whose crossing curls into a spiral wave.

The excitation u and the recovery w start as orthogonal stripes of zeros and ones, so the excited stripe runs into recovering tissue on one side only and its free end rotates. With (a, b, d, epsilon) = (0.6, 0.05, 1, 0.02) the spiral is stable for gamma = 1. With (0.7, 0.05, 1, 0.08) and gamma = 3 it breaks up into spatiotemporal chaos. The source runs a domain of 300 by 300, and the default lattice covers 192 by 192 with a spacing of 3, wide enough for the stable spiral.

Parameters:

Name Type Description Default
shape Sequence[int]

The number of rows and columns of the lattice with zero-flux edges.

 (64, 64)
spacing float

The distance between neighbouring sites.

 3.0
d float

The diffusion of the excitation.

 1.0
a float

The inverse slope of the threshold.

 0.6
b float

The offset of the threshold.

 0.05
epsilon float

The ratio of the time scales of the recovery and the excitation.

 0.02
gamma float

The exponent of the excitation in the recovery.

 1.0

Returns:

Name Type Description
system System

The system, with the excitation medium.u and the recovery medium.w.

Examples:

>>> from skcomplex.systems.biology.physiology import barkley
>>> barkley((20, 20), gamma=3.0).elements[0].values['gamma']
3.0
Source code in src/skcomplex/systems/biology/physiology/_cardiac.py 
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
@register_system(
    field='biology/physiology',
    tags=('field', 'excitable media', 'spirals', 'spatiotemporal chaos', 'pde'),
    references=('Barkley 1991', 'Datseris and Parlitz 2022'),
    defaults={'until': 400.0, 'dt': 4.0, 'solver': RK4Solver(dt=0.2)},
)
def barkley(
    shape: Sequence[int] = (64, 64),
    *,
    spacing: float = 3.0,
    d: float = 1.0,
    a: float = 0.6,
    b: float = 0.05,
    epsilon: float = 0.02,
    gamma: float = 1.0,
) -> System:
    """Return Barkley's excitable medium, started from orthogonal stripes whose crossing curls into a spiral wave.

    The excitation `u` and the recovery `w` start as orthogonal stripes of zeros and ones, so the excited stripe runs
    into recovering tissue on one side only and its free end rotates. With `(a, b, d, epsilon) = (0.6, 0.05, 1, 0.02)`
    the spiral is stable for `gamma = 1`. With `(0.7, 0.05, 1, 0.08)` and `gamma = 3` it breaks up into spatiotemporal
    chaos. The source runs a domain of 300 by 300, and the default lattice covers 192 by 192 with a spacing of 3, wide
    enough for the stable spiral.

    Args:
        shape:
            The number of rows and columns of the lattice with zero-flux edges.

        spacing:
            The distance between neighbouring sites.

        d:
            The diffusion of the excitation.

        a:
            The inverse slope of the threshold.

        b:
            The offset of the threshold.

        epsilon:
            The ratio of the time scales of the recovery and the excitation.

        gamma:
            The exponent of the excitation in the recovery.

    Returns:
        system:
            The system, with the excitation `medium.u` and the recovery `medium.w`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import barkley
        >>> barkley((20, 20), gamma=3.0).elements[0].values['gamma']
        3.0
    """
    rows, columns = shape
    u, w = np.zeros(shape), np.zeros(shape)
    u[:, : columns // 2] = 1.0
    w[: rows // 2, :] = 1.0
    values = {'u': u, 'w': w, 'a': a, 'b': b, 'epsilon': epsilon, 'gamma': gamma}
    space = LatticeSpace((rows, columns), boundary='reflecting', spacing=spacing)
    medium = FieldElement('medium', element=BarkleyElement, space=space, values=values)
    diffusion = DiffusionInteraction('diffusion', target='medium', state='u', coefficient=d)
    return System(elements=[medium], interactions=[diffusion])

cardiovascular(*, heart_rate=75.0, E_max=2.5, R_systemic=1.0, C_arterial=1.3, C_venous=50.0, volume=1500.0)

Return a closed lumped system of the systemic circulation driven by a left ventricle of varying elastance.

Blood leaves the ventricle through the aortic valve into the compliant arteries, flows through the systemic resistance into the veins, and returns through the mitral valve. The ventricle stiffens in systole, so the system beats, the valves open and close, and the total volume of blood stays constant. With the defaults the arterial pressure swings between about 80 and 120 mmHg.

Parameters:

Name Type Description Default
heart_rate float

The number of beats per minute.

 75.0
E_max float

The largest elastance of the ventricle in mmHg/ml.

 2.5
R_systemic float

The systemic resistance in mmHg s/ml.

 1.0
C_arterial float

The arterial compliance in ml/mmHg.

 1.3
C_venous float

The venous compliance in ml/mmHg.

 50.0
volume float

The total volume of blood in ml, split at the start among the ventricle, arteries, and veins.

 1500.0

Returns:

Name Type Description
system System

The system, with the volumes ventricle.V, arteries.V, and veins.V.

Examples:

>>> from skcomplex.systems.biology.physiology import cardiovascular
>>> sum(chamber.V for chamber in cardiovascular().elements)
1500.0
Source code in src/skcomplex/systems/biology/physiology/_circulation.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
@register_system(
    field='biology/physiology',
    tags=('cardiovascular', 'circulation', 'lumped system', 'hybrid'),
    references=('Suga and Sagawa 1974', 'Trzaska 2023'),
    defaults={'until': 8.0, 'dt': 0.02},
)
def cardiovascular(
    *,
    heart_rate: float = 75.0,
    E_max: float = 2.5,
    R_systemic: float = 1.0,
    C_arterial: float = 1.3,
    C_venous: float = 50.0,
    volume: float = 1500.0,
) -> System:
    """Return a closed lumped system of the systemic circulation driven by a left ventricle of varying elastance.

    Blood leaves the ventricle through the aortic valve into the compliant arteries, flows through the systemic
    resistance into the veins, and returns through the mitral valve. The ventricle stiffens in systole, so the system
    beats, the valves open and close, and the total volume of blood stays constant. With the defaults the arterial
    pressure swings between about 80 and 120 mmHg.

    Args:
        heart_rate:
            The number of beats per minute.

        E_max:
            The largest elastance of the ventricle in mmHg/ml.

        R_systemic:
            The systemic resistance in mmHg s/ml.

        C_arterial:
            The arterial compliance in ml/mmHg.

        C_venous:
            The venous compliance in ml/mmHg.

        volume:
            The total volume of blood in ml, split at the start among the ventricle, arteries, and veins.

    Returns:
        system:
            The system, with the volumes `ventricle.V`, `arteries.V`, and `veins.V`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import cardiovascular
        >>> sum(chamber.V for chamber in cardiovascular().elements)
        1500.0
    """
    period = 60.0 / heart_rate
    systole = 0.375 * period
    ventricle = ChamberElement(
        'ventricle',
        V=0.1 * volume,
        V0=10.0,
        E_min=0.08,
        E_max=E_max,
        period=period,
        systole=systole,
    )
    arteries = ChamberElement(
        'arteries',
        V=0.3 * volume,
        V0=250.0,
        E_min=1 / C_arterial,
        E_max=1 / C_arterial,
        period=period,
        systole=systole,
    )
    veins = ChamberElement(
        'veins',
        V=0.6 * volume,
        V0=500.0,
        E_min=1 / C_venous,
        E_max=1 / C_venous,
        period=period,
        systole=systole,
    )
    flows = [
        ResistiveFlowInteraction('aortic_valve', source='ventricle', target='arteries', R=0.02, valve=True),
        ResistiveFlowInteraction('systemic', source='arteries', target='veins', R=R_systemic),
        ResistiveFlowInteraction('mitral_valve', source='veins', target='ventricle', R=0.02, valve=True),
    ]
    return System(elements=[ventricle, arteries, veins], interactions=flows)

cellular_potts(shape=(40, 40), *, lattice='square', cell_shape=(5, 8), adhesion=DEFAULT_ADHESION, temperature=2.0, area_strength=1.0, perimeter_strength=0.0, target_perimeter=0.0, hits=1, random_state=0)

Return the cellular Potts model of a confluent tissue whose cells of several kinds sort by differential adhesion.

Every site of a periodic lattice belongs to one cell, and the cells start as rectangles of cell_shape sites, whose area is the target area, with kinds in equal shares at random. Label copies sample the energy sum_<i,j> w_ij J(t_i, t_j) [s_i != s_j] + lambda sum_s (a_s - A)^2 + mu sum_s (p_s - P)^2 of Nemati and de Graaf (2024), Eq. 3, extended to the adhesion J between kinds of Graner and Glazier (1992). The square lattice reads the Moore neighbourhood and the triangular lattice the six nearest neighbours, both with unit weights, and the disordered lattice is the periodic Voronoi tessellation of one jittered point per unit area, whose sites have the area of their Voronoi cell and touch along sides of weight w_ij, the length of the side over the mean length. A heterotypic adhesion above the mean of the homotypic ones raises the share of contacts between cells of one kind.

Parameters:

Name Type Description Default
shape Sequence[int]

The number of rows and columns of the lattice, or of points of the disordered lattice, multiples of the cell shape.

 (40, 40)
lattice str

'square', 'triangular', or 'disordered'.

 'square'
cell_shape Sequence[int]

The rows and columns of the initial rectangle of every cell, whose area is the target area.

 (5, 8)
adhesion ArrayLike

The adhesion energy per unit contact between kinds, a symmetric matrix with one row per kind.

 DEFAULT_ADHESION
temperature float

The temperature of the label copies.

 2.0
area_strength float

The strength lambda of the area constraint.

 1.0
perimeter_strength float

The strength mu of the perimeter constraint.

 0.0
target_perimeter float

The target perimeter P, the weighted number of contacts of a cell with other cells.

 0.0
hits int

The number of copies proposed per site and sweep.

 1
random_state int | Generator | None

The seed or generator of the kinds of the cells and of the disordered lattice.

 0

Returns:

Name Type Description
system System

The system, with the states tissue.cell.

Raises:

Type Description
ValidationError

The lattice is unknown, or the shape is not a multiple of the cell shape.

Examples:

>>> from skcomplex.systems.biology.physiology import cellular_potts
>>> system = cellular_potts((10, 16))
>>> system.interactions[0].kinds
(0, 0, 1, 1)
>>> cellular_potts((10, 16), lattice='disordered').elements[0].space.graph.number_of_nodes()
160
Source code in src/skcomplex/systems/biology/physiology/_cellular_potts.py 
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
@register_system(
    field='biology/physiology',
    tags=('cellular potts', 'tissue', 'cell sorting', 'differential adhesion', 'lattice', 'monte carlo'),
    references=('Nemati and de Graaf 2024, arXiv 2404.09055', 'Graner and Glazier 1992'),
    defaults={'until': 2000, 'dt': 1},
    observables=_cellular_potts_observables,
)
def cellular_potts(
    shape: Sequence[int] = (40, 40),
    *,
    lattice: str = 'square',
    cell_shape: Sequence[int] = (5, 8),
    adhesion: npt.ArrayLike = DEFAULT_ADHESION,
    temperature: float = 2.0,
    area_strength: float = 1.0,
    perimeter_strength: float = 0.0,
    target_perimeter: float = 0.0,
    hits: int = 1,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the cellular Potts model of a confluent tissue whose cells of several kinds sort by differential adhesion.

    Every site of a periodic lattice belongs to one cell, and the cells start as rectangles of `cell_shape` sites, whose
    area is the target area, with kinds in equal shares at random. Label copies sample the energy
    `sum_<i,j> w_ij J(t_i, t_j) [s_i != s_j] + lambda sum_s (a_s - A)^2 + mu sum_s (p_s - P)^2` of Nemati and de Graaf
    (2024), Eq. 3, extended to the adhesion `J` between kinds of Graner and Glazier (1992). The square lattice reads the
    Moore neighbourhood and the triangular lattice the six nearest neighbours, both with unit weights, and the
    disordered lattice is the periodic Voronoi tessellation of one jittered point per unit area, whose sites have the
    area of their Voronoi cell and touch along sides of weight `w_ij`, the length of the side over the mean length. A
    heterotypic adhesion above the mean of the homotypic ones raises the share of contacts between cells of one kind.

    Args:
        shape:
            The number of rows and columns of the lattice, or of points of the disordered lattice, multiples of the
            cell shape.

        lattice:
            `'square'`, `'triangular'`, or `'disordered'`.

        cell_shape:
            The rows and columns of the initial rectangle of every cell, whose area is the target area.

        adhesion:
            The adhesion energy per unit contact between kinds, a symmetric matrix with one row per kind.

        temperature:
            The temperature of the label copies.

        area_strength:
            The strength `lambda` of the area constraint.

        perimeter_strength:
            The strength `mu` of the perimeter constraint.

        target_perimeter:
            The target perimeter `P`, the weighted number of contacts of a cell with other cells.

        hits:
            The number of copies proposed per site and sweep.

        random_state:
            The seed or generator of the kinds of the cells and of the disordered lattice.

    Returns:
        system:
            The system, with the states `tissue.cell`.

    Raises:
        ValidationError:
            The lattice is unknown, or the shape is not a multiple of the cell shape.

    Examples:
        >>> from skcomplex.systems.biology.physiology import cellular_potts
        >>> system = cellular_potts((10, 16))
        >>> system.interactions[0].kinds
        (0, 0, 1, 1)
        >>> cellular_potts((10, 16), lattice='disordered').elements[0].space.graph.number_of_nodes()
        160
    """
    rows, columns = (int(size) for size in shape)
    height, width = (int(size) for size in cell_shape)
    if lattice not in LATTICES:
        error_msg = f'The lattice of the cellular Potts model is one of {", ".join(LATTICES)}, but got {lattice!r}.'
        raise ValidationError(error_msg)
    if rows % height or columns % width:
        error_msg = f'The shape {(rows, columns)} of the tissue is not a multiple of the cell shape {(height, width)}.'
        raise ValidationError(error_msg)
    rng = np.random.default_rng(random_state)
    count = (rows // height) * (columns // width)
    kinds = np.asarray(adhesion, dtype=float).shape[0]
    kind_of = tuple(int(kind) for kind in rng.permutation(np.arange(count) % kinds))
    if lattice == 'disordered':
        graph = _disordered_lattice((rows, columns), rng)
        x, y = np.asarray([point for _, point in graph.nodes(data='pos')]).T
        space: LatticeSpace | NetworkSpace = NetworkSpace(graph)
    else:
        y, x = np.indices((rows, columns)).reshape(2, -1).astype(float)
        space = LatticeSpace((rows, columns), geometry=lattice)
    cells = (y // height).astype(int) * (columns // width) + (x // width).astype(int)
    tissue = LatticeElement('tissue', element=TissueSiteElement, space=space, values={'cell': cells.astype(float)})
    sampling = CellAdhesionInteraction(
        'adhesion',
        target='tissue',
        kinds=kind_of,
        adhesion=adhesion,
        temperature=temperature,
        area_strength=area_strength,
        target_area=float(height * width),
        perimeter_strength=perimeter_strength,
        target_perimeter=target_perimeter,
        hits=hits,
    )
    return System(elements=[tissue], interactions=[sampling], time=DiscreteTime())

closed_contour_cell_exchange(n_cells=8, *, r=2.5, c=0.1, p=0.5, random_state=0)

Return the exchange of signalling molecules between logistic cells arranged in a closed contour.

Every cell runs the logistic map on its concentration and takes up a share c_i of its next state from the next cell through the power map x_{i+1}^{p_i}, x_i <- (1 - c_i) r x_i (1 - x_i) + c_i x_{i+1}^{p_i}. Mihailović et al. show that an equilibrium is asymptotically stable for every coupling when every x_i lies in ((r - 1) / 2r, (r + 1) / 2r) and above p_{i-1}^(1 / (1 - p_{i-1})). In this domain the cells synchronize to a common concentration, and a chaotic environment, r near 4, desynchronizes them.

Parameters:

Name Type Description Default
n_cells int

The number of cells of the contour.

 8
r float

The logistic parameter of every cell.

 2.5
c float | Sequence[float]

The coupling of every cell to the next, one value or one per cell.

 0.1
p float | Sequence[float]

The affinity of every cell for the uptake, one value or one per cell.

 0.5
random_state int | Generator | None

The seed or generator of the initial concentrations, uniform between 0.2 and 0.8.

 0

Returns:

Name Type Description
system System

The system, with the concentrations cells.x, cell i + 1 next to cell i.

Examples:

>>> from skcomplex.systems.biology.physiology import closed_contour_cell_exchange
>>> closed_contour_cell_exchange(5).elements[0].prepare().size_
5
Source code in src/skcomplex/systems/biology/physiology/_exchange.py 
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
@register_system(
    field='biology/physiology',
    tags=('map', 'logistic', 'coupled maps', 'ring', 'synchronization', 'cell communication'),
    references=('Mihailović et al. 2023',),
    defaults={'until': 200, 'dt': 1},
)
def closed_contour_cell_exchange(
    n_cells: int = 8,
    *,
    r: float = 2.5,
    c: float | Sequence[float] = 0.1,
    p: float | Sequence[float] = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the exchange of signalling molecules between logistic cells arranged in a closed contour.

    Every cell runs the logistic map on its concentration and takes up a share `c_i` of its next state from the next
    cell through the power map `x_{i+1}^{p_i}`, `x_i <- (1 - c_i) r x_i (1 - x_i) + c_i x_{i+1}^{p_i}`. Mihailović et
    al. show that an equilibrium is asymptotically stable for every coupling when every `x_i` lies in
    `((r - 1) / 2r, (r + 1) / 2r)` and above `p_{i-1}^(1 / (1 - p_{i-1}))`. In this domain the cells synchronize to a
    common concentration, and a chaotic environment, `r` near 4, desynchronizes them.

    Args:
        n_cells:
            The number of cells of the contour.

        r:
            The logistic parameter of every cell.

        c:
            The coupling of every cell to the next, one value or one per cell.

        p:
            The affinity of every cell for the uptake, one value or one per cell.

        random_state:
            The seed or generator of the initial concentrations, uniform between 0.2 and 0.8.

    Returns:
        system:
            The system, with the concentrations `cells.x`, cell `i + 1` next to cell `i`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import closed_contour_cell_exchange
        >>> closed_contour_cell_exchange(5).elements[0].prepare().size_
        5
    """
    rng = np.random.default_rng(random_state)
    cells = GroupElement(
        'cells',
        element=ContourCellElement,
        size=n_cells,
        values={
            'x': rng.uniform(0.2, 0.8, n_cells),
            'r': r,
            'c': np.broadcast_to(np.asarray(c, dtype=float), n_cells).copy(),
            'p': np.broadcast_to(np.asarray(p, dtype=float), n_cells).copy(),
        },
    )
    uptake = ContourUptakeInteraction('uptake', target='cells')
    return System(elements=[cells], interactions=[uptake], time=DiscreteTime())

coagulation_positive_feedback(*, injury=1.0, injury_time=1.0, duration=1.0, ka=3.0)

Return the coagulation cascade, an all-or-nothing burst of thrombin bounded by protein C and fibrinolysis.

Thrombin activates its own production through the factors of the cascade, a positive feedback that amplifies an injury above a threshold into a burst of thrombin and a clot of fibrin, while a weaker injury fades. With the default rates, the burst needs a thrombin of (ka - sqrt(ka^2 - 4)) / 2, the unstable state of the feedback. Thrombin also activates protein C, which switches the burst off, and the fibrin activates plasmin, which dissolves the clot, so the clot stays at the site and in time. The positive feedback, the off switch, and the fibrinolysis follow Konieczny et al.

Parameters:

Name Type Description Default
injury float

The rate of thrombin activation by the injury.

 1.0
injury_time float

The start of the injury.

 1.0
duration float

The duration of the injury.

 1.0
ka float

The maximal rate of the positive feedback, which is bistable above 2.

 3.0

Returns:

Name Type Description
system System

The system, with the thrombin blood.T, the protein C blood.C, the fibrin blood.F, and the plasmin blood.P.

Examples:

>>> from skcomplex.systems.biology.physiology import coagulation_positive_feedback
>>> coagulation_positive_feedback(injury=0.2).elements[0].injury
0.2
Source code in src/skcomplex/systems/biology/physiology/_regulation.py 
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
@register_system(
    field='biology/physiology',
    tags=('coagulation', 'positive feedback', 'threshold', 'excitable', 'fibrinolysis', 'ode'),
    references=('Konieczny et al. 2023',),
    defaults={'until': 120.0, 'dt': 0.1},
)
def coagulation_positive_feedback(
    *,
    injury: float = 1.0,
    injury_time: float = 1.0,
    duration: float = 1.0,
    ka: float = 3.0,
) -> System:
    """Return the coagulation cascade, an all-or-nothing burst of thrombin bounded by protein C and fibrinolysis.

    Thrombin activates its own production through the factors of the cascade, a positive feedback that amplifies an
    injury above a threshold into a burst of thrombin and a clot of fibrin, while a weaker injury fades. With the
    default rates, the burst needs a thrombin of `(ka - sqrt(ka^2 - 4)) / 2`, the unstable state of the feedback.
    Thrombin also activates protein C, which switches the burst off, and the fibrin activates plasmin, which dissolves
    the clot, so the clot stays at the site and in time. The positive feedback, the off switch, and the fibrinolysis
    follow Konieczny et al.

    Args:
        injury:
            The rate of thrombin activation by the injury.

        injury_time:
            The start of the injury.

        duration:
            The duration of the injury.

        ka:
            The maximal rate of the positive feedback, which is bistable above 2.

    Returns:
        system:
            The system, with the thrombin `blood.T`, the protein C `blood.C`, the fibrin `blood.F`, and the plasmin
            `blood.P`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import coagulation_positive_feedback
        >>> coagulation_positive_feedback(injury=0.2).elements[0].injury
        0.2
    """
    blood = CoagulationElement('blood', injury=injury, injury_time=injury_time, duration=duration, ka=ka)
    interactions = [
        RateInteraction(
            'trigger',
            target='blood',
            state='T',
            rate=lambda s: s.target.injury * pulse(s.t, blood.injury_time, blood.duration),
        ),
        ReactionInteraction(
            'thrombin_feedback',
            target='blood',
            products={'T': 1},
            propensity=lambda s: s.target.ka * hill(s.target.T, s.target.k_half, 2),
        ),
        ReactionInteraction(
            'thrombin_decay',
            target='blood',
            reactants={'T': 1},
            propensity=lambda s: s.target.d_t * s.target.T,
        ),
        ReactionInteraction(
            'thrombin_inhibition',
            target='blood',
            reactants={'T': 1},
            modifiers={'C': 1},
            propensity=lambda s: s.target.k_c * s.target.C * s.target.T,
        ),
        ReactionInteraction(
            'inhibitor_production',
            target='blood',
            products={'C': 1},
            propensity=lambda s: s.target.s_c * s.target.T,
        ),
        ReactionInteraction(
            'inhibitor_decay',
            target='blood',
            reactants={'C': 1},
            propensity=lambda s: s.target.d_c * s.target.C,
        ),
        ReactionInteraction(
            'fibrin_production',
            target='blood',
            products={'F': 1},
            propensity=lambda s: s.target.s_f * s.target.T,
        ),
        ReactionInteraction(
            'fibrinolysis',
            target='blood',
            reactants={'F': 1},
            modifiers={'P': 1},
            propensity=lambda s: s.target.k_p * s.target.P * s.target.F,
        ),
        ReactionInteraction(
            'plasmin_production',
            target='blood',
            products={'P': 1},
            propensity=lambda s: s.target.s_p * s.target.F,
        ),
        ReactionInteraction(
            'plasmin_decay',
            target='blood',
            reactants={'P': 1},
            propensity=lambda s: s.target.d_p * s.target.P,
        ),
    ]
    return System(elements=[blood], interactions=interactions)

drug_toxicity(*, k=0.5, eta=0.05, dose=1.0, period=7.0, width=1.0, doses=6)

Return the toxicity a weekly drug accumulates in the body, which the metabolism clears once the doses stop.

Every dose raises the plasma concentration, which the body eliminates within days, and the toxicity integrates the concentration, dT/dt = c - eta T. With eta = 0 the toxicity is the cumulative exposure and never falls.

Parameters:

Name Type Description Default
k float

The elimination rate of the drug, per day.

 0.5
eta float

The rate at which the metabolism lowers the toxicity, per day.

 0.05
dose float

The infusion rate during a dose.

 1.0
period float

The time between the starts of doses, in days.

 7.0
width float

The duration of a dose, in days.

 1.0
doses int

The number of doses.

 6

Returns:

Name Type Description
system System

The system, with the concentration patient.c and the toxicity patient.T.

Examples:

>>> from skcomplex.systems.biology.physiology import drug_toxicity
>>> drug_toxicity(eta=0.0).elements[0].eta
0.0
Source code in src/skcomplex/systems/biology/physiology/_toxicity.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
@register_system(
    field='biology/physiology',
    tags=('pharmacology', 'toxicity', 'ode', 'dosing'),
    references=('Padmanabhan, Meskin, and Al Moustafa 2021', 'Martin and Teo 1994'),
    defaults={'until': 84.0, 'dt': 0.2},
)
def drug_toxicity(
    *,
    k: float = 0.5,
    eta: float = 0.05,
    dose: float = 1.0,
    period: float = 7.0,
    width: float = 1.0,
    doses: int = 6,
) -> System:
    """Return the toxicity a weekly drug accumulates in the body, which the metabolism clears once the doses stop.

    Every dose raises the plasma concentration, which the body eliminates within days, and the toxicity integrates
    the concentration, `dT/dt = c - eta T`. With `eta = 0` the toxicity is the cumulative exposure and never falls.

    Args:
        k:
            The elimination rate of the drug, per day.

        eta:
            The rate at which the metabolism lowers the toxicity, per day.

        dose:
            The infusion rate during a dose.

        period:
            The time between the starts of doses, in days.

        width:
            The duration of a dose, in days.

        doses:
            The number of doses.

    Returns:
        system:
            The system, with the concentration `patient.c` and the toxicity `patient.T`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import drug_toxicity
        >>> drug_toxicity(eta=0.0).elements[0].eta
        0.0
    """
    patient = DrugToxicityElement('patient', k=k, eta=eta, dose=dose, period=period, width=width, doses=doses)
    interactions = [
        DosingInteraction(
            'dosing',
            target='patient',
            state='c',
            dose=lambda s: s.target.dose,
            period=period,
            width=width,
            count=doses,
        ),
        ReactionInteraction(
            'clearance',
            target='patient',
            reactants={'c': 1},
            propensity=lambda s: s.target.k * s.target.c,
        ),
        RateInteraction('toxicity', target='patient', state='T', rate=lambda s: s.target.c),
        ReactionInteraction(
            'toxicity_recovery',
            target='patient',
            reactants={'T': 1},
            propensity=lambda s: s.target.eta * s.target.T,
        ),
    ]
    return System(elements=[patient], interactions=interactions)

fitzhugh_nagumo(shape=(100, 100), *, D=1.0, a=0.7, b=0.8, epsilon=0.02)

Return the FitzHugh-Nagumo excitable medium, started from a broken wave that curls into a spiral.

The upper half starts excited and the left half refractory, so the free end of the wave rotates.

Parameters:

Name Type Description Default
shape Sequence[int]

The number of rows and columns of the lattice with zero-flux edges.

 (100, 100)
D float

The diffusion of the potential u.

 1.0
a float

The offset of the recovery.

 0.7
b float

The damping of the recovery.

 0.8
epsilon float

The time scale of the recovery.

 0.02

Returns:

Name Type Description
system System

The system, with the potential medium.u and the recovery medium.v.

Examples:

>>> from skcomplex.systems.biology.physiology import fitzhugh_nagumo
>>> fitzhugh_nagumo((20, 20)).interactions[0].epsilon
0.02
Source code in src/skcomplex/systems/biology/physiology/_physiology.py 
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
@register_system(
    field='biology/physiology',
    tags=('field', 'excitable media', 'spirals', 'pde'),
    references=('FitzHugh 1961', 'Nagumo, Arimoto, and Yoshizawa 1962', 'Gros 2015'),
    defaults={'until': 200.0, 'dt': 2.0, 'solver': RK4Solver(dt=0.05)},
)
def fitzhugh_nagumo(
    shape: Sequence[int] = (100, 100),
    *,
    D: float = 1.0,
    a: float = 0.7,
    b: float = 0.8,
    epsilon: float = 0.02,
) -> System:
    """Return the FitzHugh-Nagumo excitable medium, started from a broken wave that curls into a spiral.

    The upper half starts excited and the left half refractory, so the free end of the wave rotates.

    Args:
        shape:
            The number of rows and columns of the lattice with zero-flux edges.

        D:
            The diffusion of the potential `u`.

        a:
            The offset of the recovery.

        b:
            The damping of the recovery.

        epsilon:
            The time scale of the recovery.

    Returns:
        system:
            The system, with the potential `medium.u` and the recovery `medium.v`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import fitzhugh_nagumo
        >>> fitzhugh_nagumo((20, 20)).interactions[0].epsilon
        0.02
    """
    rows, columns = shape
    u, v = np.full(shape, -1.2), np.full(shape, -0.62)
    u[: rows // 2, :] = 1.5
    v[:, : columns // 2] = 1.0
    space = LatticeSpace((rows, columns), boundary='reflecting')
    medium = FieldElement('medium', element=FitzHughNagumoElement, space=space, values={'u': u, 'v': v})
    kinetics = FitzHughNagumoInteraction('kinetics', target='medium', fast='u', slow='v', a=a, b=b, epsilon=epsilon)
    diffusion = DiffusionInteraction('diffusion', target='medium', state='u', coefficient=D)
    return System(elements=[medium], interactions=[kinetics, diffusion])

glucose_dual_hormone_control(*, glucagon=True, meal=1.0, exercise=0.8, meal_time=2.0, exercise_time=12.0, duration=2.0)

Return blood glucose stabilised by antagonistic hormone loops, insulin against highs and glucagon against lows.

A meal raises the glucose and secretes insulin, which raises the uptake of glucose, and an exercise lowers it and secretes glucagon, which releases glucose from the liver. Each receptor acts in one direction only, so one loop alone corrects excursions to one side of the norm. Two opposing loops on the same quantity attenuate the excursions to both sides, after Konieczny et al. The time unit is about one hour, and the glucose is relative to its norm.

Parameters:

Name Type Description Default
glucagon bool

Whether the glucagon loop acts, False for insulin alone.

 True
meal float

The rate of glucose intake during the meal.

 1.0
exercise float

The rate of extra glucose uptake during the exercise.

 0.8
meal_time float

The start of the meal.

 2.0
exercise_time float

The start of the exercise.

 12.0
duration float

The duration of the meal and of the exercise.

 2.0

Returns:

Name Type Description
system System

The system, with the glucose organism.G, the insulin organism.I, and the glucagon organism.N.

Examples:

>>> from skcomplex.systems.biology.physiology import glucose_dual_hormone_control
>>> glucose_dual_hormone_control(glucagon=False).elements[0].k_n
0.0
Source code in src/skcomplex/systems/biology/physiology/_regulation.py 
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
@register_system(
    field='biology/physiology',
    tags=('homeostasis', 'hormones', 'insulin', 'glucagon', 'negative feedback', 'ode'),
    references=('Konieczny et al. 2023',),
    defaults={'until': 24.0, 'dt': 0.05},
)
def glucose_dual_hormone_control(
    *,
    glucagon: bool = True,
    meal: float = 1.0,
    exercise: float = 0.8,
    meal_time: float = 2.0,
    exercise_time: float = 12.0,
    duration: float = 2.0,
) -> System:
    """Return blood glucose stabilised by antagonistic hormone loops, insulin against highs and glucagon against lows.

    A meal raises the glucose and secretes insulin, which raises the uptake of glucose, and an exercise lowers it and
    secretes glucagon, which releases glucose from the liver. Each receptor acts in one direction only, so one loop
    alone corrects excursions to one side of the norm. Two opposing loops on the same quantity attenuate the
    excursions to both sides, after Konieczny et al. The time unit is about one hour, and the glucose is relative to
    its norm.

    Args:
        glucagon:
            Whether the glucagon loop acts, `False` for insulin alone.

        meal:
            The rate of glucose intake during the meal.

        exercise:
            The rate of extra glucose uptake during the exercise.

        meal_time:
            The start of the meal.

        exercise_time:
            The start of the exercise.

        duration:
            The duration of the meal and of the exercise.

    Returns:
        system:
            The system, with the glucose `organism.G`, the insulin `organism.I`, and the glucagon `organism.N`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import glucose_dual_hormone_control
        >>> glucose_dual_hormone_control(glucagon=False).elements[0].k_n
        0.0
    """
    k_n = 1.0 if glucagon else 0.0
    insulin, glucagon_level = 2.0 * 0.5, 1.0 * 0.5
    organism = GlucoseControlElement(
        'organism',
        I=insulin,
        N=glucagon_level,
        R=(0.5 + 0.5 * insulin) - k_n * glucagon_level,
        k_n=k_n,
        meal=meal,
        exercise=exercise,
        meal_time=meal_time,
        exercise_time=exercise_time,
        duration=duration,
    )
    interactions = [
        ReactionInteraction('release', target='organism', products={'G': 1}, propensity=lambda s: s.target.R),
        ReactionInteraction(
            'basal_uptake',
            target='organism',
            reactants={'G': 1},
            propensity=lambda s: s.target.k0 * s.target.G,
        ),
        ReactionInteraction(
            'insulin_uptake',
            target='organism',
            reactants={'G': 1},
            modifiers={'I': 1},
            propensity=lambda s: s.target.k_i * s.target.I * s.target.G,
        ),
        ReactionInteraction(
            'glucagon_release',
            target='organism',
            products={'G': 1},
            propensity=lambda s: s.target.k_n * s.target.N,
        ),
        RateInteraction(
            'meal',
            target='organism',
            state='G',
            rate=lambda s: s.target.meal * pulse(s.t, organism.meal_time, organism.duration),
        ),
        RateInteraction(
            'exercise',
            target='organism',
            state='G',
            rate=lambda s: -s.target.exercise * pulse(s.t, organism.exercise_time, organism.duration),
        ),
        ReactionInteraction(
            'insulin_secretion',
            target='organism',
            products={'I': 1},
            propensity=lambda s: s.target.s_i * hill(s.target.G, s.target.K, s.target.n),
        ),
        ReactionInteraction(
            'insulin_clearance',
            target='organism',
            reactants={'I': 1},
            propensity=lambda s: s.target.d_i * s.target.I,
        ),
        ReactionInteraction(
            'glucagon_secretion',
            target='organism',
            products={'N': 1},
            propensity=lambda s: s.target.s_n * (1 - hill(s.target.G, s.target.K, s.target.n)),
        ),
        ReactionInteraction(
            'glucagon_clearance',
            target='organism',
            reactants={'N': 1},
            propensity=lambda s: s.target.d_n * s.target.N,
        ),
    ]
    return System(elements=[organism], interactions=interactions)

hyperbolic_chemotaxis(sites=40, *, s=100.0, D_n=1.0, D_c=0.001, alpha=0.33, tau=1.0, n0=0.25, width=0.03)

Return cells and a chemical moving with finite speed, whose two bumps of cells spread and attract each other.

The model is the hyperbolic limit of a kinetic model of cells and a chemical that change velocity at random, with the speed s and turning rates set by the diffusions D_n = s^2 / (2 mu1) and D_c = s^2 / (2 sigma) of its parabolic limit, Keller-Segel, and mu2 = mu1. Following the test of Bellomo et al. (2017, section 6.6.2), the cells start as two Gaussian bumps of mass n0 and width width at (3 width, 3 width) and (-3 width, -3 width) in the square [-0.4, 0.4]^2, here with periodic edges, the cells secrete the chemical, and all fluxes start at rest. The fluxes relax fast, so the bumps spread out at a finite speed while the cell mass 2 n0 is conserved and the chemical accumulates where the cells are.

Parameters:

Name Type Description Default
sites int

The number of sites along each side of the square.

 40
s float

The speed of the cells and of the chemical.

 100.0
D_n float

The diffusion of the cells in the parabolic limit.

 1.0
D_c float

The diffusion of the chemical in the parabolic limit.

 0.001
alpha float

The chemotactic sensitivity of the cells.

 0.33
tau float

The time scale of the chemical.

 1.0
n0 float

The mass of each bump of cells.

 0.25
width float

The width of each bump of cells.

 0.03

Returns:

Name Type Description
system System

The system, with the cell density tissue.n, the chemical tissue.c, and their fluxes.

Examples:

>>> from skcomplex.systems.biology.physiology import hyperbolic_chemotaxis
>>> round(hyperbolic_chemotaxis(20).elements[0].values['mu1'])
5000
Source code in src/skcomplex/systems/biology/physiology/_chemotaxis.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
@register_system(
    field='biology/physiology',
    tags=('chemotaxis', 'hyperbolic', 'kinetic theory', 'field', 'pde'),
    references=('Bellomo et al. 2017',),
    defaults={'until': 0.015, 'dt': 0.0005, 'solver': ScipyODESolver(method='BDF')},
)
def hyperbolic_chemotaxis(
    sites: int = 40,
    *,
    s: float = 100.0,
    D_n: float = 1.0,
    D_c: float = 0.001,
    alpha: float = 0.33,
    tau: float = 1.0,
    n0: float = 0.25,
    width: float = 0.03,
) -> System:
    """Return cells and a chemical moving with finite speed, whose two bumps of cells spread and attract each other.

    The model is the hyperbolic limit of a kinetic model of cells and a chemical that change velocity at random, with
    the speed `s` and turning rates set by the diffusions `D_n = s^2 / (2 mu1)` and `D_c = s^2 / (2 sigma)` of its
    parabolic limit, Keller-Segel, and `mu2 = mu1`. Following the test of Bellomo et al. (2017, section 6.6.2), the
    cells start as two Gaussian bumps of mass `n0` and width `width` at `(3 width, 3 width)` and `(-3 width, -3 width)`
    in the square `[-0.4, 0.4]^2`, here with periodic edges, the cells secrete the chemical, and all fluxes start at
    rest. The fluxes relax fast, so the bumps spread out at a finite speed while the cell mass `2 n0` is conserved and
    the chemical accumulates where the cells are.

    Args:
        sites:
            The number of sites along each side of the square.

        s:
            The speed of the cells and of the chemical.

        D_n:
            The diffusion of the cells in the parabolic limit.

        D_c:
            The diffusion of the chemical in the parabolic limit.

        alpha:
            The chemotactic sensitivity of the cells.

        tau:
            The time scale of the chemical.

        n0:
            The mass of each bump of cells.

        width:
            The width of each bump of cells.

    Returns:
        system:
            The system, with the cell density `tissue.n`, the chemical `tissue.c`, and their fluxes.

    Examples:
        >>> from skcomplex.systems.biology.physiology import hyperbolic_chemotaxis
        >>> round(hyperbolic_chemotaxis(20).elements[0].values['mu1'])
        5000
    """
    spacing = 0.8 / sites
    axis = -0.4 + spacing * (np.arange(sites) + 0.5)
    x, y = np.meshgrid(axis, axis, indexing='ij')
    offset = 3 * width
    bumps = np.exp(-((x - offset) ** 2 + (y - offset) ** 2) / (2 * width**2))
    bumps += np.exp(-((x + offset) ** 2 + (y + offset) ** 2) / (2 * width**2))
    turning = s**2 / (2 * D_n)
    values = {
        'n': n0 / (2 * np.pi * width**2) * bumps,
        's': s,
        'mu1': turning,
        'mu2': turning,
        'sigma': s**2 / (2 * D_c),
        'tau': tau,
        'alpha': alpha,
    }
    space = LatticeSpace((sites, sites), boundary='periodic', spacing=spacing)
    return System(elements=[FieldElement('tissue', element=ChemotacticTissueElement, space=space, values=values)])

mackey_glass(x=0.5, *, beta=0.2, gamma=0.1, n=10.0, tau=17.0)

Return the Mackey-Glass equation x' = beta x(t - tau) / (1 + x(t - tau)^n) - gamma x of blood cell production.

It is periodic for tau = 6 and chaotic for tau = 17. The history before the start is the initial value.

Parameters:

Name Type Description Default
x float

The initial concentration, also its history.

 0.5
beta float

The production rate.

 0.2
gamma float

The decay rate.

 0.1
n float

The steepness of the production feedback.

 10.0
tau float

The delay of production.

 17.0

Returns:

Name Type Description
system System

The system, with the concentration cells.x.

Examples:

>>> from skcomplex.systems.biology.physiology import mackey_glass
>>> mackey_glass(tau=6.0).elements[0].tau
6.0
Source code in src/skcomplex/systems/biology/physiology/_physiology.py 
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
@register_system(
    field='biology/physiology',
    tags=('chaos', 'delay', 'physiology'),
    references=('Mackey and Glass 1977', 'Erneux 2009'),
    defaults={'until': 500.0, 'dt': 0.5},
)
def mackey_glass(
    x: float = 0.5,
    *,
    beta: float = 0.2,
    gamma: float = 0.1,
    n: float = 10.0,
    tau: float = 17.0,
) -> System:
    """Return the Mackey-Glass equation `x' = beta x(t - tau) / (1 + x(t - tau)^n) - gamma x` of blood cell production.

    It is periodic for `tau = 6` and chaotic for `tau = 17`. The history before the start is the initial value.

    Args:
        x:
            The initial concentration, also its history.

        beta:
            The production rate.

        gamma:
            The decay rate.

        n:
            The steepness of the production feedback.

        tau:
            The delay of production.

    Returns:
        system:
            The system, with the concentration `cells.x`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import mackey_glass
        >>> mackey_glass(tau=6.0).elements[0].tau
        6.0
    """
    return System(elements=[MackeyGlassElement('cells', x=x, beta=beta, gamma=gamma, n=n, tau=tau)])

nanotube_gap_junction_exchange(graph=None, *, alpha=1.0, beta=10.0, g=0.1, xi=0.0, x=None, random_state=0)

Return the intercellular exchange of a substance through gap junctions and tunnelling nanotubes on a network.

Every cell metabolizes the substance at the Michaelis-Menten rate alpha x / (beta + x) and exchanges it with its neighbours, by diffusion g (x_j - x_i) through gap junctions and by the transfer xi x_j through nanotubes. The state without substance is stable while the largest eigenvalue of xi A - g L stays below alpha / beta, where A is the adjacency and L the Laplacian of the network. Past the threshold xi_crit the nanotubes destabilize the communication and the concentrations grow. On a regular graph of degree k, xi_crit = alpha / (beta k).

Parameters:

Name Type Description Default
graph Graph | None

The NetworkX graph of the gap junctions, or None for a modular Erdős-Rényi network of 100 cells in ten modules, with an attachment probability of 0.03 overall and 90% of the links within modules.

 None
alpha float

The maximal rate of metabolism.

 1.0
beta float

The half-saturation concentration of metabolism.

 10.0
g float

The exchange coefficient of the gap junctions.

 0.1
xi float

The strength of the nanotubes.

 0.0
x float | Sequence[float] | None

The initial concentrations, one value or one per node in node order, or None for values drawn uniformly between 0 and 1.

 None
random_state int | Generator | None

The seed or generator of the default network and of the initial concentrations.

 0

Returns:

Name Type Description
system System

The system, with the concentrations cells.x, one value per node.

Examples:

>>> import networkx as nx
>>> from skcomplex.systems.biology.physiology import nanotube_gap_junction_exchange
>>> nanotube_gap_junction_exchange(nx.cycle_graph(6), xi=0.01).elements[0].prepare().size_
6
Source code in src/skcomplex/systems/biology/physiology/_exchange.py 
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
@register_system(
    field='biology/physiology',
    tags=('network', 'cell communication', 'gap junctions', 'nanotubes', 'stability', 'ode'),
    references=('Mihailović et al. 2023',),
    defaults={'until': 100.0, 'dt': 0.5},
)
def nanotube_gap_junction_exchange(
    graph: nx.Graph | None = None,
    *,
    alpha: float = 1.0,
    beta: float = 10.0,
    g: float = 0.1,
    xi: float = 0.0,
    x: float | Sequence[float] | None = None,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the intercellular exchange of a substance through gap junctions and tunnelling nanotubes on a network.

    Every cell metabolizes the substance at the Michaelis-Menten rate `alpha x / (beta + x)` and exchanges it with
    its neighbours, by diffusion `g (x_j - x_i)` through gap junctions and by the transfer `xi x_j` through nanotubes.
    The state without substance is stable while the largest eigenvalue of `xi A - g L` stays below `alpha / beta`,
    where `A` is the adjacency and `L` the Laplacian of the network. Past the threshold `xi_crit` the nanotubes
    destabilize the communication and the concentrations grow. On a regular graph of degree `k`,
    `xi_crit = alpha / (beta k)`.

    Args:
        graph:
            The NetworkX graph of the gap junctions, or `None` for a modular Erdős-Rényi network of 100 cells in ten
            modules, with an attachment probability of 0.03 overall and 90% of the links within modules.

        alpha:
            The maximal rate of metabolism.

        beta:
            The half-saturation concentration of metabolism.

        g:
            The exchange coefficient of the gap junctions.

        xi:
            The strength of the nanotubes.

        x:
            The initial concentrations, one value or one per node in node order, or `None` for values drawn uniformly
            between 0 and 1.

        random_state:
            The seed or generator of the default network and of the initial concentrations.

    Returns:
        system:
            The system, with the concentrations `cells.x`, one value per node.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.systems.biology.physiology import nanotube_gap_junction_exchange
        >>> nanotube_gap_junction_exchange(nx.cycle_graph(6), xi=0.01).elements[0].prepare().size_
        6
    """
    rng = np.random.default_rng(random_state)
    if graph is None:
        graph = nx.random_partition_graph([10] * 10, 0.297, 0.0033, seed=int(rng.integers(2**31)))
    size = graph.number_of_nodes()
    initial = rng.uniform(0.0, 1.0, size) if x is None else np.broadcast_to(np.asarray(x, dtype=float), size).copy()
    values = {
        'x': initial,
        'alpha': alpha,
        'beta': beta,
        'g': g,
        'xi': xi,
    }
    cells = FieldElement('cells', element=ExchangingCellElement, space=NetworkSpace(graph), values=values)
    exchange = CouplingInteraction('exchange', target='cells', read='x', write='x', function=_exchange)
    return System(elements=[cells], interactions=[exchange], time=ContinuousTime())

steerable_feedback_loop(*, v=2.0, d=1.0, K=1.0, n=4.0, T=1.0, tau=0.5, steering=0.5, start=20.0, duration=20.0)

Return a receptor-effector loop that holds a product at its set point until a steering signal shifts it.

The loop of the negative feedback system of Konieczny et al. (2023) settles at the set point where synthesis v / (1 + (x / K)^n) balances consumption d x. A steering signal, such as a bacterial toxin acting on the temperature receptors in fever, raises the sensitivity of the receptor by the factor 1 + steering, so the loop climbs to a higher set point while it acts and returns once it ends. A slow effector or a long transmission delay turns the regulation into sustained oscillations around the set point.

Parameters:

Name Type Description Default
v float

The largest rate of synthesis.

 2.0
d float

The consumption rate of the product.

 1.0
K float

The sensitivity of the receptor.

 1.0
n float

The steepness of the receptor.

 4.0
T float

The response time of the effector.

 1.0
tau float

The transmission delay between the receptor and the effector.

 0.5
steering float

The relative shift of the sensitivity by the steering signal.

 0.5
start float

The time the steering signal starts.

 20.0
duration float

The duration of the steering signal.

 20.0

Returns:

Name Type Description
system System

The system, with the product loop.x and the effector activity loop.e.

Examples:

>>> from skcomplex.systems.biology.physiology import steerable_feedback_loop
>>> steerable_feedback_loop(tau=3.0).elements[0].tau
3.0
Source code in src/skcomplex/systems/biology/physiology/_regulation.py 
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
@register_system(
    field='biology/physiology',
    tags=('negative feedback', 'homeostasis', 'set point', 'fever', 'delay'),
    references=('Wach et al. 2015', 'Konieczny et al. 2023'),
    defaults={'until': 60.0, 'dt': 0.1},
)
def steerable_feedback_loop(
    *,
    v: float = 2.0,
    d: float = 1.0,
    K: float = 1.0,
    n: float = 4.0,
    T: float = 1.0,
    tau: float = 0.5,
    steering: float = 0.5,
    start: float = 20.0,
    duration: float = 20.0,
) -> System:
    """Return a receptor-effector loop that holds a product at its set point until a steering signal shifts it.

    The loop of the negative feedback system of Konieczny et al. (2023) settles at the set point where synthesis
    `v / (1 + (x / K)^n)` balances consumption `d x`. A steering signal, such as a bacterial toxin acting on the
    temperature receptors in fever, raises the sensitivity of the receptor by the factor `1 + steering`, so the loop
    climbs to a higher set point while it acts and returns once it ends. A slow effector or a long transmission delay
    turns the regulation into sustained oscillations around the set point.

    Args:
        v:
            The largest rate of synthesis.

        d:
            The consumption rate of the product.

        K:
            The sensitivity of the receptor.

        n:
            The steepness of the receptor.

        T:
            The response time of the effector.

        tau:
            The transmission delay between the receptor and the effector.

        steering:
            The relative shift of the sensitivity by the steering signal.

        start:
            The time the steering signal starts.

        duration:
            The duration of the steering signal.

    Returns:
        system:
            The system, with the product `loop.x` and the effector activity `loop.e`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import steerable_feedback_loop
        >>> steerable_feedback_loop(tau=3.0).elements[0].tau
        3.0
    """
    loop = FeedbackLoopElement(
        'loop',
        v=v,
        d=d,
        K=K,
        n=n,
        T=T,
        tau=tau,
        steering=steering,
        start=start,
        duration=duration,
    )
    return System(elements=[loop])

two_compartment_pharmacokinetics(c0=1.0, *, k10=0.082, k12=0.0475, k21=0.00125)

Return the two-compartment pharmacokinetics of a bolus dose given into the central compartment.

The central concentration follows dUc/dt = k21 Up - (k10 + k12) Uc and the peripheral one dUp/dt = k12 Uc - k21 Up, from Uc = c0 and Up = 0. The central concentration is the bi-exponential c0 ((alpha - k21) e^(-alpha t) + (k21 - beta) e^(-beta t)) / (alpha - beta), where alpha and beta are the roots of x^2 - (k10 + k12 + k21) x + k10 k21: a fast distribution phase is followed by a slow elimination phase. The total drug Uc + Up falls only by the elimination k10 Uc. The rates are those of doxorubicin, per minute.

Parameters:

Name Type Description Default
c0 float

The bolus, the initial central concentration.

 1.0
k10 float

The elimination rate from the central compartment, per minute.

 0.082
k12 float

The transfer rate from the central to the peripheral compartment, per minute.

 0.0475
k21 float

The transfer rate from the peripheral to the central compartment, per minute.

 0.00125

Returns:

Name Type Description
system System

The system, with the central concentration patient.Uc and the peripheral one patient.Up.

Examples:

>>> from skcomplex.systems.biology.physiology import two_compartment_pharmacokinetics
>>> two_compartment_pharmacokinetics(c0=5.0).elements[0].Uc
5.0
Source code in src/skcomplex/systems/biology/physiology/_pharmacokinetics.py 
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
@register_system(
    field='biology/physiology',
    tags=('pharmacology', 'pharmacokinetics', 'ode', 'compartments'),
    references=('Padmanabhan, Meskin, and Al Moustafa 2021', 'Hendriks et al. 2012'),
    defaults={'until': 3000.0, 'dt': 10.0},
)
def two_compartment_pharmacokinetics(
    c0: float = 1.0,
    *,
    k10: float = 0.082,
    k12: float = 0.0475,
    k21: float = 0.00125,
) -> System:
    """Return the two-compartment pharmacokinetics of a bolus dose given into the central compartment.

    The central concentration follows `dUc/dt = k21 Up - (k10 + k12) Uc` and the peripheral one
    `dUp/dt = k12 Uc - k21 Up`, from `Uc = c0` and `Up = 0`. The central concentration is the bi-exponential
    `c0 ((alpha - k21) e^(-alpha t) + (k21 - beta) e^(-beta t)) / (alpha - beta)`, where `alpha` and `beta` are the
    roots of `x^2 - (k10 + k12 + k21) x + k10 k21`: a fast distribution phase is followed by a slow elimination
    phase. The total drug `Uc + Up` falls only by the elimination `k10 Uc`. The rates are those of doxorubicin, per
    minute.

    Args:
        c0:
            The bolus, the initial central concentration.

        k10:
            The elimination rate from the central compartment, per minute.

        k12:
            The transfer rate from the central to the peripheral compartment, per minute.

        k21:
            The transfer rate from the peripheral to the central compartment, per minute.

    Returns:
        system:
            The system, with the central concentration `patient.Uc` and the peripheral one `patient.Up`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import two_compartment_pharmacokinetics
        >>> two_compartment_pharmacokinetics(c0=5.0).elements[0].Uc
        5.0
    """
    reactions = [
        ReactionInteraction('distribution', reactants={('patient', 'Uc'): 1}, products={('patient', 'Up'): 1}, k=k12),
        ReactionInteraction('redistribution', reactants={('patient', 'Up'): 1}, products={('patient', 'Uc'): 1}, k=k21),
        ReactionInteraction('elimination', reactants={('patient', 'Uc'): 1}, k=k10),
    ]
    return System(elements=[TwoCompartmentElement('patient', Uc=c0)], interactions=reactions)

two_layer_sepsis(n_cells=20, *, beta=0.5 * math.pi, sigma=1.0, pathological=0.0, omega_healthy=0.0, omega_pathological=1.0, omega_immune=0.0, alpha_p=0.3 * math.pi, alpha_s=0.35 * math.pi, epsilon_p=0.03, epsilon_s=0.3, connectivity=0.5, seed=0)

Return the two-layer model of tumour disease and sepsis, parenchymal and immune cells as adaptive oscillators.

The parenchymal layer of n_cells phase oscillators couples through fixed links a_ij, drawn as 1 with the probability connectivity, and adaptive cytokine links kappa_ij. The immune layer couples through adaptive links alone, which adapt ten times faster, and every cell couples to its partner in the other layer with strength sigma. Following Eqs. (40) and (41) of Berner et al. (2023), dphi_i/dt = omega_i - sum_j (a_ij + kappa_ij) sin(phi_i - phi_j + alpha) / N - sigma sin(phi_i - phi_i') and dkappa_ij/dt = -epsilon (kappa_ij + sin(phi_i - phi_j - beta)). The healthy state is the frequency-synchronized one, and the age parameter beta and a fraction pathological of cells with another natural frequency push the layers into desynchronized, pathological states.

Parameters:

Name Type Description Default
n_cells int

The number of cells of each layer.

 20
beta float

The age parameter of the adaptation.

 0.5 * pi
sigma float

The coupling between the layers.

 1.0
pathological float

The fraction of pathological parenchymal cells.

 0.0
omega_healthy float

The natural frequency of healthy parenchymal cells.

 0.0
omega_pathological float

The natural frequency of pathological parenchymal cells.

 1.0
omega_immune float

The natural frequency of the immune cells.

 0.0
alpha_p float

The phase lag of the parenchymal layer.

 0.3 * pi
alpha_s float

The phase lag of the immune layer.

 0.35 * pi
epsilon_p float

The adaptation rate of the parenchymal links.

 0.03
epsilon_s float

The adaptation rate of the immune links.

 0.3
connectivity float

The probability of a fixed parenchymal link.

 0.5
seed int | None

The seed of the fixed links, the pathological cells, and the initial phases and couplings.

 0

Returns:

Name Type Description
system System

The system, with the phases cells.phi, the parenchymal cells first and the immune cells after, and the couplings parenchymal_links.kappa and immune_links.kappa.

Examples:

>>> from skcomplex.systems.biology.physiology import two_layer_sepsis
>>> two_layer_sepsis(5).elements[0].prepare().size_
10
Source code in src/skcomplex/systems/biology/physiology/_disease.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
@register_system(
    field='biology/physiology',
    tags=('sepsis', 'adaptive network', 'multiplex', 'phase oscillators', 'ode'),
    references=('Sawicki et al. 2022', 'Berner et al. 2023'),
    defaults={'until': 200.0, 'dt': 1.0},
)
def two_layer_sepsis(
    n_cells: int = 20,
    *,
    beta: float = 0.5 * math.pi,
    sigma: float = 1.0,
    pathological: float = 0.0,
    omega_healthy: float = 0.0,
    omega_pathological: float = 1.0,
    omega_immune: float = 0.0,
    alpha_p: float = 0.3 * math.pi,
    alpha_s: float = 0.35 * math.pi,
    epsilon_p: float = 0.03,
    epsilon_s: float = 0.3,
    connectivity: float = 0.5,
    seed: int | None = 0,
) -> System:
    """Return the two-layer model of tumour disease and sepsis, parenchymal and immune cells as adaptive oscillators.

    The parenchymal layer of `n_cells` phase oscillators couples through fixed links `a_ij`, drawn as 1 with the
    probability `connectivity`, and adaptive cytokine links `kappa_ij`. The immune layer couples through adaptive links
    alone, which adapt ten times faster, and every cell couples to its partner in the other layer with strength
    `sigma`. Following Eqs. (40) and (41) of Berner et al. (2023), `dphi_i/dt = omega_i - sum_j (a_ij + kappa_ij)
    sin(phi_i - phi_j + alpha) / N - sigma sin(phi_i - phi_i')` and `dkappa_ij/dt = -epsilon (kappa_ij + sin(phi_i -
    phi_j - beta))`. The healthy state is the frequency-synchronized one, and the age parameter `beta` and a fraction
    `pathological` of cells with another natural frequency push the layers into desynchronized, pathological states.

    Args:
        n_cells:
            The number of cells of each layer.

        beta:
            The age parameter of the adaptation.

        sigma:
            The coupling between the layers.

        pathological:
            The fraction of pathological parenchymal cells.

        omega_healthy:
            The natural frequency of healthy parenchymal cells.

        omega_pathological:
            The natural frequency of pathological parenchymal cells.

        omega_immune:
            The natural frequency of the immune cells.

        alpha_p:
            The phase lag of the parenchymal layer.

        alpha_s:
            The phase lag of the immune layer.

        epsilon_p:
            The adaptation rate of the parenchymal links.

        epsilon_s:
            The adaptation rate of the immune links.

        connectivity:
            The probability of a fixed parenchymal link.

        seed:
            The seed of the fixed links, the pathological cells, and the initial phases and couplings.

    Returns:
        system:
            The system, with the phases `cells.phi`, the parenchymal cells first and the immune cells after, and the
            couplings `parenchymal_links.kappa` and `immune_links.kappa`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import two_layer_sepsis
        >>> two_layer_sepsis(5).elements[0].prepare().size_
        10
    """
    rng = np.random.default_rng(seed)
    n = n_cells
    parenchyma = nx.complete_graph(range(n), create_using=nx.DiGraph)
    parenchyma.add_nodes_from(range(n, 2 * n))
    immune = nx.DiGraph()
    immune.add_nodes_from(range(2 * n))
    immune.add_edges_from((i, j) for i in range(n, 2 * n) for j in range(n, 2 * n) if i != j)
    partners = nx.DiGraph()
    partners.add_nodes_from(range(2 * n))
    partners.add_edges_from([(i, i + n) for i in range(n)] + [(i + n, i) for i in range(n)])
    space = NetworkSpace(parenchyma, layers=[LayerSpec('immune', immune), LayerSpec('partners', partners)])
    frequencies = np.full(2 * n, omega_immune)
    sick = rng.random(n) < pathological
    frequencies[:n] = np.where(sick, omega_pathological, omega_healthy)
    cells = LatticeElement(
        'cells',
        element=CellOscillatorElement,
        space=space,
        values={'phi': rng.uniform(0, 2 * math.pi, 2 * n), 'omega': frequencies},
    )
    edges = n * (n - 1)
    parenchymal_links = GroupElement(
        'parenchymal_links',
        element=CytokineLinkElement,
        space=space,
        on='edges',
        values={'kappa': rng.uniform(-1, 1, edges), 'a': (rng.random(edges) < connectivity).astype(float)},
    )
    immune_links = GroupElement(
        'immune_links',
        element=CytokineLinkElement,
        space=space,
        on='edges',
        layer='immune',
        values={'kappa': rng.uniform(-1, 1, edges)},
    )
    return System(
        elements=[cells, parenchymal_links, immune_links],
        interactions=[
            AdaptivePhaseCouplingInteraction(
                'parenchymal_coupling',
                cells='cells',
                links='parenchymal_links',
                N=n,
                alpha=alpha_p,
                beta=beta,
                epsilon=epsilon_p,
            ),
            AdaptivePhaseCouplingInteraction(
                'immune_coupling',
                cells='cells',
                links='immune_links',
                N=n,
                alpha=alpha_s,
                beta=beta,
                epsilon=epsilon_s,
            ),
            LayerCouplingInteraction('interlayer', cells='cells', sigma=sigma, layer='partners'),
        ],
    )

west_brown_enquist(depth=6, *, n=2, r_c=1.0, l_c=1.0)

Return the distribution network of West, Brown, and Enquist, a space-filling tree of invariant capillaries.

The tree branches n ways at every level down to the capillaries at level depth. The metabolic rate is proportional to the number of capillaries, n^depth, and the body mass to the blood volume, which the area-preserving and space-filling vessels make V = pi r_c^2 l_c n^depth (n^(depth / 3) - 1) / (n^(1 / 3) - 1), which grows as n^(4 depth / 3) for deep trees. So the metabolic rate scales as the mass to the power 3/4, Kleiber's law, which the exponent between consecutive depths approaches from below. An Experiment over depth fits it.

Parameters:

Name Type Description Default
depth int

The level of the capillaries, the number of branchings.

 6
n int

The branching ratio.

 2
r_c float

The radius of a capillary.

 1.0
l_c float

The length of a capillary.

 1.0

Returns:

Name Type Description
system System

The system, with the radii vessels.r, lengths vessels.l, and volumes vessels.volume on the edges of the tree and the observables blood_volume and capillaries of the catalog, which a run passes in observables.

Examples:

>>> from skcomplex.systems.biology.physiology import west_brown_enquist
>>> west_brown_enquist(3).elements[0].prepare().size_
14
Source code in src/skcomplex/systems/biology/physiology/_distribution_networks.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
@register_system(
    field='biology/physiology',
    tags=('allometry', 'kleiber law', 'fractal network', 'static', 'metabolic scaling'),
    references=('West, Brown, and Enquist 1997', 'Ribeiro and Netto 2024'),
    observables=_network_totals,
)
def west_brown_enquist(depth: int = 6, *, n: int = 2, r_c: float = 1.0, l_c: float = 1.0) -> System:
    """Return the distribution network of West, Brown, and Enquist, a space-filling tree of invariant capillaries.

    The tree branches `n` ways at every level down to the capillaries at level `depth`. The metabolic rate is
    proportional to the number of capillaries, `n^depth`, and the body mass to the blood volume, which the
    area-preserving and space-filling vessels make `V = pi r_c^2 l_c n^depth (n^(depth / 3) - 1) / (n^(1 / 3) - 1)`,
    which grows as `n^(4 depth / 3)` for deep trees. So the metabolic rate scales as the mass to the power 3/4,
    Kleiber's law, which the exponent between consecutive depths approaches from below. An `Experiment` over `depth`
    fits it.

    Args:
        depth:
            The level of the capillaries, the number of branchings.

        n:
            The branching ratio.

        r_c:
            The radius of a capillary.

        l_c:
            The length of a capillary.

    Returns:
        system:
            The system, with the radii `vessels.r`, lengths `vessels.l`, and volumes `vessels.volume` on the edges of
            the tree and the observables `blood_volume` and `capillaries` of the catalog, which a run passes in
            `observables`.

    Examples:
        >>> from skcomplex.systems.biology.physiology import west_brown_enquist
        >>> west_brown_enquist(3).elements[0].prepare().size_
        14
    """
    graph = nx.balanced_tree(n, depth)
    space = NetworkSpace(graph).prepare()
    _, children = space.edges_[None]
    levels = nx.single_source_shortest_path_length(graph, 0)
    vessels = GroupElement(
        'vessels',
        element=BranchingVesselElement,
        space=space,
        on='edges',
        values={
            'level': np.array([levels[int(child)] for child in children], dtype=float),
            'depth': float(depth),
            'n': float(n),
            'r_c': r_c,
            'l_c': l_c,
        },
    )
    return System(elements=[vessels], space=space, time=StaticTime())