Regulatory flux balance analysis
A batch culture on glucose and acetate with Boolean catabolite repression, which needs the extra.
- Name:
regulatory_fba - Field:
biology/systems_biology - Formalism: co-simulation
- References: Raman 2021; Covert et al. 2001
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('regulatory_fba')
entry = describe('regulatory_fba')
trajectory = create_simulation('regulatory_fba').prepare().run().trajectory_
Factory
@register_system(
field='biology/systems_biology',
tags=('regulatory flux balance', 'boolean network', 'linear program', 'metabolism', 'diauxie', 'co-simulation'),
references=('Raman 2021', 'Covert et al. 2001'),
defaults={
'until': 8.0,
'dt': SYNCHRONISATION,
'solver': CoSimulationSolver(sync_dt=SYNCHRONISATION, solvers={'culture': RK4Solver(dt=ODE_STEP)}),
},
)
def regulatory_fba(
*,
biomass: float = 0.01,
glucose: float = 10.0,
acetate: float = 10.0,
threshold: float = 0.1,
lag: int = 10,
acetate_vmax: float = 10.0,
km: float = 0.015,
) -> System:
"""Return a batch culture on glucose and acetate with Boolean catabolite repression, which needs the extra.
Boolean rules alternate with the program (Raman 2021, Eq. 9.7). A node reads whether glucose lies above a threshold,
a chain of nodes delays it, and while the last node is on it holds the genes of acetate kinase off through update
laws. So the culture leaves acetate alone while glucose lasts, stops at depletion, and grows again on acetate after
the lag of the chain.
The alternation is that of Raman, but the rules are an illustration and not those of Covert et al. (2001): the genes
of acetate kinase, the threshold, and the chain of nodes come from no source, and the diauxie of Raman is on glucose
and lactose, which the E. coli core system cannot take up.
Args:
biomass:
The initial biomass in gDW per L.
glucose:
The initial glucose in mM.
acetate:
The initial acetate in mM.
threshold:
The glucose above which the sensor is on, in mM.
lag:
The number of nodes after the sensor, each a synchronisation interval of delay.
acetate_vmax:
The largest acetate uptake in mmol per gDW per h.
km:
The half-saturation concentration of both substrates in mM.
Returns:
system:
The system.
Raises:
MissingExtraError:
The extra is not installed.
Examples:
>>> from skcomplex.systems.biology.systems_biology import regulatory_fba
>>> [system.label for system in regulatory_fba(lag=1).systems]
['control', 'culture']
"""
cell = load_network('textbook', 'cell', needed_by='regulatory_fba')
cell.parsimonious, cell.infeasible = True, 'zero'
identifier = cell._identifier(ACETATE_EXCHANGE)
cell.set_params(
**{f'lower_{identifier}': -float(acetate_vmax)},
**{cell.gene_name(gene): 0.0 for gene in ACETATE_GENES},
)
nodes = [BooleanNodeElement(f'node_{position}', x=1) for position in range(lag + 1)]
chain = [
GlucoseSensorInteraction(
f'delay_{position}',
node=f'control/node_{position}',
source=f'control/node_{position - 1}',
)
for position in range(1, lag + 1)
]
return System(
systems=[
System('control', elements=[cell, *nodes], time=DiscreteTime()),
System(
'culture',
elements=[
BiomassElement('biomass', x=biomass),
SubstrateElement('glucose', c=glucose),
SubstrateElement('acetate', c=acetate),
],
),
],
interactions=[
GlucoseSensorInteraction(
'sensor',
node='control/node_0',
source='culture/glucose',
threshold=threshold,
),
*chain,
SwitchInteraction('repression', node=f'control/node_{lag}', cell='control/cell', genes=ACETATE_GENES),
ExchangeInteraction(
'glucose_uptake',
cell='control/cell',
biomass='culture/biomass',
substrate='culture/glucose',
reaction=GLUCOSE_EXCHANGE,
km=km,
interval=SYNCHRONISATION,
),
ExchangeInteraction(
'acetate_uptake',
cell='control/cell',
biomass='culture/biomass',
substrate='culture/acetate',
reaction=ACETATE_EXCHANGE,
vmax=acetate_vmax,
km=km,
interval=SYNCHRONISATION,
),
GrowthInteraction('growth', cell='control/cell', biomass='culture/biomass', reaction=GROWTH),
],
)
Components
Elements:
MetabolicNetworkElement, fromDeclaredElementBooleanNodeElement, fromElementBiomassElement, fromElementSubstrateElement, fromElement
Interactions:
GlucoseSensorInteraction, fromInteractionSwitchInteraction, fromInteractionExchangeInteraction, fromInteractionGrowthInteraction, fromInteraction
class GlucoseSensorInteraction(Interaction):
"""Boolean node that reads whether a substrate lies above a threshold, or the previous node of a chain.
Args:
label:
The label of the interaction.
node:
The path of the Boolean node.
source:
The path of what the node reads: a substrate with a threshold, or a Boolean node without one.
threshold:
The concentration above which the node turns on, or `None` for a node that copies another.
Examples:
>>> from skcomplex.systems.biology.systems_biology import GlucoseSensorInteraction
>>> GlucoseSensorInteraction('sensor', node='a', source='glucose', threshold=0.1).threshold
0.1
"""
def __init__(
self,
label: str,
*,
node: Annotated[str, RoleSpec()] = 'node',
source: Annotated[str, RoleSpec()] = 'source',
threshold: float | None = None,
) -> None:
super().__init__(label)
self.node = node
self.source = source
self.threshold = threshold
def _get_laws(self, namespace: SimpleNamespace) -> Sequence[BaseLaw]:
"""Return the update of the node.
Args:
namespace:
The namespace of the roles.
Returns:
laws:
The update law of the node.
"""
if self.threshold is None:
return [UpdateLaw(namespace.node.x, namespace.source.x)]
return [UpdateLaw(namespace.node.x, where(namespace.source.c > float(self.threshold), 1, 0))]
class SwitchInteraction(Interaction):
"""Update of the switches of genes from a Boolean node that represses them.
Args:
label:
The label of the interaction.
node:
The path of the Boolean node, whose state `x` is 1 for on.
cell:
The path of the metabolic network.
genes:
The names of the genes the node represses.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import SwitchInteraction
>>> SwitchInteraction('repression', node='node', cell='cell', genes=('b2296',)).genes
('b2296',)
"""
def __init__(
self,
label: str,
*,
node: Annotated[str, RoleSpec()] = 'node',
cell: Annotated[str, RoleSpec()] = 'cell',
genes: Sequence[str] = (),
) -> None:
super().__init__(label)
self.node = node
self.cell = cell
self.genes = genes
def _get_laws(self, namespace: SimpleNamespace) -> Sequence[BaseLaw]:
"""Return the update of every switch.
Args:
namespace:
The namespace of the roles.
Returns:
laws:
One update law per gene.
"""
return [
UpdateLaw(getattr(namespace.cell, MetabolicNetworkElement.gene_name(gene)), 1 - namespace.node.x)
for gene in self.genes
]
class ExchangeInteraction(Interaction):
"""Exchange of a metabolite between a metabolic network and its medium, bounded by Michaelis-Menten uptake.
The law bounds the uptake, the negative flux of the exchange reaction, by `vmax c / (km + c)` and, with an interval,
by what the medium holds per biomass over the interval (Raman 2021, Section 9.6). The medium changes at the rate
`dc/dt = v x`, and a depleted medium loses nothing more.
Args:
label:
The label of the interaction.
cell:
The path of the metabolic network.
biomass:
The path of the biomass.
substrate:
The path of the metabolite.
reaction:
The name of the exchange reaction.
vmax:
The largest uptake.
km:
The half-saturation concentration.
interval:
The time over which the medium must last, or `None` for no such bound.
threshold:
The biomass above which the program solves, or `None` for always.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import ExchangeInteraction
>>> ExchangeInteraction('uptake', reaction='EX_glc__D_e').vmax
10.0
"""
def __init__(
self,
label: str,
*,
cell: Annotated[str, RoleSpec()] = 'cell',
biomass: Annotated[str, RoleSpec()] = 'biomass',
substrate: Annotated[str, RoleSpec()] = 'substrate',
reaction: str = 'EX_glc__D_e',
vmax: Annotated[float, ParameterSpec(domain=POSITIVE)] = 10.0,
km: Annotated[float, ParameterSpec(domain=POSITIVE)] = 0.015,
interval: float | None = None,
threshold: float | None = None,
) -> None:
super().__init__(label)
self.cell = cell
self.biomass = biomass
self.substrate = substrate
self.reaction = reaction
self.vmax = vmax
self.km = km
self.interval = interval
self.threshold = threshold
def _get_laws(self, namespace: SimpleNamespace) -> Sequence[BaseLaw]:
"""Return the bound of the uptake and the rate of the medium.
Args:
namespace:
The namespace of the roles and parameters.
Returns:
laws:
The linear program law of the bound and the rate law of the medium.
"""
flux = getattr(namespace.cell, MetabolicNetworkElement.flux_name(self.reaction))
amount, biomass = namespace.substrate.c, namespace.biomass.x
held = where(amount > 0, amount, 0)
constraints = [flux >= -namespace.vmax * held / (namespace.km + held)]
if self.interval is not None:
constraints.append(flux >= -held / (biomass * float(self.interval)))
active = None if self.threshold is None else biomass > float(self.threshold)
return [
LinearProgramLaw([flux], constraints, active=active, infeasible='zero'),
RateLaw(amount, flux * biomass * where(either(amount > 0, flux > 0), 1, 0)),
]
class GrowthInteraction(Interaction):
"""Growth of biomass at the rate that a reaction of a metabolic network gives, `dx/dt = v x`.
Args:
label:
The label of the interaction.
cell:
The path of the metabolic network.
biomass:
The path of the biomass.
reaction:
The name of the reaction of growth.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import GrowthInteraction
>>> GrowthInteraction('growth', cell='cell', biomass='culture', reaction='grow').reaction
'grow'
"""
def __init__(
self,
label: str,
*,
cell: Annotated[str, RoleSpec()] = 'cell',
biomass: Annotated[str, RoleSpec()] = 'biomass',
reaction: str = 'Biomass_Ecoli_core',
) -> None:
super().__init__(label)
self.cell = cell
self.biomass = biomass
self.reaction = reaction
def _get_laws(self, namespace: SimpleNamespace) -> Sequence[BaseLaw]:
"""Return the rate of the biomass.
Args:
namespace:
The namespace of the roles.
Returns:
laws:
The rate law of the biomass.
"""
flux = getattr(namespace.cell, MetabolicNetworkElement.flux_name(self.reaction))
return [RateLaw(namespace.biomass.x, flux * namespace.biomass.x)]
class MetabolicNetworkElement(DeclaredElement):
"""Metabolic network whose fluxes one linear program sets, with a balance per metabolite.
The element declares one algebraic state `v_<reaction>` per reaction, one adaptive switch `g_<gene>` per gene, and
the parameters `lower_<reaction>` and `upper_<reaction>` of its bounds. The rule of a reaction gates both bounds, so
a reaction whose rule is false carries no flux.
Args:
label:
The label of the element.
visual:
How the element is drawn, or `None` for the bars of its program.
objective:
The coefficient of every reaction in the objective, or `None` for `default_objective`, or nothing for a
network whose objective another law holds.
parsimonious:
Whether to minimise the total absolute flux after the objective.
infeasible:
The infeasibility policy of the program.
name:
The name of the network, or `None` for none.
reactions:
The reactions by name.
default_objective:
The objective of the network when `objective` is `None`.
variables:
The declared states and parameters that replace the ones of the reactions, by name, or `None` for none.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import ReactionSpec, metabolic_network_type
>>> reactions = {'uptake': ReactionSpec({'A': 1.0}), 'use': ReactionSpec({'A': -1.0}, upper=4.0)}
>>> network = metabolic_network_type('TwoSteps', reactions, {'use': 1.0})('cell')
>>> network.get_value('upper_use')
4.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
objective: Mapping[str, float] | None = None,
parsimonious: bool = False,
infeasible: str = 'raise',
name: str | None = None,
reactions: Mapping[str, ReactionSpec] = NO_REACTIONS,
default_objective: Mapping[str, float] = NO_OBJECTIVE,
variables: Mapping[str, DeclaredVariableSpec] | None = None,
) -> None:
super().__init__(label, visual=visual, variables=variables)
self.objective = objective
self.parsimonious = parsimonious
self.infeasible = infeasible
self.name = name
self.reactions = reactions
self.default_objective = default_objective
def _declare_variables(self) -> Mapping[str, DeclaredVariableSpec]:
"""Return the states and parameters of the reactions and genes, with the variables that replace them."""
flux, bound = StateSpec(algebraic=True), ParameterSpec()
declared: dict[str, DeclaredVariableSpec] = {}
for reaction, spec in self.reactions.items():
identifier = self._identifier(reaction)
declared[f'v_{identifier}'] = DeclaredVariableSpec(0.0, flux)
declared[f'lower_{identifier}'] = DeclaredVariableSpec(float(spec.lower), bound)
declared[f'upper_{identifier}'] = DeclaredVariableSpec(float(spec.upper), bound)
switch = ParameterSpec(adaptive=True, bounds=(0.0, 1.0))
for gene in self.genes_of():
declared[self.gene_name(gene)] = DeclaredVariableSpec(1.0, switch)
return {**declared, **(self.variables or {})}
def _prepare(self, context: BaseSpace | None = None) -> None:
"""Set `default_visual_`, the bars of the fluxes.
Args:
context:
The space the element inherits from its system, which a network does not read.
"""
self.default_visual_ = VisualSpec('linear_program')
@staticmethod
def _identifier(name: str) -> str:
"""Return a name as a part of an identifier."""
return re.sub(r'\W', '_', name)
@classmethod
def flux_name(cls, reaction: str) -> str:
"""Return the name of the state of the flux of a reaction.
Args:
reaction:
The name of the reaction.
Returns:
name:
The name of the state.
"""
return f'v_{cls._identifier(reaction)}'
@classmethod
def gene_name(cls, gene: str) -> str:
"""Return the name of the switch of a gene.
Args:
gene:
The name of the gene.
Returns:
name:
The name of the switch.
"""
return f'g_{cls._identifier(gene)}'
def genes_of(self) -> tuple[str, ...]:
"""Return the names of the genes that the rules of the reactions read.
Returns:
genes:
The names, in the order the rules name them.
"""
names = [name for reaction in self.reactions.values() for name in _rule_genes(reaction.rule)]
return tuple(dict.fromkeys(names))
def get_number(self, name: str) -> float:
"""Return the value of a state or parameter of the network as a number.
Args:
name:
The name of the variable, such as `upper_<reaction>`.
Returns:
number:
The value.
"""
return float(cast(float, self.get_value(name)))
def _gate(self, namespace: SimpleNamespace, rule: str) -> Expression | int:
"""Return the condition of a rule on the switches of its genes, or 1 for a reaction without genes."""
if not rule.strip():
return 1
return both(self._condition(namespace, ast.parse(rule, mode='eval').body))
def _condition(self, namespace: SimpleNamespace, node: ast.expr) -> Expression:
"""Return the condition of a node of a rule."""
if isinstance(node, ast.BoolOp):
combine = both if isinstance(node.op, ast.And) else either
return combine(*(self._condition(namespace, item) for item in node.values))
if isinstance(node, ast.Name):
return getattr(namespace, self.gene_name(node.id)) > SWITCH_THRESHOLD
error_msg = f'A gene rule of {self.label} holds {ast.dump(node)}. Use gene names, `and`, `or`, and parentheses.'
raise ValidationError(error_msg)
def _balances(self, namespace: SimpleNamespace) -> list[Expression]:
"""Return the balance of every metabolite."""
terms: dict[str, list[Expression]] = {}
for reaction, spec in self.reactions.items():
flux = getattr(namespace, self.flux_name(reaction))
for metabolite, coefficient in spec.coefficients.items():
terms.setdefault(metabolite, []).append(float(coefficient) * flux)
return [sum(items[1:], items[0]) == 0 for items in terms.values()]
def _bounds(self, namespace: SimpleNamespace) -> list[Expression]:
"""Return the lower and upper bound of every reaction, gated by its rule."""
relations = []
for reaction, spec in self.reactions.items():
flux = getattr(namespace, self.flux_name(reaction))
gate = self._gate(namespace, spec.rule)
identifier = self._identifier(reaction)
relations.append(flux >= getattr(namespace, f'lower_{identifier}') * gate)
relations.append(flux <= getattr(namespace, f'upper_{identifier}') * gate)
return relations
def _objectives(self, namespace: SimpleNamespace) -> list[LinearObjectiveSpec]:
"""Return the objective of the network and, for a parsimonious network, its total absolute flux."""
weights = self.default_objective if self.objective is None else self.objective
if not weights:
return []
terms = [float(weight) * getattr(namespace, self.flux_name(name)) for name, weight in weights.items()]
objectives = [LinearObjectiveSpec(sum(terms[1:], terms[0]))]
if self.parsimonious:
fluxes = [abs(getattr(namespace, self.flux_name(name))) for name in self.reactions]
objectives.append(LinearObjectiveSpec(sum(fluxes[1:], fluxes[0]), sense='min'))
return objectives
def _get_laws(self, namespace: SimpleNamespace) -> Sequence[LinearProgramLaw]:
"""Return the linear program of the network.
Args:
namespace:
The namespace of the symbols of the element.
Returns:
laws:
The linear program law.
"""
variables = [SR(getattr(namespace, self.flux_name(name))) for name in self.reactions]
constraints = [*self._balances(namespace), *self._bounds(namespace)]
objectives = self._objectives(namespace)
return [LinearProgramLaw(variables, constraints, objectives, infeasible=self.infeasible)]
class BooleanNodeElement(Element):
"""Node of a Boolean network, a gene or protein that is on or off.
Args:
label:
The label of the node.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The initial value, 0 or 1.
Examples:
>>> from skcomplex.systems.biology.systems_biology import BooleanNodeElement
>>> BooleanNodeElement('Sic1', x=1).x
1
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
class BiomassElement(Element):
"""Biomass of a culture.
Args:
label:
The label of the element.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The biomass in gDW per L.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import BiomassElement
>>> BiomassElement('culture', x=0.01).x
0.01
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec()] = 0.01,
) -> None:
super().__init__(label, visual=visual)
self.x = x
class SubstrateElement(Element):
"""Concentration of a metabolite of a medium.
Args:
label:
The label of the element.
visual:
How the element is drawn, or `None` for `default_visual_`.
c:
The concentration in mM.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import SubstrateElement
>>> SubstrateElement('glucose', c=10.0).c
10.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
c: Annotated[float, StateSpec()] = 10.0,
) -> None:
super().__init__(label, visual=visual)
self.c = c