diff --git a/models/ecoli/processes/rna_degradation.py b/models/ecoli/processes/rna_degradation.py index 5e51cd5348..52529152d3 100755 --- a/models/ecoli/processes/rna_degradation.py +++ b/models/ecoli/processes/rna_degradation.py @@ -1,10 +1,6 @@ -#!/usr/bin/env python - """ -RnaDegradation -RNA degradation sub-model. +Submodel for RNA degradation. -@author: Javier Carrera @organization: Covert Lab, Department of Bioengineering, Stanford University @date: Created 1/26/2015 - Updated 22/3/2017 @@ -23,25 +19,33 @@ tau = doubling time kcatEndoRNase = enzymatic activity for EndoRNases kd = RNA degradation rates - Km = Michaelis-Menten constants fitted to recapitulate first-order RNA decay: - kd * r = kcatEndoRNase * EndoRNase * r / (Km + r), non-cooperative EndoRNases - kd * r = kcatEndoRNase * EndoRNase * r/Km / (1 + sum(r/Km)), cooperation - -This sub-model encodes molecular simulation of RNA degradation as two main steps guided by RNases, "endonucleolytic cleavage" and "exonucleolytic digestion": -1. Compute total counts of RNA to be degraded (D) and total capacity for endo-cleavage (C) at each time point + Km = Michaelis-Menten constants fitted to recapitulate first-order + RNA decay: + kd * r = kcatEndoRNase * EndoRNase * r / (Km + r), + non-cooperative EndoRNases + kd * r = kcatEndoRNase * EndoRNase * r/Km / (1 + sum(r/Km)), + cooperation + +This sub-model encodes molecular simulation of RNA degradation as two main +steps guided by RNases, "endonucleolytic cleavage" and "exonucleolytic +digestion": + +1. Compute total counts of RNA to be degraded (D) and total capacity for +endo-cleavage (C) at each time point 2. Evaluate C and D. If C > D, then define a fraction of active endoRNases -3. Dissect RNA degraded into different species (mRNA, tRNA, and rRNA) by accounting endoRNases specificity -4. Update RNA fragments (assumption: fragments are represented as a pull of nucleotides) because of endonucleolytic cleavage -5. Compute total capacity of exoRNases and determine fraction of nucleotides that can be diggested -6. Update pull of metabolites (H and H2O) because of exonucleolytic digestion +3. Dissect RNA degraded into different species (mRNA, tRNA, and rRNA) by +accounting endoRNases specificity +4. Update RNA fragments (assumption: fragments are represented as a pool of +nucleotides) created because of endonucleolytic cleavage +5. Compute total capacity of exoRNases and determine fraction of nucleotides +that can be digested +6. Update pool of metabolites (H and H2O) created because of exonucleolytic +digestion """ from __future__ import division import numpy as np -import math -import numpy -import random import wholecell.processes.process from wholecell.utils.constants import REQUEST_PRIORITY_DEGRADATION @@ -72,16 +76,6 @@ def initialize(self, sim, sim_data): self.KcatExoRNase = sim_data.constants.KcatExoRNase self.KcatEndoRNases = sim_data.process.rna_decay.kcats - # Load RNase specificity - self.TargetEndoRNasesFullMRNA_indexes = sim_data.process.rna_decay.TargetEndoRNasesFullMRNA - self.TargetEndoRNasesFullTRNA_indexes = sim_data.process.rna_decay.TargetEndoRNasesFullTRNA - self.TargetEndoRNasesFullRRNA_indexes = sim_data.process.rna_decay.TargetEndoRNasesFullRRNA - - self.mrna_index = sim_data.process.rna_decay.mrna_index - self.trna_index = sim_data.process.rna_decay.trna_index - self.rrna_index = sim_data.process.rna_decay.rrna_index - self.rtrna_index = sim_data.process.rna_decay.rtrna_index - # Load information about charged tRNA self.charged_trna_names = sim_data.process.transcription.charged_trna_names self.charged_trna = self.bulkMoleculesView(self.charged_trna_names) @@ -94,11 +88,11 @@ def initialize(self, sim, sim_data): shuffleIdxs = sim_data.process.transcription.rnaDegRateShuffleIdxs self.rnaDegRates = self.rnaDegRates[shuffleIdxs] - self.isMRna = sim_data.process.transcription.rnaData["isMRna"] - self.isRRna = sim_data.process.transcription.rnaData["isRRna"] - self.isTRna = sim_data.process.transcription.rnaData["isTRna"] + self.isMRna = sim_data.process.transcription.rnaData["isMRna"].astype(np.float64) + self.isRRna = sim_data.process.transcription.rnaData["isRRna"].astype(np.float64) + self.isTRna = sim_data.process.transcription.rnaData["isTRna"].astype(np.float64) - self.rnaLens = sim_data.process.transcription.rnaData['length'].asNumber() + self.rna_lengths = sim_data.process.transcription.rnaData['length'].asNumber() # Build stoichiometric matrix endCleavageMetaboliteIds = [id_ + "[c]" for id_ in sim_data.moleculeGroups.fragmentNT_IDs] @@ -117,7 +111,7 @@ def initialize(self, sim, sim_data): self.rnas = self.bulkMoleculesView(rnaIds) self.h2o = self.bulkMoleculesView(["WATER[c]"]) self.nmps = self.bulkMoleculesView(["AMP[c]", "CMP[c]", "GMP[c]", "UMP[c]"]) - self.h = self.bulkMoleculesView(["PROTON[c]"]) + self.proton = self.bulkMoleculesView(["PROTON[c]"]) self.fragmentMetabolites = self.bulkMoleculesView(endCleavageMetaboliteIds) self.fragmentBases = self.bulkMoleculesView([id_ + "[c]" for id_ in sim_data.moleculeGroups.fragmentNT_IDs]) @@ -129,7 +123,6 @@ def initialize(self, sim, sim_data): self.rrlaIdx = sim_data.process.transcription.rnaData["id"].tolist().index("RRLA-RRNA[c]") self.rrsaIdx = sim_data.process.transcription.rnaData["id"].tolist().index("RRSA-RRNA[c]") - self.endoRnases = self.bulkMoleculesView(endoRnaseIds) self.exoRnases = self.bulkMoleculesView(exoRnaseIds) self.bulkMoleculesRequestPriorityIs(REQUEST_PRIORITY_DEGRADATION) @@ -137,162 +130,173 @@ def initialize(self, sim, sim_data): # Load Michaelis-Menten constants fitted to recapitulate first-order RNA decay model self.Km = sim_data.process.transcription.rnaData["KmEndoRNase"] + # If set to True, assume cooperation in endoRNase activity self.EndoRNaseCoop = sim_data.constants.EndoRNaseCooperation + + # If set to False, assume RNAs degrade simply by first-order kinetics self.EndoRNaseFunc = sim_data.constants.EndoRNaseFunction def calculateRequest(self): - # Compute factor that convert counts into concentration, and viceversa - cellMass = (self.readFromListener("Mass", "cellMass") * units.fg) - cellVolume = cellMass / self.cellDensity - countsToMolar = 1 / (self.nAvogadro * cellVolume) + # Compute factor that convert counts into concentration, and vice versa + cell_mass = self.readFromListener("Mass", "cellMass") * units.fg + cell_volume = cell_mass / self.cellDensity + counts_to_molar = 1 / (self.nAvogadro * cell_volume) - # View RNAs - rnasTotal = self.rnas.total().copy() - rnasTotal[self.rrsaIdx] += self.ribosome30S.total() - rnasTotal[[self.rrlaIdx, self.rrfaIdx]] += self.ribosome50S.total() - rnasTotal[[self.rrlaIdx, self.rrfaIdx, self.rrsaIdx]] += self.activeRibosomes.total() - rnasTotal[self.isTRna] += self.charged_trna.total() + # Get total counts of RNAs including rRNAs and charged tRNAs + rna_counts = self.rnas.total().copy() + rna_counts[self.rrsaIdx] += self.ribosome30S.total() + rna_counts[[self.rrlaIdx, self.rrfaIdx]] += self.ribosome50S.total() + rna_counts[[self.rrlaIdx, self.rrfaIdx, self.rrsaIdx]] += self.activeRibosomes.total() + rna_counts[self.isTRna.astype(np.bool)] += self.charged_trna.total() - # Calculate endoRNase active fraction based on Michaelis-Menten kinetics - if not self.EndoRNaseCoop: - fracEndoRnaseSaturated = countsToMolar * rnasTotal / (self.Km + (countsToMolar * rnasTotal)) + # Compute RNA concentrations + rna_conc_molar = counts_to_molar * rna_counts - # Calculate endoRNase active fraction based on generalized Michaelis-Menten kinetics + # Get counts of endoRNases + endornase_counts = self.endoRnases.total().copy() + total_kcat_endornase = units.dot(self.KcatEndoRNases, endornase_counts) + + # Calculate the fraction of active endoRNases for each RNA based on + # Michaelis-Menten kinetics if self.EndoRNaseCoop: - fracEndoRnaseSaturated = (countsToMolar * rnasTotal) / self.Km / (1 + units.sum((countsToMolar * rnasTotal) / self.Km)) - - Kd = self.rnaDegRates - Kcat = self.KcatEndoRNases - EndoR = sum(self.endoRnases.total()) - RNA = rnasTotal - FractDiffRNAdecay = units.sum(units.abs(Kd * RNA - units.sum(self.KcatEndoRNases * self.endoRnases.total()) * fracEndoRnaseSaturated)) - FractEndoRRnaCounts = EndoR.astype(float) / sum(RNA.astype(float)) - - # Calculate total counts of RNAs to degrade according to - # the total counts of active endoRNases and their cleavage activity - nRNAsTotalToDegrade = np.round((units.sum(self.KcatEndoRNases * self.endoRnases.total()) * fracEndoRnaseSaturated * (units.s * self.timeStepSec())).asNumber().sum()) - - # Dissect RNAse specificity into mRNA, tRNA, and rRNA - MrnaSpec = units.sum(fracEndoRnaseSaturated * self.isMRna) - TrnaSpec = units.sum(fracEndoRnaseSaturated * self.isTRna) - RrnaSpec = units.sum(fracEndoRnaseSaturated * self.isRRna) - - TargetEndoRNasesFullMRNA = np.zeros_like(self.TargetEndoRNasesFullMRNA_indexes) - TargetEndoRNasesFullMRNA[self.TargetEndoRNasesFullMRNA_indexes == self.mrna_index] = MrnaSpec - TargetEndoRNasesFullMRNA[self.TargetEndoRNasesFullMRNA_indexes == 0] = 0 - TargetEndoRNasesFullMRNA[self.TargetEndoRNasesFullMRNA_indexes == 1] = 1 - - TargetEndoRNasesFullTRNA = np.zeros_like(self.TargetEndoRNasesFullTRNA_indexes) - TargetEndoRNasesFullTRNA[self.TargetEndoRNasesFullTRNA_indexes == self.trna_index] = TrnaSpec - TargetEndoRNasesFullTRNA[self.TargetEndoRNasesFullTRNA_indexes == 0] = 0 - TargetEndoRNasesFullTRNA[self.TargetEndoRNasesFullTRNA_indexes == 1] = 1 - - TargetEndoRNasesFullRRNA = np.zeros_like(self.TargetEndoRNasesFullRRNA_indexes) - TargetEndoRNasesFullRRNA[self.TargetEndoRNasesFullRRNA_indexes == self.rrna_index] = RrnaSpec - TargetEndoRNasesFullRRNA[self.TargetEndoRNasesFullRRNA_indexes == self.rtrna_index] = RrnaSpec + TrnaSpec - TargetEndoRNasesFullRRNA[self.TargetEndoRNasesFullRRNA_indexes == 0] = 0 - TargetEndoRNasesFullRRNA[self.TargetEndoRNasesFullRRNA_indexes == 1] = 1 - - - # Dissect total counts of RNA degraded into mRNA, tRNA, and rRNA - TargetEndoRNasesFullMRNA = MrnaSpec - TargetEndoRNasesFullTRNA = TrnaSpec - TargetEndoRNasesFullRRNA = RrnaSpec - - nMRNAsTotalToDegrade = np.round(sum(TargetEndoRNasesFullMRNA * - self.endoRnases.total() * - self.KcatEndoRNases * (units.s * self.timeStepSec())).asNumber() + frac_endornase_saturated = ( + rna_conc_molar / self.Km / (1 + units.sum(rna_conc_molar / self.Km)) + ).asNumber() + else: + frac_endornase_saturated = ( + rna_conc_molar / (self.Km + rna_conc_molar) + ).asNumber() + + # Calculate difference in degradation rates from first-order decay + # and the number of EndoRNases per one molecule of RNA + total_endornase_counts = np.sum(endornase_counts) + diff_relative_first_order_decay = units.sum( + units.abs(self.rnaDegRates * rna_counts - + total_kcat_endornase * frac_endornase_saturated) ) - nTRNAsTotalToDegrade = np.round(sum(TargetEndoRNasesFullTRNA * - self.endoRnases.total() * - self.KcatEndoRNases * (units.s * self.timeStepSec())).asNumber() + endornase_per_rna = total_endornase_counts / np.sum(rna_counts) + + self.writeToListener("RnaDegradationListener", + "FractionActiveEndoRNases", + np.sum(frac_endornase_saturated) ) - nRRNAsTotalToDegrade = np.round(sum(TargetEndoRNasesFullRRNA * - self.endoRnases.total() * - self.KcatEndoRNases * (units.s * self.timeStepSec())).asNumber() + self.writeToListener("RnaDegradationListener", + "DiffRelativeFirstOrderDecay", + diff_relative_first_order_decay.asNumber() ) - if nRNAsTotalToDegrade != nMRNAsTotalToDegrade + nTRNAsTotalToDegrade + nRRNAsTotalToDegrade: - nRNAsTotalToDegrade = nMRNAsTotalToDegrade + nTRNAsTotalToDegrade + nRRNAsTotalToDegrade - - # Compute RNAse specificity - RNAspecificity = (fracEndoRnaseSaturated / units.sum(fracEndoRnaseSaturated)).asNumber() - - nRNAsToDegrade = np.zeros(len(RNAspecificity)) - nMRNAsToDegrade = np.zeros(len(RNAspecificity)) - nTRNAsToDegrade = np.zeros(len(RNAspecificity)) - nRRNAsToDegrade = np.zeros(len(RNAspecificity)) - - # Boolean variable (nRNAs) that tracks availability of RNAs for a given gene - nRNAs = rnasTotal.astype(np.bool) - - - # Determine mRNAs to be degraded according to RNA specificities and total counts of mRNAs degraded - while nMRNAsToDegrade.sum() < nMRNAsTotalToDegrade and (rnasTotal * self.isMRna).sum() != 0: - nMRNAsToDegrade += np.fmin( - self.randomState.multinomial(nMRNAsTotalToDegrade - nMRNAsToDegrade.sum(), 1. / sum(RNAspecificity * self.isMRna * nRNAs) * RNAspecificity * self.isMRna * nRNAs), - rnasTotal * self.isMRna + self.writeToListener( + "RnaDegradationListener", + "FractEndoRRnaCounts", + endornase_per_rna) + + if self.EndoRNaseFunc: + # Dissect RNAse specificity into mRNA, tRNA, and rRNA + mrna_specificity = np.dot(frac_endornase_saturated, self.isMRna) + trna_specificity = np.dot(frac_endornase_saturated, self.isTRna) + rrna_specificity = np.dot(frac_endornase_saturated, self.isRRna) + + n_total_mrnas_to_degrade = self._calculate_total_n_to_degrade( + mrna_specificity, + total_kcat_endornase ) - - # Determine tRNAs and rRNAs to be degraded (with equal specificity) depending on total counts degraded, respectively - while nTRNAsToDegrade.sum() < nTRNAsTotalToDegrade and (rnasTotal * self.isTRna).sum() != 0: - nTRNAsToDegrade += np.fmin( - self.randomState.multinomial(nTRNAsTotalToDegrade - nTRNAsToDegrade.sum(), 1. / sum(self.isTRna * nRNAs) * self.isTRna * nRNAs), - rnasTotal * self.isTRna + n_total_trnas_to_degrade = self._calculate_total_n_to_degrade( + trna_specificity, + total_kcat_endornase ) - while nRRNAsToDegrade.sum() < nRRNAsTotalToDegrade and (rnasTotal * self.isRRna).sum() != 0: - nRRNAsToDegrade += np.fmin( - self.randomState.multinomial(nRRNAsTotalToDegrade - nRRNAsToDegrade.sum(), 1. / sum(self.isRRna * nRNAs) * self.isRRna * nRNAs), - rnasTotal * self.isRRna + n_total_rrnas_to_degrade = self._calculate_total_n_to_degrade( + rrna_specificity, + total_kcat_endornase + ) + + # Compute RNAse specificity + rna_specificity = frac_endornase_saturated / np.sum(frac_endornase_saturated) + + # Boolean variable that tracks existence of each RNA + rna_exists = rna_counts.astype(np.bool) + + # Compute degradation probabilities of each RNA: for mRNAs, this + # is based on the specificity of each mRNA. For tRNAs and rRNAs, + # this is distributed evenly. + mrna_deg_probs = 1. / np.dot(rna_specificity, self.isMRna * rna_exists) * rna_specificity * self.isMRna * rna_exists + trna_deg_probs = 1. / np.dot(self.isTRna, rna_exists) * self.isTRna * rna_exists + rrna_deg_probs = 1. / np.dot(self.isRRna, rna_exists) * self.isRRna * rna_exists + + # Mask RNA counts into each class of RNAs + mrna_counts = rna_counts * self.isMRna + trna_counts = rna_counts * self.isTRna + rrna_counts = rna_counts * self.isRRna + + # Determine number of individual RNAs to be degraded for each class + # of RNA. + n_mrnas_to_degrade = self._get_rnas_to_degrade( + n_total_mrnas_to_degrade, mrna_deg_probs, mrna_counts + ) + + n_trnas_to_degrade = self._get_rnas_to_degrade( + n_total_trnas_to_degrade, trna_deg_probs, trna_counts ) - nRNAsToDegrade = nMRNAsToDegrade + nTRNAsToDegrade + nRRNAsToDegrade + n_rrnas_to_degrade = self._get_rnas_to_degrade( + n_total_rrnas_to_degrade, rrna_deg_probs, rrna_counts + ) + + n_rnas_to_degrade = n_mrnas_to_degrade + n_trnas_to_degrade + n_rrnas_to_degrade # First order decay with non-functional EndoRNase activity - # Determine mRNAs to be degraded by sampling a Poisson distribution (Kdeg * RNA) - if not self.EndoRNaseFunc: - nRNAsToDegrade = np.fmin( - self.randomState.poisson( (self.rnaDegRates * rnasTotal).asNumber() ), + # Determine mRNAs to be degraded by sampling a Poisson distribution + # (Kdeg * RNA) + else: + n_rnas_to_degrade = np.fmin( + self.randomState.poisson( + (self.rnaDegRates * rna_counts).asNumber() + ), self.rnas.total() ) - self.rnas.requestIs(nRNAsToDegrade) + self.rnas.requestIs(n_rnas_to_degrade) self.endoRnases.requestAll() self.exoRnases.requestAll() self.fragmentBases.requestAll() - # Calculating amount of water required for total RNA hydrolysis by endo and - # exonucleases. Assuming complete hydrolysis for now. One additional water - # for each RNA to hydrolyze the 5' diphosphate. - waterForNewRnas = (nRNAsToDegrade * (self.rnaLens - 1)).sum() + nRNAsToDegrade.sum() + # Calculate the amount of water required for total RNA hydrolysis by + # endo and exonucleases. Assuming complete hydrolysis for now. Note + # that one additional water molecule is needed for each RNA to + # hydrolyze the 5' diphosphate. + waterForNewRnas = np.dot(n_rnas_to_degrade, self.rna_lengths) waterForLeftOverFragments = self.fragmentBases.total().sum() self.h2o.requestIs(waterForNewRnas + waterForLeftOverFragments) - - self.writeToListener("RnaDegradationListener", "FractionActiveEndoRNases", sum(fracEndoRnaseSaturated)) - self.writeToListener("RnaDegradationListener", "DiffRelativeFirstOrderDecay", FractDiffRNAdecay.asNumber()) - self.writeToListener("RnaDegradationListener", "FractEndoRRnaCounts", FractEndoRRnaCounts) + def evolveState(self): + n_degraded_rna = self.rnas.counts() - self.writeToListener("RnaDegradationListener", "countRnaDegraded", self.rnas.counts()) - self.writeToListener("RnaDegradationListener", "nucleotidesFromDegradation", (self.rnas.counts() * self.rnaLens).sum()) + self.writeToListener( + "RnaDegradationListener", "countRnaDegraded", n_degraded_rna + ) + self.writeToListener( + "RnaDegradationListener", "nucleotidesFromDegradation", + np.dot(n_degraded_rna, self.rna_lengths) + ) # Calculate endolytic cleavage events - # Modeling assumption: Once a RNA is cleaved by an endonuclease it's resulting nucleotides - # are lumped together as "polymerized fragments". These fragments can carry over from - # previous timesteps. We are also assuming that during endonucleolytic cleavage the 5' - # terminal phosphate is removed. This is modeled as all of the fragments being one - # long linear chain of "fragment bases". + # Modeling assumption: Once a RNA is cleaved by an endonuclease its + # resulting nucleotides are lumped together as "polymerized fragments". + # These fragments can carry over from previous timesteps. We are also + # assuming that during endonucleolytic cleavage the 5'terminal + # phosphate is removed. This is modeled as all of the fragments being + # one long linear chain of "fragment bases". # Example: - # PPi-Base-PO4(-)-Base-PO4(-)-Base-PO4(-)-Base-OH - # ==> - # Pi-FragmentBase-PO4(-)-FragmentBase-PO4(-)-FragmentBase-PO4(-)-FragmentBase + PPi + # PPi-Base-PO4(-)-Base-PO4(-)-Base-OH + # ==> + # Pi-FragmentBase-PO4(-)-FragmentBase-PO4(-)-FragmentBase + PPi # Note: Lack of -OH on 3' end of chain - metabolitesEndoCleavage = np.dot(self.endoDegradationSMatrix, self.rnas.counts()) + metabolitesEndoCleavage = np.dot( + self.endoDegradationSMatrix, n_degraded_rna) self.rnas.countsIs(0) self.fragmentMetabolites.countsInc(metabolitesEndoCleavage) @@ -302,34 +306,99 @@ def evolveState(self): # Calculate exolytic cleavage events - # Modeling assumption: We model fragments as one long fragment chain of polymerized nucleotides. - # We are also assuming that there is no sequence specificity or bias towards which nucleotides - # are hydrolyzed. + # Modeling assumption: We model fragments as one long fragment chain of + # polymerized nucleotides. We are also assuming that there is no + # sequence specificity or bias towards which nucleotides are + # hydrolyzed. # Example: - # Pi-FragmentBase-PO4(-)-FragmentBase-PO4(-)-FragmentBase-PO4(-)-FragmentBase + 4 H2O - # ==> - # 3 NMP + 3 H(+) + # Pi-FragmentBase-PO4(-)-FragmentBase-PO4(-)-FragmentBase + 3 H2O + # ==> + # 3 NMP + 3 H(+) # Note: Lack of -OH on 3' end of chain - nExoRNases = self.exoRnases.counts() - exoCapacity = nExoRNases.sum() * self.KcatExoRNase * (units.s * self.timeStepSec()) - NucleotideRecycling = self.fragmentBases.counts().sum() + n_exoRNases = self.exoRnases.counts() + n_fragment_bases = self.fragmentBases.counts() + n_fragment_bases_sum = n_fragment_bases.sum() + + exornase_capacity = n_exoRNases.sum() * self.KcatExoRNase * ( + units.s * self.timeStepSec() + ) - if exoCapacity >= self.fragmentBases.counts().sum(): - self.nmps.countsInc(self.fragmentBases.counts()) - self.h2o.countsDec(self.fragmentBases.counts().sum()) - self.h.countsInc(self.fragmentBases.counts().sum()) + if exornase_capacity >= n_fragment_bases_sum: + self.nmps.countsInc(n_fragment_bases) + self.h2o.countsDec(n_fragment_bases_sum) + self.proton.countsInc(n_fragment_bases_sum) self.fragmentBases.countsIs(0) + total_fragment_bases_digested = n_fragment_bases_sum + else: - fragmentSpecificity = self.fragmentBases.counts() / self.fragmentBases.counts().sum() - possibleBasesToDigest = self.randomState.multinomial(exoCapacity, fragmentSpecificity) - fragmentBasesDigested = self.fragmentBases.counts() - np.fmax(self.fragmentBases.counts() - possibleBasesToDigest, 0) - self.nmps.countsInc(fragmentBasesDigested) - self.h2o.countsDec(fragmentBasesDigested.sum()) - self.h.countsInc(fragmentBasesDigested.sum()) - self.fragmentBases.countsDec(fragmentBasesDigested) - NucleotideRecycling = fragmentBasesDigested.sum() - - self.writeToListener("RnaDegradationListener", "fragmentBasesDigested", NucleotideRecycling) + fragmentSpecificity = n_fragment_bases / n_fragment_bases_sum + possibleBasesToDigest = self.randomState.multinomial( + exornase_capacity, fragmentSpecificity) + n_fragment_bases_digested = n_fragment_bases - np.fmax( + n_fragment_bases - possibleBasesToDigest, 0) + + total_fragment_bases_digested = n_fragment_bases_digested.sum() + + self.nmps.countsInc(n_fragment_bases_digested) + self.h2o.countsDec(total_fragment_bases_digested) + self.proton.countsInc(total_fragment_bases_digested) + self.fragmentBases.countsDec(n_fragment_bases_digested) + + self.writeToListener("RnaDegradationListener", + "fragmentBasesDigested", total_fragment_bases_digested) + + + def _calculate_total_n_to_degrade(self, specificity, total_kcat_endornase): + """ + Calculate the total number of RNAs to degrade for a specific class of + RNAs, based on the specificity of endoRNases on that specific class and + the total kcat value of the endoRNases. + + Args: + specificity: Sum of fraction of active endoRNases for all RNAs + in a given class + total_kcat_endornase: The summed kcat of all existing endoRNases + Returns: + Total number of RNAs to degrade for the given class of RNAs + """ + return np.round( + (specificity * total_kcat_endornase + * (units.s * self.timeStepSec())).asNumber() + ) + + + def _get_rnas_to_degrade(self, n_total_rnas_to_degrade, rna_deg_probs, + rna_counts): + """ + Distributes the total count of RNAs to degrade for each class of RNAs + into individual RNAs, based on the given degradation probabilities + of individual RNAs. The upper bound is set by the current count of the + specific RNA. + + Args: + n_total_rnas_to_degrade: Total number of RNAs to degrade for the + given class of RNAs (integer, scalar) + rna_deg_probs: Degradation probabilities of each RNA (vector of + equal length to the total number of different RNAs) + rna_counts: Current counts of each RNA molecule (vector of equal + length to the total number of different RNAs) + Returns: + Vector of equal length to rna_counts, specifying the number of + molecules to degrade for each RNA + """ + n_rnas_to_degrade = np.zeros_like(rna_counts) + + if rna_counts.sum() != 0: + while n_rnas_to_degrade.sum() < n_total_rnas_to_degrade: + n_rnas_to_degrade += np.fmin( + self.randomState.multinomial( + n_total_rnas_to_degrade - n_rnas_to_degrade.sum(), + rna_deg_probs + ), + rna_counts + ) + + return n_rnas_to_degrade diff --git a/reconstruction/ecoli/dataclasses/process/rna_decay.py b/reconstruction/ecoli/dataclasses/process/rna_decay.py index 2dc70ed220..ef52b8738d 100644 --- a/reconstruction/ecoli/dataclasses/process/rna_decay.py +++ b/reconstruction/ecoli/dataclasses/process/rna_decay.py @@ -21,11 +21,6 @@ def __init__(self, raw_data, sim_data): self._buildRnaDecayData(raw_data, sim_data) def _buildRnaDecayData(self, raw_data, sim_data): - self.mrna_index = 2 - self.rrna_index = 3 - self.trna_index = 4 - self.rtrna_index = 5 - self.endoRnaseIds = [x["endoRnase"].encode("utf-8") for x in raw_data.endoRnases] self.kcats = (1 / units.s) * np.array([x["kcat"].asNumber(1 / units.s) for x in raw_data.endoRnases]) self.StatsFit = { @@ -56,40 +51,6 @@ def _buildRnaDecayData(self, raw_data, sim_data): self.KmConvergence = [] - self.TargetEndoRNasesFullMRNA = np.zeros(len(self.endoRnaseIds)) - self.TargetEndoRNasesFullTRNA = np.zeros(len(self.endoRnaseIds)) - self.TargetEndoRNasesFullRRNA = np.zeros(len(self.endoRnaseIds)) - - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("EG10856-MONOMER[p]")] = self.mrna_index - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("EG10857-MONOMER[c]")] = self.mrna_index - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("G7175-MONOMER[c]")] = 0 - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("EG10859-MONOMER[i]")] = self.mrna_index - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("EG11299-MONOMER[c]")] = self.mrna_index - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("EG10860-MONOMER[c]")] = self.mrna_index - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("EG10861-MONOMER[c]")] = self.mrna_index - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("G7365-MONOMER[c]")] = self.mrna_index - self.TargetEndoRNasesFullMRNA[self.endoRnaseIds.index("EG10862-MONOMER[c]")] = self.mrna_index - - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("EG10856-MONOMER[p]")] = self.trna_index - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("EG10857-MONOMER[c]")] = 0 - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("G7175-MONOMER[c]")] = 1 - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("EG10859-MONOMER[i]")] = self.trna_index - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("EG11299-MONOMER[c]")] = 0 - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("EG10860-MONOMER[c]")] = self.trna_index - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("EG10861-MONOMER[c]")] = self.trna_index - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("G7365-MONOMER[c]")] = self.trna_index - self.TargetEndoRNasesFullTRNA[self.endoRnaseIds.index("EG10862-MONOMER[c]")] = self.trna_index - - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("EG10856-MONOMER[p]")] = self.rrna_index - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("EG10857-MONOMER[c]")] = self.rtrna_index - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("G7175-MONOMER[c]")] = 0 - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("EG10859-MONOMER[i]")] = self.rrna_index - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("EG11299-MONOMER[c]")] = self.rtrna_index - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("EG10860-MONOMER[c]")] = self.rrna_index - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("EG10861-MONOMER[c]")] = self.rrna_index - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("G7365-MONOMER[c]")] = self.rrna_index - self.TargetEndoRNasesFullRRNA[self.endoRnaseIds.index("EG10862-MONOMER[c]")] = self.rrna_index - def kmLossFunction(self, vMax, rnaConc, kDeg, isEndoRnase, alpha): ''' Generates the functions used for estimating the per-RNA affinities (Michaelis-Menten