Skip to content

Kinetics

Chemical kinetics: reactions, oscillating reactors, and excitable media.

CatalystPelletElement(label, visual=None, *, c=1.0, h=0.1, k=10.0, K=20.0)

Bases: Element

Catalyst pellet, the concentration c of the reactant at its surface relative to the bulk.

The mass transfer and the surface reaction of catalyst_pellet change the concentration and read their constants from the pellet.

Parameters:

Name Type Description Default
label str

The label of the pellet.

 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 surface concentration.

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

The coefficient of mass transfer from the bulk.

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

The rate constant of the surface reaction.

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

The adsorption constant, which sets the inhibition.

 20.0

Examples:

>>> from skcomplex.models.chemistry.kinetics import CatalystPelletElement
>>> CatalystPelletElement('pellet').K
20.0
Source code in src/skcomplex/models/chemistry/kinetics/_catalyst_pellet.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    c: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    h: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 20.0,
) -> None:
    super().__init__(label, visual=visual)
    self.c = c
    self.h = h
    self.k = k
    self.K = K

CubicAutocatalysisElement(label, visual=None, *, a=1.0, b=0.0, F=0.04, k=0.05)

Bases: Element

Continuously fed reactor of cubic autocatalysis, the concentrations of the reactant a and the autocatalyst b.

The reactions and the flows of cubic_autocatalysis change the concentrations and read the feed rate F and the decay rate k from the reactor.

Parameters:

Name Type Description Default
label str

The label of the reactor.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

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

The initial concentration of the reactant.

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

The initial concentration of the autocatalyst.

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

The feed rate, the inverse residence time.

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

The decay rate of the autocatalyst.

 0.05

Examples:

>>> from skcomplex.models.chemistry.kinetics import CubicAutocatalysisElement
>>> CubicAutocatalysisElement('reactor').F
0.04
Source code in src/skcomplex/models/chemistry/kinetics/_cubic_autocatalysis.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    a: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    b: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    F: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.04,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
) -> None:
    super().__init__(label, visual=visual)
    self.a = a
    self.b = b
    self.F = F
    self.k = k

CyclicCellElement(label, visual=None, *, level=0, n=14.0, threshold=1.0)

Bases: Element

Cell of Griffeath's cyclic cellular automaton, with n levels that follow each other in a cycle.

A cell moves from level k to level k + 1 modulo n when at least threshold of its neighbours hold that level, by the local rule of cyclic_cellular_automaton, so every level eats its predecessor and random initial levels organize into spiral waves.

Parameters:

Name Type Description Default
label str

The label of the cell or the lattice group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

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

The initial level, from 0 to n - 1.

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

The number of levels.

 14.0
threshold Annotated[float, ParameterSpec(domain=open(0, oo))]

The number of neighbours at the following level that advance a cell.

 1.0

Examples:

>>> from skcomplex.models.chemistry.kinetics import CyclicCellElement
>>> CyclicCellElement('cell').n
14.0
Source code in src/skcomplex/models/chemistry/kinetics/_cyclic.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    level: Annotated[float, StateSpec(domain=NN)] = 0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 14.0,
    threshold: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.level = level
    self.n = n
    self.threshold = threshold

HodgepodgeCellElement(label, visual=None, *, level=0, n=100.0, k1=2.0, k2=3.0, g=34.0)

Bases: Element

Cell of the hodgepodge machine of Gerhardt and Schuster, healthy at 0, infected below n, and ill at n.

Healthy cells catch the infection from their infected and ill neighbours, infected cells worsen by the mean of their neighbourhood plus g, and ill cells heal, by the local rule of hodgepodge.

Parameters:

Name Type Description Default
label str

The label of the cell or the lattice group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

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

The initial level of infection from 0 to n.

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

The level of illness.

 100.0
k1 Annotated[float, ParameterSpec(domain=open(0, oo))]

The weight of infected neighbours for healthy cells.

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

The weight of ill neighbours for healthy cells.

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

The speed of the infection.

 34.0

Examples:

>>> from skcomplex.models.chemistry.kinetics import HodgepodgeCellElement
>>> HodgepodgeCellElement('cell').n
100.0
Source code in src/skcomplex/models/chemistry/kinetics/_hodgepodge.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    level: Annotated[float, StateSpec(domain=NN)] = 0,
    n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
    k1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    k2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 3.0,
    g: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 34.0,
) -> None:
    super().__init__(label, visual=visual)
    self.level = level
    self.n = n
    self.k1 = k1
    self.k2 = k2
    self.g = g

OregonatorElement(label, visual=None, *, x=0.1, y=0.1, z=0.1, epsilon=0.04, delta=0.0004, q=0.0008, f=1.0)

Bases: Element

Reactor of the scaled Oregonator, the scaled concentrations of bromous acid, bromide, and the oxidized catalyst.

The scaled concentrations follow epsilon dx/dt = q y - x y + x (1 - x), delta dy/dt = -q y - x y + f z, and dz/dt = x - z. The reactions of oregonator change the concentrations, with the time scales epsilon and delta dividing the stoichiometric coefficients of x and y, and read the parameters from the reactor.

Parameters:

Name Type Description Default
label str

The label of the reactor.

 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 scaled concentration of bromous acid.

 0.1
y Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial scaled concentration of bromide.

 0.1
z Annotated[float, StateSpec(domain=closed_open(0, oo))]

The initial scaled concentration of the oxidized catalyst.

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

The ratio of the time scales of bromous acid and the catalyst.

 0.04
delta Annotated[float, ParameterSpec(domain=open(0, oo))]

The ratio of the time scales of bromide and the catalyst.

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

The ratio of rate constants that sets the lowest bromous acid concentration.

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

The stoichiometric factor of bromide production from the catalyst.

 1.0

Examples:

