Dynamic flux balance analysis of a batch culture
A batch culture of the E. coli core system, which needs the extra metabolism.
- Name:
dynamic_fba - Field:
biology/systems_biology - Formalism: co-simulation
- References: Raman 2021; Mahadevan et al. 2002; Varma and Palsson 1994
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('dynamic_fba')
entry = describe('dynamic_fba')
trajectory = create_simulation('dynamic_fba').prepare().run().trajectory_
Factory
@register_system(
field='biology/systems_biology',
tags=('dynamic flux balance', 'linear program', 'metabolism', 'batch culture', 'co-simulation'),
references=('Raman 2021', 'Mahadevan et al. 2002', 'Varma and Palsson 1994'),
defaults={'until': 10.0, 'dt': SYNCHRONISATION},
settings=_settings,
)
def dynamic_fba(
approach: str = 'static',
*,
biomass: float = 0.01,
glucose: float = 10.0,
vmax: float = 10.0,
km: float = 0.015,
interval: float = SYNCHRONISATION,
) -> System:
"""Return a batch culture of the E. coli core system, which needs the extra `metabolism`.
The program solves for the fluxes with Michaelis-Menten glucose uptake and then the smallest total flux. By the
static optimisation approach it solves in a discrete partition at every synchronisation point, and the biomass and
the glucose integrate with the fluxes held (Raman 2021, Eqs. 9.8 and 9.9). By the direct approach it solves at every
evaluation of the drift. Without glucose the program is infeasible, and the policy `'zero'` stops growth, so biomass
grows in one phase and stays.
Args:
approach:
`'static'` for the static optimisation approach or `'direct'` for the direct approach.
biomass:
The initial biomass in gDW per L.
glucose:
The initial glucose in mM.
vmax:
The largest glucose uptake in mmol per gDW per h.
km:
The half-saturation concentration of glucose in mM.
interval:
The synchronisation interval of the static optimisation approach in h.
Returns:
system:
The system.
Raises:
ValidationError:
The approach is unknown.
MissingExtraError:
The extra is not installed.
Examples:
>>> from skcomplex.systems.biology.systems_biology import dynamic_fba
>>> [system.label for system in dynamic_fba().systems]
['fba', 'culture']
"""
if approach not in APPROACHES:
error_msg = f'Unknown approach {approach!r}. Use one of {", ".join(APPROACHES)}.'
raise ValidationError(error_msg)
cell = load_network('textbook', 'cell', needed_by='dynamic_fba')
cell.parsimonious, cell.infeasible = True, 'zero'
amounts = [BiomassElement('biomass', x=biomass), SubstrateElement('glucose', c=glucose)]
static = approach == 'static'
cell_path = 'fba/cell' if static else 'cell'
biomass_path = 'culture/biomass' if static else 'biomass'
interactions = [
ExchangeInteraction(
'uptake',
cell=cell_path,
biomass=biomass_path,
substrate='culture/glucose' if static else 'glucose',
reaction=GLUCOSE_EXCHANGE,
vmax=vmax,
km=km,
interval=interval if static else None,
),
GrowthInteraction('growth', cell=cell_path, biomass=biomass_path, reaction=GROWTH),
]
if static:
return System(
systems=[System('fba', elements=[cell], time=DiscreteTime()), System('culture', elements=amounts)],
interactions=interactions,
)
return System(elements=[cell, *amounts], 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