Tumour-immune model with dendritic licensing, IL-2, IL-10, and TGF-beta
A tumour against dendritic, cytotoxic, helper, and regulatory T cells with IL-2, IL-10, and TGF-beta.
- Name:
tumour_immune_suppression - Field:
biology/oncology - Formalism: ordinary differential equations
- References: Padmanabhan, Meskin, and Al Moustafa 2021; Robertson-Tessi, El-Kareh, and Goriely 2012
Factory
@register_model(
field='biology/oncology',
tags=('immunotherapy', 'ode', 'immunosuppression', 'cytokines', 'antigenicity'),
references=(BOOK_REFERENCE, 'Robertson-Tessi, El-Kareh, and Goriely 2012'),
defaults={'until': 200.0, 'dt': 0.5, 'solver': STIFF_SOLVER},
)
def tumour_immune_suppression(
*,
alpha_a: float = 1e-3,
r1: float = 0.333,
r2: float = 333.0,
EMC: float = 1e4,
EMH: float = 1e4,
EMR: float = 1e4,
) -> System:
"""Return a tumour against dendritic, cytotoxic, helper, and regulatory T cells with IL-2, IL-10, and TGF-beta.
A tumour of one cell grows exponentially, then as a power of its size, and only its accessible fraction is killed.
Its antigen matures dendritic cells, which recruit T cells from the memory pools, and the regulatory cells with
IL-10 and TGF-beta suppress the response. A weakly antigenic tumour escapes, a moderately antigenic one oscillates,
and a strongly antigenic one is held to a few cells. The book's growth rows are garbled, and the exponential rate
is read as 0.333 per day and the power coefficient as 333 cells^(1/2) per day, within the book's range 100 to
1000, so the growth turns from exponential to power at `10^6` cells. The book gives no memory pools, which are
taken as constant at `10^4` cells each. The other parameters are the book's, the initial immune cells are
`10^-5`, cells are counts, cytokines in ng/ml, and time in days.
Args:
alpha_a:
The antigenicity of the tumour, in the book's range `10^-5` to `10^5`.
r1:
The exponential growth rate of the tumour.
r2:
The power growth coefficient of the tumour.
EMC:
The memory cytotoxic T cells.
EMH:
The memory helper T cells.
EMR:
The memory regulatory T cells.
Returns:
system:
The model, with the tumour `site.A`, the dendritic cells `site.EUD` and `site.ED`, the T cells `site.EPC`,
`site.E`, `site.EPH`, `site.EH`, `site.EPR`, and `site.ER`, and the cytokines `site.C1`, `site.C2`, and
`site.C3`.
Examples:
>>> from skcomplex.models.biology.oncology import tumour_immune_suppression
>>> tumour_immune_suppression(alpha_a=1.0).elements[0].alpha_a
1.0
"""
site = SuppressedTumourElement('site', alpha_a=alpha_a, r1=r1, r2=r2, EMC=EMC, EMH=EMH, EMR=EMR)
def accessible(s: SimpleNamespace) -> Expression:
"""Return the fraction of the tumour accessible to the cytotoxic T cells."""
fraction: Expression = s.target.A / (1 + (s.target.A ** (1 - s.target.q) / s.target.g3) ** s.target.p) ** (
1 / s.target.p
)
return fraction
def maturation(s: SimpleNamespace) -> Expression:
"""Return the maturation by IL-2 of proliferating T cells, which TGF-beta suppresses."""
factor: Expression = hill(s.target.C1, s.target.EC1_50) / (1 + s.target.C3 / s.target.EC3_50)
return factor
interactions = [
RateInteraction(
'tumour_growth',
target='site',
state='A',
rate=lambda s: s.target.A
/ ((1 / s.target.r1) ** s.target.p + (s.target.A ** (1 - s.target.n) / s.target.r2) ** s.target.p)
** (1 / s.target.p),
),
ReactionInteraction(
'kill',
target='site',
reactants={'A': 1},
modifiers={'E': 1},
propensity=lambda s: s.target.a1
* accessible(s)
* s.target.E
/ (s.target.E + accessible(s) / s.target.g1 + s.target.ER / s.target.g2)
/ (1 + s.target.C3 / s.target.C3_50),
),
RateInteraction(
'antigen_uptake',
target='site',
state='EUD',
rate=lambda s: s.target.alpha_a
* accessible(s)
/ ((1 + s.target.C2 / s.target.C2_50) * (1 + s.target.ER / s.target.ER_50)),
),
ReactionInteraction(
'licensing',
target='site',
reactants={'EUD': 1},
products={'ED': 1},
propensity=lambda s: s.target.lambda_a * s.target.EUD / (1 + s.target.EUD / s.target.EMH),
),
ReactionInteraction(
'immature_death',
target='site',
reactants={'EUD': 1},
propensity=lambda s: s.target.d1 * s.target.EUD,
),
ReactionInteraction(
'licensed_death',
target='site',
reactants={'ED': 1},
propensity=lambda s: s.target.d2 * s.target.ED,
),
RateInteraction(
'cytotoxic_priming',
target='site',
state='EPC',
rate=lambda s: s.target.a2
* s.target.EMC
* s.target.ED
/ (s.target.ED + s.target.g4 * (s.target.EMC + s.target.EMH + s.target.EMR)),
),
ReactionInteraction(
'cytotoxic_precursor_death',
target='site',
reactants={'EPC': 1},
propensity=lambda s: s.target.d3 * s.target.EPC,
),
RateInteraction(
'cytotoxic_maturation',
target='site',
state='E',
rate=lambda s: s.target.a3 * s.target.EPC * maturation(s),
),
ReactionInteraction(
'cytotoxic_death',
target='site',
reactants={'E': 1},
propensity=lambda s: s.target.d4 * s.target.E,
),
RateInteraction(
'helper_priming',
target='site',
state='EPH',
rate=lambda s: s.target.a4
* s.target.EMH
* (s.target.ED + s.target.EUD)
/ (s.target.ED + s.target.EUD + s.target.g4 * (s.target.EMC + s.target.EMH + s.target.EMR)),
),
ReactionInteraction(
'helper_precursor_death',
target='site',
reactants={'EPH': 1},
propensity=lambda s: s.target.d5 * s.target.EPH,
),
RateInteraction(
'helper_maturation',
target='site',
state='EH',
rate=lambda s: s.target.a5 * s.target.EPH * maturation(s),
),
ReactionInteraction(
'conversion',
target='site',
reactants={'EH': 1},
products={'ER': 1},
modifiers={'C3': 1},
propensity=lambda s: s.target.a6 * s.target.EH * s.target.C3 / (s.target.g5 + s.target.C3),
),
ReactionInteraction(
'helper_death',
target='site',
reactants={'EH': 1},
propensity=lambda s: s.target.d6 * s.target.EH,
),
RateInteraction(
'regulatory_priming',
target='site',
state='EPR',
rate=lambda s: s.target.a7
* s.target.EMR
* s.target.ED
/ (s.target.ED + s.target.g4 * (s.target.EMC + s.target.EMH + s.target.EMR)),
),
ReactionInteraction(
'regulatory_precursor_death',
target='site',
reactants={'EPR': 1},
propensity=lambda s: s.target.d7 * s.target.EPR,
),
RateInteraction(
'regulatory_maturation',
target='site',
state='ER',
rate=lambda s: s.target.a8 * s.target.EPR * hill(s.target.C1, s.target.EC1_50),
),
ReactionInteraction(
'regulatory_death',
target='site',
reactants={'ER': 1},
propensity=lambda s: s.target.d8 * s.target.ER,
),
RateInteraction(
'il2_secretion',
target='site',
state='C1',
rate=lambda s: s.target.a9
* s.target.EPH
/ ((1 + s.target.C3 / s.target.CC3_50) * (1 + s.target.C2 / s.target.CC2_50)),
),
ReactionInteraction(
'il2_loss',
target='site',
reactants={'C1': 1},
propensity=lambda s: s.target.C1 / s.target.tau1,
),
RateInteraction(
'regulatory_il10',
target='site',
state='C2',
rate=lambda s: s.target.r3 * s.target.ER,
),
RateInteraction(
'tumour_il10',
target='site',
state='C2',
rate=lambda s: s.target.r4 * s.target.A,
),
ReactionInteraction(
'il10_loss',
target='site',
reactants={'C2': 1},
propensity=lambda s: s.target.C2 / s.target.tau2,
),
RateInteraction(
'regulatory_tgf_beta',
target='site',
state='C3',
rate=lambda s: s.target.r5 * s.target.ER,
),
RateInteraction(
'tumour_tgf_beta',
target='site',
state='C3',
rate=lambda s: s.target.r6 * s.target.A,
),
ReactionInteraction(
'tgf_beta_loss',
target='site',
reactants={'C3': 1},
propensity=lambda s: s.target.C3 / s.target.tau3,
),
]
return System(elements=[site], interactions=interactions)
Components
class SuppressedTumourElement(Element):
"""Tumour, dendritic cells, proliferating and active T cells, IL-2, IL-10, and TGF-beta.
The equations are those of Padmanabhan et al. (4.22) to (4.33). The tumour grows by the smooth minimum
`A / ((1 / r1)^p + (A^(1 - n) / r2)^p)^(1/p)` of exponential and power growth and is killed through its fraction
`A* = A / (1 + (A^(1 - q) / g3)^p)^(1/p)` accessible to the cytotoxic T cells `E`, at the rate
`a1 A* E / (E + A* / g1 + ER / g2) / (1 + C3 / C3_50)`. Antigen matures unlicensed dendritic cells, which helper
memory cells license, and dendritic cells recruit the memory pools `EMC`, `EMH`, and `EMR` into proliferating
cytotoxic, helper, and regulatory cells, which IL-2 matures. IL-10 and TGF-beta, from regulatory and tumour cells,
suppress maturation, IL-2 production, and killing, and TGF-beta converts helper into regulatory cells.
The element holds the states and parameters and has no laws of its own: the interactions of
`tumour_immune_suppression` read its parameters.
Args:
label:
The label of the tumour site.
visual:
How the element is drawn, or `None` for `default_visual_`.
A:
The initial tumour cells.
EUD:
The initial unlicensed dendritic cells.
ED:
The initial licensed dendritic cells.
EPC:
The initial proliferating cytotoxic T cells.
E:
The initial active cytotoxic T cells.
EPH:
The initial proliferating helper T cells.
EH:
The initial active helper T cells.
EPR:
The initial proliferating regulatory T cells.
ER:
The initial active regulatory T cells.
C1:
The initial IL-2 concentration.
C2:
The initial IL-10 concentration.
C3:
The initial TGF-beta concentration.
r1:
The exponential growth rate of the tumour.
r2:
The power growth coefficient of the tumour.
n:
The exponent of the power growth.
p:
The smoothness of the transitions.
q:
The exponent of the accessible fraction of a large tumour.
a1:
The maximal kill rate of tumour cells.
g1:
The ratio of tumour to cytotoxic cells at half kill.
g2:
The ratio of regulatory to cytotoxic cells at half kill.
g3:
The size scale of the accessible fraction.
C3_50:
The TGF-beta at half kill.
alpha_a:
The antigenicity of the tumour.
C2_50:
The IL-10 at half dendritic maturation.
ER_50:
The regulatory cells at half dendritic maturation.
lambda_a:
The licensing rate of dendritic cells.
d1:
The death rate of unlicensed dendritic cells.
d2:
The death rate of licensed dendritic cells.
a2:
The recruitment rate of proliferating cytotoxic cells.
g4:
The steepness of the recruitment by dendritic cells.
d3:
The death rate of proliferating cytotoxic cells.
a3:
The maturation rate of cytotoxic cells.
EC1_50:
The IL-2 at half maturation.
EC3_50:
The TGF-beta at half maturation.
d4:
The death rate of active cytotoxic cells.
a4:
The recruitment rate of proliferating helper cells.
d5:
The death rate of proliferating helper cells.
a5:
The maturation rate of helper cells.
a6:
The conversion rate of helper into regulatory cells.
g5:
The TGF-beta at half conversion.
d6:
The death rate of active helper cells.
a7:
The recruitment rate of proliferating regulatory cells.
d7:
The death rate of proliferating regulatory cells.
a8:
The maturation rate of regulatory cells.
d8:
The death rate of active regulatory cells.
a9:
The IL-2 production by helper cells.
CC2_50:
The IL-10 at half suppression of IL-2.
CC3_50:
The TGF-beta at half suppression of IL-2.
tau1:
The lifetime of IL-2.
r3:
The IL-10 production by regulatory cells.
r4:
The IL-10 production by tumour cells.
tau2:
The lifetime of IL-10.
r5:
The TGF-beta production by regulatory cells.
r6:
The TGF-beta production by tumour cells.
tau3:
The lifetime of TGF-beta.
EMC:
The memory cytotoxic T cells.
EMH:
The memory helper T cells.
EMR:
The memory regulatory T cells.
Examples:
>>> from skcomplex.models.biology.oncology import SuppressedTumourElement
>>> SuppressedTumourElement('site').lambda_a
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
A: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 1.0,
EUD: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
ED: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
EPC: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
E: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
EPH: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
EH: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
EPR: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
ER: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e-05,
C1: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
C2: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
C3: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
r1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.333,
r2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 333.0,
n: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
p: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 3.0,
q: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2 / 3,
a1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.9,
g1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.2,
g2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 11.0,
g3: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
C3_50: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 3.5,
alpha_a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
C2_50: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.4,
ER_50: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 20000000.0,
lambda_a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
d1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.14,
d2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
a2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 23.0,
g4: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.33,
d3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
a3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 16.0,
EC1_50: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.3,
EC3_50: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.9,
d4: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
a4: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 9.9,
d5: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
a5: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.9,
a6: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.022,
g5: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.7,
d6: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
a7: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.1,
d7: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
a8: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.1,
d8: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
a9: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.7e-05,
CC2_50: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.75,
CC3_50: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.9,
tau1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.08,
r3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.4e-08,
r4: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.3e-10,
tau2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.05,
r5: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.8e-08,
r6: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.1e-07,
tau3: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.07,
EMC: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10000.0,
EMH: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10000.0,
EMR: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10000.0,
) -> None:
super().__init__(label, visual=visual)
self.A = A
self.EUD = EUD
self.ED = ED
self.EPC = EPC
self.E = E
self.EPH = EPH
self.EH = EH
self.EPR = EPR
self.ER = ER
self.C1 = C1
self.C2 = C2
self.C3 = C3
self.r1 = r1
self.r2 = r2
self.n = n
self.p = p
self.q = q
self.a1 = a1
self.g1 = g1
self.g2 = g2
self.g3 = g3
self.C3_50 = C3_50
self.alpha_a = alpha_a
self.C2_50 = C2_50
self.ER_50 = ER_50
self.lambda_a = lambda_a
self.d1 = d1
self.d2 = d2
self.a2 = a2
self.g4 = g4
self.d3 = d3
self.a3 = a3
self.EC1_50 = EC1_50
self.EC3_50 = EC3_50
self.d4 = d4
self.a4 = a4
self.d5 = d5
self.a5 = a5
self.a6 = a6
self.g5 = g5
self.d6 = d6
self.a7 = a7
self.d7 = d7
self.a8 = a8
self.d8 = d8
self.a9 = a9
self.CC2_50 = CC2_50
self.CC3_50 = CC3_50
self.tau1 = tau1
self.r3 = r3
self.r4 = r4
self.tau2 = tau2
self.r5 = r5
self.r6 = r6
self.tau3 = tau3
self.EMC = EMC
self.EMH = EMH
self.EMR = EMR
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or 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)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
class ReactionInteraction(JumpInteraction):
"""Reaction that moves units from its reactants to its products at a mass-action or a given propensity.
A reactant, a product, or a modifier is a pair of the path of an element, or of a member of a group such as
`cells[0]`, and the name of its state that holds the units, mapped to its stoichiometric coefficient. A coefficient
is a positive number, or a function of the namespace of the reaction that returns an expression of its parameters,
such as `lambda s: 1 / s.target.epsilon`. So `S + I -> 2 I` has the reactants `{('S', 'n'): 1, ('I', 'n'): 1}` and
the products `{('I', 'n'): 2}`, and two states of one element are `('cell', 'u')` and `('cell', 'v')`. A species may
also be the name of a state of the `target` alone, such as `{'u': 1, 'v': 2}`. A target that is a group or a field
reacts in every member, with the species of every member its own states, as the reactions at every site of a
reaction-diffusion field.
A subclass may bind its species under roles of its own names, such as `susceptible` and `infected`: it returns its
entries from `_get_reactants`, `_get_products`, and `_get_modifiers`, pairs of the path a role binds and a state
name, and the reaction reads every pair from the namespace of the role that binds its path.
The propensity is `k` times the product of every reactant and modifier state to the power of its coefficient, or the
expression the `propensity` function returns. The reaction changes every state by its coefficient among the products
minus its coefficient among the reactants, so it conserves the total of the states whose changes sum to zero, and
modifiers enter the propensity but do not change.
In the mode `'rate'` the reaction is a rate equation. In the mode `'jump'` it is the random event of a
`JumpInteraction` at the propensity that moves the units of the coefficients. In the mode `'langevin'` it is the
chemical Langevin equation, the rate equation with one Wiener process of amplitude the square root of the propensity
driving every state the reaction changes.
In the `'agent'` form the target is an agent group and every member reacts on its own states, the species named by
state names alone, in every step from the states at its start. The propensity of a member is `k` times the product
of its reactant and modifier states to the power of their coefficients, or the number `propensity(agent, model)`
returns, and a coefficient that is a function is the number `coefficient(agent, model)` returns. In the mode
`'rate'` the member adds the changes times its propensity, one unit step of the rate equation, and in the mode
`'jump'` it adds the changes with the probability `1 - exp(-propensity)` that an event fires in the step. The
members react at once from the states at the start of the rule under the schedule `'synchronous'`, one after the
other in their order or in a new random order every step under `'sequential'` or `'random-sequential'`, or each in
its own activation under `'member'`.
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, or `None`.
reactants:
The coefficients of the reactants, or `None` for a source.
products:
The coefficients of the products, or `None` for a sink.
k:
The rate constant.
propensity:
A function of the namespace of the reaction that returns the propensity, or `None` for mass action. The
namespace holds the tuples `reactants`, `products`, and `modifiers` of the namespaces of the elements their
pairs bind, in the order of their entries, `target`, the namespace of the target or of one of its members,
and the symbol `k`. In the agent form, a function `propensity(agent, model)` that returns a number.
modifiers:
The coefficients of the modifiers, whose states the propensity reads but the reaction does not change, or
`None`.
mode:
The form of the laws, `'rate'`, `'jump'`, or `'langevin'`.
form:
`'symbolic'` for the laws of the states, or `'agent'` for a rule in which every member of an agent group
reacts on its own states.
phase:
The phase of the step the agent form runs in, or `None` for the default phase.
schedule:
The order of the members in the agent form, `'synchronous'`, `'sequential'`, `'random-sequential'`, or
`'member'`.
Raises:
ValidationError:
The mode, the form, or the schedule is unknown, the symbolic form is not synchronous, the reaction has
neither reactants nor products, an entry is neither a pair of a path and a state name nor a state name of a
given target, a coefficient is neither a positive number nor a function, or the agent form has no target,
has an entry that is a pair, or has the mode `'langevin'`.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import ReactionInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> infection = ReactionInteraction(
... 'infection', reactants={('S', 'n'): 1, ('I', 'n'): 1}, products={('I', 'n'): 2}, k=0.0005
... )
>>> recovery = ReactionInteraction('recovery', reactants={('I', 'n'): 1}, products={('R', 'n'): 1}, k=0.1)
>>> system = System(
... elements=[CompartmentElement('S', n=990), CompartmentElement('I', n=10), CompartmentElement('R')],
... interactions=[infection, recovery],
... )
>>> states = Simulation(system, dt=1.0).check().prepare().run(until=50).trajectory_.states
>>> round(float(states[-1].sum()), 6)
1000.0
The decay `n -> 0` at every site of a field, at the rate `k n` of the site:
>>> import numpy as np
>>> from skcomplex.core.elements import FieldElement
>>> from skcomplex.core.spaces import LatticeSpace
>>> decay = ReactionInteraction('decay', target='cells', reactants={'n': 1}, k=1.0)
>>> values = {'n': np.arange(3.0)}
>>> cells = FieldElement('cells', element=CompartmentElement, space=LatticeSpace((3,)), values=values)
>>> simulation = Simulation(System(elements=[cells], interactions=[decay]), dt=1.0).check().prepare()
>>> np.round(simulation.run(until=1).trajectory_.lattice('cells.n')[-1], 3).tolist()
[0.0, 0.368, 0.736]
A coefficient that is a parameter, the decay at the rate `x` that removes `1 / epsilon` units per event:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.specs import ParameterSpec, StateSpec
>>> class ScaledElement(Element):
... def __init__(
... self,
... label,
... x: Annotated[float, StateSpec()] = 1.0,
... epsilon: Annotated[float, ParameterSpec()] = 0.5,
... ):
... super().__init__(label)
... self.x = x
... self.epsilon = epsilon
>>> scaled = ReactionInteraction(
... 'decay', target='box', reactants={'x': lambda s: 1 / s.target.epsilon}, propensity=lambda s: s.target.x
... )
>>> simulation = Simulation(System(elements=[ScaledElement('box')], interactions=[scaled]), dt=0.01)
>>> round(float(simulation.check().prepare().run(until=0.5).trajectory_.values('box.x')[-1]), 3)
0.368
Every agent of a group decays on its own states, by `k n` in a step of the rate equation:
>>> from skcomplex.core.elements import GroupElement
>>> from skcomplex.core.spaces import DiscreteTime
>>> class CellElement(Element):
... def __init__(self, label, n: Annotated[float, StateSpec()] = 3.0):
... super().__init__(label)
... self.n = n
>>> cells = GroupElement('cells', element=CellElement, size=2)
>>> decay = ReactionInteraction('decay', target='cells', reactants={'n': 1}, k=0.1, form='agent')
>>> system = System(elements=[cells], interactions=[decay], time=DiscreteTime())
>>> trajectory = Simulation(system).check().prepare().run(until=2).trajectory_
>>> np.round(trajectory.agents('cells', 2).values['n'], 6).tolist()
[2.43, 2.43]
"""
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,
k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
propensity: Callable[..., object] | None = None,
modifiers: Annotated[Stoichiometry | None, RoleSpec(variadic=True)] = None,
mode: str = 'rate',
form: str = 'symbolic',
phase: str | None = None,
schedule: str = 'synchronous',
) -> None:
super().__init__(label, target='')
self.target = target
self.reactants = reactants
self.products = products
self.k = k
self.propensity = propensity
self.modifiers = modifiers
self.mode = mode
self.form = form
self.phase = phase
self.schedule = schedule
def _get_reactants(self) -> Stoichiometry | None:
"""Return the coefficients of the reactants, which a subclass with roles of its own names returns."""
return self.reactants
def _get_products(self) -> Stoichiometry | None:
"""Return the coefficients of the products, which a subclass with roles of its own names returns."""
return self.products
def _get_modifiers(self) -> Stoichiometry | None:
"""Return the coefficients of the modifiers, which a subclass with roles of its own names returns."""
return self.modifiers
def _get_mode(self) -> str:
"""Return the form of the laws, which a subclass that sets it by another argument returns."""
return self.mode
def _get_bound(self, namespace: SimpleNamespace) -> dict[str, SimpleNamespace]:
"""Return the namespace of every path the plain roles of the reaction bind, by path."""
roles = [
name
for name, spec in self._get_specs().items()
if isinstance(spec, RoleSpec) and not spec.variadic and not spec.noise
]
return {
getattr(self, name): getattr(namespace, name)
for name in roles
if isinstance(getattr(self, name), str) and hasattr(namespace, name)
}
def _get_terms(
self,
entries: Stoichiometry | None,
namespace: SimpleNamespace,
name: str,
) -> list[tuple[Expression, Expression | float]]:
"""Return the state symbol and the coefficient of every entry of a kind of species.
A state name alone is a state of the target. A pair reads the namespace of the plain role that binds its path,
or else the namespace of its variadic role `name`, in the order of the pairs. A coefficient that is a function
is its expression at the namespace.
"""
entries = entries or {}
pairs = [entry for entry in entries if not isinstance(entry, str)]
sources: dict[object, SimpleNamespace] = {}
if pairs and name in self._get_specs():
sources = dict(zip(pairs, getattr(namespace, name), strict=True))
elif pairs:
bound = self._get_bound(namespace)
sources = {pair: bound[pair[0]] for pair in pairs}
return [
(
getattr(
self._get_member(namespace, 'target') if isinstance(entry, str) else sources[entry],
self._get_state(entry),
),
SR(coefficient(namespace)) if callable(coefficient) else coefficient,
)
for entry, coefficient in entries.items()
]
def _get_state(self, entry: tuple[str, str] | str) -> str:
"""Return the name of the state of an entry, a state name alone or the second item of a pair."""
return entry if isinstance(entry, str) else entry[1]
def _get_changes(self, namespace: SimpleNamespace) -> dict[Expression, float]:
"""Return the net change of every state the reaction changes, products minus reactants."""
changes: dict[Expression, float] = {}
for sign, entries, name in (
(-1, self._get_reactants(), 'reactants'),
(1, self._get_products(), 'products'),
):
for state, coefficient in self._get_terms(entries, namespace, name):
changes[state] = changes.get(state, 0) + sign * coefficient
return {state: change for state, change in changes.items() if change}
def _get_mass_action(self, namespace: SimpleNamespace, constant: Expression) -> Expression:
"""Return the mass action of a rate constant, times every reactant and modifier to its coefficient."""
propensity = SR(constant)
for entries, name in ((self._get_reactants(), 'reactants'), (self._get_modifiers(), 'modifiers')):
for state, coefficient in self._get_terms(entries, namespace, name):
propensity *= state**coefficient
return propensity
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the propensity of the reaction, by the user's function or by mass action at the rate constant `k`."""
if self.propensity is not None:
return SR(self.propensity(namespace))
return self._get_mass_action(namespace, self._get_parameter_symbol('k', namespace))
def _get_langevin_laws(self, propensity: Expression, changes: dict[Expression, float]) -> list[BasePart]:
"""Return the drift of the reaction and the noise of its one Wiener process on every state it changes."""
amplitude = sqrt(where(propensity > 0, propensity, 0))
rates: list[BasePart] = [RateLaw(state, change * propensity) for state, change in changes.items()]
noises: list[BasePart] = [
NoiseLaw(state, change * amplitude, process=self.label) for state, change in changes.items()
]
return rates + noises
def _is_entry(self, entry: object) -> bool:
"""Check whether an entry is a pair of a path and a state name, or a state name of a given target."""
if isinstance(entry, str):
return self.target is not None and entry.isidentifier()
return (
isinstance(entry, tuple)
and len(entry) == ENTRY_SIZE
and all(isinstance(item, str) for item in entry)
and entry[1].isidentifier()
)
def _is_coefficient(self, coefficient: object) -> bool:
"""Check whether a coefficient is a positive number or a function."""
return callable(coefficient) or (isinstance(coefficient, int | float) and coefficient > 0)
def _check_entries(self, name: str, entries: Stoichiometry | None) -> None:
"""Raise an error for entries that are neither pairs nor state names of the target, or not coefficients."""
if entries is None:
return
valid = isinstance(entries, Mapping) and all(
self._is_entry(entry) and self._is_coefficient(coefficient) for entry, coefficient in entries.items()
)
if not valid:
error_msg = (
f'The {name} of reaction {self.label!r} are {entries!r}. Map pairs of a path and a state name, or '
'state names of the target, to positive numbers or functions, such as {("S", "n"): 1} or {"n": 1}.'
)
raise ValidationError(error_msg)
def _check_form(self) -> None:
"""Raise an error for an unknown form, or an agent form without a target, with pairs, or of Langevin mode."""
if self.form not in FORMS:
error_msg = f'Reaction {self.label!r} has the form {self.form!r}. Use one of {", ".join(FORMS)}.'
raise ValidationError(error_msg)
self._check_variant('schedule', self.schedule, SCHEDULES)
if self.form != 'agent' and self.schedule != 'synchronous':
error_msg = f'The schedule {self.schedule!r} of reaction {self.label!r} runs agent rules. Use form="agent".'
raise ValidationError(error_msg)
if self.form != 'agent':
return
if self.target is None:
error_msg = f'Reaction {self.label!r} has the agent form but no target. Give the agent group as `target`.'
raise ValidationError(error_msg)
mode = self._get_mode()
if mode not in AGENT_MODES:
error_msg = (
f'Reaction {self.label!r} has the agent form and the mode {mode!r}. Use {" or ".join(AGENT_MODES)}.'
)
raise ValidationError(error_msg)
entries = [*(self._get_reactants() or {}), *(self._get_products() or {}), *(self._get_modifiers() or {})]
if not all(isinstance(entry, str) for entry in entries):
error_msg = (
f'Reaction {self.label!r} has the agent form, whose species are states of every member. Name them by '
'state names alone, such as {"n": 1}.'
)
raise ValidationError(error_msg)
def _check_jump(self) -> None:
"""Raise an error for an unknown mode or form, no reactants and products, or entries the reaction cannot take.
Raises:
ValidationError:
The mode or the form is unknown, the reaction has neither reactants nor products, an entry is neither a
pair of a path and a state name nor a state name of a given target, a coefficient is neither a positive
number nor a function, or the agent form has no target, has an entry that is a pair, or has the mode
`'langevin'`.
"""
mode = self._get_mode()
if mode not in MODES:
error_msg = f'Reaction {self.label!r} has the mode {mode!r}. Use one of {", ".join(MODES)}.'
raise ValidationError(error_msg)
reactants, products = self._get_reactants(), self._get_products()
if not reactants and not products:
error_msg = f'Reaction {self.label!r} has neither reactants nor products. Give at least one of them.'
raise ValidationError(error_msg)
for name, entries in (('reactants', reactants), ('products', products), ('modifiers', self._get_modifiers())):
self._check_entries(name, entries)
self._check_form()
def _get_number(self, value: object, agent: AgentViewProtocol, model: ModelViewProtocol) -> float:
"""Return a number, or the value of a function of an agent and the model as a number."""
found = value(agent, model) if callable(value) else value
if not isinstance(found, SupportsFloat):
error_msg = f'A function of reaction {self.label!r} returns {found!r} for an agent. Return a number.'
raise ValidationError(error_msg)
return float(found)
def _get_agent_changes(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> dict[str, float]:
"""Return the net change of every state of an agent the reaction changes, products minus reactants."""
changes: dict[str, float] = {}
for sign, entries in ((-1, self._get_reactants()), (1, self._get_products())):
for entry, coefficient in (entries or {}).items():
state = self._get_state(entry)
changes[state] = changes.get(state, 0.0) + sign * self._get_number(coefficient, agent, model)
return changes
def _get_agent_propensity(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> float:
"""Return the propensity of an agent, by the user's function or by mass action at the rate constant `k`."""
if self.propensity is not None:
return self._get_number(self.propensity, agent, model)
propensity = self._get_parameter_value('k', model)
for entries in (self._get_reactants(), self._get_modifiers()):
for entry, coefficient in (entries or {}).items():
value = float(getattr(agent, self._get_state(entry)))
propensity *= value ** self._get_number(coefficient, agent, model)
return propensity
def _get_agent_scale(self, propensity: float, model: ModelViewProtocol) -> float:
"""Return the multiple of its changes an agent adds in a step, its propensity or whether an event fires."""
if self._get_mode() == 'rate':
return propensity
return float(model.rng.random() < 1.0 - math.exp(-max(propensity, 0.0)))
def _get_agent_step(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> dict[str, float]:
"""Return the change of every state of an agent in the step, its changes times its propensity or its event."""
changes, propensity = self._get_agent_changes(agent, model), self._get_agent_propensity(agent, model)
scale = self._get_agent_scale(propensity, model)
return {state: change * scale for state, change in changes.items()}
def _apply_agent_step(self, agent: AgentViewProtocol, changes: dict[str, float]) -> None:
"""Add their changes to the states of an agent."""
for state, change in changes.items():
setattr(agent, state, float(getattr(agent, state)) + change)
def _react_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""React in one agent on its own states, in its own activation."""
if agent.alive:
self._apply_agent_step(agent, self._get_agent_step(agent, model))
def react(self, model: ModelViewProtocol) -> None:
"""React in every agent of the target on its own states, by the schedule.
Args:
model:
The model view of the interaction.
"""
self._update_agents(model, self.schedule, self._get_agent_step, self._apply_agent_step)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the laws of the reaction in its mode, or the rule of the agent form.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target, the reactants, the
products, and the modifiers.
Returns:
laws:
The rate laws, the jump law, the rate and noise laws, or the rule of the reaction.
"""
if self.form == 'agent' and self.schedule == 'member':
return [RuleLaw(self._react_member, phase=self.phase, members=self._get_role('target'))]
if self.form == 'agent':
return [RuleLaw(self.react, phase=self.phase)]
mode = self._get_mode()
if mode == 'jump':
return [*super()._get_laws(namespace)]
propensity = self._get_rate(namespace)
changes = self._get_changes(namespace)
if mode == 'langevin':
return self._get_langevin_laws(propensity, changes)
return [RateLaw(state, change * propensity) for state, change in changes.items()]
Equations
site.A: diff(A_site(t), t) == A_site(t)/((A_site(t)^(-n_site + 1)/r2_site)^p_site + (1/r1_site)^p_site)^(1/p_site) - a1_site*A_site(t)*E_site(t)/(((A_site(t)^(-q_site + 1)/g3_site)^p_site + 1)^(1/p_site)*(C3_site(t)/C3_50_site + 1)*(ER_site(t)/g2_site + A_site(t)/(g1_site*((A_site(t)^(-q_site + 1)/g3_site)^p_site + 1)^(1/p_site)) + E_site(t)))
site.EUD: diff(EUD_site(t), t) == -d1_site*EUD_site(t) - lambda_a_site*EUD_site(t)/(EUD_site(t)/EMH_site + 1) + alpha_a_site*A_site(t)/(((A_site(t)^(-q_site + 1)/g3_site)^p_site + 1)^(1/p_site)*(C2_site(t)/C2_50_site + 1)*(ER_site(t)/ER_50_site + 1))
site.ED: diff(ED_site(t), t) == -d2_site*ED_site(t) + lambda_a_site*EUD_site(t)/(EUD_site(t)/EMH_site + 1)
site.EPC: diff(EPC_site(t), t) == EMC_site*a2_site*ED_site(t)/((EMC_site + EMH_site + EMR_site)*g4_site + ED_site(t)) - d3_site*EPC_site(t)
site.E: diff(E_site(t), t) == -d4_site*E_site(t) + a3_site*C1_site(t)*EPC_site(t)/((EC1_50_site + C1_site(t))*(C3_site(t)/EC3_50_site + 1))
site.EPH: diff(EPH_site(t), t) == EMH_site*a4_site*(ED_site(t) + EUD_site(t))/((EMC_site + EMH_site + EMR_site)*g4_site + ED_site(t) + EUD_site(t)) - d5_site*EPH_site(t)
site.EH: diff(EH_site(t), t) == -d6_site*EH_site(t) - a6_site*C3_site(t)*EH_site(t)/(g5_site + C3_site(t)) + a5_site*C1_site(t)*EPH_site(t)/((EC1_50_site + C1_site(t))*(C3_site(t)/EC3_50_site + 1))
site.EPR: diff(EPR_site(t), t) == EMR_site*a7_site*ED_site(t)/((EMC_site + EMH_site + EMR_site)*g4_site + ED_site(t)) - d7_site*EPR_site(t)
site.ER: diff(ER_site(t), t) == a6_site*C3_site(t)*EH_site(t)/(g5_site + C3_site(t)) + a8_site*C1_site(t)*EPR_site(t)/(EC1_50_site + C1_site(t)) - d8_site*ER_site(t)
site.C1: diff(C1_site(t), t) == -C1_site(t)/tau1_site + a9_site*EPH_site(t)/((C2_site(t)/CC2_50_site + 1)*(C3_site(t)/CC3_50_site + 1))
site.C2: diff(C2_site(t), t) == r4_site*A_site(t) + r3_site*ER_site(t) - C2_site(t)/tau2_site
site.C3: diff(C3_site(t), t) == r6_site*A_site(t) + r5_site*ER_site(t) - C3_site(t)/tau3_site