COMETS: dynamic flux balance on a lattice with diffusion
An expanding colony of the E. coli core system on a lattice, which needs the extra metabolism.
- Name:
comets - Field:
biology/ecology - Formalism: ordinary differential equations
- References: Raman 2021; Harcombe et al. 2014
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('comets')
entry = describe('comets')
trajectory = create_simulation('comets').prepare().run().trajectory_
Factory
@register_system(
field='biology/ecology',
tags=('COMETS', 'dynamic flux balance', 'lattice', 'diffusion', 'linear program', 'metabolism'),
references=('Raman 2021', 'Harcombe et al. 2014'),
defaults={'until': 14.0, 'dt': STEP, 'solver': IMEXEulerSolver(dt=STEP)},
)
def comets(
shape: tuple[int, int] = (50, 50),
*,
box: float = 0.02,
inoculum: float = 3e-7,
glucose: float = 5.0,
plate_glucose: float = 0.0088,
biomass_diffusion: float = 3e-10,
glucose_diffusion: float = 5e-6,
vmax: float = 10.0,
km: float = 0.01,
interval: float = STEP,
) -> System:
"""Return an expanding colony of the E. coli core system on a lattice, which needs the extra `metabolism`.
Run with `IMEXEulerSolver`, every step is a growth step and then a diffusion step, the two steps that COMETS
alternates (Raman 2021, Section 11.1.3.4). The growth step solves the program of every site that holds biomass and
applies the update rules of Harcombe et al. (2014). The diffusion step solves one implicit Euler system, where
Harcombe et al. use an alternating direction implicit scheme. Biomass and glucose both diffuse, as in Harcombe et
al., while the short description of Raman names metabolites only.
The states are gDW per L and mM of the medium under a box. A site whose program has no solution has zero fluxes,
which happens when its glucose cannot pay for the maintenance flux of the network.
Harcombe et al. in what follows is the main text of Harcombe et al. (2014), and the source of each default is named
under its argument. The time step of 0.1 h of the registered run is from its In silico experiments. The length of
the registered run keeps it within the budget of the gallery and has no source.
These choices have no source in that text. The periodic edges are those of `LatticeSpace`, since the text does not
state the outer boundary. The volume under a box is the amount of glucose under it over the concentration of the
medium, since the text defers the volume to its supplement, which was not read. The molar mass of glucose is 180.16
g/mol.
The network is the E. coli core system, where Harcombe et al. use iJO1366. The death rate of 1% per step, the
largest colony height of 200 um, and the oxygen of 0.25 mmol/cm^2 of Harcombe et al. are left out, so the biomass of
a box has no upper limit and oxygen is limited only by the network.
Args:
shape:
The number of boxes along each side, by default from Harcombe et al., Colony expansion comparisons.
box:
The side of a box in cm, by default from Harcombe et al., Colony expansion comparisons.
inoculum:
The initial biomass of the central box in g, by default from Harcombe et al., Colony expansion comparisons.
glucose:
The glucose of the medium in g/L, by default from Harcombe et al., In silico experiments.
plate_glucose:
The glucose of the medium per area of the plate in mmol/cm^2, by default from Harcombe et al., In silico
experiments.
biomass_diffusion:
The diffusion coefficient of the biomass in cm^2/s, by default the value for colonies from Harcombe et al.,
Table 1.
glucose_diffusion:
The diffusion coefficient of the glucose in cm^2/s, by default from Harcombe et al., Table 1.
vmax:
The largest glucose uptake in mmol/g/h, by default from Harcombe et al., Table 1.
km:
The half-saturation concentration of glucose in mM, by default from Harcombe et al., Table 1.
interval:
The time in h over which the glucose of a site must last, by default the time step, the bound of Raman
(2021, Section 9.6).
Returns:
system:
The system.
Raises:
MissingExtraError:
The extra is not installed.
Examples:
>>> from skcomplex.systems.biology.ecology import comets
>>> [element.label for element in comets(shape=(2, 2)).elements]
['cells', 'biomass', 'glucose']
"""
source = load_network('textbook', 'cell', needed_by='comets')
plain = {name: replace(spec, rule='') for name, spec in source.reactions.items()}
concentration = glucose / GLUCOSE_MOLAR_MASS * MMOL_PER_MOL
volume = plate_glucose * box**2 / concentration
initial = np.zeros(shape)
initial[shape[0] // 2, shape[1] // 2] = inoculum / volume
space = LatticeSpace(shape, spacing=box)
images = VisualSpec('image', {'states': [MetabolicNetworkElement.flux_name(GROWTH)]})
cells = {
'name': 'ColonyCell',
'reactions': plain,
'default_objective': source.default_objective,
'infeasible': 'zero',
}
return System(
elements=[
GroupElement(
'cells',
element=MetabolicNetworkElement,
values=cells,
on='sites',
space=space,
visual=images,
),
FieldElement('biomass', element=BiomassElement, space=space, values={'x': initial}),
FieldElement(
'glucose',
element=SubstrateElement,
space=space,
values={'c': np.full(shape, concentration)},
),
],
interactions=[
ExchangeInteraction(
'uptake',
cell='cells',
biomass='biomass',
substrate='glucose',
reaction=GLUCOSE_EXCHANGE,
vmax=vmax,
km=km,
interval=interval,
threshold=0.0,
),
GrowthInteraction('growth', cell='cells', biomass='biomass', reaction=GROWTH),
DiffusionInteraction(
'biomass_spread',
target='biomass',
state='x',
coefficient=biomass_diffusion * SECONDS_PER_HOUR,
),
DiffusionInteraction(
'glucose_spread',
target='glucose',
state='c',
coefficient=glucose_diffusion * SECONDS_PER_HOUR,
),
],
space=space,
)
Components
Elements:
MetabolicNetworkElement, fromDeclaredElementBiomassElement, fromElementSubstrateElement, fromElement
Interactions:
ExchangeInteraction, fromInteractionGrowthInteraction, fromInteractionDiffusionInteraction
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
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 DiffusionInteraction(ConservationInteraction):
"""Diffusion `du/dt = D laplacian(u)` of a state `u` of a field, the conservation law of the flux `-D gradient(u)`.
The coefficient `D` is uniform in space, a number or a function of the namespace that returns an expression of
parameters and time. The interaction lowers to the Laplacian of the field, which its discretization evaluates, so
the total of the state on a periodic field stays constant and a point mass spreads as the heat kernel.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state diffuses.
state:
The name of the state of the field that diffuses.
coefficient:
A function of the namespace that returns the diffusion coefficient, or a number.
Raises:
ValidationError:
The state is not a name, or the coefficient is neither callable nor a number.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.components.elements import Element, FieldElement
>>> from skcomplex.core.components.interactions import DiffusionInteraction
>>> from skcomplex.core.components.spaces import LatticeSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((40,), spacing=0.1, boundary='periodic')
>>> field = FieldElement('rod', element=CellElement, space=space, values={'u': np.eye(40)[20]})
>>> heat = DiffusionInteraction('heat', target='rod', state='u', coefficient=0.1)
>>> trajectory = Simulation(System(elements=[field], interactions=[heat]), dt=0.01).prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('rod.u')[-1].sum()), 6)
1.0
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (FluxLaw, RateLaw)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state',))],
state: str = 'u',
coefficient: Callable[[SimpleNamespace], object] | float = 1.0,
) -> None:
super().__init__(label, target=target, state=state)
self.coefficient = coefficient
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a coefficient that is neither callable nor a number.
Raises:
ValidationError:
The type declares states, the state is not a name, or the coefficient is neither callable nor a number.
"""
super()._check()
if not (callable(self.coefficient) or isinstance(self.coefficient, int | float)):
error_msg = f'The coefficient of {self.label!r} is {self.coefficient!r}. Give a function or a number.'
raise ValidationError(error_msg)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate `D laplacian(u)` of the state, the divergence of the diffusive flux.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The rate law of the state.
"""
coefficient = self.coefficient(namespace) if callable(self.coefficient) else self.coefficient
state = getattr(self._get_member(namespace, 'target'), self.state)
return [RateLaw(state, coefficient * laplacian(state))]
Equations
cells.v_ACALD[i]: diff(v_ACALD_cells(t), t) == 0
cells.v_ACALDt[i]: diff(v_ACALDt_cells(t), t) == 0
cells.v_ACKr[i]: diff(v_ACKr_cells(t), t) == 0
cells.v_ACONTa[i]: diff(v_ACONTa_cells(t), t) == 0
cells.v_ACONTb[i]: diff(v_ACONTb_cells(t), t) == 0
cells.v_ACt2r[i]: diff(v_ACt2r_cells(t), t) == 0
cells.v_ADK1[i]: diff(v_ADK1_cells(t), t) == 0
cells.v_AKGDH[i]: diff(v_AKGDH_cells(t), t) == 0
cells.v_AKGt2r[i]: diff(v_AKGt2r_cells(t), t) == 0
cells.v_ALCD2x[i]: diff(v_ALCD2x_cells(t), t) == 0
cells.v_ATPM[i]: diff(v_ATPM_cells(t), t) == 0
cells.v_ATPS4r[i]: diff(v_ATPS4r_cells(t), t) == 0
cells.v_Biomass_Ecoli_core[i]: diff(v_Biomass_Ecoli_core_cells(t), t) == 0
cells.v_CO2t[i]: diff(v_CO2t_cells(t), t) == 0
cells.v_CS[i]: diff(v_CS_cells(t), t) == 0
cells.v_CYTBD[i]: diff(v_CYTBD_cells(t), t) == 0
cells.v_D_LACt2[i]: diff(v_D_LACt2_cells(t), t) == 0
cells.v_ENO[i]: diff(v_ENO_cells(t), t) == 0
cells.v_ETOHt2r[i]: diff(v_ETOHt2r_cells(t), t) == 0
cells.v_EX_ac_e[i]: diff(v_EX_ac_e_cells(t), t) == 0
cells.v_EX_acald_e[i]: diff(v_EX_acald_e_cells(t), t) == 0
cells.v_EX_akg_e[i]: diff(v_EX_akg_e_cells(t), t) == 0
cells.v_EX_co2_e[i]: diff(v_EX_co2_e_cells(t), t) == 0
cells.v_EX_etoh_e[i]: diff(v_EX_etoh_e_cells(t), t) == 0
cells.v_EX_for_e[i]: diff(v_EX_for_e_cells(t), t) == 0
cells.v_EX_fru_e[i]: diff(v_EX_fru_e_cells(t), t) == 0
cells.v_EX_fum_e[i]: diff(v_EX_fum_e_cells(t), t) == 0
cells.v_EX_glc__D_e[i]: diff(v_EX_glc__D_e_cells(t), t) == 0
cells.v_EX_gln__L_e[i]: diff(v_EX_gln__L_e_cells(t), t) == 0
cells.v_EX_glu__L_e[i]: diff(v_EX_glu__L_e_cells(t), t) == 0
cells.v_EX_h_e[i]: diff(v_EX_h_e_cells(t), t) == 0
cells.v_EX_h2o_e[i]: diff(v_EX_h2o_e_cells(t), t) == 0
cells.v_EX_lac__D_e[i]: diff(v_EX_lac__D_e_cells(t), t) == 0
cells.v_EX_mal__L_e[i]: diff(v_EX_mal__L_e_cells(t), t) == 0
cells.v_EX_nh4_e[i]: diff(v_EX_nh4_e_cells(t), t) == 0
cells.v_EX_o2_e[i]: diff(v_EX_o2_e_cells(t), t) == 0
cells.v_EX_pi_e[i]: diff(v_EX_pi_e_cells(t), t) == 0
cells.v_EX_pyr_e[i]: diff(v_EX_pyr_e_cells(t), t) == 0
cells.v_EX_succ_e[i]: diff(v_EX_succ_e_cells(t), t) == 0
cells.v_FBA[i]: diff(v_FBA_cells(t), t) == 0
cells.v_FBP[i]: diff(v_FBP_cells(t), t) == 0
cells.v_FORt2[i]: diff(v_FORt2_cells(t), t) == 0
cells.v_FORti[i]: diff(v_FORti_cells(t), t) == 0
cells.v_FRD7[i]: diff(v_FRD7_cells(t), t) == 0
cells.v_FRUpts2[i]: diff(v_FRUpts2_cells(t), t) == 0
cells.v_FUM[i]: diff(v_FUM_cells(t), t) == 0
cells.v_FUMt2_2[i]: diff(v_FUMt2_2_cells(t), t) == 0
cells.v_G6PDH2r[i]: diff(v_G6PDH2r_cells(t), t) == 0
cells.v_GAPD[i]: diff(v_GAPD_cells(t), t) == 0
cells.v_GLCpts[i]: diff(v_GLCpts_cells(t), t) == 0
cells.v_GLNS[i]: diff(v_GLNS_cells(t), t) == 0
cells.v_GLNabc[i]: diff(v_GLNabc_cells(t), t) == 0
cells.v_GLUDy[i]: diff(v_GLUDy_cells(t), t) == 0
cells.v_GLUN[i]: diff(v_GLUN_cells(t), t) == 0
cells.v_GLUSy[i]: diff(v_GLUSy_cells(t), t) == 0
cells.v_GLUt2r[i]: diff(v_GLUt2r_cells(t), t) == 0
cells.v_GND[i]: diff(v_GND_cells(t), t) == 0
cells.v_H2Ot[i]: diff(v_H2Ot_cells(t), t) == 0
cells.v_ICDHyr[i]: diff(v_ICDHyr_cells(t), t) == 0
cells.v_ICL[i]: diff(v_ICL_cells(t), t) == 0
cells.v_LDH_D[i]: diff(v_LDH_D_cells(t), t) == 0
cells.v_MALS[i]: diff(v_MALS_cells(t), t) == 0
cells.v_MALt2_2[i]: diff(v_MALt2_2_cells(t), t) == 0
cells.v_MDH[i]: diff(v_MDH_cells(t), t) == 0
cells.v_ME1[i]: diff(v_ME1_cells(t), t) == 0
cells.v_ME2[i]: diff(v_ME2_cells(t), t) == 0
cells.v_NADH16[i]: diff(v_NADH16_cells(t), t) == 0
cells.v_NADTRHD[i]: diff(v_NADTRHD_cells(t), t) == 0
cells.v_NH4t[i]: diff(v_NH4t_cells(t), t) == 0
cells.v_O2t[i]: diff(v_O2t_cells(t), t) == 0
cells.v_PDH[i]: diff(v_PDH_cells(t), t) == 0
cells.v_PFK[i]: diff(v_PFK_cells(t), t) == 0
cells.v_PFL[i]: diff(v_PFL_cells(t), t) == 0
cells.v_PGI[i]: diff(v_PGI_cells(t), t) == 0
cells.v_PGK[i]: diff(v_PGK_cells(t), t) == 0
cells.v_PGL[i]: diff(v_PGL_cells(t), t) == 0
cells.v_PGM[i]: diff(v_PGM_cells(t), t) == 0
cells.v_PIt2r[i]: diff(v_PIt2r_cells(t), t) == 0
cells.v_PPC[i]: diff(v_PPC_cells(t), t) == 0
cells.v_PPCK[i]: diff(v_PPCK_cells(t), t) == 0
cells.v_PPS[i]: diff(v_PPS_cells(t), t) == 0
cells.v_PTAr[i]: diff(v_PTAr_cells(t), t) == 0
cells.v_PYK[i]: diff(v_PYK_cells(t), t) == 0
cells.v_PYRt2[i]: diff(v_PYRt2_cells(t), t) == 0
cells.v_RPE[i]: diff(v_RPE_cells(t), t) == 0
cells.v_RPI[i]: diff(v_RPI_cells(t), t) == 0
cells.v_SUCCt2_2[i]: diff(v_SUCCt2_2_cells(t), t) == 0
cells.v_SUCCt3[i]: diff(v_SUCCt3_cells(t), t) == 0
cells.v_SUCDi[i]: diff(v_SUCDi_cells(t), t) == 0
cells.v_SUCOAS[i]: diff(v_SUCOAS_cells(t), t) == 0
cells.v_TALA[i]: diff(v_TALA_cells(t), t) == 0
cells.v_THD2[i]: diff(v_THD2_cells(t), t) == 0
cells.v_TKT1[i]: diff(v_TKT1_cells(t), t) == 0
cells.v_TKT2[i]: diff(v_TKT2_cells(t), t) == 0
cells.v_TPI[i]: diff(v_TPI_cells(t), t) == 0
biomass.x[i]: diff(x_biomass(t), t) == v_Biomass_Ecoli_core_cells(t)*x_biomass(t) + (1.08e-06)*laplacian(x_biomass(t))
glucose.c[i]: diff(c_glucose(t), t) == v_EX_glc__D_e_cells(t)*where(either(c_glucose(t) > 0, v_EX_glc__D_e_cells(t) > 0), 1, 0)*x_biomass(t) + 0.018000000000000002*laplacian(c_glucose(t))