Skip to content

Physiology

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

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.models.biology.physiology import AdaptiveResponderElement
>>> AdaptiveResponderElement('organism').period
10.0
Source code in src/skcomplex/models/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.models.biology.physiology import AlievPanfilovElement
>>> AlievPanfilovElement('tissue').k
8.0
Source code in src/skcomplex/models/biology/physiology/_cardiac.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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

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.models.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/models/biology/physiology/_chamber.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
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

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.models.biology.physiology import CoagulationElement
>>> CoagulationElement('blood').injury_time
1.0
Source code in src/skcomplex/models/biology/physiology/_coagulation.py 
 90
 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
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.models.biology.physiology import ContourCellElement
>>> ContourCellElement('cells').r
2.5
Source code in src/skcomplex/models/biology/physiology/_contour_cell.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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.elements import GroupElement
>>> from skcomplex.core.spaces import DiscreteTime
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.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).check().prepare().run(until=1).trajectory_.values('cells.x')[-1].round(4).tolist()
[0.7344, 0.25]
Source code in src/skcomplex/models/biology/physiology/_contour_uptake.py 
40
41
42
43
44
45
46
47
48
49
50
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,
    )

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.models.biology.physiology import DrugToxicityElement
>>> DrugToxicityElement('patient', doses=4).doses
4
Source code in src/skcomplex/models/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.models.biology.physiology import ExchangingCellElement
>>> ExchangingCellElement('cells').g
0.1
Source code in src/skcomplex/models/biology/physiology/_exchanging_cell.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
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

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.models.biology.physiology import FitzHughNagumoElement
>>> FitzHughNagumoElement('cell', u=1.5).u
1.5
Source code in src/skcomplex/models/biology/physiology/_fitz_hugh_nagumo.py 
36
37
38
39
40
41
42
43
44
45
46
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

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.models.biology.physiology import GlucoseControlElement
>>> GlucoseControlElement('organism').meal_time
2.0
Source code in src/skcomplex/models/biology/physiology/_glucose_control.py 
 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
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

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.models.biology.physiology import MackeyGlassElement
>>> MackeyGlassElement('blood', tau=6.0).tau
6.0
Source code in src/skcomplex/models/biology/physiology/_mackey_glass.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
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.models.biology.physiology import ResistiveFlowInteraction
>>> ResistiveFlowInteraction('aortic_valve', source='ventricle', target='arteries', valve=True).valve
True
Source code in src/skcomplex/models/biology/physiology/_resistive_flow.py 
45
46
47
48
49
50
51
52
53
54
55
56
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

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.models.biology.physiology import TwoCompartmentElement
>>> TwoCompartmentElement('patient').Up
0.0
Source code in src/skcomplex/models/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

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 model, 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.models.biology.physiology import adaptive_response_costs
>>> adaptive_response_costs('stimulus-response').elements[0].eta
0.0
Source code in src/skcomplex/models/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_model(
    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 model, 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.models.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 model 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 models 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 model, with the potential tissue.u and the recovery tissue.v.

Examples:

>>> from skcomplex.models.biology.physiology import aliev_panfilov
>>> aliev_panfilov((20, 20)).elements[0].values['k']
8.0
Source code in src/skcomplex/models/biology/physiology/_cardiac.py 
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
165
166
167
@register_model(
    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 model 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 models
    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 model, with the potential `tissue.u` and the recovery `tissue.v`.

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

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 model 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 model 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 model, with the volumes ventricle.V, arteries.V, and veins.V.

Examples:

>>> from skcomplex.models.biology.physiology import cardiovascular
>>> sum(chamber.V for chamber in cardiovascular().elements)
1500.0
Source code in src/skcomplex/models/biology/physiology/_circulation.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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
@register_model(
    field='biology/physiology',
    tags=('cardiovascular', 'circulation', 'lumped model', '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 model 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 model
    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 model, with the volumes `ventricle.V`, `arteries.V`, and `veins.V`.

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

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 model, with the concentrations cells.x, cell i + 1 next to cell i.

Examples:

>>> from skcomplex.models.biology.physiology import closed_contour_cell_exchange
>>> closed_contour_cell_exchange(5).elements[0].check().prepare().size_
5
Source code in src/skcomplex/models/biology/physiology/_exchange.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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
@register_model(
    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 model, with the concentrations `cells.x`, cell `i + 1` next to cell `i`.

    Examples:
        >>> from skcomplex.models.biology.physiology import closed_contour_cell_exchange
        >>> closed_contour_cell_exchange(5).elements[0].check().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 model, with the thrombin blood.T, the protein C blood.C, the fibrin blood.F, and the plasmin blood.P.

Examples:

>>> from skcomplex.models.biology.physiology import coagulation_positive_feedback
>>> coagulation_positive_feedback(injury=0.2).elements[0].injury
0.2
Source code in src/skcomplex/models/biology/physiology/_regulation.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
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
@register_model(
    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 model, with the thrombin `blood.T`, the protein C `blood.C`, the fibrin `blood.F`, and the plasmin
            `blood.P`.

    Examples:
        >>> from skcomplex.models.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 model, with the concentration patient.c and the toxicity patient.T.

Examples:

>>> from skcomplex.models.biology.physiology import drug_toxicity
>>> drug_toxicity(eta=0.0).elements[0].eta
0.0
Source code in src/skcomplex/models/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_model(
    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 model, with the concentration `patient.c` and the toxicity `patient.T`.

    Examples:
        >>> from skcomplex.models.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 model, with the potential medium.u and the recovery medium.v.

Examples:

>>> from skcomplex.models.biology.physiology import fitzhugh_nagumo
>>> fitzhugh_nagumo((20, 20)).interactions[0].epsilon
0.02
Source code in src/skcomplex/models/biology/physiology/_physiology.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
@register_model(
    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 model, with the potential `medium.u` and the recovery `medium.v`.

    Examples:
        >>> from skcomplex.models.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 model, with the glucose organism.G, the insulin organism.I, and the glucagon organism.N.

Examples:

>>> from skcomplex.models.biology.physiology import glucose_dual_hormone_control
>>> glucose_dual_hormone_control(glucagon=False).elements[0].k_n
0.0
Source code in src/skcomplex/models/biology/physiology/_regulation.py 
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 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
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
@register_model(
    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 model, with the glucose `organism.G`, the insulin `organism.I`, and the glucagon `organism.N`.

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

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 model, with the concentration cells.x.

Examples:

>>> from skcomplex.models.biology.physiology import mackey_glass
>>> mackey_glass(tau=6.0).elements[0].tau
6.0
Source code in src/skcomplex/models/biology/physiology/_physiology.py 
 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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
@register_model(
    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 model, with the concentration `cells.x`.

    Examples:
        >>> from skcomplex.models.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 model, with the concentrations cells.x, one value per node.

Examples:

>>> import networkx as nx
>>> from skcomplex.models.biology.physiology import nanotube_gap_junction_exchange
>>> nanotube_gap_junction_exchange(nx.cycle_graph(6), xi=0.01).elements[0].check().prepare().size_
6
Source code in src/skcomplex/models/biology/physiology/_exchange.py 
 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
165
@register_model(
    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 model, with the concentrations `cells.x`, one value per node.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.models.biology.physiology import nanotube_gap_junction_exchange
        >>> nanotube_gap_junction_exchange(nx.cycle_graph(6), xi=0.01).elements[0].check().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())

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 model, with the central concentration patient.Uc and the peripheral one patient.Up.

Examples:

>>> from skcomplex.models.biology.physiology import two_compartment_pharmacokinetics
>>> two_compartment_pharmacokinetics(c0=5.0).elements[0].Uc
5.0
Source code in src/skcomplex/models/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_model(
    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 model, with the central concentration `patient.Uc` and the peripheral one `patient.Up`.

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