# 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):

$a \rightarrow b$

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:

1. Transcription reactions
2. Translation reactions
3. tRNA charging reactions
4. 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:
for rxn in ijo.reactions:


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.

Note: k_deg will initially be set to 0. We will apply RNA degradation later in the tutorial.
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
"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:

1. Create the ProcessData(s) associated with the reaction and populate them with the necessary information
2. Create the MEReaction and link the appropriate ProcessData
3. Execute the MEReaction’s update method

#### 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: id (str) – Identifier of the transcribed gene. As a best practice, this ID should be prefixed with ‘RNA + _’ RNA_type (str) – Type of RNA encoded by gene sequence (mRNA, rRNA, tRNA, or ncRNA) nucleotide_sequence (str) – String of base pair abbreviations for nucleotides contained in the gene
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)


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


#### 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):

1. RNA_polymerase from data.RNA_polymerase w/ coupling coefficient (if present)
2. RNA products defined in data.RNA_products
3. Nucleotide reactants defined in data.nucleotide_counts
4. If tRNA or rRNA contained in data.RNA_types, excised base products
5. Metabolites + enzymes w/ coupling coefficients defined in data.subreactions (if present)
6. Biomass cobrame.core.component.Constraint corresponding to data.RNA_products and their associated masses
7. 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

Note: the RNA_polymerase complex is not included in the reaction. This will be added later

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:

$\mu \leq v_{biomass\_dilution} \leq \mu$

A mathematical description of the biomass constraint can be found in Biomass Dilution Constraints in ME-Model Fundamentals.

Note: This is not a complete picture of transcription because the RNA polymerase is missing.

#### 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')


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:

\begin{align} v_{dilution,RNAP, j} = \frac{l_{TU,j}}{3 c_{ribo}\kappa_{\tau}} v_{transcription,j} (\mu+r_0\kappa_{\tau}), & \forall j \in TU \end{align}

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.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


#### 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):

1. Amino acids defined in data.amino_acid_sequence. Subtracting water to account for condensation reactions during polymerization
2. Ribosome w/ translation coupling coefficient (if present)
3. mRNA defined in data.mRNA w/ translation coupling coefficient
4. mRNA + nucleotides + hydrolysis ATP cost w/ degradation coupling coefficient (if kdeg (defined in model.global_info) > 0)
5. RNA_degradosome w/ degradation coupling coefficient (if present and kdeg > 0)
6. Protein product defined in data.protein
7. Subreactions defined in data.subreactions
8. protein_biomass cobrame.core.component.Constraint corresponding to the protein product’s mass
9. 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

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:

\begin{align} v_{dilution,j} = \frac{3}{\kappa_{\tau} c_{mRNA}} \cdot (\mu + \kappa_{\tau} r_0) v_{translation,j} , & & \forall j \in mRNA \end{align}

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.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:

\begin{align} v_{dilution,ribo, j} = \frac{l_{p,j}}{c_{ribo}\kappa_{\tau}} v_{translation,j} (\mu+r_0\kappa_{\tau}) , & \forall j \in mRNA \end{align}

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$$

Note: The above reactions do not provide a complete picture of translation in that it is missing charged tRNAs to facillitate tRNA addition.

Below, we’ll correct this by adding in an tRNA charging reaction.

### 3.3.3. Add tRNA Charging Reaction¶

In [16]:

# Must add tRNA metabolite first
gene = cobrame.TranscribedGene('RNA_d', 'tRNA', sequence)

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 + _ + + _ + ” template. If tRNA initiates translation, 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 + _ + + _ + ”. If tRNA initiates translation, should be replaced with START.
In [18]:

rxn = cobrame.tRNAChargingReaction('charging_tRNA_d_GUA')
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):

1. Charged tRNA product following template: “generic_tRNA + _ + <data.codon> + _ + <data.amino_acid>”
2. tRNA metabolite (defined in data.RNA) w/ charging coupling coefficient
3. Charged amino acid (defined in data.amino_acid) w/ charging coupling coefficient
1. Synthetase (defined in data.synthetase) w/ synthetase coupling coefficient found, in part, using data.synthetase_keff
2. 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:

\begin{align} v_{dilution,j} \geq \frac{1}{\kappa_{\tau} c_{tRNA,j}} (\mu + \kappa_{\tau} r_0) v_{charging,j} , & \forall j \in tRNA \end{align}

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$$
Note: This tRNA charging reaction is still missing a tRNA synthetase which catalyzes the amino acid addition to the uncharged tRNA.

#### 3.3.3.4. Incorporate tRNA Synthetases¶

.. autoclass:: cobrame.core.component.Complex :noindex:
In [20]:

synthetase = cobrame.Complex('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


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

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 process

or

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


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("-----------------------------")
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

-----------------------------
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

.. autoclass:: cobrame.core.reaction.ComplexFormation :noindex:
In [27]:

rxn = cobrame.ComplexFormation('formation_complex_ab')
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):

1. Complex product defined in self._complex_id
2. Protein subunits with stoichiometery defined in data.stoichiometry
3. Metabolites and enzymes w/ coupling coefficients defined in data.subreactions. This often includes enzyme complex modifications by coenzymes or prosthetic groups.
4. 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('-------------------------')
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>

-------------------------
2fe2s_c + 4.27350427350427e-6*mu complex_ba + protein_a + protein_b --> complex_ab + 0.17582 prosthetic_group_biomass


#### 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): “ + _ + + _ +
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')
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')
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):

1. Complex w/ coupling coefficients defined in self.complex_data.id and self.keff
2. 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

Note: the $$k_{eff}$$ and complex_ab is not included in the reaction since no complex has been associated to it yet

#### 3.3.6.4. Associate enzyme with MetabolicReaction¶

The ComplexData object created in the previous cell can be incorporated into the MetabolicReaction using code below.

Note: the update() function is required to apply the change.
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

$\frac{\mu}{k_{eff}}$

which in its entirety represents the dilution of an enzyme following a cell doubling. The coupling constraint can be summarized as followed

\begin{align} &v_{dilution,j} = \mu \sum_{i} \left( \frac{1}{k_{eff,i}} v_{usage,i} \right), & \forall j \in Enzyme \end{align}

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. Metabolite object for the RNA product cobrame.core.component.TranscribedGene
In [37]:

building.create_transcribed_gene(me, 'b','tRNA', 'ATCG')
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
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)
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.

Created <Metabolite c at 0x7f4bf6abf6d8> in <MetabolicReaction b_to_c_FWD_complex_ab at 0x7f4bf6abf7b8>