Joint, community, and SteadyCom flux balance
Two members of the E. coli core system that share a medium, which needs the extra metabolism.
- Name:
community_fba - Field:
biology/ecology - Formalism: static solution
- References: Raman 2021; Khandelwal et al. 2013; Chan et al. 2017; Diener et al. 2020
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('community_fba')
entry = describe('community_fba')
trajectory = create_simulation('community_fba').prepare().run().trajectory_
Factory
@register_system(
field='biology/ecology',
tags=('community flux balance', 'SteadyCom', 'cross-feeding', 'linear program', 'metabolism', 'static'),
references=('Raman 2021', 'Khandelwal et al. 2013', 'Chan et al. 2017', 'Diener et al. 2020'),
defaults={'until': 0.0},
)
def community_fba(
knockouts: Sequence[str] = ('ICDHyr', 'NH4t'),
*,
abundances: Sequence[float] = (0.5, 0.5),
exchanges: Mapping[str, float] | None = None,
common_rate: bool = False,
) -> System:
"""Return two members of the E. coli core system that share a medium, which needs the extra `metabolism`.
The programs of the members merge through a medium whose law holds the growth weighted by the abundances (Raman
2021, Section 11.1.3). By default one member lacks isocitrate dehydrogenase and takes up 2-oxoglutarate, and the
other lacks ammonium transport and takes up glutamate, each from what the other secretes, so neither grows alone,
and together both grow. `CommunityAnalysis` gives SteadyCom and the cooperative trade-off of the same members. The
pair is an illustration found on the E. coli core system: the consortium of Raman (Section 11.1.3.6) has four
mutants that need amino acids, whose synthesis the core system does not hold.
Args:
knockouts:
The reaction to close in each member, or the empty string for none.
abundances:
The abundance of each member.
exchanges:
The lower bound of the net exchange of the community, by reaction, or `None` for no net uptake of the
cross-fed metabolites.
common_rate:
Whether the members grow at one specific rate, as in community FBA, rather than each at its own, as in joint
FBA.
Returns:
system:
The system.
Raises:
MissingExtraError:
The extra is not installed.
Examples:
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.biology.ecology import community_fba
>>> trajectory = Simulation(community_fba()).prepare().run(0.0).trajectory_
>>> round(float(trajectory.values('medium.objective_0')[0]), 3)
0.862
"""
shared = dict(exchanges) if exchanges is not None else dict.fromkeys(CROSS_FED, 0.0)
members = []
for position, knockout in enumerate(knockouts):
member = load_network('textbook', f'member_{position}', needed_by='community_fba')
member.objective = {}
for reaction in shared:
if shared[reaction] >= 0:
member.set_params(**{f'lower_{member._identifier(reaction)}': OPEN_UPTAKE})
if knockout:
identifier = member._identifier(knockout)
member.set_params(**{f'lower_{identifier}': 0.0, f'upper_{identifier}': 0.0})
members.append(member)
return System(
elements=[*members, community_medium_type(list(shared))('medium')],
interactions=[
CommunityInteraction(
'joint',
members=[member.label for member in members],
medium='medium',
abundances=tuple(abundances),
exchanges=shared,
growth=GROWTH,
common_rate=common_rate,
),
],
time=StaticTime(),
)
Components
Elements:
MetabolicNetworkElement, fromDeclaredElementCommunityMediumElement, fromDeclaredElement
Interactions:
CommunityInteraction, fromInteraction
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 CommunityMediumElement(DeclaredElement):
"""Medium that a community shares, with one algebraic state per shared exchange.
The element declares the state `u_<reaction>` of the net exchange of every reaction.
Args:
label:
The label of the element.
visual:
How the element is drawn, or `None` for `default_visual_`.
exchanges:
The names of the reactions.
variables:
The declared states that replace the states of the net exchanges, by name, or `None` for none.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import CommunityMediumElement
>>> CommunityMediumElement.state_name('EX_glc__D_e')
'u_EX_glc__D_e'
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
exchanges: Sequence[str] = (),
variables: Mapping[str, DeclaredVariableSpec] | None = None,
) -> None:
super().__init__(label, visual=visual, variables=variables)
self.exchanges = exchanges
def _declare_variables(self) -> Mapping[str, DeclaredVariableSpec]:
"""Return the state of the net exchange of every reaction, with the variables that replace them."""
flux = StateSpec(algebraic=True)
return {
**{self.state_name(reaction): DeclaredVariableSpec(0.0, flux) for reaction in self.exchanges},
**(self.variables or {}),
}
@staticmethod
def state_name(reaction: str) -> str:
"""Return the name of the state of the net exchange of a reaction.
Args:
reaction:
The name of the reaction.
Returns:
name:
The name of the state.
"""
return f'u_{MetabolicNetworkElement._identifier(reaction)}'
class CommunityInteraction(Interaction):
"""Joint program of a community: the medium balances the exchanges, and its law holds the objective of growth.
Args:
label:
The label of the interaction.
members:
The paths of the metabolic networks.
medium:
The path of the shared medium.
abundances:
The abundance of every member, which weighs its exchanges and its growth.
exchanges:
The lower bound of the net exchange of the community, by reaction.
growth:
The name of the reaction of growth of every member.
common_rate:
Whether every member grows at one specific rate, the balanced growth of community FBA (Raman 2021, Section
11.1.3.2), rather than each at its own, as in joint FBA.
Examples:
>>> from skcomplex.systems.biology.common.metabolism import CommunityInteraction
>>> CommunityInteraction('joint', members=['a', 'b'], medium='medium', abundances=(0.5, 0.5)).abundances
(0.5, 0.5)
"""
def __init__(
self,
label: str,
*,
members: Annotated[list[str], RoleSpec(variadic=True)],
medium: Annotated[str, RoleSpec()],
abundances: Sequence[float] = (),
exchanges: Mapping[str, float] | None = None,
growth: str = 'Biomass_Ecoli_core',
common_rate: bool = False,
) -> None:
super().__init__(label)
self.members = members
self.medium = medium
self.abundances = abundances
self.exchanges = exchanges
self.growth = growth
self.common_rate = common_rate
def _weighted(self, namespace: SimpleNamespace, reaction: str) -> Expression:
"""Return the sum over the members of the flux of a reaction weighted by the abundance."""
name = MetabolicNetworkElement.flux_name(reaction)
terms = [
float(share) * getattr(member, name)
for share, member in zip(self.abundances, namespace.members, strict=True)
]
return sum(terms[1:], terms[0])
def _get_laws(self, namespace: SimpleNamespace) -> Sequence[BaseLaw]:
"""Return the program that joins the members through the medium.
Args:
namespace:
The namespace of the roles.
Returns:
laws:
The linear program law.
"""
variables = []
relations: list[Expression] = []
for reaction, floor in (self.exchanges or {}).items():
state = getattr(namespace.medium, CommunityMediumElement.state_name(reaction))
variables.append(state)
relations += [state == self._weighted(namespace, reaction), state >= float(floor)]
name = MetabolicNetworkElement.flux_name(self.growth)
growth = [getattr(member, name) for member in namespace.members]
if self.common_rate:
relations += [item == growth[0] for item in growth[1:]]
objective = LinearObjectiveSpec(self._weighted(namespace, self.growth))
return [LinearProgramLaw([*variables, *growth], relations, [objective])]