Source code for ecoli.processes.equilibrium

"""
===========
Equilibrium
===========

This process models how ligands are bound to or unbound
from their transcription factor binding partners in a fashion
that maintains equilibrium.
"""

import numpy as np

from ecoli.library.schema import numpy_schema, bulk_name_to_idx, counts, listener_schema
from ecoli.processes.registries import topology_registry
from ecoli.processes.partition import PartitionedProcess

from wholecell.utils import units


# Register default topology for this process, associating it with process name
NAME = "ecoli-equilibrium"
TOPOLOGY = {"listeners": ("listeners",), "bulk": ("bulk",), "timestep": ("timestep",)}
topology_registry.register(NAME, TOPOLOGY)


[docs] class Equilibrium(PartitionedProcess): """Equilibrium PartitionedProcess molecule_names: list of molecules that are being iterated over size:94 """ name = NAME topology = TOPOLOGY defaults = { "jit": False, "n_avogadro": 0.0, "cell_density": 0.0, "stoichMatrix": [[]], "fluxesAndMoleculesToSS": lambda counts, volume, avogadro, random, jit: ( [], [], ), "moleculeNames": [], "seed": 0, "complex_ids": [], "reaction_ids": [], } # Constructor def __init__(self, parameters=None): super().__init__(parameters) # Simulation options # utilized in the fluxes and molecules function self.jit = self.parameters["jit"] # Get constants self.n_avogadro = self.parameters["n_avogadro"] self.cell_density = self.parameters["cell_density"] # Create matrix and method # stoichMatrix: (94, 33), molecule counts are (94,). self.stoichMatrix = self.parameters["stoichMatrix"] # fluxesAndMoleculesToSS: solves ODES to get to steady state based off # of cell density, volumes and molecule counts self.fluxesAndMoleculesToSS = self.parameters["fluxesAndMoleculesToSS"] self.product_indices = [ idx for idx in np.where(np.any(self.stoichMatrix > 0, axis=1))[0] ] # Build views # moleculeNames: list of molecules that are being iterated over size: 94 self.moleculeNames = self.parameters["moleculeNames"] self.molecule_idx = None self.seed = self.parameters["seed"] self.random_state = np.random.RandomState(seed=self.seed) self.complex_ids = self.parameters["complex_ids"] self.reaction_ids = self.parameters["reaction_ids"]
[docs] def ports_schema(self): return { "bulk": numpy_schema("bulk"), "listeners": { "mass": listener_schema({"cell_mass": 0}), "equilibrium_listener": { **listener_schema( { "reaction_rates": ( [0.0] * len(self.reaction_ids), self.reaction_ids, ) } ) }, }, "timestep": {"_default": self.parameters["time_step"]}, }
[docs] def calculate_request(self, timestep, states): # At t=0, convert all strings to indices if self.molecule_idx is None: self.molecule_idx = bulk_name_to_idx( self.moleculeNames, states["bulk"]["id"] ) # Get molecule counts moleculeCounts = counts(states["bulk"], self.molecule_idx) # Get cell mass and volume cellMass = (states["listeners"]["mass"]["cell_mass"] * units.fg).asNumber( units.g ) cellVolume = cellMass / self.cell_density # Solve ODEs to steady state self.rxnFluxes, self.req = self.fluxesAndMoleculesToSS( moleculeCounts, cellVolume, self.n_avogadro, self.random_state, jit=self.jit, ) # Request counts of molecules needed requests = {"bulk": [(self.molecule_idx, self.req.astype(int))]} return requests
[docs] def evolve_state(self, timestep, states): # Get molecule counts moleculeCounts = counts(states["bulk"], self.molecule_idx) # Get counts of molecules allocated to this process rxnFluxes = self.rxnFluxes.copy() # If we didn't get allocated all the molecules we need, make do with # what we have (decrease reaction fluxes so that they make use of what # we have, but not more). Reduces at least one reaction every iteration # so the max number of iterations is the number of reactions that were # originally expected to occur + 1 to reach the break statement. max_iterations = int(np.abs(rxnFluxes).sum()) + 1 for it in range(max_iterations): # Check if any metabolites will have negative counts with current reactions negative_metabolite_idxs = np.where( np.dot(self.stoichMatrix, rxnFluxes) + moleculeCounts < 0 )[0] if len(negative_metabolite_idxs) == 0: break # Reduce reactions that consume metabolites with negative counts limited_rxn_stoich = self.stoichMatrix[negative_metabolite_idxs, :] fwd_rxn_idxs = np.where( np.logical_and(limited_rxn_stoich < 0, rxnFluxes > 0) )[1] rev_rxn_idxs = np.where( np.logical_and(limited_rxn_stoich > 0, rxnFluxes < 0) )[1] rxnFluxes[fwd_rxn_idxs] -= 1 rxnFluxes[rev_rxn_idxs] += 1 rxnFluxes[fwd_rxn_idxs] = np.fmax(0, rxnFluxes[fwd_rxn_idxs]) rxnFluxes[rev_rxn_idxs] = np.fmin(0, rxnFluxes[rev_rxn_idxs]) else: raise ValueError( "Could not get positive counts in equilibrium with" " allocated molecules." ) # Increment changes in molecule counts deltaMolecules = np.dot(self.stoichMatrix, rxnFluxes).astype(int) update = { "bulk": [(self.molecule_idx, deltaMolecules)], "listeners": { "equilibrium_listener": { "reaction_rates": deltaMolecules[self.product_indices] / states["timestep"] } }, } return update
def test_equilibrium_listener(): from ecoli.experiments.ecoli_master_sim import EcoliSim sim = EcoliSim.from_file() sim.total_time = 2 sim.raw_output = False sim.build_ecoli() sim.run() listeners = sim.query()["agents"]["0"]["listeners"] assert isinstance(listeners["equilibrium_listener"]["reaction_rates"][0], list) assert isinstance(listeners["equilibrium_listener"]["reaction_rates"][1], list) if __name__ == "__main__": test_equilibrium_listener()