"""
======================
Polypeptide Elongation
======================
This process models the polymerization of amino acids into polypeptides
by ribosomes using an mRNA transcript as a template. Elongation terminates
once a ribosome has reached the end of an mRNA transcript. Polymerization
occurs across all ribosomes simultaneously and resources are allocated to
maximize the progress of all ribosomes within the limits of the maximum ribosome
elongation rate, available amino acids and GTP, and the length of the transcript.
"""
from typing import Any, Callable, Optional, Tuple
from numba import njit
import numpy as np
import numpy.typing as npt
from scipy.integrate import solve_ivp
from unum import Unum
# wcEcoli imports
from wholecell.utils.polymerize import buildSequences, polymerize, computeMassIncrease
from wholecell.utils.random import stochasticRound
from wholecell.utils import units
# vivarium imports
from vivarium.core.composition import simulate_process
from vivarium.library.dict_utils import deep_merge
from vivarium.library.units import units as vivunits
from vivarium.plots.simulation_output import plot_variables
# vivarium-ecoli imports
from ecoli.library.schema import (
listener_schema,
numpy_schema,
counts,
attrs,
bulk_name_to_idx,
)
from ecoli.processes.registries import topology_registry
from ecoli.processes.partition import PartitionedProcess
MICROMOLAR_UNITS = units.umol / units.L
"""Units used for all concentrations."""
REMOVED_FROM_CHARGING = {"L-SELENOCYSTEINE[c]"}
"""Amino acids to remove from charging when running with
``steady_state_trna_charging``"""
# Register default topology for this process, associating it with process name
NAME = "ecoli-polypeptide-elongation"
TOPOLOGY = {
"environment": ("environment",),
"boundary": ("boundary",),
"listeners": ("listeners",),
"active_ribosome": ("unique", "active_ribosome"),
"bulk": ("bulk",),
"polypeptide_elongation": ("process_state", "polypeptide_elongation"),
# Non-partitioned counts
"bulk_total": ("bulk",),
"timestep": ("timestep",),
}
topology_registry.register(NAME, TOPOLOGY)
DEFAULT_AA_NAMES = [
"L-ALPHA-ALANINE[c]",
"ARG[c]",
"ASN[c]",
"L-ASPARTATE[c]",
"CYS[c]",
"GLT[c]",
"GLN[c]",
"GLY[c]",
"HIS[c]",
"ILE[c]",
"LEU[c]",
"LYS[c]",
"MET[c]",
"PHE[c]",
"PRO[c]",
"SER[c]",
"THR[c]",
"TRP[c]",
"TYR[c]",
"L-SELENOCYSTEINE[c]",
"VAL[c]",
]
[docs]
class PolypeptideElongation(PartitionedProcess):
"""Polypeptide Elongation PartitionedProcess
defaults:
proteinIds: array length n of protein names
"""
name = NAME
topology = TOPOLOGY
defaults = {
"time_step": 1,
"max_time_step": 2.0,
"n_avogadro": 6.02214076e23 / units.mol,
"proteinIds": np.array([]),
"proteinLengths": np.array([]),
"proteinSequences": np.array([[]]),
"aaWeightsIncorporated": np.array([]),
"endWeight": np.array([2.99146113e-08]),
"variable_elongation": False,
"make_elongation_rates": (
lambda random, rate, timestep, variable: np.array([])
),
"next_aa_pad": 1,
"ribosomeElongationRate": 17.388824902723737,
"translation_aa_supply": {"minimal": np.array([])},
"import_threshold": 1e-05,
"aa_from_trna": np.zeros(21),
"gtpPerElongation": 4.2,
"aa_supply_in_charging": False,
"adjust_timestep_for_charging": False,
"mechanistic_translation_supply": False,
"mechanistic_aa_transport": False,
"ppgpp_regulation": False,
"disable_ppgpp_elongation_inhibition": False,
"trna_charging": False,
"translation_supply": False,
"mechanistic_supply": False,
"ribosome30S": "ribosome30S",
"ribosome50S": "ribosome50S",
"amino_acids": DEFAULT_AA_NAMES,
"aa_exchange_names": DEFAULT_AA_NAMES,
"basal_elongation_rate": 22.0,
"ribosomeElongationRateDict": {
"minimal": 17.388824902723737 * units.aa / units.s
},
"uncharged_trna_names": np.array([]),
"aaNames": DEFAULT_AA_NAMES,
"aa_enzymes": [],
"proton": "PROTON",
"water": "H2O",
"cellDensity": 1100 * units.g / units.L,
"elongation_max": 22 * units.aa / units.s,
"aa_from_synthetase": np.array([[]]),
"charging_stoich_matrix": np.array([[]]),
"charged_trna_names": [],
"charging_molecule_names": [],
"synthetase_names": [],
"ppgpp_reaction_names": [],
"ppgpp_reaction_metabolites": [],
"ppgpp_reaction_stoich": np.array([[]]),
"ppgpp_synthesis_reaction": "GDPPYPHOSKIN-RXN",
"ppgpp_degradation_reaction": "PPGPPSYN-RXN",
"aa_importers": [],
"amino_acid_export": None,
"synthesis_index": 0,
"aa_exporters": [],
"get_pathway_enzyme_counts_per_aa": None,
"import_constraint_threshold": 0,
"unit_conversion": 0,
"elong_rate_by_ppgpp": 0,
"amino_acid_import": None,
"degradation_index": 1,
"amino_acid_synthesis": None,
"rela": "RELA",
"spot": "SPOT",
"ppgpp": "ppGpp",
"kS": 100.0,
"KMtf": 1.0,
"KMaa": 100.0,
"krta": 1.0,
"krtf": 500.0,
"KD_RelA": 0.26,
"k_RelA": 75.0,
"k_SpoT_syn": 2.6,
"k_SpoT_deg": 0.23,
"KI_SpoT": 20.0,
"aa_supply_scaling": lambda aa_conc, aa_in_media: 0,
"seed": 0,
"emit_unique": False,
}
def __init__(self, parameters=None):
super().__init__(parameters)
self.max_time_step = self.parameters["max_time_step"]
# Simulation options
self.aa_supply_in_charging = self.parameters["aa_supply_in_charging"]
self.adjust_timestep_for_charging = self.parameters[
"adjust_timestep_for_charging"
]
self.mechanistic_translation_supply = self.parameters[
"mechanistic_translation_supply"
]
self.mechanistic_aa_transport = self.parameters["mechanistic_aa_transport"]
self.ppgpp_regulation = self.parameters["ppgpp_regulation"]
self.disable_ppgpp_elongation_inhibition = self.parameters[
"disable_ppgpp_elongation_inhibition"
]
self.variable_elongation = self.parameters["variable_elongation"]
self.variable_polymerize = self.ppgpp_regulation or self.variable_elongation
translation_supply = self.parameters["translation_supply"]
trna_charging = self.parameters["trna_charging"]
# Load parameters
self.n_avogadro = self.parameters["n_avogadro"]
self.proteinIds = self.parameters["proteinIds"]
self.protein_lengths = self.parameters["proteinLengths"]
self.proteinSequences = self.parameters["proteinSequences"]
self.aaWeightsIncorporated = self.parameters["aaWeightsIncorporated"]
self.endWeight = self.parameters["endWeight"]
self.make_elongation_rates = self.parameters["make_elongation_rates"]
self.next_aa_pad = self.parameters["next_aa_pad"]
self.ribosome30S = self.parameters["ribosome30S"]
self.ribosome50S = self.parameters["ribosome50S"]
self.amino_acids = self.parameters["amino_acids"]
self.aa_exchange_names = self.parameters["aa_exchange_names"]
self.aa_environment_names = [aa[:-3] for aa in self.aa_exchange_names]
self.aa_enzymes = self.parameters["aa_enzymes"]
self.ribosomeElongationRate = self.parameters["ribosomeElongationRate"]
# Amino acid supply calculations
self.translation_aa_supply = self.parameters["translation_aa_supply"]
self.import_threshold = self.parameters["import_threshold"]
# Used for figure in publication
self.trpAIndex = np.where(self.proteinIds == "TRYPSYN-APROTEIN[c]")[0][0]
self.elngRateFactor = 1.0
# Data structures for charging
self.aa_from_trna = self.parameters["aa_from_trna"]
# Set modeling method
# TODO: Test that these models all work properly
if trna_charging:
self.elongation_model = SteadyStateElongationModel(self.parameters, self)
elif translation_supply:
self.elongation_model = TranslationSupplyElongationModel(
self.parameters, self
)
else:
self.elongation_model = BaseElongationModel(self.parameters, self)
# Growth associated maintenance energy requirements for elongations
self.gtpPerElongation = self.parameters["gtpPerElongation"]
# Need to account for ATP hydrolysis for charging that has been
# removed from measured GAM (ATP -> AMP is 2 hydrolysis reactions)
# if charging reactions are not explicitly modeled
if not trna_charging:
self.gtpPerElongation += 2
# basic molecule names
self.proton = self.parameters["proton"]
self.water = self.parameters["water"]
self.rela = self.parameters["rela"]
self.spot = self.parameters["spot"]
self.ppgpp = self.parameters["ppgpp"]
self.aa_importers = self.parameters["aa_importers"]
self.aa_exporters = self.parameters["aa_exporters"]
# Numpy index for bulk molecule
self.proton_idx = None
# Names of molecules associated with tRNA charging
self.ppgpp_reaction_metabolites = self.parameters["ppgpp_reaction_metabolites"]
self.uncharged_trna_names = self.parameters["uncharged_trna_names"]
self.charged_trna_names = self.parameters["charged_trna_names"]
self.charging_molecule_names = self.parameters["charging_molecule_names"]
self.synthetase_names = self.parameters["synthetase_names"]
self.seed = self.parameters["seed"]
self.random_state = np.random.RandomState(seed=self.seed)
self.zero_aa_exchange_rates = (
MICROMOLAR_UNITS / units.s * np.zeros(len(self.amino_acids))
)
[docs]
def ports_schema(self):
return {
"environment": {
"media_id": {"_default": "", "_updater": "set"},
"exchange": {"*": {"_default": 0}},
},
"boundary": {
"external": {
aa: {"_default": 0} for aa in sorted(self.aa_environment_names)
}
},
"listeners": {
"mass": listener_schema({"cell_mass": 0.0, "dry_mass": 0.0}),
"growth_limits": listener_schema(
{
"fraction_trna_charged": (
[0.0] * len(self.uncharged_trna_names),
self.uncharged_trna_names,
),
"aa_allocated": ([0] * len(self.amino_acids), self.amino_acids),
"aa_pool_size": ([0] * len(self.amino_acids), self.amino_acids),
"aa_request_size": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"active_ribosome_allocated": 0,
"net_charged": (
[0] * len(self.uncharged_trna_names),
self.uncharged_trna_names,
),
"aas_used": ([0] * len(self.amino_acids), self.amino_acids),
"aa_count_diff": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
# Below only if trna_charging enbaled
"original_aa_supply": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_in_media": (
[False] * len(self.amino_acids),
self.amino_acids,
),
"synthetase_conc": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"uncharged_trna_conc": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"charged_trna_conc": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_conc": ([0.0] * len(self.amino_acids), self.amino_acids),
"ribosome_conc": 0.0,
"fraction_aa_to_elongate": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_supply": ([0.0] * len(self.amino_acids), self.amino_acids),
"aa_synthesis": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_import": ([0.0] * len(self.amino_acids), self.amino_acids),
"aa_export": ([0.0] * len(self.amino_acids), self.amino_acids),
"aa_importers": (
[0] * len(self.aa_importers),
self.aa_importers,
),
"aa_exporters": (
[0] * len(self.aa_exporters),
self.aa_exporters,
),
"aa_supply_enzymes_fwd": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_supply_enzymes_rev": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_supply_aa_conc": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_supply_fraction_fwd": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"aa_supply_fraction_rev": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"ppgpp_conc": 0.0,
"rela_conc": 0.0,
"spot_conc": 0.0,
"rela_syn": ([0.0] * len(self.amino_acids), self.amino_acids),
"spot_syn": 0.0,
"spot_deg": 0.0,
"spot_deg_inhibited": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"trna_charged": ([0] * len(self.amino_acids), self.amino_acids),
}
),
"ribosome_data": listener_schema(
{
"translation_supply": (
[0.0] * len(self.amino_acids),
self.amino_acids,
),
"effective_elongation_rate": 0.0,
"aa_count_in_sequence": (
[0] * len(self.amino_acids),
self.amino_acids,
),
"aa_counts": ([0.0] * len(self.amino_acids), self.amino_acids),
"actual_elongations": 0,
"actual_elongation_hist": [0] * 22,
"elongations_non_terminating_hist": [0] * 22,
"did_terminate": 0,
"termination_loss": 0,
"num_trpA_terminated": 0,
"process_elongation_rate": 0.0,
}
),
},
"bulk": numpy_schema("bulk"),
"bulk_total": numpy_schema("bulk"),
"active_ribosome": numpy_schema(
"active_ribosome", emit=self.parameters["emit_unique"]
),
"polypeptide_elongation": {
"aa_count_diff": {
"_default": {},
"_emit": True,
"_updater": "set",
"_divider": "empty_dict",
},
"gtp_to_hydrolyze": {
"_default": 0.0,
"_emit": True,
"_updater": "set",
"_divider": "zero",
},
"aa_exchange_rates": {
"_default": [0.0],
"_emit": True,
"_updater": "set",
"_divider": "zero",
},
},
"timestep": {"_default": self.parameters["time_step"]},
}
[docs]
def calculate_request(self, timestep, states):
"""
Set ribosome elongation rate based on simulation medium environment and elongation rate factor
which is used to create single-cell variability in growth rate
The maximum number of amino acids that can be elongated in a single timestep is set to 22
intentionally as the minimum number of padding values on the protein sequence matrix is set to 22.
If timesteps longer than 1.0s are used, this feature will lead to errors in the effective ribosome
elongation rate.
"""
if self.proton_idx is None:
bulk_ids = states["bulk"]["id"]
self.proton_idx = bulk_name_to_idx(self.proton, bulk_ids)
self.water_idx = bulk_name_to_idx(self.water, bulk_ids)
self.rela_idx = bulk_name_to_idx(self.rela, bulk_ids)
self.spot_idx = bulk_name_to_idx(self.spot, bulk_ids)
self.ppgpp_idx = bulk_name_to_idx(self.ppgpp, bulk_ids)
self.monomer_idx = bulk_name_to_idx(self.proteinIds, bulk_ids)
self.amino_acid_idx = bulk_name_to_idx(self.amino_acids, bulk_ids)
self.aa_enzyme_idx = bulk_name_to_idx(self.aa_enzymes, bulk_ids)
self.ppgpp_rxn_metabolites_idx = bulk_name_to_idx(
self.ppgpp_reaction_metabolites, bulk_ids
)
self.uncharged_trna_idx = bulk_name_to_idx(
self.uncharged_trna_names, bulk_ids
)
self.charged_trna_idx = bulk_name_to_idx(self.charged_trna_names, bulk_ids)
self.charging_molecule_idx = bulk_name_to_idx(
self.charging_molecule_names, bulk_ids
)
self.synthetase_idx = bulk_name_to_idx(self.synthetase_names, bulk_ids)
self.ribosome30S_idx = bulk_name_to_idx(self.ribosome30S, bulk_ids)
self.ribosome50S_idx = bulk_name_to_idx(self.ribosome50S, bulk_ids)
self.aa_importer_idx = bulk_name_to_idx(self.aa_importers, bulk_ids)
self.aa_exporter_idx = bulk_name_to_idx(self.aa_exporters, bulk_ids)
# MODEL SPECIFIC: get ribosome elongation rate
self.ribosomeElongationRate = self.elongation_model.elongation_rate(states)
# If there are no active ribosomes, return immediately
if states["active_ribosome"]["_entryState"].sum() == 0:
return {"listeners": {"ribosome_data": {}, "growth_limits": {}}}
# Build sequences to request appropriate amount of amino acids to
# polymerize for next timestep
(
proteinIndexes,
peptideLengths,
) = attrs(states["active_ribosome"], ["protein_index", "peptide_length"])
self.elongation_rates = self.make_elongation_rates(
self.random_state,
self.ribosomeElongationRate,
states["timestep"],
self.variable_elongation,
)
sequences = buildSequences(
self.proteinSequences, proteinIndexes, peptideLengths, self.elongation_rates
)
sequenceHasAA = sequences != polymerize.PAD_VALUE
aasInSequences = np.bincount(sequences[sequenceHasAA], minlength=21)
# Calculate AA supply for expected doubling of protein
dryMass = states["listeners"]["mass"]["dry_mass"] * units.fg
current_media_id = states["environment"]["media_id"]
translation_supply_rate = (
self.translation_aa_supply[current_media_id] * self.elngRateFactor
)
mol_aas_supplied = (
translation_supply_rate * dryMass * states["timestep"] * units.s
)
self.aa_supply = units.strip_empty_units(mol_aas_supplied * self.n_avogadro)
# MODEL SPECIFIC: Calculate AA request
fraction_charged, aa_counts_for_translation, requests = (
self.elongation_model.request(states, aasInSequences)
)
# Write to listeners
listeners = requests.setdefault("listeners", {})
ribosome_data_listener = listeners.setdefault("ribosome_data", {})
ribosome_data_listener["translation_supply"] = (
translation_supply_rate.asNumber()
)
growth_limits_listener = requests["listeners"].setdefault("growth_limits", {})
growth_limits_listener["fraction_trna_charged"] = np.dot(
fraction_charged, self.aa_from_trna
)
growth_limits_listener["aa_pool_size"] = counts(
states["bulk_total"], self.amino_acid_idx
)
growth_limits_listener["aa_request_size"] = aa_counts_for_translation
return requests
[docs]
def evolve_state(self, timestep, states):
"""
Set ribosome elongation rate based on simulation medium environment and elongation rate factor
which is used to create single-cell variability in growth rate
The maximum number of amino acids that can be elongated in a single timestep is set to 22
intentionally as the minimum number of padding values on the protein sequence matrix is set to 22.
If timesteps longer than 1.0s are used, this feature will lead to errors in the effective ribosome
elongation rate.
"""
update = {
"listeners": {"ribosome_data": {}, "growth_limits": {}},
"polypeptide_elongation": {},
"active_ribosome": {},
"bulk": [],
}
# Begin wcEcoli evolveState()
# Set values for metabolism in case of early return
update["polypeptide_elongation"]["gtp_to_hydrolyze"] = 0
update["polypeptide_elongation"]["aa_count_diff"] = {}
# Get number of active ribosomes
n_active_ribosomes = states["active_ribosome"]["_entryState"].sum()
update["listeners"]["growth_limits"]["active_ribosome_allocated"] = (
n_active_ribosomes
)
update["listeners"]["growth_limits"]["aa_allocated"] = counts(
states["bulk"], self.amino_acid_idx
)
# If there are no active ribosomes, return immediately
if n_active_ribosomes == 0:
return update
# Polypeptide elongation requires counts to be updated in real-time
# so make a writeable copy of bulk counts to do so
states["bulk"] = counts(states["bulk"], range(len(states["bulk"])))
# Build amino acids sequences for each ribosome to polymerize
protein_indexes, peptide_lengths, positions_on_mRNA = attrs(
states["active_ribosome"],
["protein_index", "peptide_length", "pos_on_mRNA"],
)
all_sequences = buildSequences(
self.proteinSequences,
protein_indexes,
peptide_lengths,
self.elongation_rates + self.next_aa_pad,
)
sequences = all_sequences[:, : -self.next_aa_pad].copy()
if sequences.size == 0:
return update
# Calculate elongation resource capacity
aaCountInSequence = np.bincount(sequences[(sequences != polymerize.PAD_VALUE)])
total_aa_counts = counts(states["bulk"], self.amino_acid_idx)
charged_trna_counts = counts(states["bulk"], self.charged_trna_idx)
# MODEL SPECIFIC: Get amino acid counts
aa_counts_for_translation = self.elongation_model.final_amino_acids(
total_aa_counts, charged_trna_counts
)
# Using polymerization algorithm elongate each ribosome up to the limits
# of amino acids, sequence, and GTP
result = polymerize(
sequences,
aa_counts_for_translation,
10000000, # Set to a large number, the limit is now taken care of in metabolism
self.random_state,
self.elongation_rates[protein_indexes],
variable_elongation=self.variable_polymerize,
)
sequence_elongations = result.sequenceElongation
aas_used = result.monomerUsages
nElongations = result.nReactions
next_amino_acid = all_sequences[
np.arange(len(sequence_elongations)), sequence_elongations
]
next_amino_acid_count = np.bincount(
next_amino_acid[next_amino_acid != polymerize.PAD_VALUE], minlength=21
)
# Update masses of ribosomes attached to polymerizing polypeptides
added_protein_mass = computeMassIncrease(
sequences, sequence_elongations, self.aaWeightsIncorporated
)
updated_lengths = peptide_lengths + sequence_elongations
updated_positions_on_mRNA = positions_on_mRNA + 3 * sequence_elongations
didInitialize = (sequence_elongations > 0) & (peptide_lengths == 0)
added_protein_mass[didInitialize] += self.endWeight
# Write current average elongation to listener
currElongRate = (sequence_elongations.sum() / n_active_ribosomes) / states[
"timestep"
]
# Ribosomes that reach the end of their sequences are terminated and
# dissociated into 30S and 50S subunits. The polypeptide that they are
# polymerizing is converted into a protein in BulkMolecules
terminalLengths = self.protein_lengths[protein_indexes]
didTerminate = updated_lengths == terminalLengths
terminatedProteins = np.bincount(
protein_indexes[didTerminate], minlength=self.proteinSequences.shape[0]
)
(protein_mass,) = attrs(states["active_ribosome"], ["massDiff_protein"])
update["active_ribosome"].update(
{
"delete": np.where(didTerminate)[0],
"set": {
"massDiff_protein": protein_mass + added_protein_mass,
"peptide_length": updated_lengths,
"pos_on_mRNA": updated_positions_on_mRNA,
},
}
)
update["bulk"].append((self.monomer_idx, terminatedProteins))
states["bulk"][self.monomer_idx] += terminatedProteins
nTerminated = didTerminate.sum()
nInitialized = didInitialize.sum()
update["bulk"].append((self.ribosome30S_idx, nTerminated))
update["bulk"].append((self.ribosome50S_idx, nTerminated))
states["bulk"][self.ribosome30S_idx] += nTerminated
states["bulk"][self.ribosome50S_idx] += nTerminated
# MODEL SPECIFIC: evolve
net_charged, aa_count_diff, evolve_update = self.elongation_model.evolve(
states,
total_aa_counts,
aas_used,
next_amino_acid_count,
nElongations,
nInitialized,
)
evolve_bulk_update = evolve_update.pop("bulk")
update = deep_merge(update, evolve_update)
update["bulk"].extend(evolve_bulk_update)
update["polypeptide_elongation"]["aa_count_diff"] = aa_count_diff
# GTP hydrolysis is carried out in Metabolism process for growth
# associated maintenance. This is passed to metabolism.
update["polypeptide_elongation"]["gtp_to_hydrolyze"] = (
self.gtpPerElongation * nElongations
)
# Write data to listeners
update["listeners"]["growth_limits"]["net_charged"] = net_charged
update["listeners"]["growth_limits"]["aas_used"] = aas_used
update["listeners"]["growth_limits"]["aa_count_diff"] = [
aa_count_diff.get(id_, 0) for id_ in self.amino_acids
]
ribosome_data_listener = update["listeners"].setdefault("ribosome_data", {})
ribosome_data_listener["effective_elongation_rate"] = currElongRate
ribosome_data_listener["aa_count_in_sequence"] = aaCountInSequence
ribosome_data_listener["aa_counts"] = aa_counts_for_translation
ribosome_data_listener["actual_elongations"] = sequence_elongations.sum()
ribosome_data_listener["actual_elongation_hist"] = np.histogram(
sequence_elongations, bins=np.arange(0, 23)
)[0]
ribosome_data_listener["elongations_non_terminating_hist"] = np.histogram(
sequence_elongations[~didTerminate], bins=np.arange(0, 23)
)[0]
ribosome_data_listener["did_terminate"] = didTerminate.sum()
ribosome_data_listener["termination_loss"] = (
terminalLengths - peptide_lengths
)[didTerminate].sum()
ribosome_data_listener["num_trpA_terminated"] = terminatedProteins[
self.trpAIndex
]
ribosome_data_listener["process_elongation_rate"] = (
self.ribosomeElongationRate / states["timestep"]
)
return update
[docs]
def isTimeStepShortEnough(self, inputTimeStep, timeStepSafetyFraction):
model_specific = self.elongation_model.isTimeStepShortEnough(
inputTimeStep, timeStepSafetyFraction
)
max_time_step = inputTimeStep <= self.max_time_step
return model_specific and max_time_step
[docs]
class BaseElongationModel(object):
"""
Base Model: Request amino acids according to upcoming sequence, assuming
max ribosome elongation.
"""
def __init__(self, parameters, process):
self.parameters = parameters
self.process = process
self.basal_elongation_rate = self.parameters["basal_elongation_rate"]
self.ribosomeElongationRateDict = self.parameters["ribosomeElongationRateDict"]
[docs]
def elongation_rate(self, states):
"""
Sets ribosome elongation rate accordint to the media; returns
max value of 22 amino acids/second.
"""
current_media_id = states["environment"]["media_id"]
rate = self.process.elngRateFactor * self.ribosomeElongationRateDict[
current_media_id
].asNumber(units.aa / units.s)
return np.min([self.basal_elongation_rate, rate])
[docs]
def amino_acid_counts(self, aasInSequences):
return aasInSequences
[docs]
def request(self, states, aasInSequences):
aa_counts_for_translation = self.amino_acid_counts(aasInSequences)
requests = {"bulk": [(self.process.amino_acid_idx, aa_counts_for_translation)]}
# Not modeling charging so set fraction charged to 0 for all tRNA
fraction_charged = np.zeros(len(self.process.amino_acid_idx))
return fraction_charged, aa_counts_for_translation, requests
[docs]
def final_amino_acids(self, total_aa_counts, charged_trna_counts):
return total_aa_counts
[docs]
def evolve(
self,
states,
total_aa_counts,
aas_used,
next_amino_acid_count,
nElongations,
nInitialized,
):
# Update counts of amino acids and water to reflect polymerization
# reactions
net_charged = np.zeros(len(self.parameters["uncharged_trna_names"]))
return (
net_charged,
{},
{
"bulk": [
(self.process.amino_acid_idx, -aas_used),
(self.process.water_idx, nElongations - nInitialized),
]
},
)
[docs]
def isTimeStepShortEnough(self, inputTimeStep, timeStepSafetyFraction):
return True
[docs]
class TranslationSupplyElongationModel(BaseElongationModel):
"""
Translation Supply Model: Requests minimum of 1) upcoming amino acid
sequence assuming max ribosome elongation (ie. Base Model) and 2)
estimation based on doubling the proteome in one cell cycle (does not
use ribosome elongation, computed in Parca).
"""
def __init__(self, parameters, process):
super().__init__(parameters, process)
[docs]
def elongation_rate(self, states):
"""
Sets ribosome elongation rate accordint to the media; returns
max value of 22 amino acids/second.
"""
return self.basal_elongation_rate
[docs]
def amino_acid_counts(self, aasInSequences):
# Check if this is required. It is a better request but there may be
# fewer elongations.
return np.fmin(self.process.aa_supply, aasInSequences)
[docs]
class SteadyStateElongationModel(TranslationSupplyElongationModel):
"""
Steady State Charging Model: Requests amino acids based on the
Michaelis-Menten competitive inhibition model.
"""
def __init__(self, parameters, process):
super().__init__(parameters, process)
# Cell parameters
self.cellDensity = self.parameters["cellDensity"]
# Names of molecules associated with tRNA charging
self.charged_trna_names = self.parameters["charged_trna_names"]
self.charging_molecule_names = self.parameters["charging_molecule_names"]
self.synthetase_names = self.parameters["synthetase_names"]
# Data structures for charging
self.aa_from_synthetase = self.parameters["aa_from_synthetase"]
self.charging_stoich_matrix = self.parameters["charging_stoich_matrix"]
self.charging_molecules_not_aa = np.array(
[
mol not in set(self.parameters["amino_acids"])
for mol in self.charging_molecule_names
]
)
# ppGpp synthesis
self.ppgpp_reaction_metabolites = self.parameters["ppgpp_reaction_metabolites"]
self.elong_rate_by_ppgpp = self.parameters["elong_rate_by_ppgpp"]
# Parameters for tRNA charging, ribosome elongation and ppGpp reactions
self.charging_params = {
"kS": self.parameters["kS"],
"KMaa": self.parameters["KMaa"],
"KMtf": self.parameters["KMtf"],
"krta": self.parameters["krta"],
"krtf": self.parameters["krtf"],
"max_elong_rate": float(
self.parameters["elongation_max"].asNumber(units.aa / units.s)
),
"charging_mask": np.array(
[
aa not in REMOVED_FROM_CHARGING
for aa in self.parameters["amino_acids"]
]
),
"unit_conversion": self.parameters["unit_conversion"],
}
self.ppgpp_params = {
"KD_RelA": self.parameters["KD_RelA"],
"k_RelA": self.parameters["k_RelA"],
"k_SpoT_syn": self.parameters["k_SpoT_syn"],
"k_SpoT_deg": self.parameters["k_SpoT_deg"],
"KI_SpoT": self.parameters["KI_SpoT"],
"ppgpp_reaction_stoich": self.parameters["ppgpp_reaction_stoich"],
"synthesis_index": self.parameters["synthesis_index"],
"degradation_index": self.parameters["degradation_index"],
}
# Amino acid supply calculations
self.aa_supply_scaling = self.parameters["aa_supply_scaling"]
# Manage unstable charging with too long time step by setting
# time_step_short_enough to False during updates. Other variables
# manage when to trigger an adjustment and how quickly the time step
# increases after being reduced
self.time_step_short_enough = True
self.max_time_step = self.process.max_time_step
self.time_step_increase = 1.01
self.max_amino_acid_adjustment = 0.05
self.amino_acid_synthesis = self.parameters["amino_acid_synthesis"]
self.amino_acid_import = self.parameters["amino_acid_import"]
self.amino_acid_export = self.parameters["amino_acid_export"]
self.get_pathway_enzyme_counts_per_aa = self.parameters[
"get_pathway_enzyme_counts_per_aa"
]
# Comparing two values with units is faster than converting units
# and comparing magnitudes
self.import_constraint_threshold = (
self.parameters["import_constraint_threshold"] * vivunits.mM
)
[docs]
def elongation_rate(self, states):
if (
self.process.ppgpp_regulation
and not self.process.disable_ppgpp_elongation_inhibition
):
cell_mass = states["listeners"]["mass"]["cell_mass"] * units.fg
cell_volume = cell_mass / self.cellDensity
counts_to_molar = 1 / (self.process.n_avogadro * cell_volume)
ppgpp_count = counts(states["bulk"], self.process.ppgpp_idx)
ppgpp_conc = ppgpp_count * counts_to_molar
rate = self.elong_rate_by_ppgpp(
ppgpp_conc, self.basal_elongation_rate
).asNumber(units.aa / units.s)
else:
rate = super().elongation_rate(states)
return rate
[docs]
def request(self, states, aasInSequences):
self.max_time_step = min(
self.process.max_time_step, self.max_time_step * self.time_step_increase
)
# Conversion from counts to molarity
cell_mass = states["listeners"]["mass"]["cell_mass"] * units.fg
dry_mass = states["listeners"]["mass"]["dry_mass"] * units.fg
cell_volume = cell_mass / self.cellDensity
self.counts_to_molar = 1 / (self.process.n_avogadro * cell_volume)
# ppGpp related concentrations
ppgpp_conc = self.counts_to_molar * counts(
states["bulk_total"], self.process.ppgpp_idx
)
rela_conc = self.counts_to_molar * counts(
states["bulk_total"], self.process.rela_idx
)
spot_conc = self.counts_to_molar * counts(
states["bulk_total"], self.process.spot_idx
)
# Get counts and convert synthetase and tRNA to a per AA basis
synthetase_counts = np.dot(
self.aa_from_synthetase,
counts(states["bulk_total"], self.process.synthetase_idx),
)
aa_counts = counts(states["bulk_total"], self.process.amino_acid_idx)
uncharged_trna_array = counts(
states["bulk_total"], self.process.uncharged_trna_idx
)
charged_trna_array = counts(states["bulk_total"], self.process.charged_trna_idx)
uncharged_trna_counts = np.dot(self.process.aa_from_trna, uncharged_trna_array)
charged_trna_counts = np.dot(self.process.aa_from_trna, charged_trna_array)
ribosome_counts = states["active_ribosome"]["_entryState"].sum()
# Get concentration
f = aasInSequences / aasInSequences.sum()
synthetase_conc = self.counts_to_molar * synthetase_counts
aa_conc = self.counts_to_molar * aa_counts
uncharged_trna_conc = self.counts_to_molar * uncharged_trna_counts
charged_trna_conc = self.counts_to_molar * charged_trna_counts
ribosome_conc = self.counts_to_molar * ribosome_counts
# Calculate amino acid supply
aa_in_media = np.array(
[
states["boundary"]["external"][aa] > self.import_constraint_threshold
for aa in self.process.aa_environment_names
]
)
fwd_enzyme_counts, rev_enzyme_counts = self.get_pathway_enzyme_counts_per_aa(
counts(states["bulk_total"], self.process.aa_enzyme_idx)
)
importer_counts = counts(states["bulk_total"], self.process.aa_importer_idx)
exporter_counts = counts(states["bulk_total"], self.process.aa_exporter_idx)
synthesis, fwd_saturation, rev_saturation = self.amino_acid_synthesis(
fwd_enzyme_counts, rev_enzyme_counts, aa_conc
)
import_rates = self.amino_acid_import(
aa_in_media,
dry_mass,
aa_conc,
importer_counts,
self.process.mechanistic_aa_transport,
)
export_rates = self.amino_acid_export(
exporter_counts, aa_conc, self.process.mechanistic_aa_transport
)
exchange_rates = import_rates - export_rates
supply_function = get_charging_supply_function(
self.process.aa_supply_in_charging,
self.process.mechanistic_translation_supply,
self.process.mechanistic_aa_transport,
self.amino_acid_synthesis,
self.amino_acid_import,
self.amino_acid_export,
self.aa_supply_scaling,
self.counts_to_molar,
self.process.aa_supply,
fwd_enzyme_counts,
rev_enzyme_counts,
dry_mass,
importer_counts,
exporter_counts,
aa_in_media,
)
# Calculate steady state tRNA levels and resulting elongation rate
self.charging_params["max_elong_rate"] = self.elongation_rate(states)
(
fraction_charged,
v_rib,
synthesis_in_charging,
import_in_charging,
export_in_charging,
) = calculate_trna_charging(
synthetase_conc,
uncharged_trna_conc,
charged_trna_conc,
aa_conc,
ribosome_conc,
f,
self.charging_params,
supply=supply_function,
limit_v_rib=True,
time_limit=states["timestep"],
)
# Use the supply calculated from each sub timestep while solving the charging steady state
if self.process.aa_supply_in_charging:
conversion = (
1 / self.counts_to_molar.asNumber(MICROMOLAR_UNITS) / states["timestep"]
)
synthesis = conversion * synthesis_in_charging
import_rates = conversion * import_in_charging
export_rates = conversion * export_in_charging
self.process.aa_supply = synthesis + import_rates - export_rates
# Use the supply calculated from the starting amino acid concentrations only
elif self.process.mechanistic_translation_supply:
# Set supply based on mechanistic synthesis and supply
self.process.aa_supply = states["timestep"] * (synthesis + exchange_rates)
else:
# Adjust aa_supply higher if amino acid concentrations are low
# Improves stability of charging and mimics amino acid synthesis
# inhibition and export
self.process.aa_supply *= self.aa_supply_scaling(aa_conc, aa_in_media)
aa_counts_for_translation = (
v_rib
* f
* states["timestep"]
/ self.counts_to_molar.asNumber(MICROMOLAR_UNITS)
)
total_trna = charged_trna_array + uncharged_trna_array
final_charged_trna = stochasticRound(
self.process.random_state,
np.dot(fraction_charged, self.process.aa_from_trna * total_trna),
)
# Request charged tRNA that will become uncharged
charged_trna_request = charged_trna_array - final_charged_trna
charged_trna_request[charged_trna_request < 0] = 0
uncharged_trna_request = final_charged_trna - charged_trna_array
uncharged_trna_request[uncharged_trna_request < 0] = 0
self.uncharged_trna_to_charge = uncharged_trna_request
self.aa_counts_for_translation = np.array(aa_counts_for_translation)
fraction_trna_per_aa = total_trna / np.dot(
np.dot(self.process.aa_from_trna, total_trna), self.process.aa_from_trna
)
total_charging_reactions = stochasticRound(
self.process.random_state,
np.dot(aa_counts_for_translation, self.process.aa_from_trna)
* fraction_trna_per_aa
+ uncharged_trna_request,
)
# Only request molecules that will be consumed in the charging reactions
aa_from_uncharging = -self.charging_stoich_matrix @ charged_trna_request
aa_from_uncharging[self.charging_molecules_not_aa] = 0
requested_molecules = (
-np.dot(self.charging_stoich_matrix, total_charging_reactions)
- aa_from_uncharging
)
requested_molecules[requested_molecules < 0] = 0
self.uncharged_trna_to_charge = uncharged_trna_request
# ppGpp reactions based on charged tRNA
request_ppgpp_metabolites = np.zeros(
len(self.process.ppgpp_reaction_metabolites)
)
if self.process.ppgpp_regulation:
total_trna_conc = self.counts_to_molar * (
uncharged_trna_counts + charged_trna_counts
)
updated_charged_trna_conc = total_trna_conc * fraction_charged
updated_uncharged_trna_conc = total_trna_conc - updated_charged_trna_conc
delta_metabolites, *_ = ppgpp_metabolite_changes(
updated_uncharged_trna_conc,
updated_charged_trna_conc,
ribosome_conc,
f,
rela_conc,
spot_conc,
ppgpp_conc,
self.counts_to_molar,
v_rib,
self.charging_params,
self.ppgpp_params,
states["timestep"],
request=True,
random_state=self.process.random_state,
)
request_ppgpp_metabolites = -delta_metabolites
ppgpp_request = counts(states["bulk"], self.process.ppgpp_idx)
return (
fraction_charged,
aa_counts_for_translation,
{
"bulk": [
(
self.process.charging_molecule_idx,
requested_molecules.astype(int),
),
(self.process.charged_trna_idx, charged_trna_request.astype(int)),
# Request water for transfer of AA from tRNA for initial polypeptide.
# This is severe overestimate assuming the worst case that every
# elongation is initializing a polypeptide. This excess of water
# shouldn't matter though.
(self.process.water_idx, int(aa_counts_for_translation.sum())),
(self.process.ppgpp_idx, ppgpp_request),
(
self.process.ppgpp_rxn_metabolites_idx,
request_ppgpp_metabolites.astype(int),
),
],
"listeners": {
"growth_limits": {
"original_aa_supply": self.process.aa_supply,
"aa_in_media": aa_in_media,
"synthetase_conc": synthetase_conc.asNumber(MICROMOLAR_UNITS),
"uncharged_trna_conc": uncharged_trna_conc.asNumber(
MICROMOLAR_UNITS
),
"charged_trna_conc": charged_trna_conc.asNumber(
MICROMOLAR_UNITS
),
"aa_conc": aa_conc.asNumber(MICROMOLAR_UNITS),
"ribosome_conc": ribosome_conc.asNumber(MICROMOLAR_UNITS),
"fraction_aa_to_elongate": f,
"aa_supply": self.process.aa_supply,
"aa_synthesis": synthesis * states["timestep"],
"aa_import": import_rates * states["timestep"],
"aa_export": export_rates * states["timestep"],
"aa_supply_enzymes_fwd": fwd_enzyme_counts,
"aa_supply_enzymes_rev": rev_enzyme_counts,
"aa_importers": importer_counts,
"aa_exporters": exporter_counts,
"aa_supply_aa_conc": aa_conc.asNumber(units.mmol / units.L),
"aa_supply_fraction_fwd": fwd_saturation,
"aa_supply_fraction_rev": rev_saturation,
"ppgpp_conc": ppgpp_conc.asNumber(MICROMOLAR_UNITS),
"rela_conc": rela_conc.asNumber(MICROMOLAR_UNITS),
"spot_conc": spot_conc.asNumber(MICROMOLAR_UNITS),
}
},
"polypeptide_elongation": {
"aa_exchange_rates": self.counts_to_molar
/ units.s
* (import_rates - export_rates)
},
},
)
[docs]
def final_amino_acids(self, total_aa_counts, charged_trna_counts):
charged_counts_to_uncharge = self.process.aa_from_trna @ charged_trna_counts
return np.fmin(
total_aa_counts + charged_counts_to_uncharge, self.aa_counts_for_translation
)
[docs]
def evolve(
self,
states,
total_aa_counts,
aas_used,
next_amino_acid_count,
nElongations,
nInitialized,
):
update = {
"bulk": [],
"listeners": {},
}
# Get tRNA counts
uncharged_trna = counts(states["bulk"], self.process.uncharged_trna_idx)
charged_trna = counts(states["bulk"], self.process.charged_trna_idx)
total_trna = uncharged_trna + charged_trna
# Adjust molecules for number of charging reactions that occurred
## Determine limitations for charging and uncharging reactions
charged_and_elongated_per_aa = np.fmax(
0, (aas_used - self.process.aa_from_trna @ charged_trna)
)
aa_for_charging = total_aa_counts - charged_and_elongated_per_aa
n_aa_charged = np.fmin(
aa_for_charging,
np.dot(
self.process.aa_from_trna,
np.fmin(self.uncharged_trna_to_charge, uncharged_trna),
),
)
n_uncharged_per_aa = aas_used - charged_and_elongated_per_aa
## Calculate changes in tRNA based on limitations
n_trna_charged = self.distribution_from_aa(n_aa_charged, uncharged_trna, True)
n_trna_uncharged = self.distribution_from_aa(
n_uncharged_per_aa, charged_trna, True
)
## Determine reactions that are charged and elongated in same time step without changing
## charged or uncharged counts
charged_and_elongated = self.distribution_from_aa(
charged_and_elongated_per_aa, total_trna
)
## Determine total number of reactions that occur
total_uncharging_reactions = charged_and_elongated + n_trna_uncharged
total_charging_reactions = charged_and_elongated + n_trna_charged
net_charged = total_charging_reactions - total_uncharging_reactions
charging_mol_delta = np.dot(
self.charging_stoich_matrix, total_charging_reactions
).astype(int)
update["bulk"].append((self.process.charging_molecule_idx, charging_mol_delta))
states["bulk"][self.process.charging_molecule_idx] += charging_mol_delta
## Account for uncharging of tRNA during elongation
update["bulk"].append(
(self.process.charged_trna_idx, -total_uncharging_reactions)
)
update["bulk"].append(
(self.process.uncharged_trna_idx, total_uncharging_reactions)
)
states["bulk"][self.process.charged_trna_idx] += -total_uncharging_reactions
states["bulk"][self.process.uncharged_trna_idx] += total_uncharging_reactions
# Update proton counts to reflect polymerization reactions and transfer of AA from tRNA
# Peptide bond formation releases a water but transferring AA from tRNA consumes a OH-
# Net production of H+ for each elongation, consume extra water for each initialization
# since a peptide bond doesn't form
update["bulk"].append((self.process.proton_idx, nElongations))
update["bulk"].append((self.process.water_idx, -nInitialized))
states["bulk"][self.process.proton_idx] += nElongations
states["bulk"][self.process.water_idx] += -nInitialized
# Create or degrade ppGpp
# This should come after all countInc/countDec calls since it shares some molecules with
# other views and those counts should be updated to get the proper limits on ppGpp reactions
if self.process.ppgpp_regulation:
v_rib = (nElongations * self.counts_to_molar).asNumber(
MICROMOLAR_UNITS
) / states["timestep"]
ribosome_conc = (
self.counts_to_molar * states["active_ribosome"]["_entryState"].sum()
)
updated_uncharged_trna_counts = (
counts(states["bulk_total"], self.process.uncharged_trna_idx)
- net_charged
)
updated_charged_trna_counts = (
counts(states["bulk_total"], self.process.charged_trna_idx)
+ net_charged
)
uncharged_trna_conc = self.counts_to_molar * np.dot(
self.process.aa_from_trna, updated_uncharged_trna_counts
)
charged_trna_conc = self.counts_to_molar * np.dot(
self.process.aa_from_trna, updated_charged_trna_counts
)
ppgpp_conc = self.counts_to_molar * counts(
states["bulk_total"], self.process.ppgpp_idx
)
rela_conc = self.counts_to_molar * counts(
states["bulk_total"], self.process.rela_idx
)
spot_conc = self.counts_to_molar * counts(
states["bulk_total"], self.process.spot_idx
)
# Need to include the next amino acid the ribosome sees for certain
# cases where elongation does not occur, otherwise f will be NaN
aa_at_ribosome = aas_used + next_amino_acid_count
f = aa_at_ribosome / aa_at_ribosome.sum()
limits = counts(states["bulk"], self.process.ppgpp_rxn_metabolites_idx)
(
delta_metabolites,
ppgpp_syn,
ppgpp_deg,
rela_syn,
spot_syn,
spot_deg,
spot_deg_inhibited,
) = ppgpp_metabolite_changes(
uncharged_trna_conc,
charged_trna_conc,
ribosome_conc,
f,
rela_conc,
spot_conc,
ppgpp_conc,
self.counts_to_molar,
v_rib,
self.charging_params,
self.ppgpp_params,
states["timestep"],
random_state=self.process.random_state,
limits=limits,
)
update["listeners"]["growth_limits"] = {
"rela_syn": rela_syn,
"spot_syn": spot_syn,
"spot_deg": spot_deg,
"spot_deg_inhibited": spot_deg_inhibited,
}
update["bulk"].append(
(self.process.ppgpp_rxn_metabolites_idx, delta_metabolites.astype(int))
)
states["bulk"][self.process.ppgpp_rxn_metabolites_idx] += (
delta_metabolites.astype(int)
)
# Use the difference between (expected AA supply based on expected
# doubling time and current DCW) and AA used to charge tRNA to update
# the concentration target in metabolism during the next time step
aa_used_trna = np.dot(self.process.aa_from_trna, total_charging_reactions)
aa_diff = self.process.aa_supply - aa_used_trna
if np.any(
np.abs(aa_diff / counts(states["bulk_total"], self.process.amino_acid_idx))
> self.max_amino_acid_adjustment
):
self.time_step_short_enough = False
update["listeners"]["growth_limits"]["trna_charged"] = aa_used_trna.astype(int)
return (
net_charged,
{aa: diff for aa, diff in zip(self.process.amino_acids, aa_diff)},
update,
)
[docs]
def distribution_from_aa(
self,
n_aa: npt.NDArray[np.int64],
n_trna: npt.NDArray[np.int64],
limited: bool = False,
) -> npt.NDArray[np.int64]:
"""
Distributes counts of amino acids to tRNAs that are associated with
each amino acid. Uses self.process.aa_from_trna mapping to distribute
from amino acids to tRNA based on the fraction that each tRNA species
makes up for all tRNA species that code for the same amino acid.
Args:
n_aa: counts of each amino acid to distribute to each tRNA
n_trna: counts of each tRNA to determine the distribution
limited: optional, if True, limits the amino acids
distributed to each tRNA to the number of tRNA that are
available (n_trna)
Returns:
Distributed counts for each tRNA
"""
# Determine the fraction each tRNA species makes up out of all tRNA of
# the associated amino acid
with np.errstate(invalid="ignore"):
f_trna = n_trna / np.dot(
np.dot(self.process.aa_from_trna, n_trna), self.process.aa_from_trna
)
f_trna[~np.isfinite(f_trna)] = 0
trna_counts = np.zeros(f_trna.shape, np.int64)
for count, row in zip(n_aa, self.process.aa_from_trna):
idx = row == 1
frac = f_trna[idx]
counts = np.floor(frac * count)
diff = int(count - counts.sum())
# Add additional counts to get up to counts to distribute
# Prevent adding over the number of tRNA available if limited
if diff > 0:
if limited:
for _ in range(diff):
frac[(n_trna[idx] - counts) == 0] = 0
# normalize for multinomial distribution
frac /= frac.sum()
adjustment = self.process.random_state.multinomial(1, frac)
counts += adjustment
else:
adjustment = self.process.random_state.multinomial(diff, frac)
counts += adjustment
trna_counts[idx] = counts
return trna_counts
[docs]
def isTimeStepShortEnough(self, inputTimeStep, timeStepSafetyFraction):
short_enough = True
# Needs to be less than the max time step to prevent oscillatory
# behavior
if inputTimeStep > self.max_time_step:
short_enough = False
# Decrease the max time step to get more stable charging
if (not self.time_step_short_enough) and (
self.process.adjust_timestep_for_charging
):
self.max_time_step = inputTimeStep / 2
self.time_step_short_enough = True
short_enough = False
return short_enough
[docs]
def calculate_trna_charging(
synthetase_conc: Unum,
uncharged_trna_conc: Unum,
charged_trna_conc: Unum,
aa_conc: Unum,
ribosome_conc: Unum,
f: Unum,
params: dict[str, Any],
supply: Optional[Callable] = None,
time_limit: float = 1000,
limit_v_rib: bool = False,
use_disabled_aas: bool = False,
) -> tuple[Unum, float, Unum, Unum, Unum]:
"""
Calculates the steady state value of tRNA based on charging and
incorporation through polypeptide elongation. The fraction of
charged/uncharged is also used to determine how quickly the
ribosome is elongating. All concentrations are given in units of
:py:data:`~ecoli.processes.polypeptide_elongation.MICROMOLAR_UNITS`.
Args:
synthetase_conc: concentration of synthetases associated with
each amino acid
uncharged_trna_conc: concentration of uncharged tRNA associated
with each amino acid
charged_trna_conc: concentration of charged tRNA associated with
each amino acid
aa_conc: concentration of each amino acid
ribosome_conc: concentration of active ribosomes
f: fraction of each amino acid to be incorporated to total amino
acids incorporated
params: parameters used in charging equations
supply: function to get the rate of amino acid supply (synthesis
and import) based on amino acid concentrations. If None, amino
acid concentrations remain constant during charging
time_limit: time limit to reach steady state
limit_v_rib: if True, v_rib is limited to the number of amino acids
that are available
use_disabled_aas: if False, amino acids in
:py:data:`~ecoli.processes.polypeptide_elongation.REMOVED_FROM_CHARGING`
are excluded from charging
Returns:
5-element tuple containing
- **new_fraction_charged**: fraction of total tRNA that is charged for each
amino acid species
- **v_rib**: ribosomal elongation rate in units of uM/s
- **total_synthesis**: the total amount of amino acids synthesized during charging
in units of MICROMOLAR_UNITS. Will be zeros if supply function is not given.
- **total_import**: the total amount of amino acids imported during charging
in units of MICROMOLAR_UNITS. Will be zeros if supply function is not given.
- **total_export**: the total amount of amino acids exported during charging
in units of MICROMOLAR_UNITS. Will be zeros if supply function is not given.
"""
def negative_check(trna1: npt.NDArray[np.float64], trna2: npt.NDArray[np.float64]):
"""
Check for floating point precision issues that can lead to small
negative numbers instead of 0. Adjusts both species of tRNA to
bring concentration of trna1 to 0 and keep the same total concentration.
Args:
trna1: concentration of one tRNA species (charged or uncharged)
trna2: concentration of another tRNA species (charged or uncharged)
"""
mask = trna1 < 0
trna2[mask] = trna1[mask] + trna2[mask]
trna1[mask] = 0
def dcdt(t: float, c: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
"""
Function for solve_ivp to integrate
Args:
c: 1D array of concentrations of uncharged and charged tRNAs
dims: 2 * number of amino acids (uncharged tRNA come first, then charged)
t: time of integration step
Returns:
Array of dc/dt for tRNA concentrations
dims: 2 * number of amino acids (uncharged tRNA come first, then charged)
"""
v_charging, dtrna, daa = dcdt_jit(
t,
c,
n_aas_masked,
n_aas,
mask,
params["kS"],
synthetase_conc,
params["KMaa"],
params["KMtf"],
f,
params["krta"],
params["krtf"],
params["max_elong_rate"],
ribosome_conc,
limit_v_rib,
aa_rate_limit,
v_rib_max,
)
if supply is None:
v_synthesis = np.zeros(n_aas)
v_import = np.zeros(n_aas)
v_export = np.zeros(n_aas)
else:
aa_conc = c[2 * n_aas_masked : 2 * n_aas_masked + n_aas]
v_synthesis, v_import, v_export = supply(unit_conversion * aa_conc)
v_supply = v_synthesis + v_import - v_export
daa[mask] = v_supply[mask] - v_charging
return np.hstack((-dtrna, dtrna, daa, v_synthesis, v_import, v_export))
# Convert inputs for integration
synthetase_conc = synthetase_conc.asNumber(MICROMOLAR_UNITS)
uncharged_trna_conc = uncharged_trna_conc.asNumber(MICROMOLAR_UNITS)
charged_trna_conc = charged_trna_conc.asNumber(MICROMOLAR_UNITS)
aa_conc = aa_conc.asNumber(MICROMOLAR_UNITS)
ribosome_conc = ribosome_conc.asNumber(MICROMOLAR_UNITS)
unit_conversion = params["unit_conversion"]
# Remove disabled amino acids from calculations
n_total_aas = len(aa_conc)
if use_disabled_aas:
mask = np.ones(n_total_aas, bool)
else:
mask = params["charging_mask"]
synthetase_conc = synthetase_conc[mask]
original_uncharged_trna_conc = uncharged_trna_conc[mask]
original_charged_trna_conc = charged_trna_conc[mask]
original_aa_conc = aa_conc[mask]
f = f[mask]
n_aas = len(aa_conc)
n_aas_masked = len(original_aa_conc)
# Limits for integration
aa_rate_limit = original_aa_conc / time_limit
trna_rate_limit = original_charged_trna_conc / time_limit
v_rib_max = max(0, ((aa_rate_limit + trna_rate_limit) / f).min())
# Integrate rates of charging and elongation
c_init = np.hstack(
(
original_uncharged_trna_conc,
original_charged_trna_conc,
aa_conc,
np.zeros(n_aas),
np.zeros(n_aas),
np.zeros(n_aas),
)
)
sol = solve_ivp(dcdt, [0, time_limit], c_init, method="BDF")
c_sol = sol.y.T
# Determine new values from integration results
final_uncharged_trna_conc = c_sol[-1, :n_aas_masked]
final_charged_trna_conc = c_sol[-1, n_aas_masked : 2 * n_aas_masked]
total_synthesis = c_sol[-1, 2 * n_aas_masked + n_aas : 2 * n_aas_masked + 2 * n_aas]
total_import = c_sol[
-1, 2 * n_aas_masked + 2 * n_aas : 2 * n_aas_masked + 3 * n_aas
]
total_export = c_sol[
-1, 2 * n_aas_masked + 3 * n_aas : 2 * n_aas_masked + 4 * n_aas
]
negative_check(final_uncharged_trna_conc, final_charged_trna_conc)
negative_check(final_charged_trna_conc, final_uncharged_trna_conc)
fraction_charged = final_charged_trna_conc / (
final_uncharged_trna_conc + final_charged_trna_conc
)
numerator_ribosome = 1 + np.sum(
f
* (
params["krta"] / final_charged_trna_conc
+ final_uncharged_trna_conc
/ final_charged_trna_conc
* params["krta"]
/ params["krtf"]
)
)
v_rib = params["max_elong_rate"] * ribosome_conc / numerator_ribosome
if limit_v_rib:
v_rib_max = max(
0,
(
(
original_aa_conc
+ (original_charged_trna_conc - final_charged_trna_conc)
)
/ time_limit
/ f
).min(),
)
v_rib = min(v_rib, v_rib_max)
# Replace SEL fraction charged with average
new_fraction_charged = np.zeros(n_total_aas)
new_fraction_charged[mask] = fraction_charged
new_fraction_charged[~mask] = fraction_charged.mean()
return new_fraction_charged, v_rib, total_synthesis, total_import, total_export
@njit(error_model="numpy")
def dcdt_jit(
t,
c,
n_aas_masked,
n_aas,
mask,
kS,
synthetase_conc,
KMaa,
KMtf,
f,
krta,
krtf,
max_elong_rate,
ribosome_conc,
limit_v_rib,
aa_rate_limit,
v_rib_max,
):
uncharged_trna_conc = c[:n_aas_masked]
charged_trna_conc = c[n_aas_masked : 2 * n_aas_masked]
aa_conc = c[2 * n_aas_masked : 2 * n_aas_masked + n_aas]
masked_aa_conc = aa_conc[mask]
v_charging = (
kS
* synthetase_conc
* uncharged_trna_conc
* masked_aa_conc
/ (KMaa[mask] * KMtf[mask])
/ (
1
+ uncharged_trna_conc / KMtf[mask]
+ masked_aa_conc / KMaa[mask]
+ uncharged_trna_conc * masked_aa_conc / KMtf[mask] / KMaa[mask]
)
)
numerator_ribosome = 1 + np.sum(
f
* (
krta / charged_trna_conc
+ uncharged_trna_conc / charged_trna_conc * krta / krtf
)
)
v_rib = max_elong_rate * ribosome_conc / numerator_ribosome
# Handle case when f is 0 and charged_trna_conc is 0
if not np.isfinite(v_rib):
v_rib = 0
# Limit v_rib and v_charging to the amount of available amino acids
if limit_v_rib:
v_charging = np.fmin(v_charging, aa_rate_limit)
v_rib = min(v_rib, v_rib_max)
dtrna = v_charging - v_rib * f
daa = np.zeros(n_aas)
return v_charging, dtrna, daa
[docs]
def get_charging_supply_function(
supply_in_charging: bool,
mechanistic_supply: bool,
mechanistic_aa_transport: bool,
amino_acid_synthesis: Callable,
amino_acid_import: Callable,
amino_acid_export: Callable,
aa_supply_scaling: Callable,
counts_to_molar: Unum,
aa_supply: npt.NDArray[np.float64],
fwd_enzyme_counts: npt.NDArray[np.int64],
rev_enzyme_counts: npt.NDArray[np.int64],
dry_mass: Unum,
importer_counts: npt.NDArray[np.int64],
exporter_counts: npt.NDArray[np.int64],
aa_in_media: npt.NDArray[np.bool_],
) -> Optional[Callable[[npt.NDArray[np.float64]], Tuple[Unum, Unum, Unum]]]:
"""
Get a function mapping internal amino acid concentrations to the amount of
amino acid supply expected.
Args:
supply_in_charging: True if using aa_supply_in_charging option
mechanistic_supply: True if using mechanistic_translation_supply option
mechanistic_aa_transport: True if using mechanistic_aa_transport option
amino_acid_synthesis: function to provide rates of synthesis for amino
acids based on the internal state
amino_acid_import: function to provide import rates for amino
acids based on the internal and external state
amino_acid_export: function to provide export rates for amino
acids based on the internal state
aa_supply_scaling: function to scale the amino acid supply based
on the internal state
counts_to_molar: conversion factor for counts to molar
(:py:data:`~ecoli.processes.polypeptide_elongation.MICROMOLAR_UNITS`)
aa_supply: rate of amino acid supply expected
fwd_enzyme_counts: enzyme counts in forward reactions for each amino acid
rev_enzyme_counts: enzyme counts in loss reactions for each amino acid
dry_mass: dry mass of the cell with mass units
importer_counts: counts for amino acid importers
exporter_counts: counts for amino acid exporters
aa_in_media: True for each amino acid that is present in the media
Returns:
Function that provides the amount of supply (synthesis, import, export)
for each amino acid based on the internal state of the cell
"""
# Create functions that are only dependent on amino acid concentrations for more stable
# charging and amino acid concentrations. If supply_in_charging is not set, then
# setting None will maintain constant amino acid concentrations throughout charging.
supply_function = None
if supply_in_charging:
counts_to_molar = counts_to_molar.asNumber(MICROMOLAR_UNITS)
zeros = counts_to_molar * np.zeros_like(aa_supply)
if mechanistic_supply:
if mechanistic_aa_transport:
def supply_function(aa_conc):
return (
counts_to_molar
* amino_acid_synthesis(
fwd_enzyme_counts, rev_enzyme_counts, aa_conc
)[0],
counts_to_molar
* amino_acid_import(
aa_in_media,
dry_mass,
aa_conc,
importer_counts,
mechanistic_aa_transport,
),
counts_to_molar
* amino_acid_export(
exporter_counts, aa_conc, mechanistic_aa_transport
),
)
else:
def supply_function(aa_conc):
return (
counts_to_molar
* amino_acid_synthesis(
fwd_enzyme_counts, rev_enzyme_counts, aa_conc
)[0],
counts_to_molar
* amino_acid_import(
aa_in_media,
dry_mass,
aa_conc,
importer_counts,
mechanistic_aa_transport,
),
zeros,
)
else:
def supply_function(aa_conc):
return (
counts_to_molar
* aa_supply
* aa_supply_scaling(aa_conc, aa_in_media),
zeros,
zeros,
)
return supply_function
def test_polypeptide_elongation(return_data=False):
def make_elongation_rates(random, base, time_step, variable_elongation=False):
size = 1
lengths = time_step * np.full(size, base, dtype=np.int64)
lengths = stochasticRound(random, lengths) if random else np.round(lengths)
return lengths.astype(np.int64)
test_config = {
"time_step": 2,
"proteinIds": np.array(["TRYPSYN-APROTEIN[c]"]),
"ribosome30S": "CPLX0-3953[c]",
"ribosome50S": "CPLX0-3962[c]",
"make_elongation_rates": make_elongation_rates,
"proteinLengths": np.array(
[245]
), # this is the length of proteins in proteinSequences
"translation_aa_supply": {
"minimal": (units.mol / units.fg / units.min)
* np.array(
[
6.73304301e-21,
3.63835219e-21,
2.89772671e-21,
3.88086822e-21,
5.04645651e-22,
4.45295877e-21,
2.64600664e-21,
5.35711230e-21,
1.26817689e-21,
3.81168405e-21,
5.66834531e-21,
4.30576056e-21,
1.70428208e-21,
2.24878356e-21,
2.49335033e-21,
3.47019761e-21,
3.83858460e-21,
6.34564026e-22,
1.86880523e-21,
1.40959498e-27,
5.20884460e-21,
]
)
},
"proteinSequences": np.array(
[
[
12,
10,
18,
9,
13,
1,
10,
9,
9,
16,
20,
9,
18,
15,
9,
10,
20,
4,
20,
13,
7,
15,
9,
18,
4,
10,
13,
15,
14,
1,
2,
14,
11,
8,
20,
0,
16,
13,
7,
8,
12,
13,
7,
1,
10,
0,
14,
10,
13,
7,
10,
11,
20,
5,
4,
1,
11,
14,
16,
3,
0,
5,
15,
18,
7,
2,
0,
9,
18,
9,
0,
2,
8,
6,
2,
2,
18,
3,
12,
20,
16,
0,
15,
2,
9,
20,
6,
14,
14,
16,
20,
16,
20,
7,
11,
11,
15,
10,
10,
17,
9,
14,
13,
13,
7,
6,
10,
18,
17,
10,
16,
7,
2,
10,
10,
9,
3,
1,
2,
2,
1,
16,
11,
0,
8,
7,
16,
9,
0,
5,
20,
20,
2,
8,
13,
11,
11,
1,
1,
9,
15,
9,
17,
12,
13,
14,
5,
7,
16,
1,
15,
1,
7,
1,
7,
10,
10,
14,
13,
11,
16,
7,
0,
13,
8,
0,
0,
9,
0,
0,
7,
20,
14,
9,
9,
14,
20,
4,
20,
15,
16,
16,
15,
2,
11,
9,
2,
10,
2,
1,
10,
8,
2,
7,
10,
20,
9,
20,
5,
12,
10,
14,
14,
9,
3,
20,
15,
6,
18,
7,
11,
3,
6,
20,
1,
5,
10,
0,
0,
8,
4,
1,
15,
9,
12,
5,
6,
11,
9,
0,
5,
10,
3,
11,
5,
20,
0,
5,
1,
5,
0,
0,
7,
11,
20,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
-1,
]
]
).astype(np.int8),
}
polypep_elong = PolypeptideElongation(test_config)
initial_state = {
"environment": {"media_id": "minimal"},
"bulk": np.array(
[
("CPLX0-3953[c]", 100),
("CPLX0-3962[c]", 100),
("TRYPSYN-APROTEIN[c]", 0),
("RELA", 0),
("SPOT", 0),
("H2O", 0),
("PROTON", 0),
("ppGpp", 0),
]
+ [(aa, 100) for aa in DEFAULT_AA_NAMES],
dtype=[("id", "U40"), ("count", int)],
),
"unique": {
"active_ribosome": np.array(
[(1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)],
dtype=[
("_entryState", np.bool_),
("unique_index", int),
("protein_index", int),
("peptide_length", int),
("pos_on_mRNA", int),
("massDiff_DNA", "<f8"),
("massDiff_mRNA", "<f8"),
("massDiff_metabolite", "<f8"),
("massDiff_miscRNA", "<f8"),
("massDiff_nonspecific_RNA", "<f8"),
("massDiff_protein", "<f8"),
("massDiff_rRNA", "<f8"),
("massDiff_tRNA", "<f8"),
("massDiff_water", "<f8"),
],
)
},
"listeners": {"mass": {"dry_mass": 350.0}},
}
settings = {"total_time": 200, "initial_state": initial_state, "topology": TOPOLOGY}
data = simulate_process(polypep_elong, settings)
if return_data:
return data, test_config
[docs]
def run_plot(data, config):
# plot a list of variables
bulk_ids = [
"CPLX0-3953[c]",
"CPLX0-3962[c]",
"TRYPSYN-APROTEIN[c]",
"RELA",
"SPOT",
"H2O",
"PROTON",
"ppGpp",
] + [aa for aa in DEFAULT_AA_NAMES]
variables = [(bulk_id,) for bulk_id in bulk_ids]
# format data
bulk_timeseries = np.array(data["bulk"])
for i, bulk_id in enumerate(bulk_ids):
data[bulk_id] = bulk_timeseries[:, i]
plot_variables(
data,
variables=variables,
out_dir="out/processes/polypeptide_elongation",
filename="variables",
)
[docs]
def main():
data, config = test_polypeptide_elongation(return_data=True)
run_plot(data, config)
if __name__ == "__main__":
main()