Dynamic multi-species metabolic system
Two E. coli core strains that differ in maximal glucose uptake on one medium, which needs the extra.
- Name:
dynamic_multispecies_fba - Field:
biology/ecology - Formalism: co-simulation
- References: Raman 2021; Zhuang et al. 2011
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('dynamic_multispecies_fba')
entry = describe('dynamic_multispecies_fba')
trajectory = create_simulation('dynamic_multispecies_fba').prepare().run().trajectory_
Factory
@register_system(
field='biology/ecology',
tags=('dynamic flux balance', 'competition', 'linear program', 'metabolism', 'co-simulation'),
references=('Raman 2021', 'Zhuang et al. 2011'),
defaults={
'until': 8.0,
'dt': SYNCHRONISATION,
'solver': CoSimulationSolver(sync_dt=SYNCHRONISATION, solvers={'culture': RK4Solver(dt=ODE_STEP)}),
},
)
def dynamic_multispecies_fba(
*,
fast_vmax: float = 10.0,
slow_vmax: float = 8.0,
biomass: float = 0.01,
glucose: float = 10.0,
km: float = 0.015,
) -> System:
"""Return two E. coli core strains that differ in maximal glucose uptake on one medium, which needs the extra.
Each strain has its own program, uptake kinetics, and biomass, and they read one glucose pool (Raman 2021, Section
11.1.3.4). The strain with the larger uptake grows faster, so its share rises until the glucose runs out. The two
strains and their uptakes of 10 and 8 are an illustration of the framework, whose example in the source is the
competition between Rhodoferax and Geobacter (Zhuang et al. 2011).
Args:
fast_vmax:
The largest glucose uptake of the fast strain, in mmol per gDW per h.
slow_vmax:
The largest glucose uptake of the slow strain, in mmol per gDW per h.
biomass:
The initial biomass of each strain in gDW per L.
glucose:
The initial glucose in mM.
km:
The half-saturation concentration of glucose in mM.
Returns:
system:
The system.
Raises:
MissingExtraError:
The extra is not installed.
Examples:
>>> from skcomplex.systems.biology.ecology import dynamic_multispecies_fba
>>> [element.label for element in dynamic_multispecies_fba().systems[1].elements]
['biomass_fast', 'biomass_slow', 'glucose']
"""
cells = []
interactions = []
for name, vmax in (('fast', fast_vmax), ('slow', slow_vmax)):
cell = load_network('textbook', f'cell_{name}', needed_by='dynamic_multispecies_fba')
cell.parsimonious, cell.infeasible = True, 'zero'
cells.append(cell)
interactions += [
ExchangeInteraction(
f'uptake_{name}',
cell=f'fba/cell_{name}',
biomass=f'culture/biomass_{name}',
substrate='culture/glucose',
reaction=GLUCOSE_EXCHANGE,
vmax=vmax,
km=km,
interval=SYNCHRONISATION,
),
GrowthInteraction(
f'growth_{name}',
cell=f'fba/cell_{name}',
biomass=f'culture/biomass_{name}',
reaction=GROWTH,
),
]
return System(
systems=[
System('fba', elements=cells, time=DiscreteTime()),
System(
'culture',
elements=[
BiomassElement('biomass_fast', x=biomass),
BiomassElement('biomass_slow', x=biomass),
SubstrateElement('glucose', c=glucose),
],
),
],
interactions=interactions,
)
Components
Elements:
MetabolicNetworkElement, fromDeclaredElementBiomassElement, fromElementSubstrateElement, fromElement
Interactions:
ExchangeInteraction, fromInteractionGrowthInteraction, fromInteraction
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 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