3. Building a ME-model¶
In [1]:
from __future__ import print_function
import cobrame
from cobrame.util import dogma, building
import cobrame.util.building
import cobra
import cobra.test
from collections import defaultdict
#import warnings
#warnings.filterwarnings('ignore')
/home/sbrg-cjlloyd/cobrapy/cobra/io/sbml3.py:24: UserWarning: Install lxml for faster SBML I/O
warn("Install lxml for faster SBML I/O")
/home/sbrg-cjlloyd/cobrapy/cobra/io/__init__.py:12: UserWarning: cobra.io.sbml requires libsbml
warn("cobra.io.sbml requires libsbml")
3.1. Overview¶
COBRAme is constructed entirely over COBRApy. This means that ME-model reactions will have all of the same properties, methods, and functions as a COBRApy reaction. However, one key difference between M and ME models is that many reactions involved in gene expression are effecively templates that are constructed identically but vary based on characteristics of the gene being expressed. For example, a gene with a given nucleotide sequence is always translated following the same rules provided by the codon table for that organism.
In order to facilliate the template nature of many gene expression
reactions, COBRAme reactions are constructed and their components are
manipulated through the use of ProcessData
classes. These act as
information vessels for holding the information assocatied with a
cellular process in simple, standard datatypes such as dictionaries and
strings.
This tutorial will go step-by-step through the process of creating a generic enzyme catalyzed reaction (i.e. metabolic reaction):
which requires the formation and coupling of complex_ab in order to proceed.
In order for this reaction to carry flux in the model we will additionally need to first add the corresponding:
- Transcription reactions
- Translation reactions
- tRNA charging reactions
- Complex formation reactions
Once these are added we will add in the synthesis of key macromolecular components (ribosome, RNA polymerase, etc.) and show how they are coupled to their respective reactions. The derived coupling coefficients will also be described. For more on the derivation of the coupling coefficients, reference the supplemental text of O’brien et. al. 2013
3.2. Initializing new ME-Models¶
When applying some constraints in the ME-model, metabolite properties are required. For instance, to calculate the total biomass (by molecular weight) produced by a particular macromolecule, the amino acid, nucleotide, etc. molecular weights are required. To enable these calculations, all metabolites from iJO1366, along with their metabolite attributes are added to the newly initialized ME-model.
Further the reactions from iJO1366 will be added to the ME-model to demonstrate ME-model solving procedures.
In [2]:
# create empty ME-model
me = cobrame.MEModel('test')
ijo = cobra.test.create_test_model('ecoli')
In [3]:
# Add all metabolites and reactions from iJO1366 to the new ME-model
for met in ijo.metabolites:
me.add_metabolites(met)
for rxn in ijo.reactions:
me.add_reaction(rxn)
The ME-model contains a “global_info” attribute which stores information used to calculate coupling constraints, along with other functions. The specifics of each of these constraints will be discussed when they are implemented.
In [4]:
# "Translational capacity" of organism
me.global_info['kt'] = 4.5 # (in h-1)scott 2010, RNA-to-protein curve fit
me.global_info['r0'] = 0.087 # scott 2010, RNA-to-protein curve fit
me.global_info['k_deg'] = 1.0/5. * 60.0 # 1/5 1/min 60 min/h # h-1
# Molecular mass of RNA component of ribosome
me.global_info['m_rr'] = 1453. # in kDa
# Average molecular mass of an amino acid
me.global_info['m_aa'] = 109. / 1000. # in kDa
# Proportion of RNA that is rRNA
me.global_info['f_rRNA'] = .86
me.global_info['m_nt'] = 324. / 1000. # in kDa
me.global_info['f_mRNA'] = .02
# tRNA associated global information
me.global_info['m_tRNA'] = 25000. / 1000. # in kDA
me.global_info['f_tRNA'] = .12
# Define the types of biomass that will be synthesized in the model
me.add_biomass_constraints_to_model(["protein_biomass", "mRNA_biomass", "tRNA_biomass", "rRNA_biomass",
"ncRNA_biomass", "DNA_biomass", "lipid_biomass", "constituent_biomass",
"prosthetic_group_biomass", "peptidoglycan_biomass"])
Define sequence of gene that will be expressed in tutorial
In [5]:
sequence = ("ATG" + "TTT" * 12 + "TAT" * 12 +
"ACG" * 12 + "GAT" * 12 + "AGT" * 12 + "TGA")
3.3. Adding Reactions without utility functions¶
We’ll first demonstrate how transcription, translation, tRNA charging,
complex formation, and metabolic reactions can be added to a model
without using any of the utility functions provided in
cobrame.util.building.py
. The second half of the tutorial will show
how these utility functions can be used to add these reactions.
The basic workflow for adding any reaction to a ME-model using COBRAme occurs in three steps:
- Create the ProcessData(s) associated with the reaction and populate them with the necessary information
- Create the MEReaction and link the appropriate ProcessData
- Execute the MEReaction’s update method
3.3.1. Add Transcription Reaction¶
3.3.1.1. Add TranscribedGene metabolite to model¶
Transcription reactions is unique in that they occur at a transcription unit level and can code for multiple transcript products. Therefore the nucleotide sequence of both the transcription unit and the RNA transcripts must be defined in order to correctly construct a transcription reaction.
-
class
cobrame.core.component.
TranscribedGene
(id, rna_type, nucleotide_sequence)[source] Metabolite class for gene created from
cobrame.core.reaction.TranscriptionReaction
Parameters: -
left_pos
¶ int – Left position of gene on the sequence of the (+) strain
-
right_pos
¶ int – Right position of gene on the sequence of the (+) strain
-
strand
¶ str –
- (+) if the RNA product is on the leading strand
- (-) if the RNA product is on the comple(mentary strand
-
In [6]:
gene = cobrame.TranscribedGene('RNA_a', 'mRNA', sequence)
me.add_metabolites([gene])
When adding the TranscribedGene
above, the RNA_type
and
nucleotide_sequence
was assigned to the gene. This sequence cannot
be determined from the transcription unit (TU) sequence because a single
TU often contains several different RNAs.
3.3.1.2. Add TranscriptionData to model¶
-
class
cobrame.core.processdata.
TranscriptionData
(id, model, rna_products=set([]))[source] Class for storing information needed to define a transcription reaction
Parameters: - id (str) – Identifier of the transcription unit, typically beginning with ‘TU’
- model (
cobrame.core.model.MEModel
) – ME-model that the TranscriptionData is associated with
-
nucleotide_sequence
¶ str – String of base pair abbreviations for nucleotides contained in the transcription unit
-
RNA_products
¶ set – IDs of
cobrame.core.component.TranscribedGene
that the transcription unit encodes. Each member should be prefixed with “RNA + _”
-
RNA_polymerase
¶ str – ID of the
cobrame.core.component.RNAP
that transcribes the transcription unit. Different IDs are used for different sigma factors
-
subreactions
¶ collections.DefaultDict(int)
– Dictionary of {cobrame.core.processdata.SubreactionData
ID: num_usages} required for the transcription unit to be transcribed
In [7]:
transcription_data = cobrame.TranscriptionData('TU_a',me,rna_products={'RNA_a'})
transcription_data.nucleotide_sequence = sequence
3.3.1.3. Add TranscriptionReaction to model¶
And point TranscriptionReaction
to TranscriptionData
-
class
cobrame.core.reaction.
TranscriptionReaction
(id)[source] Transcription of a TU to produced TranscribedGene.
RNA is transcribed on a transcription unit (TU) level. This type of reaction produces all of the RNAs contained within a TU, as well as accounts for the splicing/excision of RNA between tRNAs and rRNAs. The appropriate RNA_biomass constrain is produced based on the molecular weight of the RNAs being transcribed
Parameters: id (str) – Identifier of the transcription reaction. As a best practice, this ID should be prefixed with ‘transcription + _’
In [8]:
transcription_rxn = cobrame.TranscriptionReaction('transcription_TU_a')
transcription_rxn.transcription_data = transcription_data
me.add_reactions([transcription_rxn])
3.3.1.4. Update TranscriptionReaction¶
-
TranscriptionReaction.
update
(verbose=True)[source] Creates reaction using the associated transcription data and adds chemical formula to RNA products
This function adds the following components to the reaction stoichiometry (using ‘data’ as shorthand for
cobrame.core.processdata.TranscriptionData
):- RNA_polymerase from data.RNA_polymerase w/ coupling coefficient (if present)
- RNA products defined in data.RNA_products
- Nucleotide reactants defined in data.nucleotide_counts
- If tRNA or rRNA contained in data.RNA_types, excised base products
- Metabolites + enzymes w/ coupling coefficients defined in data.subreactions (if present)
- Biomass
cobrame.core.component.Constraint
corresponding to data.RNA_products and their associated masses - Demand reactions for each transcript product of this reaction
Parameters: verbose (bool) – Prints when new metabolites are added to the model when executing update()
In [9]:
transcription_rxn.update()
print(transcription_rxn.reaction)
86 atp_c + 38 ctp_c + 12 gtp_c + 50 utp_c --> RNA_a + 59.172286 mRNA_biomass + 186 ppi_c
/home/sbrg-cjlloyd/cobrame/cobrame/core/reaction.py:813 UserWarning: RNA Polymerase () not found
This reaction now produces a small amount of the a mRNA_biomass
metabolite (constraint). This term has a coefficient corresponding to
the molecular weight (in \(kDA\)) of the RNA being transcribed. This
constraint will be implemented into a \(v_{biomasss\_dilution}\)
reaction with the form:
A mathematical description of the biomass
constraint can be found in
Biomass Dilution Constraints in ME-Model Fundamentals.
3.3.1.5. Incorporate RNA Polymerase¶
For the purposes of this tutorial, we’ll skip the steps required to add the reactions to form the RNA_polymerase. The steps however are identical to those outlined in add enzyme complexes below
-
class
cobrame.core.component.
RNAP
(id)[source] Metabolite class for RNA polymerase complexes. Inherits from
cobrame.core.component.Complex
Parameters: id (str) – Identifier of the RNA Polymerase.
In [10]:
RNAP = cobrame.RNAP('RNA_polymerase')
me.add_metabolites(RNAP)
Associate RNA_polymerase with all TranscriptionData
and update
In [11]:
for data in me.transcription_data:
data.RNA_polymerase = RNAP.id
me.reactions.transcription_TU_a.update()
print(me.reactions.transcription_TU_a.reaction)
0.00088887053605567*mu + 0.000347992814865795 RNA_polymerase + 86 atp_c + 38 ctp_c + 12 gtp_c + 50 utp_c --> RNA_a + 59.172286 mRNA_biomass + 186 ppi_c
The coefficient for RNA_polymerase is the first instance in this tutorial where a coupling constraint is imposed. In this case the constraint couples the formation of a RNA_polymerase metabolite to its transcription flux. This constraint is formulated as in O’brien et. al. 2013, with assumption that \(k_{rnap} = 3 \cdot k_{ribosome}\) based on data from Proshkin et al. 2010:
where:
- \(\kappa_{\tau}\) and \(r_0\) are phenomenological parameters from Scott et. al. 2010 that describe the linear relationship between the observed RNA/protein ratio of E. coli and its growth rate (\(\mu\))
- \(c_{ribo} = \frac{m_{rr}}{f_{rRNA}\cdot m_{aa}}\) where: \(m_{rr}\) is the mass of rRNA per ribosome. \(f_{rRNA}\) is the fraction of total RNA that is rRNA \(m_{aa}\) is the molecular weight of an average amino acid
- \(v_{transcription, j}\) is the rate of transcription for \(TU_j\)
- \(l_{TU, j}\) is number of nucleotides in \(TU_j\)
3.3.2. Add Translation Reaction¶
3.3.2.1. Add TranslationData to model¶
In order to add a TranslationData object to a ME-model the user must additionally specifify the mRNA id and protein id of the translation reaction that will be added. This information as well as a nucleotide sequence is the only information required to add a translation reaction.
-
class
cobrame.core.processdata.
TranslationData
(id, model, mrna, protein)[source] Class for storing information about a translation reaction.
Parameters: - id (str) – Identifier of the gene being translated, typically the locus tag
- model (
cobrame.core.model.MEModel
) – ME-model that the TranslationData is associated with - mrna (str) – ID of the mRNA that is being translated
- protein (str) – ID of the protein product.
-
mRNA
¶ str – ID of the mRNA that is being translated
-
protein
¶ str – ID of the protein product.
-
subreactions
¶ collections.DefaultDict(int)
– Dictionary of {cobrame.core.processdata.SubreactionData.id
: num_usages} required for the mRNA to be translated
-
nucleotide_sequence
¶ str – String of base pair abbreviations for nucleotides contained in the gene being translated
In [12]:
data = cobrame.TranslationData('a', me, 'RNA_a', 'protein_a')
data.nucleotide_sequence = sequence
3.3.2.2. Add TranslationReaction to model¶
By associating the TranslationReaction with its corresponding TranslationData object and running the update function, COBRAme will create a reaction reaction for the nucleotide sequence given based on the organisms codon table and prespecified translation machinery.
-
class
cobrame.core.reaction.
TranslationReaction
(id)[source] Reaction class for the translation of a TranscribedGene to a TranslatedGene
Parameters: id (str) – Identifier of the translation reaction. As a best practice, this ID should be prefixed with ‘translation + _’
In [13]:
rxn = cobrame.TranslationReaction('translation_a')
rxn.translation_data = data
me.add_reaction(rxn)
3.3.2.3. Update TranslationReaction¶
-
TranslationReaction.
update
(verbose=True)[source] Creates reaction using the associated translation data and adds chemical formula to protein product
This function adds the following components to the reaction stoichiometry (using ‘data’ as shorthand for
cobrame.core.processdata.TranslationData
):- Amino acids defined in data.amino_acid_sequence. Subtracting water to account for condensation reactions during polymerization
- Ribosome w/ translation coupling coefficient (if present)
- mRNA defined in data.mRNA w/ translation coupling coefficient
- mRNA + nucleotides + hydrolysis ATP cost w/ degradation coupling coefficient (if kdeg (defined in model.global_info) > 0)
- RNA_degradosome w/ degradation coupling coefficient (if present and kdeg > 0)
- Protein product defined in data.protein
- Subreactions defined in data.subreactions
- protein_biomass
cobrame.core.component.Constraint
corresponding to the protein product’s mass - Subtract mRNA_biomass
cobrame.core.component.Constraint
defined by mRNA degradation coupling coefficinet (if kdeg > 0)
Parameters: verbose (bool) – Prints when new metabolites are added to the model when executing update()
In [14]:
rxn.update()
print(rxn.reaction)
/home/sbrg-cjlloyd/cobrame/cobrame/core/reaction.py:1051 UserWarning: ribosome not found
/home/sbrg-cjlloyd/cobrame/cobrame/core/reaction.py:1094 UserWarning: RNA_degradosome not found
0.000498399634202103*mu + 0.000195123456790123 + 0.00598079561042524*(mu + 0.3915)/mu RNA_a + 12 asp__L_c + 0.276611796982167*(mu + 0.3915)/mu atp_c + 0.353897348367627*(mu + 0.3915)/mu mRNA_biomass + met__L_c + 12 phe__L_c + 12 ser__L_c + 12 thr__L_c + 12 tyr__L_c --> 0.276611796982167*(mu + 0.3915)/mu adp_c + 0.514348422496571*(mu + 0.3915)/mu amp_c + 0.227270233196159*(mu + 0.3915)/mu cmp_c + 0.0717695473251029*(mu + 0.3915)/mu gmp_c + 60.0 - 0.276611796982167*(mu + 0.3915)/mu h2o_c + 0.276611796982167*(mu + 0.3915)/mu h_c + 0.276611796982167*(mu + 0.3915)/mu pi_c + protein_a + 7.500606 protein_biomass + 0.299039780521262*(mu + 0.3915)/mu ump_c
In this case the constraint couples the formation of a mRNA metabolite to its translation flux. This constraint is formulated as in O’brien et. al. 2013:
where:
- \(\kappa_{\tau}\) and \(r_0\) are phenomenological parameters from Scott et. al. 2010 that describe the linear relationship between the observed RNA/protein ratio of E. coli and its growth rate (\(\mu\))
- \(c_{mRNA} = \frac{m_{nt}}{f_{mRNA}\cdot m_{aa}}\) where: \(m_{nt}\) is the molecular weight of an average mRNA nucleotide. \(f_{mRNA}\) is the fraction of total RNA that is mRNA \(m_{aa}\) is the molecular weight of an average amino acid
- \(v_{translation, j}\) is the rate of translation for \(mRNA_j\)
3.3.2.4. Incorporate Ribosome¶
-
class
cobrame.core.component.
Ribosome
(id)[source] Metabolite class for Ribosome complexes. Inherits from
cobrame.core.component.Complex
Parameters: id (str) – Identifier of the Ribosome.
In [15]:
ribosome = cobrame.Ribosome('ribosome')
me.add_metabolites([ribosome])
me.reactions.translation_a.update()
print(me.reactions.translation_a.reaction)
0.000498399634202103*mu + 0.000195123456790123 + 0.00598079561042524*(mu + 0.3915)/mu RNA_a + 12 asp__L_c + 0.276611796982167*(mu + 0.3915)/mu atp_c + 0.353897348367627*(mu + 0.3915)/mu mRNA_biomass + met__L_c + 12 phe__L_c + 0.000874533914506385*mu + 0.00034238002752925 ribosome + 12 ser__L_c + 12 thr__L_c + 12 tyr__L_c --> 0.276611796982167*(mu + 0.3915)/mu adp_c + 0.514348422496571*(mu + 0.3915)/mu amp_c + 0.227270233196159*(mu + 0.3915)/mu cmp_c + 0.0717695473251029*(mu + 0.3915)/mu gmp_c + 60.0 - 0.276611796982167*(mu + 0.3915)/mu h2o_c + 0.276611796982167*(mu + 0.3915)/mu h_c + 0.276611796982167*(mu + 0.3915)/mu pi_c + protein_a + 7.500606 protein_biomass + 0.299039780521262*(mu + 0.3915)/mu ump_c
/home/sbrg-cjlloyd/cobrame/cobrame/core/reaction.py:1094 UserWarning: RNA_degradosome not found
This imposes a new coupling constraint for the ribosome. In this case the constraint couples the formation of a ribosome to its translation flux. This constraint is formulated as in O’brien et. al. 2013:
where:
- \(\kappa_{\tau}\) and \(r_0\) are phenomenological parameters from Scott et. al. 2010 that describe the linear relationship between the observed RNA/protein ratio of E. coli and its growth rate (\(\mu\))
- \(c_{ribo} = \frac{m_{rr}}{f_{rRNA}\cdot m_{aa}}\) where: \(m_{nt}\) is the mass of rRNA per ribosome. \(f_{rRNA}\) is the fraction of total RNA that is rRNA \(m_{aa}\) is the molecular weight of an average amino acid
- \(v_{translation, j}\) is the rate of translation for \(mRNA_j\)
- \(l_{p, j}\) is number of amino acids in peptide translated from \(mRNA_j\)
Below, we’ll correct this by adding in an tRNA charging reaction.
3.3.3. Add tRNA Charging Reaction¶
3.3.3.1. Add tRNAData to model¶
In [16]:
# Must add tRNA metabolite first
gene = cobrame.TranscribedGene('RNA_d', 'tRNA', sequence)
me.add_metabolites([gene])
-
class
cobrame.core.processdata.
tRNAData
(id, model, amino_acid, rna, codon)[source] Class for storing information about a tRNA charging reaction.
Parameters: - id (str) – Identifier for tRNA charging process. As best practice, this should be follow “tRNA + _ + <tRNA_locus> + _ + <codon>” template. If tRNA initiates translation, <codon> should be replaced with START.
- model (
cobrame.core.model.MEModel
) – ME-model that the tRNAData is associated with - amino_acid (str) – Amino acid that the tRNA transfers to an peptide
- rna (str) – ID of the uncharged tRNA metabolite. As a best practice, this ID should be prefixed with ‘RNA + _’
-
subreactions
¶ collections.DefaultDict(int)
– Dictionary of {cobrame.core.processdata.SubreactionData.id
: num_usages} required for the tRNA to be charged
-
synthetase
¶ str – ID of the tRNA synthetase required to charge the tRNA with an amino acid
-
synthetase_keff
¶ float – Effective turnover rate of the tRNA synthetase
In [17]:
data = cobrame.tRNAData('tRNA_d_GUA', me, 'val__L_c', 'RNA_d', 'GUA')
3.3.3.2. Add tRNAChargingReaction to model¶
And point tRNAChargingReaction
to tRNAData
-
class
cobrame.core.reaction.
tRNAChargingReaction
(id)[source] Reaction class for the charging of a tRNA with an amino acid
Parameters: id (str) – Identifier for the charging reaction. As a best practice, ID should follow the template “charging_tRNA + _ + <tRNA_locus> + _ + <codon>”. If tRNA initiates translation, <codon> should be replaced with START.
In [18]:
rxn = cobrame.tRNAChargingReaction('charging_tRNA_d_GUA')
me.add_reaction(rxn)
rxn.tRNA_data = data
3.3.3.3. Update tRNAChargingReaction¶
-
tRNAChargingReaction.
update
(verbose=True)[source] Creates reaction using the associated tRNA data
This function adds the following components to the reaction stoichiometry (using ‘data’ as shorthand for
cobrame.core.processdata.tRNAData
):- Charged tRNA product following template: “generic_tRNA + _ + <data.codon> + _ + <data.amino_acid>”
- tRNA metabolite (defined in data.RNA) w/ charging coupling coefficient
- Charged amino acid (defined in data.amino_acid) w/ charging coupling coefficient
- Synthetase (defined in data.synthetase) w/ synthetase coupling coefficient found, in part, using data.synthetase_keff
- Post transcriptional modifications defined in data.subreactions
Parameters: verbose (bool) – Prints when new metabolites are added to the model when executing update()
In [19]:
#Setting verbose=False suppresses print statements indicating that new metabolites were created
rxn.update(verbose=False)
print(rxn.reaction)
0.000116266666666667*mu + 4.55184e-5 RNA_d + 0.000116266666666667*mu + 4.55184e-5 val__L_c --> generic_tRNA_GUA_val__L_c
This reaction creates one generic_charged_tRNA
equivalement that can
then be used in a translation reaction
The coefficient for RNA_d
and lys__L_c
are defined by:
where:
- \(\kappa_{\tau}\) and \(r_0\) are phenomenological parameters from Scott et. al. 2010 that describe the linear relationship between the observed RNA/protein ratio of E. coli and its growth rate (\(\mu\))
- \(c_{tRNA, j} = \frac{m_{tRNA}}{f_{tRNA}\cdot m_{aa}}\) where: \(m_{tRNA}\) is molecular weight of an average tRNA. \(f_{tRNA}\) is the fraction of total RNA that is tRNA \(m_{aa}\) is the molecular weight of an average amino acid
- \(v_{charging, j}\) is the rate of charging for \(tRNA_j\)
3.3.3.4. Incorporate tRNA Synthetases¶
.. autoclass:: cobrame.core.component.Complex :noindex:
In [20]:
synthetase = cobrame.Complex('synthetase')
me.add_metabolites(synthetase)
Associate synthetase with tRNAData
and update
In [21]:
data.synthetase = synthetase.id
rxn.update()
print(rxn.reaction)
0.000116266666666667*mu + 4.55184e-5 RNA_d + 4.27350427350427e-6*mu*(0.000116266666666667*mu + 1.0000455184) synthetase + 0.000116266666666667*mu + 4.55184e-5 val__L_c --> generic_tRNA_GUA_val__L_c
The synthetase coupling was reformulated from O’brien et. al. 2013 enable more modularity in the ME-model. A more complete mathematical description of the tRNA synthetase coupling constraints can be found in the tRNA.ipynb
3.3.4. Add tRNAs to Translation¶
Here we take advantage of an additional subclass of ProcessData
,
called a SubreactionData
object. This class is used to lump together
processeses that occur as a result of many individual reactions,
including translation elongation, ribosome formation, tRNA modification,
etc. Since each of these steps often involve an enzyme that requires its
own coupling constraint, this process allows these processes to be
lumped into one reaction while still enabling each subprocess to be
modified.
TranslationData
objects have an subreaction_from_sequence
method
that returns any subreactions that have been added to the model and are
part of translation elongation (i.e. tRNA). Since no tRNA-mediated amino
acid addition subreactions have been added to the model, the below call
returns nothing.
In [22]:
print(me.process_data.a.subreactions_from_sequence)
{}
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction phe_addition_at_UUU not in model
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction tyr_addition_at_UAU not in model
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction thr_addition_at_ACG not in model
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction asp_addition_at_GAU not in model
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction ser_addition_at_AGU not in model
UserWarnings
are returned to indicate that tRNA subreactions have
not been added for each codon.
Below, we add the SubreactionData (excluding enzymes) for the addition
of an amino acid using information from the E. coli codon table. The
charge tRNA does not act as an enzyme in this case because it’s coupling
is handled in the tRNAChargingReaction
3.3.4.1. Add Subreactions for tRNA addition to model¶
-
class
cobrame.core.processdata.
SubreactionData
(id, model)[source] Parameters: - id (str) – Identifier of the subreaction data. As a best practice, if the subreaction data details a modification, the ID should be prefixed with “mod + _”
- model (
cobrame.core.model.MEModel
) – ME-model that the SubreactionData is associated with
-
enzyme
¶ list or str or None – List of
cobrame.core.component.Complex.id
s for enzymes that catalyze this processor
String of single
cobrame.core.component.Complex.id
for enzyme that catalyzes this process
-
keff
¶ float – Effective turnover rate of enzyme(s) in subreaction process
-
_element_contribution
¶ dict – If subreaction adds a chemical moiety to a macromolecules via a modification or other means, net element contribution of the modification process should be accounted for. This can be used to mass balance check each of the individual processes.
Dictionary of {element: net_number_of_contributions}
In [23]:
data = cobrame.SubreactionData('asp_addition_at_GAU', me)
data.stoichiometry = {'generic_tRNA_GAU_asp__L_c': -1,
'gtp_c': -1, 'gdp_c': 1, 'h_c': 1,
'pi_c': 1}
Now calling subreactions_from_sequence
returns the number of tRNA
subreactions that should be added to the TranslationData
In [24]:
translation_subreactions = me.process_data.a.subreactions_from_sequence
print(translation_subreactions)
{'asp_addition_at_GAU': 12}
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction phe_addition_at_UUU not in model
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction tyr_addition_at_UAU not in model
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction thr_addition_at_ACG not in model
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:826 UserWarning: tRNA addition subreaction ser_addition_at_AGU not in model
Updating TranslationData.subreactions
with the tRNA subreactions
incorporates this information into the TranslationReaction
In [25]:
print("Before adding tRNA subreaction")
print("------------------------------")
print(me.reactions.translation_a.reaction)
print("")
# Link tranlation_data to subreactions and update
for subreaction, value in translation_subreactions.items():
me.process_data.a.subreactions[subreaction] = value
me.reactions.translation_a.update(verbose=False)
print("After adding tRNA subreaction")
print("-----------------------------")
print(me.reactions.translation_a.reaction)
Before adding tRNA subreaction
------------------------------
0.000498399634202103*mu + 0.000195123456790123 + 0.00598079561042524*(mu + 0.3915)/mu RNA_a + 12 asp__L_c + 0.276611796982167*(mu + 0.3915)/mu atp_c + 0.353897348367627*(mu + 0.3915)/mu mRNA_biomass + met__L_c + 12 phe__L_c + 0.000874533914506385*mu + 0.00034238002752925 ribosome + 12 ser__L_c + 12 thr__L_c + 12 tyr__L_c --> 0.276611796982167*(mu + 0.3915)/mu adp_c + 0.514348422496571*(mu + 0.3915)/mu amp_c + 0.227270233196159*(mu + 0.3915)/mu cmp_c + 0.0717695473251029*(mu + 0.3915)/mu gmp_c + 60.0 - 0.276611796982167*(mu + 0.3915)/mu h2o_c + 0.276611796982167*(mu + 0.3915)/mu h_c + 0.276611796982167*(mu + 0.3915)/mu pi_c + protein_a + 7.500606 protein_biomass + 0.299039780521262*(mu + 0.3915)/mu ump_c
After adding tRNA subreaction
-----------------------------
0.000498399634202103*mu + 0.000195123456790123 + 0.00598079561042524*(mu + 0.3915)/mu RNA_a + 12 asp__L_c + 0.276611796982167*(mu + 0.3915)/mu atp_c + 12.0 generic_tRNA_GAU_asp__L_c + 12.0 gtp_c + 0.353897348367627*(mu + 0.3915)/mu mRNA_biomass + met__L_c + 12 phe__L_c + 0.000874533914506385*mu + 0.00034238002752925 ribosome + 12 ser__L_c + 12 thr__L_c + 12 tyr__L_c --> 0.276611796982167*(mu + 0.3915)/mu adp_c + 0.514348422496571*(mu + 0.3915)/mu amp_c + 0.227270233196159*(mu + 0.3915)/mu cmp_c + 12.0 gdp_c + 0.0717695473251029*(mu + 0.3915)/mu gmp_c + 60.0 - 0.276611796982167*(mu + 0.3915)/mu h2o_c + 12.0 + 0.276611796982167*(mu + 0.3915)/mu h_c + 12.0 + 0.276611796982167*(mu + 0.3915)/mu pi_c + protein_a + 7.500606 protein_biomass + 0.299039780521262*(mu + 0.3915)/mu ump_c
/home/sbrg-cjlloyd/cobrame/cobrame/core/reaction.py:1094 UserWarning: RNA_degradosome not found
/home/sbrg-cjlloyd/cobrame/cobrame/core/processdata.py:229 UserWarning: No element contribution input for subreaction (asp_addition_at_GAU), calculating based on stoichiometry instead
3.3.5. Add Complex Formation Reaction¶
3.3.5.1. Add ComplexData to model¶
For COBRAme models, the reaction gene-protein-reaction rule (GPR) is replaced with a metabolite representing the synthesis of the enzyme(s) catalyzing a reaction. This metabolite is formed explicitly in a ME model by seperate reaction to transcribe the gene(s) and translate the protein(s) the compose the complex.
-
class
cobrame.core.processdata.
ComplexData
(id, model)[source] Contains all information associated with the formation of an functional enzyme complex.
This can include any enzyme complex modifications required for the enzyme to become active.
Parameters: - id (str) – Identifier of the complex data. As a best practice, this should typically use the same ID as the complex being formed. In cases with multiple ways to form complex ‘_ + alt’ or similar suffixes can be used.
- model (
cobrame.core.model.MEModel
) – ME-model that the ComplexData is associated with
-
stoichiometry
¶ collections.DefaultDict(int)
– Dictionary containing {protein_id: count} for all protein subunits comprising enzyme complex
-
subreactions
¶ dict – Dictionary of {subreaction_data_id: count} for all complex formation subreactions/modifications. This can include cofactor/prosthetic group binding or enzyme side group addition.
In [26]:
data = cobrame.ComplexData('complex_ab', me)
data.stoichiometry = {'protein_a': 1, 'protein_b': 1}
3.3.5.2. Add ComplexFormation reaction to model¶
And point ComplexFormation
to ComplexData
In [27]:
rxn = cobrame.ComplexFormation('formation_complex_ab')
me.add_reaction(rxn)
rxn.complex_data_id = data.id
rxn._complex_id = data.id
3.3.5.3. Update ComplexFormation reaction¶
-
ComplexFormation.
update
(verbose=True)[source] Creates reaction using the associated complex data and adds chemical formula to complex metabolite product.
This function adds the following components to the reaction stoichiometry (using ‘data’ as shorthand for
cobrame.core.processdata.ComplexData
):- Complex product defined in self._complex_id
- Protein subunits with stoichiometery defined in data.stoichiometry
- Metabolites and enzymes w/ coupling coefficients defined in data.subreactions. This often includes enzyme complex modifications by coenzymes or prosthetic groups.
- Biomass
cobrame.core.component.Constraint
corresponding to modifications detailed in data.subreactions, if any
Parameters: verbose (bool) – Prints when new metabolites are added to the model when executing update()
In [28]:
rxn.update(verbose=False)
print(me.reactions.formation_complex_ab.reaction)
protein_a + protein_b --> complex_ab
3.3.5.4. Apply modification to complex formation reaction¶
Many enzyme complexes in an ME-model require cofactors or prosthetic groups in order to properly function. Information about such processes are stored as ModificationData.
For instance, we can add the modification of an iron-sulfur cluster, a common prosthetic group, by doing the following:
In [29]:
# Define the stoichiometry of the modification
mod_data = cobrame.SubreactionData('mod_2fe2s_c', me)
mod_data.stoichiometry = {'2fe2s_c': -1}
# this process can also be catalyzed by a chaperone
mod_data.enzyme = 'complex_ba'
mod_data.keff = 65. # default value
Associate modification to complex and update()
its formation
In [30]:
complex_data = me.process_data.complex_ab
complex_data.subreactions['mod_2fe2s_c'] = 1
Update ComplexFormation reaction
In [31]:
print('Before adding modification')
print('--------------------------')
print(me.reactions.formation_complex_ab.reaction)
me.reactions.formation_complex_ab.update()
print('\n')
print('After adding modification')
print('-------------------------')
print(me.reactions.formation_complex_ab.reaction)
Before adding modification
--------------------------
protein_a + protein_b --> complex_ab
Created <Complex complex_ba at 0x7f4bf69a8b38> in <ComplexFormation formation_complex_ab at 0x7f4bf5f13ef0>
After adding modification
-------------------------
2fe2s_c + 4.27350427350427e-6*mu complex_ba + protein_a + protein_b --> complex_ab + 0.17582 prosthetic_group_biomass
3.3.6. Add Metabolic Reaction¶
3.3.6.1. Add StoichiometricData to model¶
MetabolicReactions require, at a minimum, one corresponding StoichiometricData. StoichiometricData essentially holds the information contained in an M-model reaction. This includes the metabolite stoichiometry and the upper and lower bound of the reaction. As a best practice, StoichiometricData typically uses an ID equivalent to the M-model reaction ID.
So first, we will create a StoichiometricData object to define the stoichiometry of the conversion of a to b. Only one StoichiometricData object should be created for both reversible and irreversible reactions
-
class
cobrame.core.processdata.
StoichiometricData
(id, model)[source] Encodes the stoichiometry for a metabolic reaction.
StoichiometricData defines the metabolite stoichiometry and upper/lower bounds of metabolic reaction
Parameters: - id (str) – Identifier of the metabolic reaction. Should be identical to the M-model reactions in most cases.
- model (
cobrame.core.model.MEModel
) – ME-model that the StoichiometricData is associated with
-
_stoichiometry
¶ dict – Dictionary of {metabolite_id: stoichiometry} for reaction
-
subreactions
¶ collections.DefaultDict(int)
– Cases where multiple enzymes (often carriers ie. Acyl Carrier Protein) are involved in a metabolic reactions.
-
upper_bound
¶ int – Upper reaction bound of metabolic reaction. Should be identical to the M-model reactions in most cases.
-
lower_bound
¶ int – Lower reaction bound of metabolic reaction. Should be identical to the M-model reactions in most cases.
In [32]:
# unique to COBRAme, construct a stoichiometric data object with the reaction information
data = cobrame.StoichiometricData('a_to_b', me)
stoichiometry = {'a':-1, 'b': 1}
data._stoichiometry = stoichiometry
data.lower_bound = -1000
data.upper_bound = 1000
3.3.6.2. Add MetabolicReaction to model¶
The StoichiometricData for this reversible reaction is then assigned to two different MetabolicReactions (Due to the enzyme dilution constraint, all enzyme catalyzed reactions must be reverisble; more on this later). The MetabolicReactions require: - The associated StoichiometricData - The reverse flag set to True for reverse reactions, False for forward reactions - Enzyme \(K_{eff}\) for reaction (discussed later, dafault=65)
These fields are then processed and the actual model reaction is created using the MetabolicReaction’s update() function
-
class
cobrame.core.reaction.
MetabolicReaction
(id)[source] Irreversible metabolic reaction including required enzymatic complex
This reaction class’s update function processes the information contained in the complex data for the enzyme that catalyzes this reaction as well as the stoichiometric data which contains the stoichiometry of the metabolic conversion being performed (i.e. the stoichiometry of the M-model reaction analog)
Parameters: id (str) – Identifier of the metabolic reaction. As a best practice, this ID should use the following template (FWD=forward, REV=reverse): “<StoichiometricData.id> + _ + <FWD or REV> + _ + <Complex.id>” -
keff
¶ float – The turnover rete (keff) couples enzymatic dilution to metabolic flux
-
reverse
¶ boolean – If True, the reaction corresponds to the reverse direction of the reaction. This is necessary since all reversible enzymatic reactions in an ME-model are broken into two irreversible reactions
-
In [33]:
# Create a forward ME Metabolic Reaction and associate the stoichiometric data to it
rxn_fwd = cobrame.MetabolicReaction('a_to_b_FWD_complex_ab')
me.add_reaction(rxn_fwd)
rxn_fwd.stoichiometric_data = data
rxn_fwd.reverse = False
rxn_fwd.keff = 65.
# Create a reverse ME Metabolic Reaction and associate the stoichiometric data to it
rxn_rev = cobrame.MetabolicReaction('a_to_b_REV_complex_ab')
me.add_reaction(rxn_rev)
rxn_rev.stoichiometric_data = data
rxn_rev.reverse = True
rxn_rev.keff = 65.
3.3.6.3. Update MetabolicReactions¶
-
MetabolicReaction.
update
(verbose=True)[source] Creates reaction using the associated stoichiometric data and complex data.
This function adds the following components to the reaction stoichiometry (using ‘data’ as shorthand for
cobrame.core.processdata.StoichiometricData
):- Complex w/ coupling coefficients defined in self.complex_data.id and self.keff
- Metabolite stoichiometry defined in data.stoichiometry. Sign is flipped if self.reverse == True
Also sets the lower and upper bounds based on self.reverse and data.upper_bound and data.lower_bound.
Parameters: verbose (bool) – Prints when new metabolites are added to the model when executing update()
In [34]:
rxn_fwd.update(verbose=False)
rxn_rev.update(verbose=False)
print(me.reactions.a_to_b_FWD_complex_ab.reaction)
print(me.reactions.a_to_b_REV_complex_ab.reaction)
a --> b
b --> a
complex_ab
is not included in the
reaction since no complex has been associated to it yet3.3.6.4. Associate enzyme with MetabolicReaction¶
The ComplexData object created in the previous cell can be incorporated into the MetabolicReaction using code below.
In [35]:
data = me.process_data.complex_ab
me.reactions.a_to_b_FWD_complex_ab.complex_data = data
print('Forward reaction (before update): %s' %
(me.reactions.a_to_b_FWD_complex_ab.reaction))
me.reactions.a_to_b_FWD_complex_ab.update()
print('Forward reaction (after update): %s' %
(me.reactions.a_to_b_FWD_complex_ab.reaction))
print('')
me.reactions.a_to_b_REV_complex_ab.complex_data = data
print('Reverse reaction (before update): %s' %
(me.reactions.a_to_b_REV_complex_ab.reaction))
me.reactions.a_to_b_REV_complex_ab.update()
print('Reverse reaction (after update): %s' %
(me.reactions.a_to_b_REV_complex_ab.reaction))
Forward reaction (before update): a --> b
Forward reaction (after update): a + 4.27350427350427e-6*mu complex_ab --> b
Reverse reaction (before update): b --> a
Reverse reaction (after update): b + 4.27350427350427e-6*mu complex_ab --> a
The coefficient for complex_ab is determined by the expression
which in its entirety represents the dilution of an enzyme following a cell doubling. The coupling constraint can be summarized as followed
Where
- \(v_{usage,i}\) is the flux through the metabolic reaction
- \(k_{eff}\) is the turnover rate for the process and conveys the productivity of the enzyme complex. Physically, it can be thought of as the number of reactions the enzyme can catalyze per cell division.
By default the \(k_{eff}\) for a MetabolicReaction is set to 65 but this can be changed using the code below.
3.3.6.5. Different Keff for forward reaction¶
In [36]:
me.reactions.a_to_b_FWD_complex_ab.keff = .00001
me.reactions.a_to_b_FWD_complex_ab.update()
# The forward and reverse direction can have differing keffs
print('Forward reaction')
print('----------------')
print(me.reactions.a_to_b_FWD_complex_ab.reaction)
print('')
print('Reverse reaction')
print('----------------')
print(me.reactions.a_to_b_REV_complex_ab.reaction)
Forward reaction
----------------
a + 27.7777777777778*mu complex_ab --> b
Reverse reaction
----------------
b + 4.27350427350427e-6*mu complex_ab --> a
3.4. Adding Reactions using utility functions¶
Add reactions using some of the utility functions provided in
cobrame.util.building.py
3.4.1. Transcription¶
Using the utility functions to create the TranscribedGene metabolite has the advantage of forcing the assignment of sequence, strand and RNA_type.
-
cobrame.util.building.
create_transcribed_gene
(me_model, locus_id, rna_type, seq, left_pos=None, right_pos=None, strand=None)[source] - Creates a TranscribedGene metabolite object and adds it to the ME-model
Parameters: - me_model (
cobrame.core.model.MEModel
) – The MEModel object to which the reaction will be added - locus_id (str) – Locus ID of RNA product. The TranscribedGene will be added as “RNA + _ + locus_id”
- left_pos (int or None) – Left position of gene on the sequence of the (+) strain
- right_pos (int or None) – Right position of gene on the sequence of the (+) strain
- seq (str) – Nucleotide sequence of RNA product. Amino acid sequence, codon counts, etc. will be calculated based on this string.
- strand (str or None) –
- (+) if the RNA product is on the leading strand
- (-) if the RNA product is on the complementary strand
- rna_type (str) – Type of RNA of the product. tRNA, rRNA, or mRNA Used for determining how RNA product will be processed.
Returns: Metabolite object for the RNA product
Return type: - me_model (
In [37]:
building.create_transcribed_gene(me, 'b','tRNA', 'ATCG')
building.add_transcription_reaction(me, 'TU_b', {'b'}, sequence)
print(me.reactions.transcription_TU_b.reaction)
me.reactions.transcription_TU_b.update()
86 atp_c + 38 ctp_c + 12 gtp_c + 182 h2o_c + 50 utp_c --> RNA_b + 85 amp_c + 37 cmp_c + 11 gmp_c + 182 h_c + 186 ppi_c + 1.2817349999999998 tRNA_biomass + 49 ump_c
/home/sbrg-cjlloyd/cobrame/cobrame/core/reaction.py:813 UserWarning: RNA Polymerase () not found
3.4.2. Translation¶
add_translation_reaction
assumes that the RNA and protein have the
same locus_id. It creates the appropriate TranslationData
and
TranslationReaction
instance, links the two together and updates the
TranslationReaction
.
-
cobrame.util.building.
add_translation_reaction
(me_model, locus_id, dna_sequence, update=False)[source] Creates and adds a TranslationReaction to the ME-model as well as the associated TranslationData
A dna_sequence is required in order to add a TranslationReaction to the ME-model
Parameters: - me_model (
cobra.core.model.MEModel
) – The MEModel object to which the reaction will be added - locus_id (str) – Locus ID of RNA product. The TranslationReaction will be added as “translation + _ + locus_id” The TranslationData will be added as “locus_id”
- dna_sequence (str) – DNA sequence of the RNA product. This string should be reverse transcribed if it originates on the complement strand.
- update (bool) – If True, use TranslationReaction’s update function to update and add reaction stoichiometry
- me_model (
In [38]:
building.add_translation_reaction(me, 'b', dna_sequence=sequence, update=True)
print(me.reactions.translation_b.reaction)
/home/sbrg-cjlloyd/cobrame/cobrame/core/reaction.py:1094 UserWarning: RNA_degradosome not found
0.000498399634202103*mu + 0.000195123456790123 + 0.00598079561042524*(mu + 0.3915)/mu RNA_b + 12 asp__L_c + 0.00448559670781893*(mu + 0.3915)/mu atp_c + 0.0076657950617284*(mu + 0.3915)/mu mRNA_biomass + met__L_c + 12 phe__L_c + 0.000874533914506385*mu + 0.00034238002752925 ribosome + 12 ser__L_c + 12 thr__L_c + 12 tyr__L_c --> 0.00448559670781893*(mu + 0.3915)/mu adp_c + 0.00598079561042524*(mu + 0.3915)/mu amp_c + 0.00598079561042524*(mu + 0.3915)/mu cmp_c + 0.00598079561042524*(mu + 0.3915)/mu gmp_c + 60.0 - 0.00448559670781893*(mu + 0.3915)/mu h2o_c + 0.00448559670781893*(mu + 0.3915)/mu h_c + 0.00448559670781893*(mu + 0.3915)/mu pi_c + protein_b + 7.500606 protein_biomass + 0.00598079561042524*(mu + 0.3915)/mu ump_c
3.4.3. Complex Formation¶
Alternatively, ComplexData has a create_complex_formation()
function
to create the sythesis reaction following the naming conventions. It
contains an update()
function which incorporates changes in the
ComplexData
-
ComplexData.
create_complex_formation
(verbose=True)[source] creates a complex formation reaction
This assumes none exists already. Will create a reaction (prefixed by “formation”) which forms the complex
Parameters: verbose (bool) – If True, print if a metabolite is added to model during update
In [39]:
data = cobrame.ComplexData('complex_ba', me)
data.stoichiometry = {'protein_a': 1, 'protein_b': 1}
data.create_complex_formation()
print(me.reactions.formation_complex_ba.reaction)
protein_a + protein_b --> complex_ba
3.4.4. Metabolic Reaction¶
-
cobrame.util.building.
add_metabolic_reaction_to_model
(me_model, stoichiometric_data_id, directionality, complex_id=None, spontaneous=False, update=False, keff=65)[source] Creates and add a MetabolicReaction to a MEModel.
Parameters: - me_model (
cobrame.core.model.MEModel
) – MEModel that the MetabolicReaction will be added to - stoichiometric_data_id (str) – ID of the StoichiometricData for the reaction being added
- directionality (str) –
- Forward: Add reaction that occurs in the forward direction
- Reverse: Add reaction that occurs in the reverse direction
- complex_id (str or None) – ID of the ComplexData for the enzyme that catalyze the reaction being added.
- spontaneous (bool) –
- If True and complex_id=’’ add reaction as spontaneous reaction
- If False and complex_id=’’ add reaction as orphan (CPLX_dummy catalyzed)
- me_model (
In [40]:
stoich_data = cobrame.StoichiometricData('b_to_c', me)
stoich_data._stoichiometry = {'b': -1, 'c': 1}
stoich_data.lower_bound = 0
stoich_data.upper_bound = 1000.
building.add_metabolic_reaction_to_model(me, stoich_data.id, 'forward', complex_id='complex_ab',
update=True)
print('Reaction b_to_c')
print('---------------')
print(me.reactions.b_to_c_FWD_complex_ab.reaction)
Created <Metabolite c at 0x7f4bf6abf6d8> in <MetabolicReaction b_to_c_FWD_complex_ab at 0x7f4bf6abf7b8>
Reaction b_to_c
---------------
b + 4.27350427350427e-6*mu complex_ab --> c