>>> from skcomplex.models.chemistry.kinetics import OregonatorElement
>>> OregonatorElement('reactor').epsilon
0.04
Source code in src/skcomplex/models/chemistry/kinetics/_oregonator.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    y: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    z: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.04,
    delta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.0004,
    q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0008,
    f: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.z = z
    self.epsilon = epsilon
    self.delta = delta
    self.q = q
    self.f = f

PulsedElectrodeElement(label, visual=None, *, eta=0.0, c=1.0, C=1.0, alpha=0.5, tau=10.0, kappa=0.05, current=5.0, period=10.0, width=2.0)

Bases: Element

Electrode of a stirred electrochemical reactor, charged by a train of current pulses.

In units of the exchange current and of the thermal voltage, the overpotential follows C deta/dt = i(t) - i_F, with the Butler-Volmer current i_F = c exp(alpha eta) - exp(-(1 - alpha) eta), and the reactant concentration follows dc/dt = (1 - c) / tau - kappa i_F.

Parameters:

Name Type Description Default
label str

The label of the electrode.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
eta Annotated[float, StateSpec()]

The initial overpotential.

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

The initial reactant concentration, relative to the feed.

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

The capacitance of the double layer.

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

The transfer coefficient.

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

The residence time of the reactor.

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

The consumption of reactant per unit of Faradaic current.

 0.05
current Annotated[float, ParameterSpec()]

The current of a pulse.

 5.0
period float

The time between the starts of pulses.

 10.0
width float

The duration of a pulse.

 2.0

Examples:

>>> from skcomplex.models.chemistry.kinetics import PulsedElectrodeElement
>>> PulsedElectrodeElement('electrode').period
10.0
Source code in src/skcomplex/models/chemistry/kinetics/_pulsed_electrode.py 
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
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    eta: Annotated[float, StateSpec()] = 0.0,
    c: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    C: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    current: Annotated[float, ParameterSpec()] = 5.0,
    period: float = 10.0,
    width: float = 2.0,
) -> None:
    super().__init__(label, visual=visual)
    self.eta = eta
    self.c = c
    self.C = C
    self.alpha = alpha
    self.tau = tau
    self.kappa = kappa
    self.current = current
    self.period = period
    self.width = width

SchloglElement(label, visual=None, *, x=50.0, c1=1.5e-06, c2=0.0025, c3=0.0002, c4=3.5, A=100000.0, B=200000.0)

Bases: Element

Reactor of Schlögl's reaction, the number of molecules of X with the buffered species A and B.

The reactions of chemical_langevin_equation change x and read the rate constants and the constant numbers of molecules of A and B from the reactor.

Parameters:

Name Type Description Default
label str

The label of the reactor.

 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 number of molecules of X.

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

The rate constant of the autocatalytic production.

 1.5e-06
c2 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The rate constant of the reverse of the autocatalytic production.

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

The rate constant of the production from B.

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

The rate constant of the conversion into B.

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

The constant number of molecules of A.

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

The constant number of molecules of B.

 200000.0

Examples:

>>> from skcomplex.models.chemistry.kinetics import SchloglElement
>>> SchloglElement('reactor').x
50.0
Source code in src/skcomplex/models/chemistry/kinetics/_reactions.py 
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 50.0,
    c1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.5e-06,
    c2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0025,
    c3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0002,
    c4: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.5,
    A: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 100000.0,
    B: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 200000.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.c1 = c1
    self.c2 = c2
    self.c3 = c3
    self.c4 = c4
    self.A = A
    self.B = B

brusselator_ode(u=1.0, v=1.0, *, a=1.0, b=3.0)

Return the well-mixed Brusselator, which oscillates on a limit cycle when b > 1 + a**2.

The reactions A -> U, B + U -> V + D, 2U + V -> 3U, and U -> E with the buffered reactants A and B give du/dt = a - (b + 1) u + u^2 v and dv/dt = b u - u^2 v, with the fixed point (a, b / a).

Parameters:

Name Type Description Default
u float

The initial concentration of the autocatalytic intermediate.

 1.0
v float

The initial concentration of the second intermediate.

 1.0
a float

The constant supply of the first reactant.

 1.0
b float

The constant supply of the second reactant.

 3.0

Returns:

Name Type Description
system System

The model, with the concentrations reactor.u and reactor.v.

Examples:

>>> from skcomplex.models.chemistry.kinetics import brusselator_ode
>>> brusselator_ode(b=1.5).elements[0].b
1.5
Source code in src/skcomplex/models/chemistry/kinetics/_oscillators.py 
 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
