Three-node phosphorylation oscillator
Three proteins whose phosphorylation cycles regulate each other in a loop that oscillates.
- Name:
phosphorylation_oscillator - Field:
biology/systems_biology - Formalism: ordinary differential equations
- References: Raman 2021; Mogilner, Wollman, and Marshall 2006
Factory
@register_model(
field='biology/systems_biology',
tags=('oscillator', 'phosphorylation', 'Michaelis-Menten', 'negative feedback', 'ode'),
references=(RAMAN, 'Mogilner, Wollman, and Marshall 2006'),
defaults={'until': 40.0, 'dt': 0.05},
)
def phosphorylation_oscillator(
A: float = 1.0,
B: float = 0.0,
C: float = 0.0,
*,
kp1: float = 1.0,
kd1: float = 2.0,
kp2: float = 2.0,
kd2: float = 1.0,
kp3: float = 2.0,
kd3: float = 1.0,
K: float = 0.01,
) -> System:
"""Return three proteins whose phosphorylation cycles regulate each other in a loop that oscillates.
C converts A, A converts B, and B converts C, each reaction with Michaelis-Menten kinetics and each total
conserved, so the loop is a negative feedback of three steps. The book prints the equations and no parameter
values. The defaults are a choice: totals of 1, Michaelis constants `K` far below the totals, which make every
cycle an ultrasensitive switch, and rates that put every switch at half the total of its enzyme. The proteins then
switch on and off in turn on a limit cycle of period about 4.2. With large `K` the switches become graded and the
loop settles on a fixed point. The amounts are fractions of the totals.
Args:
A:
The initial amount of unphosphorylated A.
B:
The initial amount of unphosphorylated B.
C:
The initial amount of unphosphorylated C.
kp1:
The maximal rate that restores A.
kd1:
The catalytic rate of the conversion of A by C.
kp2:
The catalytic rate of the conversion of B by A.
kd2:
The maximal rate that removes B.
kp3:
The catalytic rate of the conversion of C by B.
kd3:
The maximal rate that removes C.
K:
The Michaelis constant of every reaction.
Returns:
system:
The model, with the amounts `loop.A`, `loop.B`, and `loop.C`.
Examples:
>>> from skcomplex.models.biology.systems_biology import phosphorylation_oscillator
>>> phosphorylation_oscillator(K=0.1).elements[0].K2
0.1
"""
loop = PhosphorylationLoopElement(
'loop',
A=A,
B=B,
C=C,
kp1=kp1,
kd1=kd1,
kp2=kp2,
kd2=kd2,
kp3=kp3,
kd3=kd3,
K1=K,
K2=K,
K3=K,
)
reactions = [
MichaelisMentenInteraction(
'a_restoration',
target='loop',
products={'A': 1},
substrate=lambda s: s.target.A_T - s.target.A,
vmax=lambda s: s.target.kp1,
km=lambda s: s.target.K1,
),
MichaelisMentenInteraction(
'a_conversion',
target='loop',
reactants={'A': 1},
substrate='A',
vmax=lambda s: s.target.kd1,
km=lambda s: s.target.K1,
modifiers={'C': 1},
),
MichaelisMentenInteraction(
'b_conversion',
target='loop',
products={'B': 1},
substrate=lambda s: s.target.B_T - s.target.B,
vmax=lambda s: s.target.kp2,
km=lambda s: s.target.K2,
modifiers={'A': 1},
),
MichaelisMentenInteraction(
'b_removal',
target='loop',
reactants={'B': 1},
substrate='B',
vmax=lambda s: s.target.kd2,
km=lambda s: s.target.K2,
),
MichaelisMentenInteraction(
'c_conversion',
target='loop',
products={'C': 1},
substrate=lambda s: s.target.C_T - s.target.C,
vmax=lambda s: s.target.kp3,
km=lambda s: s.target.K3,
modifiers={'B': 1},
),
MichaelisMentenInteraction(
'c_removal',
target='loop',
reactants={'C': 1},
substrate='C',
vmax=lambda s: s.target.kd3,
km=lambda s: s.target.K3,
),
]
return System(elements=[loop], interactions=reactions)
Components
class PhosphorylationLoopElement(Element):
"""Three proteins, each phosphorylated and dephosphorylated by Michaelis-Menten kinetics, that regulate each other.
The states are the unphosphorylated amounts, and each total is conserved. With `u_T = A_T - A`, and so on:
- `dA/dt = kp1 u_T / (K1 + u_T) - kd1 C A / (K1 + A)`,
- `dB/dt = kp2 A (B_T - B) / (K2 + B_T - B) - kd2 B / (K2 + B)`,
- `dC/dt = kp3 B (C_T - C) / (K3 + C_T - C) - kd3 C / (K3 + C)`.
The element holds the states and parameters and has no laws of its own: the interactions of
`phosphorylation_oscillator` read its parameters.
Args:
label:
The label of the loop.
visual:
How the element is drawn, or `None` for `default_visual_`.
A:
The initial amount of unphosphorylated A.
B:
The initial amount of unphosphorylated B.
C:
The initial amount of unphosphorylated C.
kp1:
The maximal rate that restores A.
kd1:
The catalytic rate of the conversion of A by C.
kp2:
The catalytic rate of the conversion of B by A.
kd2:
The maximal rate that removes B.
kp3:
The catalytic rate of the conversion of C by B.
kd3:
The maximal rate that removes C.
K1:
The Michaelis constant of both reactions of A.
K2:
The Michaelis constant of both reactions of B.
K3:
The Michaelis constant of both reactions of C.
A_T:
The total amount of A.
B_T:
The total amount of B.
C_T:
The total amount of C.
Examples:
>>> from skcomplex.models.biology.systems_biology import PhosphorylationLoopElement
>>> PhosphorylationLoopElement('loop', A=0.5).A
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
A: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
B: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
C: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
kp1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
kd1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
kp2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
kd2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
kp3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
kd3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
K1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
K2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
K3: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
A_T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
B_T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
C_T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.A = A
self.B = B
self.C = C
self.kp1 = kp1
self.kd1 = kd1
self.kp2 = kp2
self.kd2 = kd2
self.kp3 = kp3
self.kd3 = kd3
self.K1 = K1
self.K2 = K2
self.K3 = K3
self.A_T = A_T
self.B_T = B_T
self.C_T = C_T
class MichaelisMentenInteraction(ReactionInteraction):
"""Reaction at the Michaelis-Menten propensity `vmax S / (km + S)`, times every modifier to its coefficient.
The substrate `S` is a state of the target, or a function of the namespace that returns an amount that is no state
of its own, such as the phosphorylated form `A_T - A` of a protein whose total is conserved. The modifiers are the
enzymes that catalyse the reaction, which scale the propensity but do not change, so `vmax` is a catalytic rate when
the reaction has modifiers and a maximal rate when it has none. The maximal rate and the Michaelis constant are
functions of the namespace, which reads the parameters of the target, or numbers.
Args:
label:
The label of the reaction, unique among the parts of its system.
target:
The path of the element, group, or field whose states the species named by a state alone are.
reactants:
The coefficients of the reactants, or `None` for a source.
products:
The coefficients of the products, or `None` for a sink.
substrate:
The name of the state of the target that saturates the propensity, or a function of the namespace that
returns the amount that does.
vmax:
A function of the namespace that returns the maximal or catalytic rate, or a number.
km:
A function of the namespace that returns the Michaelis constant, or a number.
modifiers:
The coefficients of the modifiers, the enzymes that scale the propensity, or `None`.
mode:
The form of the laws, `'rate'`, `'jump'`, or `'langevin'`.
Raises:
ValidationError:
The substrate is neither a name nor a function, a rate is neither a function nor a number, or the reaction
cannot take its species or mode.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.common.kinetics import MichaelisMentenInteraction
>>> uptake = MichaelisMentenInteraction(
... 'uptake', target='cell', reactants={'n': 1}, substrate='n', vmax=1.0, km=100.0
... )
>>> system = System(elements=[CompartmentElement('cell', n=500.0)], interactions=[uptake])
>>> trajectory = Simulation(system, dt=1.0).check().prepare().run(until=10.0).trajectory_
>>> round(float(trajectory.values('cell.n')[-1]), 3)
491.679
"""
def __init__(
self,
label: str,
*,
target: Annotated[str | None, RoleSpec()] = None,
reactants: Annotated[Stoichiometry | None, RoleSpec(variadic=True)] = None,
products: Annotated[Stoichiometry | None, RoleSpec(variadic=True)] = None,
substrate: str | Callable[[SimpleNamespace], object] = 'n',
vmax: Callable[[SimpleNamespace], object] | float = 1.0,
km: Callable[[SimpleNamespace], object] | float = 1.0,
modifiers: Annotated[Stoichiometry | None, RoleSpec(variadic=True)] = None,
mode: str = 'rate',
) -> None:
super().__init__(label, target=target, reactants=reactants, products=products, modifiers=modifiers, mode=mode)
self.substrate = substrate
self.vmax = vmax
self.km = km
def _check(self) -> None:
"""Raise an error for a substrate or rates the reaction cannot take.
Raises:
ValidationError:
The substrate is neither a name nor a function, a rate is neither a function nor a number, or the
reaction cannot take its species or mode.
"""
super()._check()
if not (callable(self.substrate) or (isinstance(self.substrate, str) and self.substrate.isidentifier())):
error_msg = (
f'The substrate of {self.label!r} is {self.substrate!r}. Give the name of a state of the target or a '
'function of the namespace.'
)
raise ValidationError(error_msg)
for name in ('vmax', 'km'):
value = getattr(self, name)
if not (callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
def _evaluate(self, value: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return value(namespace) if callable(value) else value
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the saturating propensity of the substrate, times every modifier to its coefficient.
Args:
namespace:
The namespace of the reaction, with the namespace of the target.
Returns:
propensity:
The propensity of the reaction.
"""
if callable(self.substrate):
substrate = SR(self.substrate(namespace))
else:
substrate = SR(getattr(namespace.target, self.substrate))
vmax, km = SR(self._evaluate(self.vmax, namespace)), SR(self._evaluate(self.km, namespace))
propensity = vmax * substrate / (km + substrate)
for state, coefficient in self._get_terms(self._get_modifiers(), namespace, 'modifiers'):
propensity *= state**coefficient
return propensity
Equations
loop.A: diff(A_loop(t), t) == -kd1_loop*A_loop(t)*C_loop(t)/(K1_loop + A_loop(t)) + (A_T_loop - A_loop(t))*kp1_loop/(A_T_loop + K1_loop - A_loop(t))
loop.B: diff(B_loop(t), t) == (B_T_loop - B_loop(t))*kp2_loop*A_loop(t)/(B_T_loop + K2_loop - B_loop(t)) - kd2_loop*B_loop(t)/(K2_loop + B_loop(t))
loop.C: diff(C_loop(t), t) == (C_T_loop - C_loop(t))*kp3_loop*B_loop(t)/(C_T_loop + K3_loop - C_loop(t)) - kd3_loop*C_loop(t)/(K3_loop + C_loop(t))