@register_model(
    field='chemistry/kinetics',
    tags=('ode', 'oscillation', 'limit cycle', 'autocatalysis'),
    references=('Prigogine and Lefever 1968', 'Datseris and Parlitz 2022', 'Trzaska 2023'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def brusselator_ode(u: float = 1.0, v: float = 1.0, *, a: float = 1.0, b: float = 3.0) -> System:
    """Return the well-mixed Brusselator, which oscillates on a limit cycle when `b > 1 + a**2`.

    The reactions `A -> U`, `B + U -> V + D`, `2U + V -> 3U`, and `U -> E` with the buffered reactants `A` and `B`
    give `du/dt = a - (b + 1) u + u^2 v` and `dv/dt = b u - u^2 v`, with the fixed point `(a, b / a)`.

    Args:
        u:
            The initial concentration of the autocatalytic intermediate.

        v:
            The initial concentration of the second intermediate.

        a:
            The constant supply of the first reactant.

        b:
            The constant supply of the second reactant.

    Returns:
        system:
            The model, with the concentrations `reactor.u` and `reactor.v`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import brusselator_ode
        >>> brusselator_ode(b=1.5).elements[0].b
        1.5
    """
    return System(
        elements=[BrusselatorReactionElement('reactor', u=u, v=v, a=a, b=b)],
        interactions=[
            ReactionInteraction('feed', target='reactor', products={'u': 1}, modifiers={'a': 1}),
            ReactionInteraction(
                'conversion',
                target='reactor',
                reactants={'u': 1},
                products={'v': 1},
                modifiers={'b': 1},
            ),
            CubicAutocatalysisInteraction('autocatalysis', target='reactor', substrate='v', autocatalyst='u'),
            ReactionInteraction('decay', target='reactor', reactants={'u': 1}),
        ],
    )

catalyst_pellet(c=1.0, *, h=0.1, k=10.0, K=20.0)

Return the surface concentration of a catalyst pellet, bistable through Langmuir-Hinshelwood kinetics.

The concentration at the surface, relative to the bulk, follows dc/dt = h (1 - c) - k c / (1 + K c)^2, where mass transfer from the bulk balances a reaction that the adsorbed reactant inhibits. For the default parameters a reacting state of low concentration and an inhibited state of high concentration are both stable.

Parameters:

Name Type Description Default
c float

The initial surface concentration.

 1.0
h float

The coefficient of mass transfer from the bulk.

 0.1
k float

The rate constant of the surface reaction.

 10.0
K float

The adsorption constant, which sets the inhibition.

 20.0

Returns:

Name Type Description
system System

The model, with the concentration pellet.c.

Examples:

>>> from skcomplex.models.chemistry.kinetics import catalyst_pellet
>>> catalyst_pellet(k=5.0).elements[0].k
5.0
Source code in src/skcomplex/models/chemistry/kinetics/_reactors.py 
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
@register_model(
    field='chemistry/kinetics',
    tags=('ode', 'catalysis', 'bistability', 'Langmuir-Hinshelwood'),
    references=('Trzaska 2023', 'Aris 1975'),
    defaults={'until': 100.0, 'dt': 1.0},
)
def catalyst_pellet(c: float = 1.0, *, h: float = 0.1, k: float = 10.0, K: float = 20.0) -> System:
    """Return the surface concentration of a catalyst pellet, bistable through Langmuir-Hinshelwood kinetics.

    The concentration at the surface, relative to the bulk, follows `dc/dt = h (1 - c) - k c / (1 + K c)^2`, where
    mass transfer from the bulk balances a reaction that the adsorbed reactant inhibits. For the default parameters a
    reacting state of low concentration and an inhibited state of high concentration are both stable.

    Args:
        c:
            The initial surface concentration.

        h:
            The coefficient of mass transfer from the bulk.

        k:
            The rate constant of the surface reaction.

        K:
            The adsorption constant, which sets the inhibition.

    Returns:
        system:
            The model, with the concentration `pellet.c`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import catalyst_pellet
        >>> catalyst_pellet(k=5.0).elements[0].k
        5.0
    """
    return System(
        elements=[CatalystPelletElement('pellet', c=c, h=h, k=k, K=K)],
        interactions=[
            RelaxationInteraction('transfer', target='pellet', state='c', goal=1.0, time=_get_transfer_time),
            ReactionInteraction('reaction', reactants={('pellet', 'c'): 1}, propensity=_get_surface_reaction),
        ],
    )

chemical_langevin_equation(x=50.0, *, c1=1.5e-06, c2=0.0025, c3=0.0002, c4=3.5, A=100000.0, B=200000.0)

Return the chemical Langevin equation of Schlögl's bistable reaction, whose noise switches between two states.

The channels A + 2X -> 3X, 3X -> A + 2X, B -> X, and X -> B have the propensities c1 A x (x - 1) / 2, c2 x (x - 1) (x - 2) / 6, c3 B, and c4 x. Every channel adds its propensity times its change of x to the drift and drives x with its own Wiener process of amplitude the square root of its propensity. The default rate constants are those of Gillespie 2000 scaled to a volume five times smaller, so the stable states lie near 17 and 113 molecules and the noise switches between them within 200 time units.

Parameters:

Name Type Description Default
x float

The initial number of molecules of X.

 50.0
c1 float

The rate constant of the autocatalytic production.

 1.5e-06
c2 float

The rate constant of the reverse of the autocatalytic production.

 0.0025
c3 float

The rate constant of the production from B.

 0.0002
c4 float

The rate constant of the conversion into B.

 3.5
A float

The constant number of molecules of A.

 100000.0
B float

The constant number of molecules of B.

 200000.0

Returns:

Name Type Description
system System

The model, with the number of molecules reactor.x.

Examples:

>>> from skcomplex.models.chemistry.kinetics import chemical_langevin_equation
>>> chemical_langevin_equation(c4=3.0).elements[0].c4
3.0
Source code in src/skcomplex/models/chemistry/kinetics/_reactions.py 
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
@register_model(
    field='chemistry/kinetics',
    tags=('sde', 'chemical Langevin equation', 'bistability', 'noise'),
    references=('Schlögl 1972', 'Gillespie 2000', 'Trzaska 2023'),
    defaults={'until': 200.0, 'dt': 0.05, 'solver': EulerMaruyamaSolver(dt=1e-3)},
)
def chemical_langevin_equation(
    x: float = 50.0,
    *,
    c1: float = 1.5e-6,
    c2: float = 2.5e-3,
    c3: float = 2e-4,
    c4: float = 3.5,
    A: float = 1e5,
    B: float = 2e5,
) -> System:
    """Return the chemical Langevin equation of Schlögl's bistable reaction, whose noise switches between two states.

    The channels `A + 2X -> 3X`, `3X -> A + 2X`, `B -> X`, and `X -> B` have the propensities `c1 A x (x - 1) / 2`,
    `c2 x (x - 1) (x - 2) / 6`, `c3 B`, and `c4 x`. Every channel adds its propensity times its change of `x` to the
    drift and drives `x` with its own Wiener process of amplitude the square root of its propensity. The default rate
    constants are those of Gillespie 2000 scaled to a volume five times smaller, so the stable states lie near 17 and
    113 molecules and the noise switches between them within 200 time units.

    Args:
        x:
            The initial number of molecules of `X`.

        c1:
            The rate constant of the autocatalytic production.

        c2:
            The rate constant of the reverse of the autocatalytic production.

        c3:
            The rate constant of the production from `B`.

        c4:
            The rate constant of the conversion into `B`.

        A:
            The constant number of molecules of `A`.

        B:
            The constant number of molecules of `B`.

    Returns:
        system:
            The model, with the number of molecules `reactor.x`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import chemical_langevin_equation
        >>> chemical_langevin_equation(c4=3.0).elements[0].c4
        3.0
    """
    return System(
        elements=[SchloglElement('reactor', x=x, c1=c1, c2=c2, c3=c3, c4=c4, A=A, B=B)],
        interactions=[
            ReactionInteraction(
                'autocatalysis',
                reactants={X: 2},
                products={X: 3},
                propensity=_get_autocatalysis,
                mode='langevin',
            ),
            ReactionInteraction(
                'reversal',
                reactants={X: 3},
                products={X: 2},
                propensity=_get_reversal,
                mode='langevin',
            ),
            ReactionInteraction('production', products={X: 1}, propensity=_get_production, mode='langevin'),
            ReactionInteraction('conversion', reactants={X: 1}, propensity=_get_conversion, mode='langevin'),
        ],
    )

chemical_langevin_network(enzyme=100.0, substrate=300.0, *, k_binding=0.01, k_unbinding=0.1, k_catalysis=0.1)

Return the chemical Langevin equation of the enzyme reaction E + S <-> C -> E + P.

Every reaction adds its rate times its change of every species to the drift, and one Wiener process with the amplitude the square root of its rate drives every species it changes. Since one increment changes all of them, the total enzyme E + C and the total substrate S + C + P stay constant under the noise, as in the jump process.

Parameters:

Name Type Description Default
enzyme float

The initial amount of free enzyme.

 100.0
substrate float

The initial amount of substrate.

 300.0
k_binding float

The rate constant of the binding of enzyme and substrate.

 0.01
k_unbinding float

The rate constant of the unbinding of the complex.

 0.1
k_catalysis float

The rate constant of the conversion of the complex into enzyme and product.

 0.1

Returns:

Name Type Description
system System

The model, with the amounts E.n, S.n, C.n, and P.n.

Examples:

>>> from skcomplex.models.chemistry.kinetics import chemical_langevin_network
>>> [reaction.mode for reaction in chemical_langevin_network().interactions]
['langevin', 'langevin', 'langevin']
Source code in src/skcomplex/models/chemistry/kinetics/_reactions.py 
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
@register_model(
    field='chemistry/kinetics',
    tags=('enzyme kinetics', 'sde', 'chemical Langevin equation', 'noise'),
    references=('Michaelis and Menten 1913', 'Gillespie 2000'),
    defaults={'until': 50.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=1e-3)},
)
def chemical_langevin_network(
    enzyme: float = 100.0,
    substrate: float = 300.0,
    *,
    k_binding: float = 0.01,
    k_unbinding: float = 0.1,
    k_catalysis: float = 0.1,
) -> System:
    """Return the chemical Langevin equation of the enzyme reaction `E + S <-> C -> E + P`.

    Every reaction adds its rate times its change of every species to the drift, and one Wiener process with the
    amplitude the square root of its rate drives every species it changes. Since one increment changes all of them, the
    total enzyme `E + C` and the total substrate `S + C + P` stay constant under the noise, as in the jump process.

    Args:
        enzyme:
            The initial amount of free enzyme.

        substrate:
            The initial amount of substrate.

        k_binding:
            The rate constant of the binding of enzyme and substrate.

        k_unbinding:
            The rate constant of the unbinding of the complex.

        k_catalysis:
            The rate constant of the conversion of the complex into enzyme and product.

    Returns:
        system:
            The model, with the amounts `E.n`, `S.n`, `C.n`, and `P.n`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import chemical_langevin_network
        >>> [reaction.mode for reaction in chemical_langevin_network().interactions]
        ['langevin', 'langevin', 'langevin']
    """
    return System(
        elements=[
            CompartmentElement('E', n=enzyme),
            CompartmentElement('S', n=substrate),
            CompartmentElement('C'),
            CompartmentElement('P'),
        ],
        interactions=_get_enzyme_reactions(k_binding, k_unbinding, k_catalysis, 'langevin'),
    )

cubic_autocatalysis(a=0.5, b=0.5, *, F=0.04, k=0.05)

Return cubic autocatalysis A + 2B -> 3B with the decay B -> C in a continuously fed reactor.

The flow through the reactor relaxes the reactant to its feed concentration 1 and washes out the autocatalyst at the rate F, so the dimensionless concentrations follow da/dt = -a b^2 + F (1 - a) and db/dt = a b^2 - (F + k) b. The washout state (1, 0) is always stable, and a reacting state coexists with it when F + k < sqrt(F) / 2.

Parameters:

Name Type Description Default
a float

The initial concentration of the reactant.

 0.5
b float

The initial concentration of the autocatalyst.

 0.5
F float

The feed rate, the inverse residence time.

 0.04
k float

The decay rate of the autocatalyst.

 0.05

Returns:

Name Type Description
system System

The model, with the concentrations reactor.a and reactor.b.

Examples:

>>> from skcomplex.models.chemistry.kinetics import cubic_autocatalysis
>>> cubic_autocatalysis(k=0.06).elements[0].k
0.06
Source code in src/skcomplex/models/chemistry/kinetics/_oscillators.py 
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
@register_model(
    field='chemistry/kinetics',
    tags=('ode', 'autocatalysis', 'bistability', 'reactor'),
    references=('Gray and Scott 1983', 'Trzaska 2023'),
    defaults={'until': 500.0, 'dt': 1.0},
)
def cubic_autocatalysis(a: float = 0.5, b: float = 0.5, *, F: float = 0.04, k: float = 0.05) -> System:
    """Return cubic autocatalysis `A + 2B -> 3B` with the decay `B -> C` in a continuously fed reactor.

    The flow through the reactor relaxes the reactant to its feed concentration 1 and washes out the autocatalyst at
    the rate `F`, so the dimensionless concentrations follow `da/dt = -a b^2 + F (1 - a)` and
    `db/dt = a b^2 - (F + k) b`. The washout state `(1, 0)` is always stable, and a reacting state coexists with it
    when `F + k < sqrt(F) / 2`.

    Args:
        a:
            The initial concentration of the reactant.

        b:
            The initial concentration of the autocatalyst.

        F:
            The feed rate, the inverse residence time.

        k:
            The decay rate of the autocatalyst.

    Returns:
        system:
            The model, with the concentrations `reactor.a` and `reactor.b`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import cubic_autocatalysis
        >>> cubic_autocatalysis(k=0.06).elements[0].k
        0.06
    """
    return System(
        elements=[CubicAutocatalysisElement('reactor', a=a, b=b, F=F, k=k)],
        interactions=[
            CubicAutocatalysisInteraction('autocatalysis', target='reactor', substrate='a', autocatalyst='b'),
            ReactionInteraction('decay', target='reactor', reactants={'b': 1}, modifiers={'k': 1}),
            RelaxationInteraction('feed', target='reactor', state='a', goal=1.0, time=_get_residence),
            RelaxationInteraction('outflow', target='reactor', state='b', goal=0.0, time=_get_residence),
        ],
    )

cyclic_cellular_automaton(shape=(64, 64), *, n=14, threshold=1, random_state=0)

Return Griffeath's cyclic cellular automaton, whose random levels organize into spiral waves.

Parameters:

Name Type Description Default
shape Sequence[int]

The number of rows and columns.

 (64, 64)
n int

The number of levels.

 14
threshold int

The number of von Neumann neighbours at the following level that advance a cell.

 1
random_state int | Generator | None

The seed or generator of the initial levels.

 0

Returns:

Name Type Description
system System

The model, with the levels cells.level.

Examples:

>>> from skcomplex.models.chemistry.kinetics import cyclic_cellular_automaton
>>> cyclic_cellular_automaton((10, 10), n=8).elements[0].values['n']
8
Source code in src/skcomplex/models/chemistry/kinetics/_cyclic.py 
 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
@register_model(
    field='chemistry/kinetics',
    tags=('lattice', 'excitable media', 'spirals', 'cellular automaton'),
    references=('Griffeath 1988', 'Fisch, Gravner, and Griffeath 1991', 'Dilão 2023'),
    defaults={'until': 200, 'dt': 1},
)
def cyclic_cellular_automaton(
    shape: Sequence[int] = (64, 64),
    *,
    n: int = 14,
    threshold: int = 1,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Griffeath's cyclic cellular automaton, whose random levels organize into spiral waves.

    Args:
        shape:
            The number of rows and columns.

        n:
            The number of levels.

        threshold:
            The number of von Neumann neighbours at the following level that advance a cell.

        random_state:
            The seed or generator of the initial levels.

    Returns:
        system:
            The model, with the levels `cells.level`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import cyclic_cellular_automaton
        >>> cyclic_cellular_automaton((10, 10), n=8).elements[0].values['n']
        8
    """
    levels = np.random.default_rng(random_state).integers(0, n, size=tuple(shape)).astype(float)
    space = LatticeSpace(tuple(shape))
    return System(
        elements=[
            LatticeElement(
                'cells',
                element=CyclicCellElement,
                space=space,
                values={'level': levels, 'n': n, 'threshold': threshold},
            ),
        ],
        interactions=[
            LocalRuleInteraction(
                'advance',
                target='cells',
                state='level',
                rule=_advance_cyclically,
                neighborhood=VonNeumannNeighborhood(center=True),
            ),
        ],
    )

hodgepodge(shape=(100, 100), *, n=100.0, k1=2.0, k2=3.0, g=34.0, random_state=0)

Return the hodgepodge machine, a model of a chemical reaction whose waves form spirals.

Parameters:

Name Type Description Default
shape Sequence[int]

The number of rows and columns.

 (100, 100)
n float

The level of illness.

 100.0
k1 float

The weight of infected neighbours for healthy cells.

 2.0
k2 float

The weight of ill neighbours for healthy cells.

 3.0
g float

The speed of the infection.

 34.0
random_state int | Generator | None

The seed or generator of the initial levels.

 0

Returns:

Name Type Description
system System

The model, with the levels cells.level.

Examples:

>>> from skcomplex.models.chemistry.kinetics import hodgepodge
>>> hodgepodge((10, 10)).elements[0].values['g']
34.0
Source code in src/skcomplex/models/chemistry/kinetics/_hodgepodge.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
@register_model(
    field='chemistry/kinetics',
    tags=('lattice', 'excitable media', 'spirals'),
    references=('Gerhardt and Schuster 1989', 'Charbonneau 2017'),
    defaults={'until': 200, 'dt': 1},
)
def hodgepodge(
    shape: Sequence[int] = (100, 100),
    *,
    n: float = 100.0,
    k1: float = 2.0,
    k2: float = 3.0,
    g: float = 34.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the hodgepodge machine, a model of a chemical reaction whose waves form spirals.

    Args:
        shape:
            The number of rows and columns.

        n:
            The level of illness.

        k1:
            The weight of infected neighbours for healthy cells.

        k2:
            The weight of ill neighbours for healthy cells.

        g:
            The speed of the infection.

        random_state:
            The seed or generator of the initial levels.

    Returns:
        system:
            The model, with the levels `cells.level`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import hodgepodge
        >>> hodgepodge((10, 10)).elements[0].values['g']
        34.0
    """
    levels = np.random.default_rng(random_state).integers(0, int(n) + 1, size=tuple(shape)).astype(float)
    space = LatticeSpace(tuple(shape))
    return System(
        elements=[
            LatticeElement(
                'cells',
                element=HodgepodgeCellElement,
                space=space,
                values={'level': levels, 'n': n, 'k1': k1, 'k2': k2, 'g': g},
            ),
        ],
        interactions=[
            LocalRuleInteraction(
                'infection',
                target='cells',
                state='level',
                rule=_infect,
                neighborhood=MooreNeighborhood(center=True),
            ),
        ],
    )

michaelis_menten(enzyme=100.0, substrate=300.0, *, k_binding=0.01, k_unbinding=0.1, k_catalysis=0.1, stochastic=False)

Return the enzyme reaction E + S <-> C -> E + P of Michaelis and Menten, as rate equations or as jumps.

The total enzyme E + C and the total substrate S + C + P are conserved. In stochastic mode the amounts are numbers of molecules and the reactions are the random events of Gillespie's algorithm.

Parameters:

Name Type Description Default
enzyme float

The initial amount of free enzyme.

 100.0
substrate float

The initial amount of substrate.

 300.0
k_binding float

The rate constant of the binding of enzyme and substrate.

 0.01
k_unbinding float

The rate constant of the unbinding of the complex.

 0.1
k_catalysis float

The rate constant of the conversion of the complex into enzyme and product.

 0.1
stochastic bool

Whether the reactions are random events.

 False

Returns:

Name Type Description
system System

The model, with the amounts E.n, S.n, C.n, and P.n.

Examples:

>>> from skcomplex.models.chemistry.kinetics import michaelis_menten
>>> [element.label for element in michaelis_menten().elements]
['E', 'S', 'C', 'P']
Source code in src/skcomplex/models/chemistry/kinetics/_reactions.py 
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
@register_model(
    field='chemistry/kinetics',
    tags=('enzyme kinetics', 'ode', 'jump process', 'gillespie'),
    references=('Michaelis and Menten 1913', 'Gillespie 1977', 'Trzaska 2023'),
    defaults={'until': 50.0, 'dt': 0.5},
)
def michaelis_menten(
    enzyme: float = 100.0,
    substrate: float = 300.0,
    *,
    k_binding: float = 0.01,
    k_unbinding: float = 0.1,
    k_catalysis: float = 0.1,
    stochastic: bool = False,
) -> System:
    """Return the enzyme reaction `E + S <-> C -> E + P` of Michaelis and Menten, as rate equations or as jumps.

    The total enzyme `E + C` and the total substrate `S + C + P` are conserved. In stochastic mode the amounts are
    numbers of molecules and the reactions are the random events of Gillespie's algorithm.

    Args:
        enzyme:
            The initial amount of free enzyme.

        substrate:
            The initial amount of substrate.

        k_binding:
            The rate constant of the binding of enzyme and substrate.

        k_unbinding:
            The rate constant of the unbinding of the complex.

        k_catalysis:
            The rate constant of the conversion of the complex into enzyme and product.

        stochastic:
            Whether the reactions are random events.

    Returns:
        system:
            The model, with the amounts `E.n`, `S.n`, `C.n`, and `P.n`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import michaelis_menten
        >>> [element.label for element in michaelis_menten().elements]
        ['E', 'S', 'C', 'P']
    """
    return System(
        elements=[
            CompartmentElement('E', n=enzyme),
            CompartmentElement('S', n=substrate),
            CompartmentElement('C'),
            CompartmentElement('P'),
        ],
        interactions=_get_enzyme_reactions(k_binding, k_unbinding, k_catalysis, 'jump' if stochastic else 'rate'),
    )

nanoparticle_corona(nanoparticle=18.0, *, radius=35.0, proteins=(600000.0, 15000.0, 8800.0), corona=(14.0, 4.0, 0.0), protein_radii=(5.0, 4.0, 8.3), association=(2400.0, 30000.0, 2000.0), dissociation=(0.002, 3e-05, 0.002))

Return the corona that albumin, high-density lipoprotein, and fibrinogen form on a nanoparticle in plasma.

Each protein i binds a free nanoparticle at the rate n_i alpha_i U P_i and leaves it at the rate delta_i B_i, where n_i = 4 (R + r_i)^2 / r_i^2 is the number of its sites on a nanoparticle of radius R. The book gives only the balance of the complexes, and the model takes the free nanoparticles and the free proteins from mass action, so U + sum B_i and every P_i + B_i are conserved. The parameters and initial concentrations are the book's, without the binding to a target. Concentrations are in nM, association rates in 1/(M s), and time is in seconds. Albumin, the most abundant, covers the nanoparticles first, and lipoprotein, which binds more tightly, replaces it within hours, which is the Vroman effect.

Parameters:

Name Type Description Default
nanoparticle float

The initial concentration of free nanoparticles.

 18.0
radius float

The radius of a nanoparticle, in nm.

 35.0
proteins Sequence[float]

The initial concentrations of free albumin, lipoprotein, and fibrinogen.

 (600000.0, 15000.0, 8800.0)
corona Sequence[float]

The initial concentrations of nanoparticles bound to each protein.

 (14.0, 4.0, 0.0)
protein_radii Sequence[float]

The radii of the proteins, in nm.

 (5.0, 4.0, 8.3)
association Sequence[float]

The association rates of the proteins.

 (2400.0, 30000.0, 2000.0)
dissociation Sequence[float]

The dissociation rates of the proteins, per second.

 (0.002, 3e-05, 0.002)

Returns:

Name Type Description
system System

The model, with the free nanoparticles nanoparticle.n, the free proteins albumin.n, lipoprotein.n, and fibrinogen.n, and the complexes albumin_corona.n, lipoprotein_corona.n, and fibrinogen_corona.n.

Examples:

>>> from skcomplex.models.chemistry.kinetics import nanoparticle_corona
>>> round(nanoparticle_corona().interactions[0].k, 7)
0.0006144
Source code in src/skcomplex/models/chemistry/kinetics/_nanoparticles.py 
 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
@register_model(
    field='chemistry/kinetics',
    tags=('nanomedicine', 'protein corona', 'mass action', 'ode', 'Vroman effect'),
    references=('Padmanabhan, Meskin, and Al Moustafa 2021', 'Dell\'Orco et al. 2010'),
    defaults={'until': 10800.0, 'dt': 60.0, 'solver': ScipyODESolver(method='LSODA', rtol=1e-8, atol=1e-12)},
)
def nanoparticle_corona(
    nanoparticle: float = 18.0,
    *,
    radius: float = 35.0,
    proteins: Sequence[float] = (6e5, 1.5e4, 8.8e3),
    corona: Sequence[float] = (14.0, 4.0, 0.0),
    protein_radii: Sequence[float] = (5.0, 4.0, 8.3),
    association: Sequence[float] = (2.4e3, 3e4, 2e3),
    dissociation: Sequence[float] = (2e-3, 3e-5, 2e-3),
) -> System:
    """Return the corona that albumin, high-density lipoprotein, and fibrinogen form on a nanoparticle in plasma.

    Each protein `i` binds a free nanoparticle at the rate `n_i alpha_i U P_i` and leaves it at the rate
    `delta_i B_i`, where `n_i = 4 (R + r_i)^2 / r_i^2` is the number of its sites on a nanoparticle of radius `R`.
    The book gives only the balance of the complexes, and the model takes the free nanoparticles and the free proteins
    from mass action, so `U + sum B_i` and every `P_i + B_i` are conserved. The parameters and initial
    concentrations are the book's, without the binding to a target. Concentrations are in nM, association rates in
    1/(M s), and time is in seconds. Albumin, the most abundant, covers the nanoparticles first, and lipoprotein,
    which binds more tightly, replaces it within hours, which is the Vroman effect.

    Args:
        nanoparticle:
            The initial concentration of free nanoparticles.

        radius:
            The radius of a nanoparticle, in nm.

        proteins:
            The initial concentrations of free albumin, lipoprotein, and fibrinogen.

        corona:
            The initial concentrations of nanoparticles bound to each protein.

        protein_radii:
            The radii of the proteins, in nm.

        association:
            The association rates of the proteins.

        dissociation:
            The dissociation rates of the proteins, per second.

    Returns:
        system:
            The model, with the free nanoparticles `nanoparticle.n`, the free proteins `albumin.n`,
            `lipoprotein.n`, and `fibrinogen.n`, and the complexes `albumin_corona.n`, `lipoprotein_corona.n`, and
            `fibrinogen_corona.n`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import nanoparticle_corona
        >>> round(nanoparticle_corona().interactions[0].k, 7)
        0.0006144
    """
    elements = [CompartmentElement('nanoparticle', n=nanoparticle)]
    interactions: list[ReactionInteraction] = []
    rows = zip(PROTEINS, proteins, corona, protein_radii, association, dissociation, strict=True)
    for name, free, bound, protein_radius, alpha, delta in rows:
        complex_ = f'{name}_corona'
        elements += [CompartmentElement(name, n=free), CompartmentElement(complex_, n=bound)]
        k = _count_binding_sites(radius, protein_radius) * alpha * LITRES_PER_NANOMOLE
        interactions += [
            ReactionInteraction(
                f'{name}_binding',
                reactants={('nanoparticle', 'n'): 1, (name, 'n'): 1},
                products={(complex_, 'n'): 1},
                k=k,
            ),
            ReactionInteraction(
                f'{name}_release',
                reactants={(complex_, 'n'): 1},
                products={('nanoparticle', 'n'): 1, (name, 'n'): 1},
                k=delta,
            ),
        ]
    return System(elements=elements, interactions=interactions)

oregonator(x=0.1, y=0.1, z=0.1, *, epsilon=0.04, delta=0.0004, q=0.0008, f=1.0)

Return the scaled Oregonator, the reduction of the FKN mechanism of the Belousov-Zhabotinsky reaction.

The scaled concentrations of bromous acid x, bromide y, and the oxidized catalyst z follow epsilon dx/dt = q y - x y + x (1 - x), delta dy/dt = -q y - x y + f z, and dz/dt = x - z. The five reactions of the mechanism, bromate with bromide Y -> X at q y, bromous acid with bromide X + Y -> 0 at x y, the autocatalysis X -> 2X + Z at x, the disproportionation 2X -> 0 at x^2, and the reduction of the catalyst Z -> f Y at z, change x and y by stoichiometric coefficients scaled by 1 / epsilon and 1 / delta. With the standard parameters the reaction settles on a relaxation limit cycle, and the small epsilon and delta make it stiff.

Parameters:

Name Type Description Default
x float

The initial scaled concentration of bromous acid.

 0.1
y float

The initial scaled concentration of bromide.

 0.1
z float

The initial scaled concentration of the oxidized catalyst.

 0.1
epsilon float

The ratio of the time scales of bromous acid and the catalyst.

 0.04
delta float

The ratio of the time scales of bromide and the catalyst.

 0.0004
q float

The ratio of rate constants that sets the lowest bromous acid concentration.

 0.0008
f float

The stoichiometric factor of bromide production from the catalyst.

 1.0

Returns:

Name Type Description
system System

The model, with the concentrations reactor.x, reactor.y, and reactor.z.

Examples:

>>> from skcomplex.models.chemistry.kinetics import oregonator
>>> oregonator(f=0.8).elements[0].f
0.8
Source code in src/skcomplex/models/chemistry/kinetics/_oscillators.py 
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
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
@register_model(
    field='chemistry/kinetics',
    tags=('ode', 'oscillation', 'relaxation oscillation', 'stiff', 'Belousov-Zhabotinsky'),
    references=('Field and Noyes 1974', 'Field, Körös, and Noyes 1972', 'Tyson 1985', 'Trzaska 2023'),
    defaults={'until': 30.0, 'dt': 0.05, 'solver': ScipyODESolver(method='LSODA', rtol=1e-6, atol=1e-9)},
)
def oregonator(
    x: float = 0.1,
    y: float = 0.1,
    z: float = 0.1,
    *,
    epsilon: float = 0.04,
    delta: float = 4e-4,
    q: float = 8e-4,
    f: float = 1.0,
) -> System:
    """Return the scaled Oregonator, the reduction of the FKN mechanism of the Belousov-Zhabotinsky reaction.

    The scaled concentrations of bromous acid `x`, bromide `y`, and the oxidized catalyst `z` follow
    `epsilon dx/dt = q y - x y + x (1 - x)`, `delta dy/dt = -q y - x y + f z`, and `dz/dt = x - z`. The five reactions
    of the mechanism, bromate with bromide `Y -> X` at `q y`, bromous acid with bromide `X + Y -> 0` at `x y`, the
    autocatalysis `X -> 2X + Z` at `x`, the disproportionation `2X -> 0` at `x^2`, and the reduction of the catalyst
    `Z -> f Y` at `z`, change `x` and `y` by stoichiometric coefficients scaled by `1 / epsilon` and `1 / delta`. With
    the standard parameters the reaction settles on a relaxation limit cycle, and the small `epsilon` and `delta` make
    it stiff.

    Args:
        x:
            The initial scaled concentration of bromous acid.

        y:
            The initial scaled concentration of bromide.

        z:
            The initial scaled concentration of the oxidized catalyst.

        epsilon:
            The ratio of the time scales of bromous acid and the catalyst.

        delta:
            The ratio of the time scales of bromide and the catalyst.

        q:
            The ratio of rate constants that sets the lowest bromous acid concentration.

        f:
            The stoichiometric factor of bromide production from the catalyst.

    Returns:
        system:
            The model, with the concentrations `reactor.x`, `reactor.y`, and `reactor.z`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import oregonator
        >>> oregonator(f=0.8).elements[0].f
        0.8
    """
    return System(
        elements=[OregonatorElement('reactor', x=x, y=y, z=z, epsilon=epsilon, delta=delta, q=q, f=f)],
        interactions=[
            ReactionInteraction(
                'bromate_bromide',
                target='reactor',
                reactants={'y': _get_y_scale},
                products={'x': _get_x_scale},
                propensity=_get_bromate_bromide,
            ),
            ReactionInteraction(
                'bromous_bromide',
                target='reactor',
                reactants={'x': _get_x_scale, 'y': _get_y_scale},
                propensity=_get_bromous_bromide,
            ),
            ReactionInteraction(
                'autocatalysis',
                target='reactor',
                products={'x': _get_x_scale, 'z': 1},
                modifiers={'x': 1},
            ),
            ReactionInteraction(
                'disproportionation',
                target='reactor',
                reactants={'x': _get_x_scale},
                propensity=_get_disproportionation,
            ),
            ReactionInteraction(
                'catalyst_reduction',
                target='reactor',
                reactants={'z': 1},
                products={'y': _get_y_yield},
            ),
        ],
    )

pulsed_electrochemical_reactor(eta=0.0, c=1.0, *, C=1.0, alpha=0.5, tau=10.0, kappa=0.05, current=5.0, period=10.0, width=2.0)

Return a stirred electrochemical reactor whose electrode is driven by a train of current pulses.

Every pulse charges the double layer, the overpotential drives a Butler-Volmer current that consumes the reactant, and the feed restores it between pulses.

Parameters:

Name Type Description Default
eta float

The initial overpotential.

 0.0
c float

The initial reactant concentration, relative to the feed.

 1.0
C float

The capacitance of the double layer.

 1.0
alpha float

The transfer coefficient.

 0.5
tau float

The residence time of the reactor.

 10.0
kappa float

The consumption of reactant per unit of Faradaic current.

 0.05
current float

The current of a pulse.

 5.0
period float

The time between the starts of pulses.

 10.0
width float

The duration of a pulse.

 2.0

Returns:

Name Type Description
system System

The model, with the overpotential electrode.eta and the concentration electrode.c.

Examples:

>>> from skcomplex.models.chemistry.kinetics import pulsed_electrochemical_reactor
>>> pulsed_electrochemical_reactor(width=1.0).elements[0].width
1.0
Source code in src/skcomplex/models/chemistry/kinetics/_reactors.py 
 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
139
140
141
142
143
144
@register_model(
    field='chemistry/kinetics',
    tags=('ode', 'electrochemistry', 'pulses', 'driven'),
    references=('Trzaska 2023', 'Bard and Faulkner 2001'),
    defaults={'until': 60.0, 'dt': 0.1},
)
def pulsed_electrochemical_reactor(
    eta: float = 0.0,
    c: float = 1.0,
    *,
    C: float = 1.0,
    alpha: float = 0.5,
    tau: float = 10.0,
    kappa: float = 0.05,
    current: float = 5.0,
    period: float = 10.0,
    width: float = 2.0,
) -> System:
    """Return a stirred electrochemical reactor whose electrode is driven by a train of current pulses.

    Every pulse charges the double layer, the overpotential drives a Butler-Volmer current that consumes the reactant,
    and the feed restores it between pulses.

    Args:
        eta:
            The initial overpotential.

        c:
            The initial reactant concentration, relative to the feed.

        C:
            The capacitance of the double layer.

        alpha:
            The transfer coefficient.

        tau:
            The residence time of the reactor.

        kappa:
            The consumption of reactant per unit of Faradaic current.

        current:
            The current of a pulse.

        period:
            The time between the starts of pulses.

        width:
            The duration of a pulse.

    Returns:
        system:
            The model, with the overpotential `electrode.eta` and the concentration `electrode.c`.

    Examples:
        >>> from skcomplex.models.chemistry.kinetics import pulsed_electrochemical_reactor
        >>> pulsed_electrochemical_reactor(width=1.0).elements[0].width
        1.0
    """
    electrode = PulsedElectrodeElement(
        'electrode',
        eta=eta,
        c=c,
        C=C,
        alpha=alpha,
        tau=tau,
        kappa=kappa,
        current=current,
        period=period,
        width=width,
    )
    return System(elements=[electrode])