Documentation for COBRApy¶
For installation instructions, please see INSTALL.rst.
Many of the examples below are viewable as IPython notebooks, which can be viewed at nbviewer.
 {
 “cells”: [
 {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“# Getting Started”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Loading a model and inspecting it”
]
}, {
}, {
“cell_type”: “code”, “execution_count”: 1, “metadata”: {}, “outputs”: [], “source”: [
“from __future__ import print_functionn”, “n”, “import cobran”, “import cobra.testn”, “n”, “# “ecoli” and “salmonella” are also valid argumentsn”, “model = cobra.test.create_test_model(“textbook”)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The reactions, metabolites, and genes attributes of the cobrapy model are a special type of list called a cobra.DictList, and each one is made up of cobra.Reaction, cobra.Metabolite and cobra.Gene objects respectively.”
]
}, {
“cell_type”: “code”, “execution_count”: 2, “metadata”: {
“scrolled”: true
}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“95n”, “72n”, “137n”
]
}
], “source”: [
“print(len(model.reactions))n”, “print(len(model.metabolites))n”, “print(len(model.genes))”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“When using [Jupyter notebook](https://jupyternotebookbeginnerguide.readthedocs.io/en/latest/) this type of information is rendered as a table.”
]
}, {
“cell_type”: “code”, “execution_count”: 3, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/html”: [
“n”, ” <table>n”, ” <tr>n”, ” <td><strong>Name</strong></td>n”, ” <td>e_coli_core</td>n”, ” </tr><tr>n”, ” <td><strong>Memory address</strong></td>n”,
 <<<<<<< HEAD
” <td>0x01158878d0</td>n”,
” </tr><tr>n”, ” <td><strong>Number of metabolites</strong></td>n”, ” <td>72</td>n”, ” </tr><tr>n”, ” <td><strong>Number of reactions</strong></td>n”, ” <td>95</td>n”, ” </tr><tr>n”, ” <td><strong>Objective expression</strong></td>n”, ” <td>1.0*Biomass_Ecoli_core  1.0*Biomass_Ecoli_core_reverse_2cdba</td>n”, ” </tr><tr>n”, ” <td><strong>Compartments</strong></td>n”, ” <td>cytosol, extracellular</td>n”, ” </tr>n”, ” </table>”
], “text/plain”: [
 <<<<<<< HEAD
“<Model e_coli_core at 0x1158878d0>”
]
}, “execution_count”: 3, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“model”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Just like a regular list, objects in the DictList can be retrieved by index. For example, to get the 30th reaction in the model (at index 29 because of [0indexing](https://en.wikipedia.org/wiki/Zerobased_numbering)):”
]
}, {
“cell_type”: “code”, “execution_count”: 4, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/html”: [
“n”, ” <table>n”, ” <tr>n”, ” <td><strong>Reaction identifier</strong></td><td>EX_glu__L_e</td>n”, ” </tr><tr>n”, ” <td><strong>Name</strong></td><td>LGlutamate exchange</td>n”, ” </tr><tr>n”, ” <td><strong>Memory address</strong></td>n”,
 <<<<<<< HEAD
” <td>0x011615e2e8</td>n”,
” </tr><tr>n”, ” <td><strong>Stoichiometry</strong></td>n”, ” <td>n”, ” <p style=’textalign:right’>glu__L_e –> </p>n”, ” <p style=’textalign:right’>LGlutamate –> </p>n”, ” </td>n”, ” </tr><tr>n”, ” <td><strong>GPR</strong></td><td></td>n”, ” </tr><tr>n”, ” <td><strong>Lower bound</strong></td><td>0.0</td>n”, ” </tr><tr>n”, ” <td><strong>Upper bound</strong></td><td>1000.0</td>n”, ” </tr>n”, ” </table>n”, ” “
], “text/plain”: [
 <<<<<<< HEAD
“<Reaction EX_glu__L_e at 0x11615e2e8>”
]
}, “execution_count”: 4, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“model.reactions[29]”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Additionally, items can be retrieved by their id using the DictList.get_by_id() function. For example, to get the cytosolic atp metabolite object (the id is “atp_c”), we can do the following:”
]
}, {
“cell_type”: “code”, “execution_count”: 5, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/html”: [
“n”, ” <table>n”, ” <tr>n”, ” <td><strong>Metabolite identifier</strong></td><td>atp_c</td>n”, ” </tr><tr>n”, ” <td><strong>Name</strong></td><td>ATP</td>n”, ” </tr><tr>n”, ” <td><strong>Memory address</strong></td>n”,
 <<<<<<< HEAD
” <td>0x01160d4630</td>n”,
” </tr><tr>n”, ” <td><strong>Formula</strong></td><td>C10H12N5O13P3</td>n”, ” </tr><tr>n”, ” <td><strong>Compartment</strong></td><td>c</td>n”, ” </tr><tr>n”, ” <td><strong>In 13 reaction(s)</strong></td><td>n”,
 <<<<<<< HEAD
” PPS, ADK1, ATPS4r, GLNS, SUCOAS, GLNabc, PGK, ATPM, PPCK, ACKr, PFK, Biomass_Ecoli_core, PYK</td>n”,
” </tr>n”, ” </table>”
], “text/plain”: [
 <<<<<<< HEAD
“<Metabolite atp_c at 0x1160d4630>”
]
}, “execution_count”: 5, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“model.metabolites.get_by_id(“atp_c”)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“As an added bonus, users with an interactive shell such as IPython will be able to tabcomplete to list elements inside a list. While this is not recommended behavior for most code because of the possibility for characters like “” inside ids, this is very useful while in an interactive prompt:”
]
}, {
“cell_type”: “code”, “execution_count”: 6, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“(10.0, 1000.0)”
]
}, “execution_count”: 6, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“model.reactions.EX_glc__D_e.bounds”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Reactions”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We will consider the reaction glucose 6phosphate isomerase, which interconverts glucose 6phosphate and fructose 6phosphate. The reaction id for this reaction in our test model is PGI.”
]
}, {
“cell_type”: “code”, “execution_count”: 7, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/html”: [
“n”, ” <table>n”, ” <tr>n”, ” <td><strong>Reaction identifier</strong></td><td>PGI</td>n”, ” </tr><tr>n”, ” <td><strong>Name</strong></td><td>glucose6phosphate isomerase</td>n”, ” </tr><tr>n”, ” <td><strong>Memory address</strong></td>n”,
 <<<<<<< HEAD
” <td>0x0116188e48</td>n”,
” </tr><tr>n”, ” <td><strong>Stoichiometry</strong></td>n”, ” <td>n”, ” <p style=’textalign:right’>g6p_c <=> f6p_c</p>n”, ” <p style=’textalign:right’>DGlucose 6phosphate <=> DFructose 6phosphate</p>n”, ” </td>n”, ” </tr><tr>n”, ” <td><strong>GPR</strong></td><td>b4025</td>n”, ” </tr><tr>n”, ” <td><strong>Lower bound</strong></td><td>1000.0</td>n”, ” </tr><tr>n”, ” <td><strong>Upper bound</strong></td><td>1000.0</td>n”, ” </tr>n”, ” </table>n”, ” “
], “text/plain”: [
 <<<<<<< HEAD
“<Reaction PGI at 0x116188e48>”
]
}, “execution_count”: 7, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi = model.reactions.get_by_id(“PGI”)n”, “pgi”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We can view the full name and reaction catalyzed as strings”
]
}, {
“cell_type”: “code”, “execution_count”: 8, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“glucose6phosphate isomerasen”, “g6p_c <=> f6p_cn”
]
}
], “source”: [
“print(pgi.name)n”, “print(pgi.reaction)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We can also view reaction upper and lower bounds. Because the pgi.lower_bound < 0, and pgi.upper_bound > 0, pgi is reversible.”
]
}, {
“cell_type”: “code”, “execution_count”: 9, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“1000.0 < pgi < 1000.0n”, “Truen”
]
}
], “source”: [
“print(pgi.lower_bound, “< pgi <”, pgi.upper_bound)n”, “print(pgi.reversibility)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The lower and upper bound of reactions can also be modified, and the reversibility attribute will automatically be updated. The preferred method for manipulating bounds is using reaction.bounds, e.g.”
]
}, {
“cell_type”: “code”, “execution_count”: 10, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“0 < pgi < 1000.0n”, “Reversibility after modification: Falsen”, “Reversibility after resetting: Truen”
]
}
], “source”: [
“old_bounds = pgi.boundsn”, “pgi.bounds = (0, 1000.0)n”, “print(pgi.lower_bound, “< pgi <”, pgi.upper_bound)n”, “print(“Reversibility after modification:”, pgi.reversibility)n”, “pgi.bounds = old_boundsn”, “print(“Reversibility after resetting:”, pgi.reversibility)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Bounds can also be modified oneatatime using reaction.lower_bound or reaction.upper_bound. This approach is less desirable than setting both bounds simultaneously with reaction.bounds because a user might accidently set a lower bound higher than an upper bound (or vice versa). Currently, cobrapy will automatically adjust the other bound (e.g. the bound the user didn’t manually adjust) so that this violation doesn’t occur, but this feature may be removed in the near future. “
]
}, {
“cell_type”: “code”, “execution_count”: 11, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“Upper bound prior to setting new lower bound: 1000.0n”, “Upper bound after setting new lower bound: 1100n”
]
}
], “source”: [
“old_bounds = pgi.boundsn”, “print(‘Upper bound prior to setting new lower bound:’, pgi.upper_bound)n”, “pgi.lower_bound = 1100n”, “print(‘Upper bound after setting new lower bound:’, pgi.upper_bound)n”, “pgi.bounds = old_bounds”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We can also ensure the reaction is mass balanced. This function will return elements which violate mass balance. If it comes back empty, then the reaction is mass balanced.”
]
}, {
“cell_type”: “code”, “execution_count”: 12, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“{}”
]
}, “execution_count”: 12, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi.check_mass_balance()”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“In order to add a metabolite, we pass in a dict with the metabolite object and its coefficient”
]
}, {
“cell_type”: “code”, “execution_count”: 13, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“‘g6p_c + h_c <=> f6p_c’”
]
}, “execution_count”: 13, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi.add_metabolites({model.metabolites.get_by_id(“h_c”): 1})n”, “pgi.reaction”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The reaction is no longer mass balanced”
]
}, {
“cell_type”: “code”, “execution_count”: 14, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“{‘charge’: 1.0, ‘H’: 1.0}”
]
}, “execution_count”: 14, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi.check_mass_balance()”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We can remove the metabolite, and the reaction will be balanced once again.”
]
}, {
“cell_type”: “code”, “execution_count”: 15, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“g6p_c <=> f6p_cn”, “{}n”
]
}
], “source”: [
“pgi.subtract_metabolites({model.metabolites.get_by_id(“h_c”): 1})n”, “print(pgi.reaction)n”, “print(pgi.check_mass_balance())”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“It is also possible to build the reaction from a string. However, care must be taken when doing this to ensure reaction id’s match those in the model. The direction of the arrow is also used to update the upper and lower bounds.”
]
}, {
“cell_type”: “code”, “execution_count”: 16, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“unknown metabolite ‘green_eggs’ createdn”, “unknown metabolite ‘ham’ createdn”
]
}
], “source”: [
“pgi.reaction = “g6p_c –> f6p_c + h_c + green_eggs + ham””
]
}, {
“cell_type”: “code”, “execution_count”: 17, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“‘g6p_c –> f6p_c + green_eggs + h_c + ham’”
]
}, “execution_count”: 17, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi.reaction”
]
}, {
“cell_type”: “code”, “execution_count”: 18, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“‘g6p_c <=> f6p_c’”
]
}, “execution_count”: 18, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi.reaction = “g6p_c <=> f6p_c”n”, “pgi.reaction”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Metabolites”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We will consider cytosolic atp as our metabolite, which has the id “atp_c” in our test model.”
]
}, {
“cell_type”: “code”, “execution_count”: 19, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/html”: [
“n”, ” <table>n”, ” <tr>n”, ” <td><strong>Metabolite identifier</strong></td><td>atp_c</td>n”, ” </tr><tr>n”, ” <td><strong>Name</strong></td><td>ATP</td>n”, ” </tr><tr>n”, ” <td><strong>Memory address</strong></td>n”,
 <<<<<<< HEAD
” <td>0x01160d4630</td>n”,
” </tr><tr>n”, ” <td><strong>Formula</strong></td><td>C10H12N5O13P3</td>n”, ” </tr><tr>n”, ” <td><strong>Compartment</strong></td><td>c</td>n”, ” </tr><tr>n”, ” <td><strong>In 13 reaction(s)</strong></td><td>n”,
 <<<<<<< HEAD
” PPS, ADK1, ATPS4r, GLNS, SUCOAS, GLNabc, PGK, ATPM, PPCK, ACKr, PFK, Biomass_Ecoli_core, PYK</td>n”,
” </tr>n”, ” </table>”
], “text/plain”: [
 <<<<<<< HEAD
“<Metabolite atp_c at 0x1160d4630>”
]
}, “execution_count”: 19, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“atp = model.metabolites.get_by_id(“atp_c”)n”, “atp”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We can print out the metabolite name and compartment (cytosol in this case) directly as string.”
]
}, {
“cell_type”: “code”, “execution_count”: 20, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“ATPn”, “cn”
]
}
], “source”: [
“print(atp.name)n”, “print(atp.compartment)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We can see that ATP is a charged molecule in our model.”
]
}, {
“cell_type”: “code”, “execution_count”: 21, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“4”
]
}, “execution_count”: 21, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“atp.charge”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We can see the chemical formula for the metabolite as well.”
]
}, {
“cell_type”: “code”, “execution_count”: 22, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“C10H12N5O13P3n”
]
}
], “source”: [
“print(atp.formula)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The reactions attribute gives a frozenset of all reactions using the given metabolite. We can use this to count the number of reactions which use atp.”
]
}, {
“cell_type”: “code”, “execution_count”: 23, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“13”
]
}, “execution_count”: 23, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“len(atp.reactions)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“A metabolite like glucose 6phosphate will participate in fewer reactions.”
]
}, {
“cell_type”: “code”, “execution_count”: 24, “metadata”: {}, “outputs”: [
 {
 “data”: {
“text/plain”: [
 <<<<<<< HEAD
“frozenset({<Reaction Biomass_Ecoli_core at 0x1161337b8>,n”, ” <Reaction G6PDH2r at 0x1160a3a20>,n”, ” <Reaction GLCpts at 0x1160a3da0>,n”, ” <Reaction PGI at 0x116188e48>})”
” <Reaction GLCpts at 0x117a9d0f0>,n”, ” <Reaction PGI at 0x117afacc0>})”
 >>>>>>> origin/devel
]
}, “execution_count”: 24, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“model.metabolites.get_by_id(“g6p_c”).reactions”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Genes”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The gene_reaction_rule is a boolean representation of the gene requirements for this reaction to be active as described in [Schellenberger et al 2011 Nature Protocols 6(9):1290307](http://dx.doi.org/doi:10.1038/nprot.2011.308).n”, “n”, “The GPR is stored as the gene_reaction_rule for a Reaction object as a string.”
]
}, {
“cell_type”: “code”, “execution_count”: 25, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“‘b4025’”
]
}, “execution_count”: 25, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“gpr = pgi.gene_reaction_rulen”, “gpr”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Corresponding gene objects also exist. These objects are tracked by the reactions itself, as well as by the model”
]
}, {
“cell_type”: “code”, “execution_count”: 26, “metadata”: {}, “outputs”: [
 {
 “data”: {
“text/plain”: [
 <<<<<<< HEAD
“frozenset({<Gene b4025 at 0x11610a2b0>})”
]
}, “execution_count”: 26, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi.genes”
]
}, {
“cell_type”: “code”, “execution_count”: 27, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/html”: [
“n”, ” <table>n”, ” <tr>n”, ” <td><strong>Gene identifier</strong></td><td>b4025</td>n”, ” </tr><tr>n”, ” <td><strong>Name</strong></td><td>pgi</td>n”, ” </tr><tr>n”, ” <td><strong>Memory address</strong></td>n”,
 <<<<<<< HEAD
” <td>0x011610a2b0</td>n”,
” </tr><tr>n”, ” <td><strong>Functional</strong></td><td>True</td>n”, ” </tr><tr>n”, ” <td><strong>In 1 reaction(s)</strong></td><td>n”, ” PGI</td>n”, ” </tr>n”, ” </table>”
], “text/plain”: [
 <<<<<<< HEAD
“<Gene b4025 at 0x11610a2b0>”
]
}, “execution_count”: 27, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi_gene = model.genes.get_by_id(“b4025”)n”, “pgi_gene”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Each gene keeps track of the reactions it catalyzes”
]
}, {
“cell_type”: “code”, “execution_count”: 28, “metadata”: {}, “outputs”: [
 {
 “data”: {
“text/plain”: [
 <<<<<<< HEAD
“frozenset({<Reaction PGI at 0x116188e48>})”
]
}, “execution_count”: 28, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi_gene.reactions”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Altering the gene_reaction_rule will create new gene objects if necessary and update all relationships.”
]
}, {
“cell_type”: “code”, “execution_count”: 29, “metadata”: {}, “outputs”: [
 {
 “data”: {
“text/plain”: [
 <<<<<<< HEAD
“frozenset({<Gene eggs at 0x1160245c0>, <Gene spam at 0x116024080>})”
]
}, “execution_count”: 29, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi.gene_reaction_rule = “(spam or eggs)”n”, “pgi.genes”
]
}, {
“cell_type”: “code”, “execution_count”: 30, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“frozenset()”
]
}, “execution_count”: 30, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“pgi_gene.reactions”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Newly created genes are also added to the model”
]
}, {
“cell_type”: “code”, “execution_count”: 31, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/html”: [
“n”, ” <table>n”, ” <tr>n”, ” <td><strong>Gene identifier</strong></td><td>spam</td>n”, ” </tr><tr>n”, ” <td><strong>Name</strong></td><td></td>n”, ” </tr><tr>n”, ” <td><strong>Memory address</strong></td>n”,
 <<<<<<< HEAD
” <td>0x0116024080</td>n”,
” </tr><tr>n”, ” <td><strong>Functional</strong></td><td>True</td>n”, ” </tr><tr>n”, ” <td><strong>In 1 reaction(s)</strong></td><td>n”, ” PGI</td>n”, ” </tr>n”, ” </table>”
], “text/plain”: [
 <<<<<<< HEAD
“<Gene spam at 0x116024080>”
]
}, “execution_count”: 31, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“model.genes.get_by_id(“spam”)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The delete_model_genes function will evaluate the GPR and set the upper and lower bounds to 0 if the reaction is knocked out. This function can preserve existing deletions or reset them using the cumulative_deletions flag.”
]
}, {
“cell_type”: “code”, “execution_count”: 32, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“after 1 KO: 1000 < flux_PGI < 1000n”, “after 2 KO: 0 < flux_PGI < 0n”
]
}
], “source”: [
“cobra.manipulation.delete_model_genes(n”, ” model, [“spam”], cumulative_deletions=True)n”, “print(“after 1 KO: %4d < flux_PGI < %4d” % (pgi.lower_bound, pgi.upper_bound))n”, “n”, “cobra.manipulation.delete_model_genes(n”, ” model, [“eggs”], cumulative_deletions=True)n”, “print(“after 2 KO: %4d < flux_PGI < %4d” % (pgi.lower_bound, pgi.upper_bound))”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The undelete_model_genes can be used to reset a gene deletion”
]
}, {
“cell_type”: “code”, “execution_count”: 33, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“1000 < pgi < 1000n”
]
}
], “source”: [
“cobra.manipulation.undelete_model_genes(model)n”, “print(pgi.lower_bound, “< pgi <”, pgi.upper_bound)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Making changes reversibly using models as contexts”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Quite often, one wants to make small changes to a model and evaluate the impacts of these. For example, we may want to knockout all reactions sequentially, and see what the impact of this is on the objective function. One way of doing this would be to create a new copy of the model before each knockout with model.copy(). However, even with small models, this is a very slow approach as models are quite complex objects. Better then would be to do the knockout, optimizing and then manually resetting the reaction bounds before proceeding with the next reaction. Since this is such a common scenario however, cobrapy allows us to use the model as a context, to have changes reverted automatically.”
]
}, {
“cell_type”: “code”, “execution_count”: 34, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“ACALD blocked (bounds: (0, 0)), new growth rate 0.873922n”, “ACALDt blocked (bounds: (0, 0)), new growth rate 0.873922n”, “ACKr blocked (bounds: (0, 0)), new growth rate 0.873922n”, “ACONTa blocked (bounds: (0, 0)), new growth rate 0.000000n”, “ACONTb blocked (bounds: (0, 0)), new growth rate 0.000000n”
]
}
], “source”: [
“model = cobra.test.create_test_model(‘textbook’)n”, “for reaction in model.reactions[:5]:n”, ” with model as model:n”, ” reaction.knock_out()n”, ” model.optimize()n”, ” print(‘%s blocked (bounds: %s), new growth rate %f’ %n”, ” (reaction.id, str(reaction.bounds), model.objective.value))”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“If we look at those knocked reactions, see that their bounds have all been reverted.”
]
}, {
“cell_type”: “code”, “execution_count”: 35, “metadata”: {}, “outputs”: [
 {
 “data”: {
 “text/plain”: [
“[(1000.0, 1000.0),n”, ” (1000.0, 1000.0),n”, ” (1000.0, 1000.0),n”, ” (1000.0, 1000.0),n”, ” (1000.0, 1000.0)]”
]
}, “execution_count”: 35, “metadata”: {}, “output_type”: “execute_result”
}
], “source”: [
“[reaction.bounds for reaction in model.reactions[:5]]”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Nested contexts are also supported”
]
}, {
“cell_type”: “code”, “execution_count”: 36, “metadata”: {}, “outputs”: [
 {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“original objective: 1.0*Biomass_Ecoli_core  1.0*Biomass_Ecoli_core_reverse_2cdban”, “print objective in first context: 1.0*ATPM_reverse_5b752 + 1.0*ATPMn”, “print objective in second context: 1.0*ACALD  1.0*ACALD_reverse_fda2bn”, “objective after exiting second context: 1.0*ATPM_reverse_5b752 + 1.0*ATPMn”, “back to original objective: 1.0*Biomass_Ecoli_core  1.0*Biomass_Ecoli_core_reverse_2cdban”
]
}
], “source”: [
“print(‘original objective: ‘, model.objective.expression)n”, “with model:n”, ” model.objective = ‘ATPM’n”, ” print(‘print objective in first context:’, model.objective.expression)n”, ” with model:n”, ” model.objective = ‘ACALD’n”, ” print(‘print objective in second context:’, model.objective.expression)n”, ” print(‘objective after exiting second context:’,n”, ” model.objective.expression)n”, “print(‘back to original objective:’, model.objective.expression)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Most methods that modify the model are supported like this including adding and removing reactions and metabolites and setting the objective. Supported methods and functions mention this in the corresponding documentation.”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“While it does not have any actual effect, for syntactic convenience it is also possible to refer to the model by a different name than outside the context. Such as”
]
}, {
“cell_type”: “code”, “execution_count”: 37, “metadata”: {}, “outputs”: [], “source”: [
“with model as inner:n”, ” inner.reactions.PFK.knock_out”
]
}
], “metadata”: {
 “kernelspec”: {
“display_name”: “Python 3”, “language”: “python”, “name”: “python3”
}, “language_info”: {
 “codemirror_mode”: {
“name”: “ipython”, “version”: 3
}, “file_extension”: “.py”, “mimetype”: “text/xpython”, “name”: “python”, “nbconvert_exporter”: “python”, “pygments_lexer”: “ipython3”, “version”: “3.6.5”
}
}, “nbformat”: 4, “nbformat_minor”: 1
}
Global Configuration¶
With cobra > 0.13.4, we introduce a global configuration object. For now, you can configure default reaction bounds and optimization solver which will be respected by newly created reactions and models.
The configuration object¶
You can get a configuration object1 in the following way:
[1]:
import cobra
[2]:
cobra_config = cobra.Configuration()
1The configuration object is a singleton. That means only one instance can exist and it is respected everywhere in COBRApy.
Reaction bounds¶
The object has the following attributes which you can inspect but also change as desired.
[3]:
cobra_config.lower_bound
[3]:
1000.0
[4]:
cobra_config.upper_bound
[4]:
1000.0
[5]:
cobra_config.bounds
[5]:
(1000.0, 1000.0)
Changing bounds¶
If you modify the above values before creating a reaction they will be used.
[6]:
cobra_config.bounds = 10, 20
[7]:
cobra.Reaction("R1")
[7]:
Reaction identifier  R1 
Name  
Memory address  0x07f0426135fd0 
Stoichiometry 
> > 
GPR  
Lower bound  0.0 
Upper bound  20 
Please note that by default reactions are irreversible. You can change this behavior by unsetting the lower bound argument.
[8]:
cobra.Reaction("R2", lower_bound=None)
[8]:
Reaction identifier  R2 
Name  
Memory address  0x07f04260d4438 
Stoichiometry 
<=> <=> 
GPR  
Lower bound  10 
Upper bound  20 
N.B.: Most models define reaction bounds explicitly which takes precedence over the configured values.
[9]:
from cobra.test import create_test_model
[10]:
model = create_test_model("textbook")
[11]:
model.reactions.ACt2r
[11]:
Reaction identifier  ACt2r 
Name  R acetate reversible transport via proton  symport 
Memory address  0x07f042607c780 
Stoichiometry 
ac_e + h_e <=> ac_c + h_c Acetate + H+ <=> Acetate + H+ 
GPR  
Lower bound  1000.0 
Upper bound  1000.0 
Solver¶
You can define the default solver used by newly instantiated models. The default solver depends on your environment. In order we test for the availability of Gurobi, CPLEX, and GLPK. GLPK is assumed to always be present in the environment.
[12]:
model.solver
[12]:
<optlang.cplex_interface.Model at 0x7f04260d4b00>
Changing solver¶
[13]:
cobra_config.solver = "glpk_exact"
[14]:
new_model = create_test_model("textbook")
[15]:
new_model.solver
[15]:
<optlang.glpk_exact_interface.Model at 0x7f04260d47b8>
Changing global configuration values is mostly useful at the beginning of a work session.
Building a Model¶
Model, Reactions and Metabolites¶
This simple example demonstrates how to create a model, create a reaction, and then add the reaction to the model.
We’ll use the ‘3OAS140’ reaction from the STM_1.0 model:
1.0 malACP[c] + 1.0 h[c] + 1.0 ddcaACP[c] \(\rightarrow\) 1.0 co2[c] + 1.0 ACP[c] + 1.0 3omrsACP[c]
First, create the model and reaction.
[1]:
from cobra import Model, Reaction, Metabolite
[2]:
model = Model('example_model')
reaction = Reaction('R_3OAS140')
reaction.name = '3 oxoacyl acyl carrier protein synthase n C140 '
reaction.subsystem = 'Cell Envelope Biosynthesis'
reaction.lower_bound = 0. # This is the default
reaction.upper_bound = 1000. # This is the default
We need to create metabolites as well. If we were using an existing model, we could use Model.get_by_id
to get the appropriate Metabolite objects instead.
[3]:
ACP_c = Metabolite(
'ACP_c',
formula='C11H21N2O7PRS',
name='acylcarrierprotein',
compartment='c')
omrsACP_c = Metabolite(
'M3omrsACP_c',
formula='C25H45N2O9PRS',
name='3Oxotetradecanoylacylcarrierprotein',
compartment='c')
co2_c = Metabolite('co2_c', formula='CO2', name='CO2', compartment='c')
malACP_c = Metabolite(
'malACP_c',
formula='C14H22N2O10PRS',
name='Malonylacylcarrierprotein',
compartment='c')
h_c = Metabolite('h_c', formula='H', name='H', compartment='c')
ddcaACP_c = Metabolite(
'ddcaACP_c',
formula='C23H43N2O8PRS',
name='DodecanoylACPnC120ACP',
compartment='c')
Side note: SId
It is highly recommended that the ids for reactions, metabolites and genes are valid SBML identifiers (SId
). SId
is a data type derived from the basic XML typestring, but with restrictions about the characters permitted and the sequences in which those characters may appear.
letter ::= ’a’..’z’,’A’..’Z’
digit ::= ’0’..’9’
idChar ::= letter  digit  ’_’
SId ::= ( letter  ’_’ ) idChar*
The main limitation is that ids cannot start with numbers. Using SId
s allows serialization to SBML. In addition features such as code completion and object access via the dot syntax will work in cobrapy
.
Adding metabolites to a reaction uses a dictionary of the metabolites and their stoichiometric coefficients. A group of metabolites can be added all at once, or they can be added one at a time.
[4]:
reaction.add_metabolites({
malACP_c: 1.0,
h_c: 1.0,
ddcaACP_c: 1.0,
co2_c: 1.0,
ACP_c: 1.0,
omrsACP_c: 1.0
})
reaction.reaction # This gives a string representation of the reaction
[4]:
'ddcaACP_c + h_c + malACP_c > ACP_c + M3omrsACP_c + co2_c'
The gene_reaction_rule is a boolean representation of the gene requirements for this reaction to be active as described in Schellenberger et al 2011 Nature Protocols 6(9):1290307. We will assign the gene reaction rule string, which will automatically create the corresponding gene objects.
[5]:
reaction.gene_reaction_rule = '( STM2378 or STM1197 )'
reaction.genes
[5]:
frozenset({<Gene STM1197 at 0x7fef047474c0>, <Gene STM2378 at 0x7fef04747dc0>})
At this point in time, the model is still empty
[6]:
print(f'{len(model.reactions)} reactions initially')
print(f'{len(model.metabolites)} metabolites initially')
print(f'{len(model.genes)} genes initially')
0 reactions initially
0 metabolites initially
0 genes initially
We will add the reaction to the model, which will also add all associated metabolites and genes
[7]:
model.add_reactions([reaction])
# The objects have been added to the model
print(f'{len(model.reactions)} reactions')
print(f'{len(model.metabolites)} metabolites')
print(f'{len(model.genes)} genes')
1 reactions
6 metabolites
2 genes
We can iterate through the model objects to observe the contents
[8]:
# Iterate through the the objects in the model
print("Reactions")
print("")
for x in model.reactions:
print("%s : %s" % (x.id, x.reaction))
print("")
print("Metabolites")
print("")
for x in model.metabolites:
print('%9s : %s' % (x.id, x.formula))
print("")
print("Genes")
print("")
for x in model.genes:
associated_ids = (i.id for i in x.reactions)
print("%s is associated with reactions: %s" %
(x.id, "{" + ", ".join(associated_ids) + "}"))
Reactions

R_3OAS140 : ddcaACP_c + h_c + malACP_c > ACP_c + M3omrsACP_c + co2_c
Metabolites

malACP_c : C14H22N2O10PRS
h_c : H
ddcaACP_c : C23H43N2O8PRS
co2_c : CO2
ACP_c : C11H21N2O7PRS
M3omrsACP_c : C25H45N2O9PRS
Genes

STM2378 is associated with reactions: {R_3OAS140}
STM1197 is associated with reactions: {R_3OAS140}
Objective¶
Last we need to set the objective of the model. Here, we just want this to be the maximization of the flux in the single reaction we added and we do this by assigning the reaction’s identifier to the objective
property of the model.
[9]:
model.objective = 'R_3OAS140'
The created objective is a symbolic algebraic expression and we can examine it by printing it
[10]:
print(model.objective.expression)
print(model.objective.direction)
1.0*R_3OAS140  1.0*R_3OAS140_reverse_60acb
max
which here shows that the solver will maximize the flux in the forward direction.
Model Validation¶
For exchange with other tools you can validate and export the model to SBML. For more information on serialization and available formats see the section “Reading and Writing Models”
[11]:
import tempfile
from pprint import pprint
from cobra.io import write_sbml_model, validate_sbml_model
with tempfile.NamedTemporaryFile(suffix='.xml') as f_sbml:
write_sbml_model(model, filename=f_sbml.name)
report = validate_sbml_model(filename=f_sbml.name)
pprint(report)
(<Model example_model at 0x7fef046eb160>,
{'COBRA_CHECK': [],
'COBRA_ERROR': [],
'COBRA_FATAL': [],
'COBRA_WARNING': [],
'SBML_ERROR': [],
'SBML_FATAL': [],
'SBML_SCHEMA_ERROR': [],
'SBML_WARNING': []})
The model is valid with no COBRA or SBML errors or warnings.
Exchanges, Sinks and Demands¶
Boundary reactions can be added using the model’s method add_boundary
. There are three different types of predefined boundary reactions: exchange, demand, and sink reactions. All of them are unbalanced pseudo reactions, that means they fulfill a function for modeling by adding to or removing metabolites from the model system but are not based on real biology. An exchange reaction is a reversible reaction that adds to or removes an extracellular metabolite from the extracellular compartment.
A demand reaction is an irreversible reaction that consumes an intracellular metabolite. A sink is similar to an exchange but specifically for intracellular metabolites, i.e., a reversible reaction that adds or removes an intracellular metabolite.
[12]:
print("exchanges", model.exchanges)
print("demands", model.demands)
print("sinks", model.sinks)
There are no boundary reactions in this model. Therefore specific types of boundary reactions such as 'exchanges', 'demands' or 'sinks' cannot be identified.
There are no boundary reactions in this model. Therefore specific types of boundary reactions such as 'exchanges', 'demands' or 'sinks' cannot be identified.
There are no boundary reactions in this model. Therefore specific types of boundary reactions such as 'exchanges', 'demands' or 'sinks' cannot be identified.
exchanges []
demands []
sinks []
Boundary reactions are defined on metabolites. First we add two metabolites to the model then we define the boundary reactions. We add glycogen to the cytosolic compartment c
and CO2 to the external compartment e
.
[13]:
model.add_metabolites([
Metabolite(
'glycogen_c',
name='glycogen',
compartment='c'
),
Metabolite(
'co2_e',
name='CO2',
compartment='e'
),
])
[14]:
# create exchange reaction
model.add_boundary(model.metabolites.get_by_id("co2_e"), type="exchange")
[14]:
Reaction identifier  EX_co2_e 
Name  CO2 exchange 
Memory address  0x07fef04703d90 
Stoichiometry 
co2_e <=> CO2 <=> 
GPR  
Lower bound  1000.0 
Upper bound  1000.0 
[15]:
# create exchange reaction
model.add_boundary(model.metabolites.get_by_id("glycogen_c"), type="sink")
[15]:
Reaction identifier  SK_glycogen_c 
Name  glycogen sink 
Memory address  0x07fef046eb2e0 
Stoichiometry 
glycogen_c <=> glycogen <=> 
GPR  
Lower bound  1000.0 
Upper bound  1000.0 
[16]:
# Now we have an additional exchange and sink reaction in the model
print("exchanges", model.exchanges)
print("sinks", model.sinks)
print("demands", model.demands)
exchanges [<Reaction EX_co2_e at 0x7fef04703d90>]
sinks [<Reaction SK_glycogen_c at 0x7fef046eb2e0>]
demands []
To create a demand reaction instead of a sink use type demand
instead of sink
.
Information on all boundary reactions is available via the model’s property boundary
.
[17]:
# boundary reactions
model.boundary
[17]:
[<Reaction EX_co2_e at 0x7fef04703d90>,
<Reaction SK_glycogen_c at 0x7fef046eb2e0>]
A neat trick to get all metabolic reactions is
[18]:
# metabolic reactions
set(model.reactions)  set(model.boundary)
[18]:
{<Reaction R_3OAS140 at 0x7fef04737d00>}
Reading and Writing Models¶
Cobrapy supports reading and writing models in SBML (with and without FBC), JSON, YAML, MAT, and pickle formats. Generally, SBML with FBC version 2 is the preferred format for general use. The JSON format may be more useful for cobrapyspecific functionality.
The package also ships with test models in various formats for testing purposes.
[1]:
import cobra.test
import os
from os.path import join
data_dir = cobra.test.data_dir
print("mini test files: ")
print(", ".join(i for i in os.listdir(data_dir) if i.startswith("mini")))
textbook_model = cobra.test.create_test_model("textbook")
ecoli_model = cobra.test.create_test_model("ecoli")
salmonella_model = cobra.test.create_test_model("salmonella")
mini test files:
mini.json, mini.mat, mini.pickle, mini.yml, mini_cobra.xml, mini_fbc1.xml, mini_fbc2.xml, mini_fbc2.xml.bz2, mini_fbc2.xml.gz
SBML¶
The Systems Biology Markup Language is an XMLbased standard format for distributing models which has support for COBRA models through the FBC extension version 2.
Cobrapy has native support for reading and writing SBML with FBCv2. Please note that all id’s in the model must conform to the SBML SID requirements in order to generate a valid SBML file.
[2]:
cobra.io.read_sbml_model(join(data_dir, "mini_fbc2.xml"))
[2]:
Name  mini_textbook 
Memory address  0x01074fd080 
Number of metabolites  23 
Number of reactions  18 
Objective expression  1.0*ATPM_reverse_5b752  1.0*PFK_reverse_d24a6 + 1.0*PFK + 1.0*ATPM 
Compartments  cytosol, extracellular 
[3]:
cobra.io.write_sbml_model(textbook_model, "test_fbc2.xml")
There are other dialects of SBML prior to FBC 2 which have previously been use to encode COBRA models. The primary ones is the “COBRA” dialect which used the “notes” fields in SBML files.
Cobrapy can use libsbml, which must be installed separately (see installation instructions) to read and write these files. When reading in a model, it will automatically detect whether FBC was used or not. When writing a model, the use_fbc_package flag can be used can be used to write files in this legacy “cobra” format.
Consider having the lxml package installed as it can speed up parsing considerably.
[4]:
cobra.io.read_sbml_model(join(data_dir, "mini_cobra.xml"))
[4]:
Name  mini_textbook 
Memory address  0x0112fa6b38 
Number of metabolites  23 
Number of reactions  18 
Objective expression  1.0*ATPM_reverse_5b752  1.0*PFK_reverse_d24a6 + 1.0*PFK + 1.0*ATPM 
Compartments  cytosol, extracellular 
[5]:
cobra.io.write_sbml_model(
textbook_model, "test_cobra.xml", use_fbc_package=False)
JSON¶
Cobrapy models have a JSON (JavaScript Object Notation) representation. This format was created for interoperability with escher.
[6]:
cobra.io.load_json_model(join(data_dir, "mini.json"))
[6]:
Name  mini_textbook 
Memory address  0x0113061080 
Number of metabolites  23 
Number of reactions  18 
Objective expression  1.0*ATPM_reverse_5b752  1.0*PFK_reverse_d24a6 + 1.0*PFK + 1.0*ATPM 
Compartments  cytosol, extracellular 
[7]:
cobra.io.save_json_model(textbook_model, "test.json")
YAML¶
Cobrapy models have a YAML (YAML Ain’t Markup Language) representation. This format was created for more human readable model representations and automatic diffs between models.
[8]:
cobra.io.load_yaml_model(join(data_dir, "mini.yml"))
[8]:
Name  mini_textbook 
Memory address  0x0113013390 
Number of metabolites  23 
Number of reactions  18 
Objective expression  1.0*ATPM_reverse_5b752  1.0*PFK_reverse_d24a6 + 1.0*PFK + 1.0*ATPM 
Compartments  extracellular, cytosol 
[9]:
cobra.io.save_yaml_model(textbook_model, "test.yml")
MATLAB¶
Often, models may be imported and exported solely for the purposes of working with the same models in cobrapy and the MATLAB cobra toolbox. MATLAB has its own “.mat” format for storing variables. Reading and writing to these mat files from python requires scipy.
A mat file can contain multiple MATLAB variables. Therefore, the variable name of the model in the MATLAB file can be passed into the reading function:
[10]:
cobra.io.load_matlab_model(
join(data_dir, "mini.mat"), variable_name="mini_textbook")
[10]:
Name  mini_textbook 
Memory address  0x0113000b70 
Number of metabolites  23 
Number of reactions  18 
Objective expression  1.0*ATPM_reverse_5b752  1.0*PFK_reverse_d24a6 + 1.0*PFK + 1.0*ATPM 
Compartments  c, e 
If the mat file contains only a single model, cobra can figure out which variable to read from, and the variable_name parameter is unnecessary.
[11]:
cobra.io.load_matlab_model(join(data_dir, "mini.mat"))
[11]:
Name  mini_textbook 
Memory address  0x0113758438 
Number of metabolites  23 
Number of reactions  18 
Objective expression  1.0*ATPM_reverse_5b752  1.0*PFK_reverse_d24a6 + 1.0*PFK + 1.0*ATPM 
Compartments  c, e 
Saving models to mat files is also relatively straightforward
[12]:
cobra.io.save_matlab_model(textbook_model, "test.mat")
Simulating with FBA¶
Simulations using flux balance analysis can be solved using Model.optimize()
. This will maximize or minimize (maximizing is the default) flux through the objective reactions.
[1]:
import cobra.test
model = cobra.test.create_test_model("textbook")
Running FBA¶
[2]:
solution = model.optimize()
print(solution)
<Solution 0.874 at 0x112eb3d30>
The Model.optimize() function will return a Solution object. A solution object has several attributes:
objective_value
: the objective valuestatus
: the status from the linear programming solverfluxes
: a pandas series with flux indexed by reaction identifier. The flux for a reaction variable is the difference of the primal values for the forward and reverse reaction variables.shadow_prices
: a pandas series with shadow price indexed by the metabolite identifier.
For example, after the last call to model.optimize()
, if the optimization succeeds it’s status will be optimal. In case the model is infeasible an error is raised.
[3]:
solution.objective_value
[3]:
0.8739215069684307
The solvers that can be used with cobrapy are so fast that for many small to midsize models computing the solution can be even faster than it takes to collect the values from the solver and convert to them python objects. With model.optimize
, we gather values for all reactions and metabolites and that can take a significant amount of time if done repeatedly. If we are only interested in the flux value of a single reaction or the objective, it is faster to instead use model.slim_optimize
which only does the optimization and returns the objective value leaving it up to you to fetch other values that you may need.
[4]:
%%time
model.optimize().objective_value
CPU times: user 3.84 ms, sys: 672 µs, total: 4.51 ms
Wall time: 6.16 ms
[4]:
0.8739215069684307
[5]:
%%time
model.slim_optimize()
CPU times: user 229 µs, sys: 19 µs, total: 248 µs
Wall time: 257 µs
[5]:
0.8739215069684307
Analyzing FBA solutions¶
Models solved using FBA can be further analyzed by using summary methods, which output printed text to give a quick representation of model behavior. Calling the summary method on the entire model displays information on the input and output behavior of the model, along with the optimized objective.
[6]:
model.summary()
IN FLUXES OUT FLUXES OBJECTIVES
  
o2_e 21.8 h2o_e 29.2 Biomass_Ecol... 0.874
glc__D_e 10 co2_e 22.8
nh4_e 4.77 h_e 17.5
pi_e 3.21
In addition, the inputoutput behavior of individual metabolites can also be inspected using summary methods. For instance, the following commands can be used to examine the overall redox balance of the model
[7]:
model.metabolites.nadh_c.summary()
PRODUCING REACTIONS  Nicotinamide adenine dinucleotide  reduced (nadh_c)

% FLUX RXN ID REACTION
   
42% 16 GAPD g3p_c + nad_c + pi_c <=> 13dpg_c + h_c + nadh_c
24% 9.28 PDH coa_c + nad_c + pyr_c > accoa_c + co2_c + nadh_c
13% 5.06 AKGDH akg_c + coa_c + nad_c > co2_c + nadh_c + succ...
13% 5.06 MDH mal__L_c + nad_c <=> h_c + nadh_c + oaa_c
8% 3.1 Biomass... 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0....
CONSUMING REACTIONS  Nicotinamide adenine dinucleotide  reduced (nadh_c)

% FLUX RXN ID REACTION
   
100% 38.5 NADH16 4.0 h_c + nadh_c + q8_c > 3.0 h_e + nad_c + q...
Or to get a sense of the main energy production and consumption reactions
[8]:
model.metabolites.atp_c.summary()
PRODUCING REACTIONS  ATP (atp_c)

% FLUX RXN ID REACTION
   
67% 45.5 ATPS4r adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c
23% 16 PGK 3pg_c + atp_c <=> 13dpg_c + adp_c
7% 5.06 SUCOAS atp_c + coa_c + succ_c <=> adp_c + pi_c + succoa_c
3% 1.76 PYK adp_c + h_c + pep_c > atp_c + pyr_c
CONSUMING REACTIONS  ATP (atp_c)

% FLUX RXN ID REACTION
   
76% 52.3 Biomass... 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0....
12% 8.39 ATPM atp_c + h2o_c > adp_c + h_c + pi_c
11% 7.48 PFK atp_c + f6p_c > adp_c + fdp_c + h_c
0% 0.223 GLNS atp_c + glu__L_c + nh4_c > adp_c + gln__L_c +...
Changing the Objectives¶
The objective function is determined from the objective_coefficient attribute of the objective reaction(s). Generally, a “biomass” function which describes the composition of metabolites which make up a cell is used.
[9]:
biomass_rxn = model.reactions.get_by_id("Biomass_Ecoli_core")
Currently in the model, there is only one reaction in the objective (the biomass reaction), with an linear coefficient of 1.
[10]:
from cobra.util.solver import linear_reaction_coefficients
linear_reaction_coefficients(model)
[10]:
{<Reaction Biomass_Ecoli_core at 0x112eab4a8>: 1.0}
The objective function can be changed by assigning Model.objective, which can be a reaction object (or just it’s name), or a dict
of {Reaction: objective_coefficient}
.
[11]:
# change the objective to ATPM
model.objective = "ATPM"
# The upper bound should be 1000, so that we get
# the actual optimal value
model.reactions.get_by_id("ATPM").upper_bound = 1000.
linear_reaction_coefficients(model)
[11]:
{<Reaction ATPM at 0x112eab470>: 1.0}
[12]:
model.optimize().objective_value
[12]:
174.99999999999966
We can also have more complicated objectives including quadratic terms.
Running FVA¶
FBA will not give always give unique solution, because multiple flux states can achieve the same optimum. FVA (or flux variability analysis) finds the ranges of each metabolic flux at the optimum.
[13]:
from cobra.flux_analysis import flux_variability_analysis
[14]:
flux_variability_analysis(model, model.reactions[:10])
[14]:
maximum  minimum  

ACALD  2.208811e30  5.247085e14 
ACALDt  0.000000e+00  5.247085e14 
ACKr  0.000000e+00  8.024953e14 
ACONTa  2.000000e+01  2.000000e+01 
ACONTb  2.000000e+01  2.000000e+01 
ACt2r  0.000000e+00  8.024953e14 
ADK1  3.410605e13  0.000000e+00 
AKGDH  2.000000e+01  2.000000e+01 
AKGt2r  0.000000e+00  2.902643e14 
ALCD2x  0.000000e+00  4.547474e14 
Setting parameter fraction_of_optimium=0.90
would give the flux ranges for reactions at 90% optimality.
[15]:
cobra.flux_analysis.flux_variability_analysis(
model, model.reactions[:10], fraction_of_optimum=0.9)
[15]:
maximum  minimum  

ACALD  0.000000e+00  2.692308 
ACALDt  0.000000e+00  2.692308 
ACKr  6.635712e30  4.117647 
ACONTa  2.000000e+01  8.461538 
ACONTb  2.000000e+01  8.461538 
ACt2r  0.000000e+00  4.117647 
ADK1  1.750000e+01  0.000000 
AKGDH  2.000000e+01  2.500000 
AKGt2r  2.651196e16  1.489362 
ALCD2x  0.000000e+00  2.333333 
The standard FVA may contain loops, i.e. high absolute flux values that only can be high if they are allowed to participate in loops (a mathematical artifact that cannot happen in vivo). Use the loopless
argument to avoid such loops. Below, we can see that FRD7 and SUCDi reactions can participate in loops but that this is avoided when using the looplesss FVA.
[16]:
loop_reactions = [model.reactions.FRD7, model.reactions.SUCDi]
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=False)
[16]:
maximum  minimum  

FRD7  980.0  0.0 
SUCDi  1000.0  20.0 
[17]:
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=True)
[17]:
maximum  minimum  

FRD7  0.0  0.0 
SUCDi  20.0  20.0 
Running FVA in summary methods¶
Flux variability analysis can also be embedded in calls to summary methods. For instance, the expected variability in substrate consumption and product formation can be quickly found by
[18]:
model.optimize()
model.summary(fva=0.95)
IN FLUXES OUT FLUXES OBJECTIVES
  
id Flux Range id Flux Range ATPM 175
     
o2_e 60 [55.9, 60] co2_e 60 [54.2, 60]
glc__D_e 10 [9.5, 10] h2o_e 60 [54.2, 60]
nh4_e 0 [0, 0.673] for_e 0 [0, 5.83]
pi_e 0 [0, 0.171] h_e 0 [0, 5.83]
ac_e 0 [0, 2.06]
acald_e 0 [0, 1.35]
pyr_e 0 [0, 1.35]
etoh_e 0 [0, 1.17]
lac__D_e 0 [0, 1.13]
succ_e 0 [0, 0.875]
akg_e 0 [0, 0.745]
glu__L_e 0 [0, 0.673]
Similarly, variability in metabolite mass balances can also be checked with flux variability analysis.
[19]:
model.metabolites.pyr_c.summary(fva=0.95)
PRODUCING REACTIONS  Pyruvate (pyr_c)

% FLUX RANGE RXN ID REACTION
    
50% 10 [1.25, 18.8] PYK adp_c + h_c + pep_c > atp_c + pyr_c
50% 10 [9.5, 10] GLCpts glc__D_e + pep_c > g6p_c + pyr_c
0% 0 [0, 8.75] ME1 mal__L_c + nad_c > co2_c + nadh_c +...
0% 0 [0, 8.75] ME2 mal__L_c + nadp_c > co2_c + nadph_c...
CONSUMING REACTIONS  Pyruvate (pyr_c)

% FLUX RANGE RXN ID REACTION
    
100% 20 [13, 28.8] PDH coa_c + nad_c + pyr_c > accoa_c + c...
0% 0 [0, 8.75] PPS atp_c + h2o_c + pyr_c > amp_c + 2.0...
0% 0 [0, 5.83] PFL coa_c + pyr_c > accoa_c + for_c
0% 0 [0, 1.35] PYRt2 h_e + pyr_e <=> h_c + pyr_c
0% 0 [0, 1.13] LDH_D lac__D_c + nad_c <=> h_c + nadh_c + p...
0% 0 [0, 0.132] Biomass... 1.496 3pg_c + 3.7478 accoa_c + 59.81 ...
In these summary methods, the values are reported as a the center point +/ the range of the FVA solution, calculated from the maximum and minimum values.
Running pFBA¶
Parsimonious FBA (often written pFBA) finds a flux distribution which gives the optimal growth rate, but minimizes the total sum of flux. This involves solving two sequential linear programs, but is handled transparently by cobrapy. For more details on pFBA, please see Lewis et al. (2010).
[20]:
model.objective = 'Biomass_Ecoli_core'
fba_solution = model.optimize()
pfba_solution = cobra.flux_analysis.pfba(model)
These functions should give approximately the same objective value.
[21]:
abs(fba_solution.fluxes["Biomass_Ecoli_core"]  pfba_solution.fluxes[
"Biomass_Ecoli_core"])
[21]:
7.7715611723760958e16
Running geometric FBA¶
Geometric FBA finds a unique optimal flux distribution which is central to the range of possible fluxes. For more details on geometric FBA, please see K Smallbone, E Simeonidis (2009).
[22]:
geometric_fba_sol = cobra.flux_analysis.geometric_fba(model)
geometric_fba_sol
[22]:
fluxes  reduced_costs  

ACALD  1.785214e14  0.0 
ACALDt  1.785214e14  0.0 
ACKr  0.000000e+00  0.0 
ACONTa  6.007250e+00  0.0 
ACONTb  6.007250e+00  0.0 
...  ...  ... 
TALA  1.496984e+00  0.0 
THD2  1.522652e14  0.0 
TKT1  1.496984e+00  0.0 
TKT2  1.181498e+00  0.0 
TPI  7.477382e+00  0.0 
95 rows × 2 columns
Simulating Deletions¶
[2]:
import pandas
from time import time
import cobra.test
from cobra.flux_analysis import (
single_gene_deletion, single_reaction_deletion, double_gene_deletion,
double_reaction_deletion)
cobra_model = cobra.test.create_test_model("textbook")
ecoli_model = cobra.test.create_test_model("ecoli")
Knocking out single genes and reactions¶
A commonly asked question when analyzing metabolic models is what will happen if a certain reaction was not allowed to have any flux at all. This can tested using cobrapy by
[3]:
print('complete model: ', cobra_model.optimize())
with cobra_model:
cobra_model.reactions.PFK.knock_out()
print('pfk knocked out: ', cobra_model.optimize())
complete model: <Solution 0.874 at 0x7f8130c3c760>
pfk knocked out: <Solution 0.704 at 0x7f8130c3c070>
For evaluating genetic manipulation strategies, it is more interesting to examine what happens if given genes are knocked out as doing so can affect no reactions in case of redundancy, or more reactions if gene when is participating in more than one reaction.
[4]:
print('complete model: ', cobra_model.optimize())
with cobra_model:
cobra_model.genes.b1723.knock_out()
print('pfkA knocked out: ', cobra_model.optimize())
cobra_model.genes.b3916.knock_out()
print('pfkB knocked out: ', cobra_model.optimize())
complete model: <Solution 0.874 at 0x7f8131dbc7c0>
pfkA knocked out: <Solution 0.874 at 0x7f816ddee820>
pfkB knocked out: <Solution 0.704 at 0x7f8130c3cf40>
Single Deletions¶
Perform all single gene deletions on a model
[5]:
deletion_results = single_gene_deletion(cobra_model)
These can also be done for only a subset of genes
[6]:
single_gene_deletion(cobra_model, cobra_model.genes[:20])
[6]:
ids  growth  status  

0  {b0351}  0.873922  optimal 
1  {b0727}  0.858307  optimal 
2  {b2587}  0.873922  optimal 
3  {b0474}  0.873922  optimal 
4  {b3734}  0.374230  optimal 
5  {b0726}  0.858307  optimal 
6  {b3735}  0.374230  optimal 
7  {b1241}  0.873922  optimal 
8  {b1276}  0.873922  optimal 
9  {b1849}  0.873922  optimal 
10  {b3732}  0.374230  optimal 
11  {s0001}  0.211141  optimal 
12  {b0118}  0.873922  optimal 
13  {b0356}  0.873922  optimal 
14  {b2296}  0.873922  optimal 
15  {b0116}  0.782351  optimal 
16  {b1478}  0.873922  optimal 
17  {b3115}  0.873922  optimal 
18  {b3736}  0.374230  optimal 
19  {b3733}  0.374230  optimal 
This can also be done for reactions
[7]:
single_reaction_deletion(cobra_model, cobra_model.reactions[:20])
[7]:
ids  growth  status  

0  {ACALDt}  0.873922  optimal 
1  {ATPS4r}  0.374230  optimal 
2  {EX_ac_e}  0.873922  optimal 
3  {AKGDH}  0.858307  optimal 
4  {ACKr}  0.873922  optimal 
5  {ETOHt2r}  0.873922  optimal 
6  {CS}  0.000000  optimal 
7  {ALCD2x}  0.873922  optimal 
8  {ACONTa}  0.000000  optimal 
9  {D_LACt2}  0.873922  optimal 
10  {ADK1}  0.873922  optimal 
11  {Biomass_Ecoli_core}  0.000000  optimal 
12  {ATPM}  0.916647  optimal 
13  {CYTBD}  0.211663  optimal 
14  {ACt2r}  0.873922  optimal 
15  {CO2t}  0.461670  optimal 
16  {ACALD}  0.873922  optimal 
17  {ENO}  0.000000  optimal 
18  {AKGt2r}  0.873922  optimal 
19  {ACONTb}  0.000000  optimal 
Double Deletions¶
Double deletions run in a similar way.
[8]:
double_gene_deletion(
cobra_model, cobra_model.genes[5:]).round(4)
[8]:
ids  growth  status  

0  {b2935}  0.8739  optimal 
1  {b2465, b2464}  0.8739  optimal 
2  {b0008, b2464}  0.8739  optimal 
3  {b0008, b3919}  0.7040  optimal 
4  {b3919, b2465}  0.7040  optimal 
5  {b0008, b2935}  0.8739  optimal 
6  {b3919}  0.7040  optimal 
7  {b2465, b2935}  0.8739  optimal 
8  {b0008, b2465}  0.8739  optimal 
9  {b3919, b2464}  0.7040  optimal 
10  {b2464}  0.8739  optimal 
11  {b0008}  0.8739  optimal 
12  {b3919, b2935}  0.7040  optimal 
13  {b2935, b2464}  0.8739  optimal 
14  {b2465}  0.8739  optimal 
By default, the double deletion function will automatically use multiprocessing, splitting the task over up to 4 cores if they are available. The number of cores can be manually specified as well. Setting use of a single core will disable use of the multiprocessing library, which often aids debugging.
[9]:
start = time() # start timer()
double_gene_deletion(
ecoli_model, ecoli_model.genes[:25], processes=2)
t1 = time()  start
print("Double gene deletions for 200 genes completed in "
"%.2f sec with 2 cores" % t1)
start = time() # start timer()
double_gene_deletion(
ecoli_model, ecoli_model.genes[:25], processes=1)
t2 = time()  start
print("Double gene deletions for 200 genes completed in "
"%.2f sec with 1 core" % t2)
print("Speedup of %.2fx" % (t2 / t1))
Double gene deletions for 200 genes completed in 9.02 sec with 2 cores
Double gene deletions for 200 genes completed in 15.48 sec with 1 core
Speedup of 1.72x
Double deletions can also be run for reactions.
[10]:
double_reaction_deletion(
cobra_model, cobra_model.reactions[2:7]).round(4)
[10]:
ids  growth  status  

0  {ACt2r, ACKr}  0.8739  optimal 
1  {ACONTa, ACKr}  0.0000  optimal 
2  {ACONTa, ACONTb}  0.0000  optimal 
3  {ACONTa}  0.0000  optimal 
4  {ACONTb, ACt2r}  0.0000  optimal 
5  {ACt2r}  0.8739  optimal 
6  {ACONTa, ADK1}  0.0000  optimal 
7  {ACONTa, ACt2r}  0.0000  optimal 
8  {ADK1, ACt2r}  0.8739  optimal 
9  {ACONTb, ACKr}  0.0000  optimal 
10  {ADK1}  0.8739  optimal 
11  {ACONTb, ADK1}  0.0000  optimal 
12  {ADK1, ACKr}  0.8739  optimal 
13  {ACKr}  0.8739  optimal 
14  {ACONTb}  0.0000  optimal 
Accessing individual deletion results¶
Note that the indices for deletions are python set objects. This is the appropriate type since the order of deletions does not matter. Deleting reaction 1 and reaction 2 will have the same effect as deleting reaction 2 and reaction 1.
To make it easier to access results all DataFrames returned by COBRAPpy deletion functions have a knockout
indexer that makes that a bit simpler. Each entry in the indexer is treated as a single deletion entry. So you need to pass sets for double deletions.
[11]:
single = single_reaction_deletion(cobra_model)
double = double_reaction_deletion(cobra_model)
print(single.knockout["ATPM"])
print(double.knockout[{"ATPM", "TKT1"}])
ids growth status
89 {ATPM} 0.916647 optimal
ids growth status
2238 {ATPM, TKT1} 0.90584 optimal
This can be used to get several deletions at once and will also work for Reaction or Gene objects (depending on what you deleted) directly.
[12]:
atpm = cobra_model.reactions.ATPM
tkt1 = cobra_model.reactions.TKT1
pfk = cobra_model.reactions.PFK
print(single.knockout[atpm, tkt1, pfk])
print(double.knockout[{atpm, tkt1}, {atpm, pfk}, {atpm}])
ids growth status
15 {PFK} 0.704037 optimal
17 {TKT1} 0.864759 optimal
89 {ATPM} 0.916647 optimal
ids growth status
762 {ATPM} 0.916647 optimal
2238 {ATPM, TKT1} 0.905840 optimal
2533 {PFK, ATPM} 0.704037 optimal
[ ]:
Production envelopes¶
Production envelopes (aka phenotype phase planes) will show distinct phases of optimal growth with different use of two different substrates. For more information, see Edwards et al.
Cobrapy supports calculating these production envelopes and they can easily be plotted using your favorite plotting package. Here, we will make one for the “textbook” E. coli core model and demonstrate plotting using matplotlib.
[1]:
import cobra.test
from cobra.flux_analysis import production_envelope
model = cobra.test.create_test_model("textbook")
We want to make a phenotype phase plane to evaluate uptakes of Glucose and Oxygen.
[2]:
prod_env = production_envelope(model, ["EX_glc__D_e", "EX_o2_e"])
[3]:
prod_env.head()
[3]:
carbon_source  carbon_yield_maximum  carbon_yield_minimum  flux_maximum  flux_minimum  mass_yield_maximum  mass_yield_minimum  EX_glc__D_e  EX_o2_e  

0  EX_glc__D_e  1.442300e13  0.0  0.000000  0.0  NaN  NaN  10.0  60.000000 
1  EX_glc__D_e  1.310050e+00  0.0  0.072244  0.0  NaN  NaN  10.0  56.842105 
2  EX_glc__D_e  2.620100e+00  0.0  0.144488  0.0  NaN  NaN  10.0  53.684211 
3  EX_glc__D_e  3.930150e+00  0.0  0.216732  0.0  NaN  NaN  10.0  50.526316 
4  EX_glc__D_e  5.240200e+00  0.0  0.288975  0.0  NaN  NaN  10.0  47.368421 
If we specify the carbon source, we can also get the carbon and mass yield. For example, temporarily setting the objective to produce acetate instead we could get production envelope as follows and pandas to quickly plot the results.
[4]:
prod_env = production_envelope(
model, ["EX_o2_e"], objective="EX_ac_e", carbon_sources="EX_glc__D_e")
[5]:
prod_env.head()
[5]:
carbon_source  carbon_yield_maximum  carbon_yield_minimum  flux_maximum  flux_minimum  mass_yield_maximum  mass_yield_minimum  EX_o2_e  

0  EX_glc__D_e  2.385536e15  0.0  0.000000  0.0  2.345496e15  0.0  60.000000 
1  EX_glc__D_e  5.263158e02  0.0  1.578947  0.0  5.174819e02  0.0  56.842105 
2  EX_glc__D_e  1.052632e01  0.0  3.157895  0.0  1.034964e01  0.0  53.684211 
3  EX_glc__D_e  1.578947e01  0.0  4.736842  0.0  1.552446e01  0.0  50.526316 
4  EX_glc__D_e  2.105263e01  0.0  6.315789  0.0  2.069927e01  0.0  47.368421 
[6]:
%matplotlib inline
[7]:
prod_env.plot(
kind='line', x='EX_o2_e', y='carbon_yield_maximum');
Previous versions of cobrapy included more tailored plots for phase planes which have now been dropped in order to improve maintainability and enhance the focus of cobrapy. Plotting for cobra models is intended for another package.
Flux sampling¶
Basic usage¶
The easiest way to get started with flux sampling is using the sample
function in the flux_analysis
submodule. sample
takes at least two arguments: a cobra model and the number of samples you want to generate.
[1]:
from cobra.test import create_test_model
from cobra.sampling import sample
model = create_test_model("textbook")
s = sample(model, 100)
s.head()
[1]:
ACALD  ACALDt  ACKr  ACONTa  ACONTb  ACt2r  ADK1  AKGDH  AKGt2r  ALCD2x  ...  RPI  SUCCt2_2  SUCCt3  SUCDi  SUCOAS  TALA  THD2  TKT1  TKT2  TPI  

0  2.060626  0.766231  1.746726  6.136642  6.136642  1.746726  13.915541  2.174506  0.242290  1.294395  ...  6.117270  33.457990  34.319917  704.483302  2.174506  6.109618  0.230408  6.109618  6.106540  3.122076 
1  1.518217  1.265778  0.253608  9.081331  9.081331  0.253608  7.194475  5.979050  0.225992  0.252439  ...  5.072733  39.902893  40.343192  718.488475  5.979050  4.991843  0.137019  4.991843  4.959315  4.172389 
2  3.790368  1.292543  0.457502  9.340755  9.340755  0.457502  23.435794  1.652395  0.333891  2.497825  ...  0.674220  0.153276  1.506968  844.889698  1.652395  0.673601  9.198001  0.673601  0.673352  7.770955 
3  5.173189  4.511308  2.333962  7.364836  7.364836  2.333962  11.725401  2.504044  0.051420  0.661881  ...  0.681200  7.506732  9.110446  885.755585  2.504044  0.656561  7.514520  0.656561  0.646653  8.450394 
4  6.787036  5.645414  1.521566  6.373250  6.373250  1.521566  4.823373  3.452123  0.126943  1.141621  ...  0.510598  9.307459  10.941500  749.854462  3.452123  0.474878  6.235982  0.474878  0.460514  8.908012 
5 rows × 95 columns
By default sample uses the optgp
method based on the method presented here as it is suited for larger models and can run in parallel. By default the sampler uses a single process. This can be changed by using the processes
argument.
[2]:
print("One process:")
%time s = sample(model, 1000)
print("Two processes:")
%time s = sample(model, 1000, processes=2)
One process:
CPU times: user 19.7 s, sys: 918 ms, total: 20.6 s
Wall time: 16.1 s
Two processes:
CPU times: user 1.31 s, sys: 154 ms, total: 1.46 s
Wall time: 8.76 s
Alternatively you can also user Artificial Centering HitandRun for sampling by setting the method to achr
. achr
does not support parallel execution but has good convergence and is almost Markovian.
[3]:
s = sample(model, 100, method="achr")
In general setting up the sampler is expensive since initial search directions are generated by solving many linear programming problems. Thus, we recommend to generate as many samples as possible in one go. However, this might require finer control over the sampling procedure as described in the following section.
Advanced usage¶
Sampler objects¶
The sampling process can be controlled on a lower level by using the sampler classes directly.
[4]:
from cobra.sampling import OptGPSampler, ACHRSampler
Both sampler classes have standardized interfaces and take some additional argument. For instance the thinning
factor. “Thinning” means only recording samples every n iterations. A higher thinning factors mean less correlated samples but also larger computation times. By default the samplers use a thinning factor of 100 which creates roughly uncorrelated samples. If you want less samples but better mixing feel free to increase this parameter. If you want to study convergence for your own
model you might want to set it to 1 to obtain all iterates.
[5]:
achr = ACHRSampler(model, thinning=10)
OptGPSampler
has an additional processes
argument specifying how many processes are used to create parallel sampling chains. This should be in the order of your CPU cores for maximum efficiency. As noted before class initialization can take up to a few minutes due to generation of initial search directions. Sampling on the other hand is quick.
[6]:
optgp = OptGPSampler(model, processes=4)
Sampling and validation¶
Both samplers have a sample function that generates samples from the initialized object and act like the sample
function described above, only that this time it will only accept a single argument, the number of samples. For OptGPSampler
the number of samples should be a multiple of the number of processes, otherwise it will be increased to the nearest multiple automatically.
[7]:
s1 = achr.sample(100)
s2 = optgp.sample(100)
You can call sample
repeatedly and both samplers are optimized to generate large amount of samples without falling into “numerical traps”. All sampler objects have a validate
function in order to check if a set of points are feasible and give detailed information about feasibility violations in a form of a short code denoting feasibility. Here the short code is a combination of any of the following letters:
“v”  valid point
“l”  lower bound violation
“u”  upper bound violation
“e”  equality violation (meaning the point is not a steady state)
For instance for a random flux distribution (should not be feasible):
[8]:
import numpy as np
bad = np.random.uniform(1000, 1000, size=len(model.reactions))
achr.validate(np.atleast_2d(bad))
[8]:
array(['le'], dtype='<U3')
And for our generated samples:
[9]:
achr.validate(s1)
[9]:
array(['v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v'], dtype='<U3')
Even though most models are numerically stable enought that the sampler should only generate valid samples we still urge to check this. validate
is pretty fast and works quickly even for large models and many samples. If you find invalid samples you do not necessarily have to rerun the entire sampling but can exclude them from the sample DataFrame.
[10]:
s1_valid = s1[achr.validate(s1) == "v"]
len(s1_valid)
[10]:
100
Batch sampling¶
Sampler objects are made for generating billions of samples, however using the sample
function might quickly fill up your RAM when working with genomescale models. Here, the batch
method of the sampler objects might come in handy. batch
takes two arguments, the number of samples in each batch and the number of batches. This will make sense with a small example.
Let’s assume we want to quantify what proportion of our samples will grow. For that we might want to generate 10 batches of 50 samples each and measure what percentage of the individual 100 samples show a growth rate larger than 0.1. Finally, we want to calculate the mean and standard deviation of those individual percentages.
[11]:
counts = [np.mean(s.Biomass_Ecoli_core > 0.1) for s in optgp.batch(100, 10)]
print("Usually {:.2f}% + {:.2f}% grow...".format(
np.mean(counts) * 100.0, np.std(counts) * 100.0))
Usually 14.50% + 2.16% grow...
Adding constraints¶
Flux sampling will respect additional contraints defined in the model. For instance we can add a constraint enforcing growth in asimilar manner as the section before.
[12]:
co = model.problem.Constraint(model.reactions.Biomass_Ecoli_core.flux_expression, lb=0.1)
model.add_cons_vars([co])
Note that this is only for demonstration purposes. usually you could set the lower bound of the reaction directly instead of creating a new constraint.
[13]:
s = sample(model, 10)
print(s.Biomass_Ecoli_core)
0 0.124471
1 0.151331
2 0.108145
3 0.144076
4 0.110480
5 0.109024
6 0.111399
7 0.139682
8 0.103511
9 0.116880
Name: Biomass_Ecoli_core, dtype: float64
As we can see our new constraint was respected.
Loopless FBA¶
The goal of this procedure is identification of a thermodynamically consistent flux state without loops, as implied by the name. You can find a more detailed description in the method section at the end of the notebook.
[1]:
%matplotlib inline
import plot_helper
import cobra.test
from cobra import Reaction, Metabolite, Model
from cobra.flux_analysis.loopless import add_loopless, loopless_solution
from cobra.flux_analysis import pfba
Loopless solution¶
Classical loopless approaches as described below are computationally expensive to solve due to the added mixedinteger constraints. A much faster, and pragmatic approach is instead to postprocess flux distributions to simply set fluxes to zero wherever they can be zero without changing the fluxes of any exchange reactions in the model. CycleFreeFlux is an algorithm that can be used to achieve this and in cobrapy it is implemented in the
cobra.flux_analysis.loopless_solution
function. loopless_solution
will identify the closest flux distribution (using only loopless elementary flux modes) to the original one. Note that this will not remove loops which you explicitly requested, for instance by forcing a loop reaction to carry nonzero flux.
Using a larger model than the simple example above, this can be demonstrated as follows
[2]:
salmonella = cobra.test.create_test_model('salmonella')
nominal = salmonella.optimize()
loopless = loopless_solution(salmonella)
[3]:
import pandas
df = pandas.DataFrame(dict(loopless=loopless.fluxes, nominal=nominal.fluxes))
[4]:
df.plot.scatter(x='loopless', y='nominal')
[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x10f7cb3c8>
This functionality can also be used in FVA by using the loopless=True
argument to avoid getting high flux ranges for reactions that essentially only can reach high fluxes if they are allowed to participate in loops (see the simulation notebook) leading to much narrower flux ranges.
Loopless model¶
Cobrapy also includes the “classical” loopless formulation by Schellenberger et. al. implemented in cobra.flux_analysis.add_loopless
modify the model with additional mixedinteger constraints that make thermodynamically infeasible loops impossible. This is much slower than the strategy provided above and should only be used if one of the two following cases applies:
You want to combine a nonlinear (e.g. quadratic) objective with the loopless condition
You want to force the model to be infeasible in the presence of loops independent of the set reaction bounds.
We will demonstrate this with a toy model which has a simple loop cycling A \(\rightarrow\) B \(\rightarrow\) C \(\rightarrow\) A, with A allowed to enter the system and C allowed to leave. A graphical view of the system is drawn below:
[5]:
plot_helper.plot_loop()
[6]:
model = Model()
model.add_metabolites([Metabolite(i) for i in "ABC"])
model.add_reactions([Reaction(i) for i in ["EX_A", "DM_C", "v1", "v2", "v3"]])
model.reactions.EX_A.add_metabolites({"A": 1})
model.reactions.DM_C.add_metabolites({"C": 1})
model.reactions.v1.add_metabolites({"A": 1, "B": 1})
model.reactions.v2.add_metabolites({"B": 1, "C": 1})
model.reactions.v3.add_metabolites({"C": 1, "A": 1})
model.objective = 'DM_C'
While this model contains a loop, a flux state exists which has no flux through reaction v\(_3\), and is identified by loopless FBA.
[7]:
with model:
add_loopless(model)
solution = model.optimize()
print("loopless solution: status = " + solution.status)
print("loopless solution flux: v3 = %.1f" % solution.fluxes["v3"])
loopless solution: status = optimal
loopless solution flux: v3 = 0.0
If there is no forced flux through a loopless reaction, parsimonious FBA will also have no flux through the loop.
[8]:
solution = pfba(model)
print("parsimonious solution: status = " + solution.status)
print("loopless solution flux: v3 = %.1f" % solution.fluxes["v3"])
parsimonious solution: status = optimal
loopless solution flux: v3 = 0.0
However, if flux is forced through v\(_3\), then there is no longer a feasible loopless solution, but the parsimonious solution will still exist.
[9]:
model.reactions.v3.lower_bound = 1
with model:
add_loopless(model)
try:
solution = model.optimize()
except:
print('model is infeasible')
model is infeasible
cobra/util/solver.py:398 UserWarning: solver status is 'infeasible'
[10]:
solution = pfba(model)
print("parsimonious solution: status = " + solution.status)
print("loopless solution flux: v3 = %.1f" % solution.fluxes["v3"])
parsimonious solution: status = optimal
loopless solution flux: v3 = 1.0
Method¶
loopless_solution
is based on a given reference flux distribution. It will look for a new flux distribution with the following requirements:
The objective value is the same as in the reference fluxes.
All exchange fluxes have the same value as in the reference distribution.
All nonexchange fluxes have the same sign (flow in the same direction) as the reference fluxes.
The sum of absolute nonexchange fluxes is minimized.
As proven in the original publication this will identify the “leastloopy” solution closest to the reference fluxes.
If you are using add_loopless
this will use the method described here. In summary, it will add \(G \approx \Delta G\) proxy variables and make loops thermodynamically infeasible. This is achieved by the following formulation.
Here the index j runs over all reactions and the index i only over internal ones. \(a_i\) are indicator variables which equal one if the reaction flux flows in hte forward direction and 0 otherwise. They are used to force the G proxies to always carry the opposite sign of the flux (as it is the case for the “real” \(\Delta G\) values). \(N_{int}\) is the nullspace matrix for internal reactions and is used to find thermodinamically “correct” values for G.
Consistency testing¶
For most problems, multiple flux states can achieve the same optimum and thus we try to obtain a consistent network. By this, we mean that there will be mulitple blocked reactions in the network, which gives rise to this inconsistency. To solve this problem, we use algorithms which can detect all the blocked reactions and also give us consistent networks.
Let us take a toy network, like so:
Here, \(v_{x}\), where \(x \in \{1, 2, \ldots, 6\}\) represent the flux carried by the reactions as shown above.
[1]:
import cobra
[2]:
test_model = cobra.Model("test_model")
v1 = cobra.Reaction("v1")
v2 = cobra.Reaction("v2")
v3 = cobra.Reaction("v3")
v4 = cobra.Reaction("v4")
v5 = cobra.Reaction("v5")
v6 = cobra.Reaction("v6")
test_model.add_reactions([v1, v2, v3, v4, v5, v6])
v1.reaction = "> 2 A"
v2.reaction = "A <> B"
v3.reaction = "A > D"
v4.reaction = "A > C"
v5.reaction = "C > D"
v6.reaction = "D >"
v1.bounds = (0.0, 3.0)
v2.bounds = (3.0, 3.0)
v3.bounds = (0.0, 3.0)
v4.bounds = (0.0, 3.0)
v5.bounds = (0.0, 3.0)
v6.bounds = (0.0, 3.0)
test_model.objective = v6
unknown metabolite 'A' created
unknown metabolite 'B' created
unknown metabolite 'D' created
unknown metabolite 'C' created
Using FVA¶
The first approach we can follow is to use FVA (Flux Variability Analysis) which among many other applications, is used to detect blocked reactions. The cobra.flux_analysis.find_blocked_reactions()
function will return a list of all the blocked reactions obtained using FVA.
[3]:
cobra.flux_analysis.find_blocked_reactions(test_model)
[3]:
['v2']
As we see above, we are able to obtain the blocked reaction, which in this case is \(v_2\).
Using FASTCC¶
The second approach to obtaining consistent network in cobrapy
is to use FASTCC. Using this method, you can expect to efficiently obtain an accurate consistent network. For more details regarding the algorithm, please see Vlassis N, Pacheco MP, Sauter T (2014).
[4]:
consistent_model = cobra.flux_analysis.fastcc(test_model)
consistent_model.reactions
[4]:
[<Reaction v1 at 0x7fc71ddea5c0>,
<Reaction v3 at 0x7fc71ddea630>,
<Reaction v4 at 0x7fc71ddea668>,
<Reaction v5 at 0x7fc71ddea6a0>,
<Reaction v6 at 0x7fc71ddea6d8>]
Similar to the FVA approach, we are able to identify that \(v_2\) is indeed the blocked reaction.
Gapfillling¶
Model gap filling is the task of figuring out which reactions have to be added to a model to make it feasible. Several such algorithms have been reported e.g. Kumar et al. 2009 and Reed et al. 2006. Cobrapy has a gap filling implementation that is very similar to that of Reed et al. where we use a mixedinteger linear program to figure out the smallest number of reactions that need to be added for a userdefined collection of reactions, i.e. a universal model. Briefly, the problem that we try to solve is
Minimize:
subject to
Where l, u are lower and upper bounds for reaction i and z is an indicator variable that is zero if the reaction is not used and otherwise 1, c is a userdefined cost associated with using the ith reaction, \(v^\star\) is the flux of the objective and t a lower bound for that objective. To demonstrate, let’s take a model and remove some essential reactions from it.
[1]:
import cobra.test
from cobra.flux_analysis import gapfill
model = cobra.test.create_test_model("salmonella")
In this model DFructose6phosphate is an essential metabolite. We will remove all the reactions using it, and at them to a separate model.
[2]:
universal = cobra.Model("universal_reactions")
for i in [i.id for i in model.metabolites.f6p_c.reactions]:
reaction = model.reactions.get_by_id(i)
universal.add_reaction(reaction.copy())
model.remove_reactions([reaction])
Now, because of these gaps, the model won’t grow.
[3]:
model.optimize().objective_value
[3]:
0.0
We will use can use the model’s original objective, growth, to figure out which of the removed reactions are required for the model be feasible again. This is very similar to making the ‘nogrowth but growth (NGG)’ predictions of Kumar et al. 2009.
[4]:
solution = gapfill(model, universal, demand_reactions=False)
for reaction in solution[0]:
print(reaction.id)
GF6PTA
F6PP
TKT2
FBP
MAN6PI
We can obtain multiple possible reaction sets by having the algorithm go through multiple iterations.
[5]:
result = gapfill(model, universal, demand_reactions=False, iterations=4)
for i, entries in enumerate(result):
print(" Run %d " % (i + 1))
for e in entries:
print(e.id)
 Run 1 
GF6PTA
F6PP
TKT2
FBP
MAN6PI
 Run 2 
GF6PTA
TALA
PGI
F6PA
MAN6PI
 Run 3 
GF6PTA
F6PP
TKT2
FBP
MAN6PI
 Run 4 
GF6PTA
TALA
PGI
F6PA
MAN6PI
We can also instead of using the original objective, specify a given metabolite that we want the model to be able to produce.
[6]:
with model:
model.objective = model.add_boundary(model.metabolites.f6p_c, type='demand')
solution = gapfill(model, universal)
for reaction in solution[0]:
print(reaction.id)
FBP
Finally, note that using mixedinteger linear programming is computationally quite expensive and for larger models you may want to consider alternative gap filling methods and reconstruction methods.
Growth media¶
The availability of nutrients has a major impact on metabolic fluxes and cobrapy
provides some helpers to manage the exchanges between the external environment and your metabolic model. In experimental settings the “environment” is usually constituted by the growth medium, ergo the concentrations of all metabolites and cofactors available to the modeled organism. However, constraintbased metabolic models only consider fluxes. Thus, you can not simply use concentrations since fluxes have
the unit mmol / [gDW h]
(concentration per gram dry weight of cells and hour).
Also, you are setting an upper bound for the particular import flux and not the flux itself. There are some crude approximations. For instance, if you supply 1 mol of glucose every 24h to 1 gram of bacteria you might set the upper exchange flux for glucose to 1 mol / [1 gDW * 24 h]
since that is the nominal maximum that can be imported. There is no guarantee however that glucose will be consumed with that flux. Thus, the preferred data for exchange fluxes are direct flux measurements as the
ones obtained from timecourse exametabolome measurements for instance.
So how does that look in COBRApy? The current growth medium of a model is managed by the medium
attribute.
[1]:
from cobra.test import create_test_model
model = create_test_model("textbook")
model.medium
[1]:
{'EX_co2_e': 1000.0,
'EX_glc__D_e': 10.0,
'EX_h_e': 1000.0,
'EX_h2o_e': 1000.0,
'EX_nh4_e': 1000.0,
'EX_o2_e': 1000.0,
'EX_pi_e': 1000.0}
This will return a dictionary that contains the upper flux bounds for all active exchange fluxes (the ones having nonzero flux bounds). Right now we see that we have enabled aerobic growth. You can modify a growth medium of a model by assigning a dictionary to model.medium
that maps exchange reactions to their respective upper import bounds. For now let us enforce anaerobic growth by shutting off the oxygen import.
[2]:
medium = model.medium
medium["EX_o2_e"] = 0.0
model.medium = medium
model.medium
[2]:
{'EX_co2_e': 1000.0,
'EX_glc__D_e': 10.0,
'EX_h_e': 1000.0,
'EX_h2o_e': 1000.0,
'EX_nh4_e': 1000.0,
'EX_pi_e': 1000.0}
As we can see oxygen import is now removed from the list of active exchanges and we can verify that this also leads to a lower growth rate.
[3]:
model.slim_optimize()
[3]:
0.21166294973530736
There is a small trap here. model.medium
can not be assigned to directly. So the following will not work:
[4]:
model.medium["EX_co2_e"] = 0.0
model.medium
[4]:
{'EX_co2_e': 1000.0,
'EX_glc__D_e': 10.0,
'EX_h_e': 1000.0,
'EX_h2o_e': 1000.0,
'EX_nh4_e': 1000.0,
'EX_pi_e': 1000.0}
As you can see EX_co2_e
is not set to zero. This is because model.medium is just a copy of the current exchange fluxes. Assigning to it directly with model.medium[...] = ...
will not change the model. You have to assign an entire dictionary with the changed import flux upper bounds:
[5]:
medium = model.medium
medium["EX_co2_e"] = 0.0
model.medium = medium
model.medium # now it worked
[5]:
{'EX_glc__D_e': 10.0,
'EX_h_e': 1000.0,
'EX_h2o_e': 1000.0,
'EX_nh4_e': 1000.0,
'EX_pi_e': 1000.0}
Setting the growth medium also connects to the context manager, so you can set a specific growth medium in a reversible manner.
[6]:
model = create_test_model("textbook")
with model:
medium = model.medium
medium["EX_o2_e"] = 0.0
model.medium = medium
print(model.slim_optimize())
print(model.slim_optimize())
model.medium
0.21166294973530736
0.8739215069684102
[6]:
{'EX_co2_e': 1000.0,
'EX_glc__D_e': 10.0,
'EX_h_e': 1000.0,
'EX_h2o_e': 1000.0,
'EX_nh4_e': 1000.0,
'EX_o2_e': 1000.0,
'EX_pi_e': 1000.0}
So the medium change is only applied within the with
block and reverted automatically.
Minimal media¶
In some cases you might be interested in the smallest growth medium that can maintain a specific growth rate, the so called “minimal medium”. For this we provide the function minimal_medium
which by default obtains the medium with the lowest total import flux. This function needs two arguments: the model and the minimum growth rate (or other objective) the model has to achieve.
[7]:
from cobra.medium import minimal_medium
max_growth = model.slim_optimize()
minimal_medium(model, max_growth)
[7]:
EX_glc__D_e 10.000000
EX_nh4_e 4.765319
EX_o2_e 21.799493
EX_pi_e 3.214895
dtype: float64
So we see that growth is actually limited by glucose import.
Alternatively you might be interested in a minimal medium with the smallest number of active imports. This can be achieved by using the minimize_components
argument (note that this uses a MIP formulation and will therefore be much slower).
[8]:
minimal_medium(model, 0.1, minimize_components=True)
[8]:
EX_glc__D_e 10.000000
EX_nh4_e 1.042503
EX_pi_e 0.703318
dtype: float64
When minimizing the number of import fluxes there may be many alternative solutions. To obtain several of those you can also pass a positive integer to minimize_components
which will give you at most that many alternative solutions. Let us try that with our model and also use the open_exchanges
argument which will assign a large upper bound to all import reactions in the model. The return type will be a pandas.DataFrame
.
[9]:
minimal_medium(model, 0.8, minimize_components=8, open_exchanges=True)
[9]:
0  1  2  3  

EX_fru_e  0.000000  521.357767  0.000000  0.000000 
EX_glc__D_e  0.000000  0.000000  0.000000  519.750758 
EX_gln__L_e  0.000000  40.698058  18.848678  0.000000 
EX_glu__L_e  348.101944  0.000000  0.000000  0.000000 
EX_mal__L_e  0.000000  0.000000  1000.000000  0.000000 
EX_nh4_e  0.000000  0.000000  0.000000  81.026921 
EX_o2_e  500.000000  0.000000  0.000000  0.000000 
EX_pi_e  66.431529  54.913419  12.583458  54.664344 
So there are 4 alternative solutions in total. One aerobic and three anaerobic ones using different carbon sources.
Boundary reactions¶
Apart from exchange reactions there are other types of boundary reactions such as demand or sink reactions. cobrapy
uses various heuristics to identify those and they can be accessed by using the appropriate attribute.
For exchange reactions:
[10]:
ecoli = create_test_model("ecoli")
ecoli.exchanges[0:5]
[10]:
[<Reaction EX_12ppd__R_e at 0x131b4a58d0>,
<Reaction EX_12ppd__S_e at 0x131b471c50>,
<Reaction EX_14glucan_e at 0x131b471e10>,
<Reaction EX_15dap_e at 0x131b471e48>,
<Reaction EX_23camp_e at 0x131b471f98>]
For demand reactions:
[11]:
ecoli.demands
[11]:
[<Reaction DM_4CRSOL at 0x131b3162b0>,
<Reaction DM_5DRIB at 0x131b4712e8>,
<Reaction DM_AACALD at 0x131b471400>,
<Reaction DM_AMOB at 0x131b4714e0>,
<Reaction DM_MTHTHF at 0x131b4715f8>,
<Reaction DM_OXAM at 0x131b4716d8>]
For sink reactions:
[12]:
ecoli.sinks
[12]:
[]
All boundary reactions (any reaction that consumes or introduces mass into the system) can be obtained with the boundary
attribute:
[13]:
ecoli.boundary[0:10]
[13]:
[<Reaction DM_4CRSOL at 0x131b3162b0>,
<Reaction DM_5DRIB at 0x131b4712e8>,
<Reaction DM_AACALD at 0x131b471400>,
<Reaction DM_AMOB at 0x131b4714e0>,
<Reaction DM_MTHTHF at 0x131b4715f8>,
<Reaction DM_OXAM at 0x131b4716d8>,
<Reaction EX_12ppd__R_e at 0x131b4a58d0>,
<Reaction EX_12ppd__S_e at 0x131b471c50>,
<Reaction EX_14glucan_e at 0x131b471e10>,
<Reaction EX_15dap_e at 0x131b471e48>]
Solvers¶
A constraintbased reconstruction and analysis model for biological systems is actually just an application of a class of discrete optimization problems typically solved with linear, mixed integer or quadratic programming techniques. Cobrapy does not implement any algorithm to find solutions to such problems but rather creates a biologically motivated abstraction to these techniques to make it easier to think of how metabolic systems work without paying much attention to how that formulates to an optimization problem.
The actual solving is instead done by tools such as the free software glpk or commercial tools gurobi and cplex which are all made available as a common programmers interface via the optlang package.
When you have defined your model, you can switch solver backend by simply assigning to the model.solver
property.
[1]:
import cobra.test
model = cobra.test.create_test_model('textbook')
[2]:
model.solver = 'glpk'
# or if you have cplex installed
model.solver = 'cplex'
For information on how to configure and tune the solver, please see the documentation for optlang project and note that model.solver
is simply an optlang object of class Model
.
[3]:
type(model.solver)
[3]:
optlang.cplex_interface.Model
Internal solver interfaces¶
Cobrapy also contains its own solver interfaces but these are now deprecated and will be removed completely in the near future. For documentation of how to use these, please refer to older documentation.
Tailored constraints, variables and objectives¶
Thanks to the use of symbolic expressions via the optlang mathematical modeling package, it is relatively straightforward to add new variables, constraints and advanced objectives that cannot be easily formulated as a combination of different reaction and their corresponding upper and lower bounds. Here we demonstrate this optlang functionality which is exposed via the model.solver.interface
.
Constraints¶
Suppose we want to ensure that two reactions have the same flux in our model. We can add this criteria as constraint to our model using the optlang solver interface by simply defining the relevant expression as follows.
[1]:
import cobra.test
model = cobra.test.create_test_model('textbook')
[2]:
same_flux = model.problem.Constraint(
model.reactions.FBA.flux_expression  model.reactions.NH4t.flux_expression,
lb=0,
ub=0)
model.add_cons_vars(same_flux)
The flux for our reaction of interest is obtained by the model.reactions.FBA.flux_expression
which is simply the sum of the forward and reverse flux, i.e.,
[3]:
model.reactions.FBA.flux_expression
[3]:
1.0*FBA  1.0*FBA_reverse_84806
Now I can maximize growth rate whilst the fluxes of reactions ‘FBA’ and ‘NH4t’ are constrained to be (near) identical.
[4]:
solution = model.optimize()
print(solution.fluxes['FBA'], solution.fluxes['NH4t'],
solution.objective_value)
4.66274904774 4.66274904774 0.855110960926157
It is also possible to add many constraints at once. For large models, with constraints involving many reactions, the efficient way to do this is to first build a dictionary of the linear coefficients for every flux, and then add the constraint at once. For example, suppose we want to add a constrain on the sum of the absolute values of every flux in the network to be less than 100:
[5]:
coefficients = dict()
for rxn in model.reactions:
coefficients[rxn.forward_variable] = 1.
coefficients[rxn.reverse_variable] = 1.
constraint = model.problem.Constraint(0, lb=0, ub=100)
model.add_cons_vars(constraint)
model.solver.update()
constraint.set_linear_coefficients(coefficients=coefficients)
Objectives¶
Simple objective such as the maximization of the flux through one or more reactions can conveniently be done by simply assigning to the model.objective
property as we have seen in previous chapters, e.g.,
[5]:
model = cobra.test.create_test_model('textbook')
with model:
model.objective = {model.reactions.Biomass_Ecoli_core: 1}
model.optimize()
print(model.reactions.Biomass_Ecoli_core.flux)
0.8739215069684307
The objectives mathematical expression is seen by
[6]:
model.objective.expression
[6]:
1.0*Biomass_Ecoli_core_reverse_2cdba + 1.0*Biomass_Ecoli_core
But suppose we need a more complicated objective, such as minimizing the Euclidean distance of the solution to the origin minus another variable, while subject to additional linear constraints. This is an objective function with both linear and quadratic components.
Consider the example problem:
min \(\frac{1}{2}\left(x^2 + y^2 \right)  y\)
subject to
\(x + y = 2\)
\(x \ge 0\)
\(y \ge 0\)
This (admittedly very artificial) problem can be visualized graphically where the optimum is indicated by the blue dot on the line of feasible solutions.
[7]:
%matplotlib inline
import plot_helper
plot_helper.plot_qp2()
We return to the textbook model and set the solver to one that can handle quadratic objectives such as cplex. We then add the linear constraint that the sum of our x and y reactions, that we set to FBA and NH4t, must equal 2.
[8]:
model.solver = 'cplex'
sum_two = model.problem.Constraint(
model.reactions.FBA.flux_expression + model.reactions.NH4t.flux_expression,
lb=2,
ub=2)
model.add_cons_vars(sum_two)
Next we add the quadratic objective
[9]:
quadratic_objective = model.problem.Objective(
0.5 * model.reactions.NH4t.flux_expression**2 + 0.5 *
model.reactions.FBA.flux_expression**2 
model.reactions.FBA.flux_expression,
direction='min')
model.objective = quadratic_objective
solution = model.optimize(objective_sense=None)
[10]:
print(solution.fluxes['NH4t'], solution.fluxes['FBA'])
0.5 1.5
Variables¶
We can also create additional variables to facilitate studying the effects of new constraints and variables. Suppose we want to study the difference in flux between nitrogen and carbon uptake whilst we block other reactions. For this it will may help to add another variable representing this difference.
[11]:
model = cobra.test.create_test_model('textbook')
difference = model.problem.Variable('difference')
We use constraints to define what values this variable shall take
[12]:
constraint = model.problem.Constraint(
model.reactions.EX_glc__D_e.flux_expression 
model.reactions.EX_nh4_e.flux_expression  difference,
lb=0,
ub=0)
model.add_cons_vars([difference, constraint])
Now we can access that difference directly during our knockout exploration by looking at its primal value.
[13]:
for reaction in model.reactions[:5]:
with model:
reaction.knock_out()
model.optimize()
print(model.solver.variables.difference.primal)
5.234680806802543
5.2346808068025386
5.234680806802525
1.8644444444444337
1.8644444444444466
Dynamic Flux Balance Analysis (dFBA) in COBRApy¶
The following notebook shows a simple, but slow example of implementing dFBA using COBRApy and scipy.integrate.solve_ivp. This notebook shows a static optimization approach (SOA) implementation and should not be considered production ready.
The model considers only basic MichaelisMenten limited growth on glucose.
[1]:
import numpy as np
from tqdm import tqdm
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
%matplotlib inline
Create or load a cobrapy model. Here, we use the ‘textbook’ ecoli core model.
[2]:
import cobra
from cobra.test import create_test_model
model = create_test_model('textbook')
Set up the dynamic system¶
Dynamic flux balance analysis couples a dynamic system in external cellular concentrations to a pseudosteady state metabolic model.
In this notebook, we define the function add_dynamic_bounds(model, y)
to convert the external metabolite concentrations into bounds on the boundary fluxes in the metabolic model.
[7]:
def add_dynamic_bounds(model, y):
"""Use external concentrations to bound the uptake flux of glucose."""
biomass, glucose = y # expand the boundary species
glucose_max_import = 10 * glucose / (5 + glucose)
model.reactions.EX_glc__D_e.lower_bound = glucose_max_import
def dynamic_system(t, y):
"""Calculate the time derivative of external species."""
biomass, glucose = y # expand the boundary species
# Calculate the specific exchanges fluxes at the given external concentrations.
with model:
add_dynamic_bounds(model, y)
cobra.util.add_lp_feasibility(model)
feasibility = cobra.util.fix_objective_as_constraint(model)
lex_constraints = cobra.util.add_lexicographic_constraints(
model, ['Biomass_Ecoli_core', 'EX_glc__D_e'], ['max', 'max'])
# Since the calculated fluxes are specific rates, we multiply them by the
# biomass concentration to get the bulk exchange rates.
fluxes = lex_constraints.values
fluxes *= biomass
# This implementation is **not** efficient, so I display the current
# simulation time using a progress bar.
if dynamic_system.pbar is not None:
dynamic_system.pbar.update(1)
dynamic_system.pbar.set_description('t = {:.3f}'.format(t))
return fluxes
dynamic_system.pbar = None
def infeasible_event(t, y):
"""
Determine solution feasibility.
Avoiding infeasible solutions is handled by solve_ivp's builtin event detection.
This function resolves the LP to determine whether or not the solution is feasible
(and if not, how far it is from feasibility). When the sign of this function changes
from epsilon to positive, we know the solution is no longer feasible.
"""
with model:
add_dynamic_bounds(model, y)
cobra.util.add_lp_feasibility(model)
feasibility = cobra.util.fix_objective_as_constraint(model)
return feasibility  infeasible_event.epsilon
infeasible_event.epsilon = 1E6
infeasible_event.direction = 1
infeasible_event.terminal = True
Run the dynamic FBA simulation¶
[4]:
ts = np.linspace(0, 15, 100) # Desired integration resolution and interval
y0 = [0.1, 10]
with tqdm() as pbar:
dynamic_system.pbar = pbar
sol = solve_ivp(
fun=dynamic_system,
events=[infeasible_event],
t_span=(ts.min(), ts.max()),
y0=y0,
t_eval=ts,
rtol=1e6,
atol=1e8,
method='BDF'
)
t = 5.804: : 185it [00:16, 11.27it/s]
Because the culture runs out of glucose, the simulation terminates early. The exact time of this ‘cell death’ is recorded in sol.t_events
.
[5]:
sol
[5]:
message: 'A termination event occurred.'
nfev: 179
njev: 2
nlu: 14
sol: None
status: 1
success: True
t: array([0. , 0.15151515, 0.3030303 , 0.45454545, 0.60606061,
0.75757576, 0.90909091, 1.06060606, 1.21212121, 1.36363636,
1.51515152, 1.66666667, 1.81818182, 1.96969697, 2.12121212,
2.27272727, 2.42424242, 2.57575758, 2.72727273, 2.87878788,
3.03030303, 3.18181818, 3.33333333, 3.48484848, 3.63636364,
3.78787879, 3.93939394, 4.09090909, 4.24242424, 4.39393939,
4.54545455, 4.6969697 , 4.84848485, 5. , 5.15151515,
5.3030303 , 5.45454545, 5.60606061, 5.75757576])
t_events: [array([5.80191035])]
y: array([[ 0.1 , 0.10897602, 0.11871674, 0.12927916, 0.14072254,
0.15310825, 0.16649936, 0.18095988, 0.19655403, 0.21334507,
0.23139394, 0.25075753, 0.27148649, 0.29362257, 0.31719545,
0.34221886, 0.36868605, 0.3965646 , 0.42579062, 0.4562623 ,
0.48783322, 0.52030582, 0.55342574, 0.58687742, 0.62028461,
0.65321433, 0.685188 , 0.71570065, 0.74425054, 0.77037369,
0.79368263, 0.81390289, 0.83089676, 0.84467165, 0.85535715,
0.8631722 , 0.86843813, 0.8715096 , 0.8727423 ],
[10. , 9.8947027 , 9.78040248, 9.65642157, 9.52205334,
9.37656372, 9.21919615, 9.04917892, 8.86573366, 8.6680879 ,
8.45549026, 8.22722915, 7.98265735, 7.72122137, 7.442497 ,
7.14623236, 6.83239879, 6.50124888, 6.15338213, 5.78981735,
5.41206877, 5.02222068, 4.62299297, 4.21779303, 3.81071525,
3.40650104, 3.01042208, 2.6280723 , 2.26504645, 1.92656158,
1.61703023, 1.33965598, 1.09616507, 0.88670502, 0.70995892,
0.56344028, 0.44387781, 0.34762375, 0.27100065]])
Plot timelines of biomass and glucose¶
[6]:
ax = plt.subplot(111)
ax.plot(sol.t, sol.y.T[:, 0])
ax2 = plt.twinx(ax)
ax2.plot(sol.t, sol.y.T[:, 1], color='r')
ax.set_ylabel('Biomass', color='b')
ax2.set_ylabel('Glucose', color='r')
[6]:
Text(0, 0.5, 'Glucose')
Using the COBRA toolbox with cobrapy¶
This example demonstrates using COBRA toolbox commands in MATLAB from python through pymatbridge.
[1]:
%load_ext pymatbridge
Starting MATLAB on ZMQ socket ipc:///tmp/pymatbridge57ff542902d94e1a8ed044e391fb0df7
Send 'exit' command to kill the server
....MATLAB started and connected!
[2]:
import cobra.test
m = cobra.test.create_test_model("textbook")
The model_to_pymatbridge function will send the model to the workspace with the given variable name.
[3]:
from cobra.io.mat import model_to_pymatbridge
model_to_pymatbridge(m, variable_name="model")
Now in the MATLAB workspace, the variable name ‘model’ holds a COBRA toolbox struct encoding the model.
[4]:
%%matlab
model
model =
rev: [95x1 double]
metNames: {72x1 cell}
b: [72x1 double]
metCharge: [72x1 double]
c: [95x1 double]
csense: [72x1 char]
genes: {137x1 cell}
metFormulas: {72x1 cell}
rxns: {95x1 cell}
grRules: {95x1 cell}
rxnNames: {95x1 cell}
description: [11x1 char]
S: [72x95 double]
ub: [95x1 double]
lb: [95x1 double]
mets: {72x1 cell}
subSystems: {95x1 cell}
First, we have to initialize the COBRA toolbox in MATLAB.
[5]:
%%matlab silent
warning('off'); % this works around a pymatbridge bug
addpath(genpath('~/cobratoolbox/'));
initCobraToolbox();
Commands from the COBRA toolbox can now be run on the model
[6]:
%%matlab
optimizeCbModel(model)
ans =
x: [95x1 double]
f: 0.8739
y: [71x1 double]
w: [95x1 double]
stat: 1
origStat: 5
solver: 'glpk'
time: 3.2911
FBA in the COBRA toolbox should give the same result as cobrapy (but maybe just a little bit slower :))
[7]:
%time
m.optimize().f
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 5.48 µs
[7]:
0.8739215069684909
FAQ¶
This document will address frequently asked questions not addressed in other pages of the documentation.
How do I install cobrapy?¶
Please see the INSTALL.rst file.
How do I cite cobrapy?¶
Please cite the 2013 publication: 10.1186/17520509774
How do I rename reactions or metabolites?¶
TL;DR Use Model.repair
afterwards
When renaming metabolites or reactions, there are issues because cobra indexes based off of ID’s, which can cause errors. For example:
[1]:
from __future__ import print_function
import cobra.test
model = cobra.test.create_test_model()
for metabolite in model.metabolites:
metabolite.id = "test_" + metabolite.id
try:
model.metabolites.get_by_id(model.metabolites[0].id)
except KeyError as e:
print(repr(e))
The Model.repair function will rebuild the necessary indexes
[2]:
model.repair()
model.metabolites.get_by_id(model.metabolites[0].id)
[2]:
Metabolite identifier  test_dcaACP_c 
Name  DecanoylACPnC100ACP 
Memory address  0x0110f09630 
Formula  C21H39N2O8PRS 
How do I delete a gene?¶
That depends on what precisely you mean by delete a gene.
If you want to simulate the model with a gene knockout, use the cobra.manipulation.delete_model_genes
function. The effects of this function are reversed by cobra.manipulation.undelete_model_genes
.
[3]:
model = cobra.test.create_test_model()
PGI = model.reactions.get_by_id("PGI")
print("bounds before knockout:", (PGI.lower_bound, PGI.upper_bound))
cobra.manipulation.delete_model_genes(model, ["STM4221"])
print("bounds after knockouts", (PGI.lower_bound, PGI.upper_bound))
bounds before knockout: (1000.0, 1000.0)
bounds after knockouts (0.0, 0.0)
If you want to actually remove all traces of a gene from a model, this is more difficult because this will require changing all the gene_reaction_rule
strings for reactions involving the gene.
How do I change the reversibility of a Reaction?¶
Reaction.reversibility
is a property in cobra which is computed when it is requested from the lower and upper bounds.
[4]:
model = cobra.test.create_test_model()
model.reactions.get_by_id("PGI").reversibility
[4]:
True
Trying to set it directly will result in an error or warning:
[5]:
try:
model.reactions.get_by_id("PGI").reversibility = False
except Exception as e:
print(repr(e))
cobra/core/reaction.py:501 UserWarning: Setting reaction reversibility is ignored
The way to change the reversibility is to change the bounds to make the reaction irreversible.
[6]:
model.reactions.get_by_id("PGI").lower_bound = 10
model.reactions.get_by_id("PGI").reversibility
[6]:
False
How do I generate an LP file from a COBRA model?¶
For optlang based solvers¶
With optlang solvers, the LP formulation of a model is obtained by it’s string representation. All solvers behave the same way.
[7]:
with open('test.lp', 'w') as out:
out.write(str(model.solver))
For cobrapy’s internal solvers¶
With the internal solvers, we first create the problem and use functions bundled with the solver.
Please note that unlike the LP file format, the MPS file format does not specify objective direction and is always a minimization. Some (but not all) solvers will rewrite the maximization as a minimization.
[8]:
model = cobra.test.create_test_model()
# glpk through cglpk
glpk = cobra.solvers.cglpk.create_problem(model)
glpk.write("test.lp")
glpk.write("test.mps") # will not rewrite objective
# cplex
cplex = cobra.solvers.cplex_solver.create_problem(model)
cplex.write("test.lp")
cplex.write("test.mps") # rewrites objective
How do I visualize my flux solutions?¶
Please browse the visualization packages on our website for the most recent list of tools.
API Reference¶
This page contains autogenerated API reference documentation 1.
cobra
¶
Subpackages¶
cobra.core
¶
Submodules¶
cobra.core.configuration
¶Provide a global configuration object.
Define a global configuration object. 

class
cobra.core.configuration.
Configuration
(**kwargs)[source]¶ Define a global configuration object.
The attributes of this singleton object are used as default values by cobra functions.

solver
[source]¶ The default solver for new models. The solver choices are the ones provided by optlang and depend on solvers installed in your environment.
 Type
{“glpk”, “cplex”, “gurobi”, “glpk_exact”}

lower_bound
¶ The standard lower bound for reversible reactions (default 1000).
 Type
float, optional

bounds
[source]¶ The default reaction bounds for newly created reactions. The bounds are in the form of lower_bound, upper_bound (default 1000.0, 1000.0).
 Type
tuple of floats

processes
¶ A default number of processes to use where multiprocessing is possible. The default number corresponds to the number of available cores (hyperthreads) minus one.
 Type

cache_directory
[source]¶ A path where the model cache should reside if caching is desired. The default directory depends on the operating system.
 Type
pathlib.Path or str, optional

max_cache_size
¶ The allowed maximum size of the model cache in bytes (default 1 GB).
 Type
int, optional

cache_expiration
¶ The expiration time in seconds for the model cache if any (default None).
 Type
int, optional

_set_default_cache_directory
(self) → None[source]¶ Set the platformdependent default cache directory.

property
solver
(self) → types.ModuleType[source] Return the optlang solver interface.

property
bounds
(self) → Tuple[Optional[Number], Optional[Number]][source] Return the lower, upper reaction bound pair.

property
cache_directory
(self) → pathlib.Path[source] Return the model cache directory.

cobra.core.dictlist
¶A combined dict and list 

class
cobra.core.dictlist.
DictList
(*args)[source]¶ Bases:
list
A combined dict and list
This object behaves like a list, but has the O(1) speed benefits of a dict when looking up elements by their id.

_check
(self, id)[source]¶ make sure duplicate id’s are not added. This function is called before adding in elements.

query
(self, search_function, attribute=None)[source]¶ Query the list
 Parameters
search_function (a string, regular expression or function) – Used to find the matching elements in the list.  a regular expression (possibly compiled), in which case the given attribute of the object should match the regular expression.  a function which takes one argument and returns True for desired values
attribute (string or None) – the name attribute of the object to passed as argument to the search_function. If this is None, the object itself is used.
 Returns
a new list of objects which match the query
 Return type
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model('textbook') >>> model.reactions.query(lambda x: x.boundary) >>> import re >>> regex = re.compile('^g', flags=re.IGNORECASE) >>> model.metabolites.query(regex, attribute='name')

_extend_nocheck
(self, iterable)[source]¶ extends without checking for uniqueness
This function should only be used internally by DictList when it can guarantee elements are already unique (as in when coming from self or other DictList). It will be faster because it skips these checks.

__sub__
(self, other)[source]¶ x.__sub__(y) <==> x  y
 Parameters
other (iterable) – other must contain only unique id’s present in the list

__isub__
(self, other)[source]¶ x.__sub__(y) <==> x = y
 Parameters
other (iterable) – other must contain only unique id’s present in the list

__add__
(self, other)[source]¶ x.__add__(y) <==> x + y
 Parameters
other (iterable) – other must contain only unique id’s which do not intersect with self

__iadd__
(self, other)[source]¶ x.__iadd__(y) <==> x += y
 Parameters
other (iterable) – other must contain only unique id’s whcih do not intersect with self

__getstate__
(self)[source]¶ gets internal state
This is only provided for backwards compatibility so older versions of cobrapy can load pickles generated with cobrapy. In reality, the “_dict” state is ignored when loading a pickle

__setstate__
(self, state)[source]¶ sets internal state
Ignore the passed in state and recalculate it. This is only for compatibility with older pickles which did not correctly specify the initialization class

cobra.core.formula
¶Describes a Chemical Formula 

class
cobra.core.formula.
Formula
(formula=None)[source]¶ Bases:
cobra.core.object.Object
Describes a Chemical Formula
 Parameters
formula (string) – A legal formula string contains only letters and numbers.
cobra.core.gene
¶Parses compiled ast of a gene_reaction_rule and identifies genes 

A Gene in a cobra model 

convert compiled ast to gene_reaction_rule str 

evaluate compiled ast of gene_reaction_rule with knockouts 

parse gpr into AST 

cobra.core.gene.
replacements
= [['.', '__COBRA_DOT__'], ["'", '__COBRA_SQUOTE__'], ['"', '__COBRA_DQUOTE__'], [':', '__COBRA_COLON__'], ['/', '__COBRA_FSLASH__'], ['\\', '__COBRA_BSLASH'], ['', '__COBRA_DASH__'], ['=', '__COBRA_EQ__']][source]¶

cobra.core.gene.
ast2str
(expr, level=0, names=None)[source]¶ convert compiled ast to gene_reaction_rule str
 Parameters
expr (str) – string for a gene reaction rule, e.g “a and b”
level (int) – internal use only
names (dict) – Dict where each element id a gene identifier and the value is the gene name. Use this to get a rule str which uses names instead. This should be done for display purposes only. All gene_reaction_rule strings which are computed with should use the id.
 Returns
The gene reaction rule
 Return type
string

cobra.core.gene.
eval_gpr
(expr, knockouts)[source]¶ evaluate compiled ast of gene_reaction_rule with knockouts

class
cobra.core.gene.
GPRCleaner
[source]¶ Bases:
ast.NodeTransformer
Parses compiled ast of a gene_reaction_rule and identifies genes
Parts of the tree are rewritten to allow periods in gene ID’s and bitwise boolean operations

cobra.core.gene.
parse_gpr
(str_expr)[source]¶ parse gpr into AST
 Parameters
str_expr (string) – string with the gene reaction rule to parse
 Returns
elements ast_tree and gene_ids as a set
 Return type

class
cobra.core.gene.
Gene
(id=None, name='', functional=True)[source]¶ Bases:
cobra.core.species.Species
A Gene in a cobra model
 Parameters
id (string) – The identifier to associate the gene with
name (string) – A longer human readable name for the gene
functional (bool) – Indicates whether the gene is functional. If it is not functional then it cannot be used in an enzyme complex nor can its products be used.

property
functional
(self)[source]¶ A flag indicating if the gene is functional.
Changing the flag is reverted upon exit if executed within the model as context.

knock_out
(self)[source]¶ Knockout gene by marking it as nonfunctional and setting all associated reactions bounds to zero.
The change is reverted upon exit if executed within the model as context.

remove_from_model
(self, model=None, make_dependent_reactions_nonfunctional=True)[source]¶ Removes the association
 Parameters
model (cobra model) – The model to remove the gene from
make_dependent_reactions_nonfunctional (bool) – If True then replace the gene with ‘False’ in the gene association, else replace the gene with ‘True’
Deprecated since version 0.4: Use cobra.manipulation.delete_model_genes to simulate knockouts and cobra.manipulation.remove_genes to remove genes from the model.
cobra.core.group
¶Define the group class.
Manage groups via this implementation of the SBML group specification. 

class
cobra.core.group.
Group
(id, name='', members=None, kind=None)[source]¶ Bases:
cobra.core.object.Object
Manage groups via this implementation of the SBML group specification.
Group is a class for holding information regarding a pathways, subsystems, or other custom groupings of objects within a cobra.Model object.
 Parameters
id (str) – The identifier to associate with this group
name (str, optional) – A human readable name for the group
members (iterable, optional) – A DictList containing references to cobra.Modelassociated objects that belong to the group.
kind ({"collection", "classification", "partonomy"}, optional) – The kind of group, as specified for the Groups feature in the SBML level 3 package specification. Can be any of “classification”, “partonomy”, or “collection”. The default is “collection”. Please consult the SBML level 3 package specification to ensure you are using the proper value for kind. In short, members of a “classification” group should have an “isa” relationship to the group (e.g. member isa polar compound, or member isa transporter). Members of a “partonomy” group should have a “partof” relationship (e.g. member is partof glycolysis). Members of a “collection” group do not have an implied relationship between the members, so use this value for kind when in doubt (e.g. member is a gapfilled reaction, or member is involved in a disease phenotype).
cobra.core.metabolite
¶Define the Metabolite class.
Metabolite is a class for holding information regarding 

class
cobra.core.metabolite.
Metabolite
(id=None, formula=None, name='', charge=None, compartment=None)[source]¶ Bases:
cobra.core.species.Species
Metabolite is a class for holding information regarding a metabolite in a cobra.Reaction object.
 Parameters

property
constraint
(self)[source]¶ Get the constraints associated with this metabolite from the solve
 Returns
the optlang constraint for this metabolite
 Return type
optlang.<interface>.Constraint

property
elements
(self)[source]¶ Dictionary of elements as keys and their count in the metabolite as integer. When set, the formula property is update accordingly

property
y
(self)[source]¶ The shadow price for the metabolite in the most recent solution
Shadow prices are computed from the dual values of the bounds in the solution.

property
shadow_price
(self)[source]¶ The shadow price in the most recent solution.
Shadow price is the dual value of the corresponding constraint in the model.
Warning
Accessing shadow prices through a Solution object is the safer, preferred, and only guaranteed to be correct way. You can see how to do so easily in the examples.
Shadow price is retrieved from the currently defined self._model.solver. The solver status is checked but there are no guarantees that the current solver state is the one you are looking for.
If you modify the underlying model after an optimization, you will retrieve the old optimization values.
 Raises
RuntimeError – If the underlying model was never optimized beforehand or the metabolite is not part of a model.
OptimizationError – If the solver status is anything other than ‘optimal’.
Examples
>>> import cobra >>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> solution = model.optimize() >>> model.metabolites.glc__D_e.shadow_price 0.09166474637510488 >>> solution.shadow_prices.glc__D_e 0.091664746375104883

remove_from_model
(self, destructive=False)[source]¶ Removes the association from self.model
The change is reverted upon exit when using the model as a context.
 Parameters
destructive (bool) – If False then the metabolite is removed from all associated reactions. If True then all associated reactions are removed from the Model.

summary
(self, solution=None, fva=None)[source]¶ Create a summary of the producing and consuming fluxes.
 Parameters
solution (cobra.Solution, optional) – A previous model solution to use for generating the summary. If
None
, the summary method will generate a parsimonious flux distribution (default None).fva (pandas.DataFrame or float, optional) – Whether or not to include flux variability analysis in the output. If given, fva should either be a previous FVA solution matching the model or a float between 0 and 1 representing the fraction of the optimum objective to be searched (default None).
 Returns
 Return type
See also
Reaction.summary()
,Model.summary()
cobra.core.model
¶Define the Model class.
Class representation for a cobra model 

class
cobra.core.model.
Model
(id_or_model=None, name=None)[source]¶ Bases:
cobra.core.object.Object
Class representation for a cobra model
 Parameters
id_or_model (Model, string) – Either an existing Model object in which case a new model object is instantiated with the same properties as the original model, or an identifier to associate with the model as a string.
name (string) – Human readable name for the model

reactions
¶ A DictList where the key is the reaction identifier and the value a Reaction
 Type

metabolites
¶ A DictList where the key is the metabolite identifier and the value a Metabolite
 Type

__getstate__
(self)[source]¶ Get state for serialization.
Ensures that the context stack is cleared prior to serialization, since partial functions cannot be pickled reliably.

property
solver
(self)[source]¶ Get or set the attached solver instance.
The associated the solver object, which manages the interaction with the associated solver, e.g. glpk.
This property is useful for accessing the optimization problem directly and to define additional nonmetabolic constraints.
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> new = model.problem.Constraint(model.objective.expression, >>> lb=0.99) >>> model.solver.add(new)

__add__
(self, other_model)[source]¶ Add the content of another model to this model (+).
The model is copied as a new object, with a new model identifier, and copies of all the reactions in the other model are added to this model. The objective is the sum of the objective expressions for the two models.

__iadd__
(self, other_model)[source]¶ Incrementally add the content of another model to this model (+=).
Copies of all the reactions in the other model are added to this model. The objective is the sum of the objective expressions for the two models.

copy
(self)[source]¶ Provides a partial ‘deepcopy’ of the Model. All of the Metabolite, Gene, and Reaction objects are created anew but in a faster fashion than deepcopy

add_metabolites
(self, metabolite_list)[source]¶ Will add a list of metabolites to the model object and add new constraints accordingly.
The change is reverted upon exit when using the model as a context.
 Parameters
metabolite_list (A list of cobra.core.Metabolite objects) –

remove_metabolites
(self, metabolite_list, destructive=False)[source]¶ Remove a list of metabolites from the the object.
The change is reverted upon exit when using the model as a context.

add_reaction
(self, reaction)[source]¶ Will add a cobra.Reaction object to the model, if reaction.id is not in self.reactions.
 Parameters
reaction (cobra.Reaction) – The reaction to add
(0.6) Use ~cobra.Model.add_reactions instead (Deprecated) –

add_boundary
(self, metabolite, type='exchange', reaction_id=None, lb=None, ub=None, sbo_term=None)[source]¶ Add a boundary reaction for a given metabolite.
There are three different types of predefined boundary reactions: exchange, demand, and sink reactions. An exchange reaction is a reversible, unbalanced reaction that adds to or removes an extracellular metabolite from the extracellular compartment. A demand reaction is an irreversible reaction that consumes an intracellular metabolite. A sink is similar to an exchange but specifically for intracellular metabolites, i.e., a reversible reaction that adds or removes an intracellular metabolite.
If you set the reaction type to something else, you must specify the desired identifier of the created reaction along with its upper and lower bound. The name will be given by the metabolite name and the given type.
 Parameters
metabolite (cobra.Metabolite) – Any given metabolite. The compartment is not checked but you are encouraged to stick to the definition of exchanges and sinks.
type (str, {"exchange", "demand", "sink"}) – Using one of the predefined reaction types is easiest. If you want to create your own kind of boundary reaction choose any other string, e.g., ‘myboundary’.
reaction_id (str, optional) – The ID of the resulting reaction. This takes precedence over the autogenerated identifiers but beware that it might make boundary reactions harder to identify afterwards when using model.boundary or specifically model.exchanges etc.
lb (float, optional) – The lower bound of the resulting reaction.
ub (float, optional) – The upper bound of the resulting reaction.
sbo_term (str, optional) – A correct SBO term is set for the available types. If a custom type is chosen, a suitable SBO term should also be set.
 Returns
The created boundary reaction.
 Return type
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> demand = model.add_boundary(model.metabolites.atp_c, type="demand") >>> demand.id 'DM_atp_c' >>> demand.name 'ATP demand' >>> demand.bounds (0, 1000.0) >>> demand.build_reaction_string() 'atp_c > '

add_reactions
(self, reaction_list)[source]¶ Add reactions to the model.
Reactions with identifiers identical to a reaction already in the model are ignored.
The change is reverted upon exit when using the model as a context.
 Parameters
reaction_list (list) – A list of cobra.Reaction objects

remove_reactions
(self, reactions, remove_orphans=False)[source]¶ Remove reactions from the model.
The change is reverted upon exit when using the model as a context.

add_groups
(self, group_list)[source]¶ Add groups to the model.
Groups with identifiers identical to a group already in the model are ignored.
If any group contains members that are not in the model, these members are added to the model as well. Only metabolites, reactions, and genes can have groups.
 Parameters
group_list (list) – A list of cobra.Group objects to add to the model.

remove_groups
(self, group_list)[source]¶ Remove groups from the model.
Members of each group are not removed from the model (i.e. metabolites, reactions, and genes in the group stay in the model after any groups containing them are removed).
 Parameters
group_list (list) – A list of cobra.Group objects to remove from the model.

get_associated_groups
(self, element)[source]¶ Returns a list of groups that an element (reaction, metabolite, gene) is associated with.
 Parameters
element (cobra.Reaction, cobra.Metabolite, or cobra.Gene) –
 Returns
All groups that the provided object is a member of
 Return type
list of cobra.Group

add_cons_vars
(self, what, **kwargs)[source]¶ Add constraints and variables to the model’s mathematical problem.
Useful for variables and constraints that can not be expressed with reactions and simple lower and upper bounds.
Additions are reversed upon exit if the model itself is used as context.
 Parameters
what (list or tuple of optlang variables or constraints.) – The variables or constraints to add to the model. Must be of class optlang.interface.Variable or optlang.interface.Constraint.
**kwargs (keyword arguments) – Passed to solver.add()

remove_cons_vars
(self, what)[source]¶ Remove variables and constraints from the model’s mathematical problem.
Remove variables and constraints that were added directly to the model’s underlying mathematical problem. Removals are reversed upon exit if the model itself is used as context.
 Parameters
what (list or tuple of optlang variables or constraints.) – The variables or constraints to add to the model. Must be of class optlang.interface.Variable or optlang.interface.Constraint.

property
problem
(self)[source]¶ The interface to the model’s underlying mathematical problem.
Solutions to cobra models are obtained by formulating a mathematical problem and solving it. Cobrapy uses the optlang package to accomplish that and with this property you can get access to the problem interface directly.
 Returns
The problem interface that defines methods for interacting with the problem and associated solver directly.
 Return type
optlang.interface

property
variables
(self)[source]¶ The mathematical variables in the cobra model.
In a cobra model, most variables are reactions. However, for specific use cases, it may also be useful to have other types of variables. This property defines all variables currently associated with the model’s problem.
 Returns
A container with all associated variables.
 Return type
optlang.container.Container

property
constraints
(self)[source]¶ The constraints in the cobra model.
In a cobra model, most constraints are metabolites and their stoichiometries. However, for specific use cases, it may also be useful to have other types of constraints. This property defines all constraints currently associated with the model’s problem.
 Returns
A container with all associated constraints.
 Return type
optlang.container.Container

property
boundary
(self)[source]¶ Boundary reactions in the model. Reactions that either have no substrate or product.

property
exchanges
(self)[source]¶ Exchange reactions in model. Reactions that exchange mass with the exterior. Uses annotations and heuristics to exclude nonexchanges such as sink reactions.

property
demands
(self)[source]¶ Demand reactions in model. Irreversible reactions that accumulate or consume a metabolite in the inside of the model.

property
sinks
(self)[source]¶ Sink reactions in model. Reversible reactions that accumulate or consume a metabolite in the inside of the model.

_populate_solver
(self, reaction_list, metabolite_list=None)[source]¶ Populate attached solver with constraints and variables that model the provided reactions.

slim_optimize
(self, error_value=float('nan'), message=None)[source]¶ Optimize model without creating a solution object.
Creating a full solution object implies fetching shadow prices and flux values for all reactions and metabolites from the solver object. This necessarily takes some time and in cases where only one or two values are of interest, it is recommended to instead use this function which does not create a solution object returning only the value of the objective. Note however that the optimize() function uses efficient means to fetch values so if you need fluxes/shadow prices for more than say 4 reactions/metabolites, then the total speed increase of slim_optimize versus optimize is expected to be small or even negative depending on how you fetch the values after optimization.

optimize
(self, objective_sense=None, raise_error=False)[source]¶ Optimize the model using flux balance analysis.
 Parameters
objective_sense ({None, 'maximize' 'minimize'}, optional) – Whether fluxes should be maximized or minimized. In case of None, the previous direction is used.
raise_error (bool) –
 If true, raise an OptimizationError if solver status is not
optimal.
Notes
Only the most commonly used parameters are presented here. Additional parameters for cobra.solvers may be available and specified with the appropriate keyword argument.

repair
(self, rebuild_index=True, rebuild_relationships=True)[source]¶ Update all indexes and pointers in a model

property
objective
(self)[source]¶ Get or set the solver objective
Before introduction of the optlang based problems, this function returned the objective reactions as a list. With optlang, the objective is not limited a simple linear summation of individual reaction fluxes, making that return value ambiguous. Henceforth, use cobra.util.solver.linear_reaction_coefficients to get a dictionary of reactions with their linear coefficients (empty if there are none)
The set value can be dictionary (reactions as keys, linear coefficients as values), string (reaction identifier), int (reaction index), Reaction or problem.Objective or sympy expression directly interpreted as objectives.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
objective_direction
(self)[source]¶ Get or set the objective direction.
When using a HistoryManager context, this attribute can be set temporarily, reversed when exiting the context.

summary
(self, solution=None, fva=None)[source]¶ Create a summary of the exchange fluxes of the model.
 Parameters
solution (cobra.Solution, optional) – A previous model solution to use for generating the summary. If
None
, the summary method will generate a parsimonious flux distribution (default None).fva (pandas.DataFrame or float, optional) – Whether or not to include flux variability analysis in the output. If given, fva should either be a previous FVA solution matching the model or a float between 0 and 1 representing the fraction of the optimum objective to be searched (default None).
 Returns
 Return type
cobra.ModelSummary
See also
Reaction.summary()
,Metabolite.summary()

__enter__
(self)[source]¶ Record all future changes to the model, undoing them when a call to __exit__ is received

__exit__
(self, type, value, traceback)[source]¶ Pop the top context manager and trigger the undo functions

merge
(self, right, prefix_existing=None, inplace=True, objective='left')[source]¶ Merge two models to create a model with the reactions from both models.
Custom constraints and variables from right models are also copied to left model, however note that, constraints and variables are assumed to be the same if they have the same name.
 rightcobra.Model
The model to add reactions from
 prefix_existingstring
Prefix the reaction identifier in the right that already exist in the left model with this string.
 inplacebool
Add reactions from right directly to left model object. Otherwise, create a new model leaving the left model untouched. When done within the model as context, changes to the models are reverted upon exit.
 objectivestring
One of ‘left’, ‘right’ or ‘sum’ for setting the objective of the resulting model to that of the corresponding model or the sum of both.
cobra.core.object
¶cobra.core.reaction
¶Define the Reaction class.
Reaction is a class for holding information regarding 

class
cobra.core.reaction.
Reaction
(id=None, name='', subsystem='', lower_bound=0.0, upper_bound=None)[source]¶ Bases:
cobra.core.object.Object
Reaction is a class for holding information regarding a biochemical reaction in a cobra.Model object.
Reactions are by default irreversible with bounds (0.0, cobra.Configuration().upper_bound) if no bounds are provided on creation. To create an irreversible reaction use lower_bound=None, resulting in reaction bounds of (cobra.Configuration().lower_bound, cobra.Configuration().upper_bound).
 Parameters

property
flux_expression
(self)[source]¶ Forward flux expression
 Returns
The expression representing the the forward flux (if associated with model), otherwise None. Representing the net flux if model.reversible_encoding == ‘unsplit’ or None if reaction is not associated with a model
 Return type
sympy expression

property
forward_variable
(self)[source]¶ An optlang variable representing the forward flux
 Returns
An optlang variable for the forward flux or None if reaction is not associated with a model.
 Return type
optlang.interface.Variable

property
reverse_variable
(self)[source]¶ An optlang variable representing the reverse flux
 Returns
An optlang variable for the reverse flux or None if reaction is not associated with a model.
 Return type
optlang.interface.Variable

property
objective_coefficient
(self)[source]¶ Get the coefficient for this reaction in a linear objective (float)
Assuming that the objective of the associated model is summation of fluxes from a set of reactions, the coefficient for each reaction can be obtained individually using this property. A more general way is to use the model.objective property directly.

property
lower_bound
(self)[source]¶ Get or set the lower bound
Setting the lower bound (float) will also adjust the associated optlang variables associated with the reaction. Infeasible combinations, such as a lower bound higher than the current upper bound will update the other bound.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
upper_bound
(self)[source]¶ Get or set the upper bound
Setting the upper bound (float) will also adjust the associated optlang variables associated with the reaction. Infeasible combinations, such as a upper bound lower than the current lower bound will update the other bound.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
bounds
(self)[source]¶ Get or set the bounds directly from a tuple
Convenience method for setting upper and lower bounds in one line using a tuple of lower and upper bound. Invalid bounds will raise an AssertionError.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
flux
(self)[source]¶ The flux value in the most recent solution.
Flux is the primal value of the corresponding variable in the model.
Warning
Accessing reaction fluxes through a Solution object is the safer, preferred, and only guaranteed to be correct way. You can see how to do so easily in the examples.
Reaction flux is retrieved from the currently defined self._model.solver. The solver status is checked but there are no guarantees that the current solver state is the one you are looking for.
If you modify the underlying model after an optimization, you will retrieve the old optimization values.
 Raises
RuntimeError – If the underlying model was never optimized beforehand or the reaction is not part of a model.
OptimizationError – If the solver status is anything other than ‘optimal’.
AssertionError – If the flux value is not within the bounds.
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> solution = model.optimize() >>> model.reactions.PFK.flux 7.477381962160283 >>> solution.fluxes.PFK 7.4773819621602833

property
reduced_cost
(self)[source]¶ The reduced cost in the most recent solution.
Reduced cost is the dual value of the corresponding variable in the model.
Warning
Accessing reduced costs through a Solution object is the safer, preferred, and only guaranteed to be correct way. You can see how to do so easily in the examples.
Reduced cost is retrieved from the currently defined self._model.solver. The solver status is checked but there are no guarantees that the current solver state is the one you are looking for.
If you modify the underlying model after an optimization, you will retrieve the old optimization values.
 Raises
RuntimeError – If the underlying model was never optimized beforehand or the reaction is not part of a model.
OptimizationError – If the solver status is anything other than ‘optimal’.
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> solution = model.optimize() >>> model.reactions.PFK.reduced_cost 8.673617379884035e18 >>> solution.reduced_costs.PFK 8.6736173798840355e18

property
gene_name_reaction_rule
(self)[source]¶ Display gene_reaction_rule with names intead.
Do NOT use this string for computation. It is intended to give a representation of the rule using more familiar gene names instead of the often cryptic ids.

property
functional
(self)[source]¶ All required enzymes for reaction are functional.
 Returns
True if the geneproteinreaction (GPR) rule is fulfilled for this reaction, or if reaction is not associated to a model, otherwise False.
 Return type

property
x
(self)[source]¶ The flux through the reaction in the most recent solution.
Flux values are computed from the primal values of the variables in the solution.

property
y
(self)[source]¶ The reduced cost of the reaction in the most recent solution.
Reduced costs are computed from the dual values of the variables in the solution.

property
reversibility
(self)[source]¶ Whether the reaction can proceed in both directions (reversible)
This is computed from the current upper and lower bounds.

property
boundary
(self)[source]¶ Whether or not this reaction is an exchange reaction.
Returns True if the reaction has either no products or reactants.

_update_awareness
(self)[source]¶ Make sure all metabolites and genes that are associated with this reaction are aware of it.

remove_from_model
(self, remove_orphans=False)[source]¶ Removes the reaction from a model.
This removes all associations between a reaction the associated model, metabolites and genes.
The change is reverted upon exit when using the model as a context.
 Parameters
remove_orphans (bool) – Remove orphaned genes and metabolites from the model as well

delete
(self, remove_orphans=False)[source]¶ Removes the reaction from a model.
This removes all associations between a reaction the associated model, metabolites and genes.
The change is reverted upon exit when using the model as a context.
Deprecated, use reaction.remove_from_model instead.
 Parameters
remove_orphans (bool) – Remove orphaned genes and metabolites from the model as well

__setstate__
(self, state)[source]¶ Probably not necessary to set _model as the cobra.Model that contains self sets the _model attribute for all metabolites and genes in the reaction.
However, to increase performance speed we do want to let the metabolite and gene know that they are employed in this reaction

__add__
(self, other)[source]¶ Add two reactions
The stoichiometry will be the combined stoichiometry of the two reactions, and the gene reaction rule will be both rules combined by an and. All other attributes (i.e. reaction bounds) will match those of the first reaction

__imul__
(self, coefficient)[source]¶ Scale coefficients in a reaction by a given value
E.g. A > B becomes 2A > 2B.
If coefficient is less than zero, the reaction is reversed and the bounds are swapped.

get_coefficient
(self, metabolite_id)[source]¶ Return the stoichiometric coefficient of a metabolite.
 Parameters
metabolite_id (str or cobra.Metabolite) –

get_coefficients
(self, metabolite_ids)[source]¶ Return the stoichiometric coefficients for a list of metabolites.
 Parameters
metabolite_ids (iterable) – Containing
str
or ``cobra.Metabolite``s.

add_metabolites
(self, metabolites_to_add, combine=True, reversibly=True)[source]¶ Add metabolites and stoichiometric coefficients to the reaction. If the final coefficient for a metabolite is 0 then it is removed from the reaction.
The change is reverted upon exit when using the model as a context.
 Parameters
metabolites_to_add (dict) – Dictionary with metabolite objects or metabolite identifiers as keys and coefficients as values. If keys are strings (name of a metabolite) the reaction must already be part of a model and a metabolite with the given name must exist in the model.
combine (bool) – Describes behavior a metabolite already exists in the reaction. True causes the coefficients to be added. False causes the coefficient to be replaced.
reversibly (bool) – Whether to add the change to the context to make the change reversibly or not (primarily intended for internal use).

subtract_metabolites
(self, metabolites, combine=True, reversibly=True)[source]¶ Subtract metabolites from a reaction.
That means add the metabolites with 1*coefficient. If the final coefficient for a metabolite is 0 then the metabolite is removed from the reaction.
Notes
A final coefficient < 0 implies a reactant.
The change is reverted upon exit when using the model as a context.
 Parameters
metabolites (dict) – Dictionary where the keys are of class Metabolite and the values are the coefficients. These metabolites will be added to the reaction.
combine (bool) – Describes behavior a metabolite already exists in the reaction. True causes the coefficients to be added. False causes the coefficient to be replaced.
reversibly (bool) – Whether to add the change to the context to make the change reversibly or not (primarily intended for internal use).

build_reaction_string
(self, use_metabolite_names=False)[source]¶ Generate a human readable reaction string

check_mass_balance
(self)[source]¶ Compute mass and charge balance for the reaction
returns a dict of {element: amount} for unbalanced elements. “charge” is treated as an element in this dict This should be empty for balanced reactions.

_associate_gene
(self, cobra_gene)[source]¶ Associates a cobra.Gene object with a cobra.Reaction.
 Parameters
cobra_gene (cobra.core.Gene.Gene) –

_dissociate_gene
(self, cobra_gene)[source]¶ Dissociates a cobra.Gene object with a cobra.Reaction.
 Parameters
cobra_gene (cobra.core.Gene.Gene) –

build_reaction_from_string
(self, reaction_str, verbose=True, fwd_arrow=None, rev_arrow=None, reversible_arrow=None, term_split='+')[source]¶ Builds reaction from reaction equation reaction_str using parser
Takes a string and using the specifications supplied in the optional arguments infers a set of metabolites, metabolite compartments and stoichiometries for the reaction. It also infers the reversibility of the reaction from the reaction arrow.
Changes to the associated model are reverted upon exit when using the model as a context.
 Parameters
reaction_str (string) – a string containing a reaction formula (equation)
verbose (bool) – setting verbosity of function
fwd_arrow (re.compile) – for forward irreversible reaction arrows
rev_arrow (re.compile) – for backward irreversible reaction arrows
reversible_arrow (re.compile) – for reversible reaction arrows
term_split (string) – dividing individual metabolite entries

summary
(self, solution=None, fva=None)[source]¶ Create a summary of the reaction flux.
 Parameters
solution (cobra.Solution, optional) – A previous model solution to use for generating the summary. If
None
, the summary method will generate a parsimonious flux distribution (default None).fva (pandas.DataFrame or float, optional) – Whether or not to include flux variability analysis in the output. If given, fva should either be a previous FVA solution matching the model or a float between 0 and 1 representing the fraction of the optimum objective to be searched (default None).
 Returns
 Return type
cobra.summary.ReactionSummary
See also
Metabolite.summary()
,Model.summary()
cobra.core.solution
¶Provide unified interfaces to optimization solutions.
A unified interface to a cobra.Model optimization solution. 

Legacy support for an interface to a cobra.Model optimization solution. 

Generate a solution representation of the current solver state. 

class
cobra.core.solution.
Solution
(objective_value, status, fluxes, reduced_costs=None, shadow_prices=None, **kwargs)[source]¶ Bases:
object
A unified interface to a cobra.Model optimization solution.
Notes
Solution is meant to be constructed by get_solution please look at that function to fully understand the Solution class.

fluxes
¶ Contains the reaction fluxes (primal values of variables).
 Type
pandas.Series

reduced_costs
¶ Contains reaction reduced costs (dual values of variables).
 Type
pandas.Series

shadow_prices
¶ Contains metabolite shadow prices (dual values of constraints).
 Type
pandas.Series


class
cobra.core.solution.
LegacySolution
(f, x=None, x_dict=None, y=None, y_dict=None, solver=None, the_time=0, status='NA', **kwargs)[source]¶ Bases:
object
Legacy support for an interface to a cobra.Model optimization solution.

x
¶ List or Array of the fluxes (primal values).
 Type
iterable

y
¶ List or Array of the dual values.
 Type
iterable
Warning
The LegacySolution class and its interface is deprecated.


cobra.core.solution.
get_solution
(model, reactions=None, metabolites=None, raise_error=False)[source]¶ Generate a solution representation of the current solver state.
 Parameters
model (cobra.Model) – The model whose reactions to retrieve values for.
reactions (list, optional) – An iterable of cobra.Reaction objects. Uses model.reactions by default.
metabolites (list, optional) – An iterable of cobra.Metabolite objects. Uses model.metabolites by default.
raise_error (bool) – If true, raise an OptimizationError if solver status is not optimal.
 Returns
 Return type
Note
This is only intended for the optlang solver interfaces and not the legacy solvers.
cobra.core.species
¶Species is a class for holding information regarding 

class
cobra.core.species.
Species
(id=None, name=None)[source]¶ Bases:
cobra.core.object.Object
Species is a class for holding information regarding a chemical Species
 Parameters
id (string) – An identifier for the chemical species
name (string) – A human readable name.

__getstate__
(self)[source]¶ Remove the references to container reactions when serializing to avoid problems associated with recursion.
Package Contents¶
Define a global configuration object. 

A combined dict and list 

A Gene in a cobra model 

Metabolite is a class for holding information regarding 

Class representation for a cobra model 

Defines common behavior of object in cobra.core 

Reaction is a class for holding information regarding 

Manage groups via this implementation of the SBML group specification. 

A unified interface to a cobra.Model optimization solution. 

Legacy support for an interface to a cobra.Model optimization solution. 

Species is a class for holding information regarding 

Generate a solution representation of the current solver state. 

class
cobra.core.
Configuration
(**kwargs)[source]¶ Define a global configuration object.
The attributes of this singleton object are used as default values by cobra functions.

solver
¶ The default solver for new models. The solver choices are the ones provided by optlang and depend on solvers installed in your environment.
 Type
{“glpk”, “cplex”, “gurobi”, “glpk_exact”}

lower_bound
¶ The standard lower bound for reversible reactions (default 1000).
 Type
float, optional

bounds
¶ The default reaction bounds for newly created reactions. The bounds are in the form of lower_bound, upper_bound (default 1000.0, 1000.0).
 Type
tuple of floats

processes
¶ A default number of processes to use where multiprocessing is possible. The default number corresponds to the number of available cores (hyperthreads) minus one.
 Type

cache_directory
¶ A path where the model cache should reside if caching is desired. The default directory depends on the operating system.
 Type
pathlib.Path or str, optional

max_cache_size
¶ The allowed maximum size of the model cache in bytes (default 1 GB).
 Type
int, optional

cache_expiration
¶ The expiration time in seconds for the model cache if any (default None).
 Type
int, optional

_set_default_solver
(self) → None¶ Set the default solver from a preferred order.

_set_default_processes
(self) → None¶ Set the default number of processes.

_set_default_cache_directory
(self) → None¶ Set the platformdependent default cache directory.

property
solver
(self) → types.ModuleType Return the optlang solver interface.

property
bounds
(self) → Tuple[Optional[Number], Optional[Number]] Return the lower, upper reaction bound pair.

property
cache_directory
(self) → pathlib.Path Return the model cache directory.

__repr__
(self) → str¶ Return a string representation of the current configuration values.

_repr_html_
(self) → str¶ Return a rich HTML representation of the current configuration values.
Notes
This special method is used automatically in Jupyter notebooks to display a result from a cell.


class
cobra.core.
DictList
(*args)[source]¶ Bases:
list
A combined dict and list
This object behaves like a list, but has the O(1) speed benefits of a dict when looking up elements by their id.

has_id
(self, id)¶

_check
(self, id)¶ make sure duplicate id’s are not added. This function is called before adding in elements.

_generate_index
(self)¶ rebuild the _dict index

get_by_id
(self, id)¶ return the element with a matching id

list_attr
(self, attribute)¶ return a list of the given attribute for every object

get_by_any
(self, iterable)¶ Get a list of members using several different ways of indexing

query
(self, search_function, attribute=None)¶ Query the list
 Parameters
search_function (a string, regular expression or function) – Used to find the matching elements in the list.  a regular expression (possibly compiled), in which case the given attribute of the object should match the regular expression.  a function which takes one argument and returns True for desired values
attribute (string or None) – the name attribute of the object to passed as argument to the search_function. If this is None, the object itself is used.
 Returns
a new list of objects which match the query
 Return type
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model('textbook') >>> model.reactions.query(lambda x: x.boundary) >>> import re >>> regex = re.compile('^g', flags=re.IGNORECASE) >>> model.metabolites.query(regex, attribute='name')

_replace_on_id
(self, new_object)¶ Replace an object by another with the same id.

append
(self, object)¶ append object to end

union
(self, iterable)¶ adds elements with id’s not already in the model

extend
(self, iterable)¶ extend list by appending elements from the iterable

_extend_nocheck
(self, iterable)¶ extends without checking for uniqueness
This function should only be used internally by DictList when it can guarantee elements are already unique (as in when coming from self or other DictList). It will be faster because it skips these checks.

__sub__
(self, other)¶ x.__sub__(y) <==> x  y
 Parameters
other (iterable) – other must contain only unique id’s present in the list

__isub__
(self, other)¶ x.__sub__(y) <==> x = y
 Parameters
other (iterable) – other must contain only unique id’s present in the list

__add__
(self, other)¶ x.__add__(y) <==> x + y
 Parameters
other (iterable) – other must contain only unique id’s which do not intersect with self

__iadd__
(self, other)¶ x.__iadd__(y) <==> x += y
 Parameters
other (iterable) – other must contain only unique id’s whcih do not intersect with self

__reduce__
(self)¶ Helper for pickle.

__getstate__
(self)¶ gets internal state
This is only provided for backwards compatibility so older versions of cobrapy can load pickles generated with cobrapy. In reality, the “_dict” state is ignored when loading a pickle

__setstate__
(self, state)¶ sets internal state
Ignore the passed in state and recalculate it. This is only for compatibility with older pickles which did not correctly specify the initialization class

index
(self, id, *args)¶ Determine the position in the list
id: A string or a
Object

__contains__
(self, object)¶ DictList.__contains__(object) <==> object in DictList
object: str or
Object

__copy__
(self)¶

insert
(self, index, object)¶ insert object before index

pop
(self, *args)¶ remove and return item at index (default last).

add
(self, x)¶ Opposite of remove. Mirrors set.add

remove
(self, x)¶ Warning
Internal use only

reverse
(self)¶ reverse IN PLACE

sort
(self, cmp=None, key=None, reverse=False)¶ stable sort IN PLACE
cmp(x, y) > 1, 0, 1

__getitem__
(self, i)¶ x.__getitem__(y) <==> x[y]

__setitem__
(self, i, y)¶ Set self[key] to value.

__delitem__
(self, index)¶ Delete self[key].

__getslice__
(self, i, j)¶

__setslice__
(self, i, j, y)¶

__delslice__
(self, i, j)¶

__getattr__
(self, attr)¶

__dir__
(self)¶ Default dir() implementation.


class
cobra.core.
Gene
(id=None, name='', functional=True)[source]¶ Bases:
cobra.core.species.Species
A Gene in a cobra model
 Parameters
id (string) – The identifier to associate the gene with
name (string) – A longer human readable name for the gene
functional (bool) – Indicates whether the gene is functional. If it is not functional then it cannot be used in an enzyme complex nor can its products be used.

property
functional
(self)¶ A flag indicating if the gene is functional.
Changing the flag is reverted upon exit if executed within the model as context.

knock_out
(self)¶ Knockout gene by marking it as nonfunctional and setting all associated reactions bounds to zero.
The change is reverted upon exit if executed within the model as context.

remove_from_model
(self, model=None, make_dependent_reactions_nonfunctional=True)¶ Removes the association
 Parameters
model (cobra model) – The model to remove the gene from
make_dependent_reactions_nonfunctional (bool) – If True then replace the gene with ‘False’ in the gene association, else replace the gene with ‘True’
Deprecated since version 0.4: Use cobra.manipulation.delete_model_genes to simulate knockouts and cobra.manipulation.remove_genes to remove genes from the model.

_repr_html_
(self)¶

class
cobra.core.
Metabolite
(id=None, formula=None, name='', charge=None, compartment=None)[source]¶ Bases:
cobra.core.species.Species
Metabolite is a class for holding information regarding a metabolite in a cobra.Reaction object.
 Parameters

_set_id_with_model
(self, value)¶

property
constraint
(self)¶ Get the constraints associated with this metabolite from the solve
 Returns
the optlang constraint for this metabolite
 Return type
optlang.<interface>.Constraint

property
elements
(self)¶ Dictionary of elements as keys and their count in the metabolite as integer. When set, the formula property is update accordingly

property
formula_weight
(self)¶ Calculate the formula weight

property
y
(self)¶ The shadow price for the metabolite in the most recent solution
Shadow prices are computed from the dual values of the bounds in the solution.

property
shadow_price
(self)¶ The shadow price in the most recent solution.
Shadow price is the dual value of the corresponding constraint in the model.
Warning
Accessing shadow prices through a Solution object is the safer, preferred, and only guaranteed to be correct way. You can see how to do so easily in the examples.
Shadow price is retrieved from the currently defined self._model.solver. The solver status is checked but there are no guarantees that the current solver state is the one you are looking for.
If you modify the underlying model after an optimization, you will retrieve the old optimization values.
 Raises
RuntimeError – If the underlying model was never optimized beforehand or the metabolite is not part of a model.
OptimizationError – If the solver status is anything other than ‘optimal’.
Examples
>>> import cobra >>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> solution = model.optimize() >>> model.metabolites.glc__D_e.shadow_price 0.09166474637510488 >>> solution.shadow_prices.glc__D_e 0.091664746375104883

remove_from_model
(self, destructive=False)¶ Removes the association from self.model
The change is reverted upon exit when using the model as a context.
 Parameters
destructive (bool) – If False then the metabolite is removed from all associated reactions. If True then all associated reactions are removed from the Model.

summary
(self, solution=None, fva=None)¶ Create a summary of the producing and consuming fluxes.
 Parameters
solution (cobra.Solution, optional) – A previous model solution to use for generating the summary. If
None
, the summary method will generate a parsimonious flux distribution (default None).fva (pandas.DataFrame or float, optional) – Whether or not to include flux variability analysis in the output. If given, fva should either be a previous FVA solution matching the model or a float between 0 and 1 representing the fraction of the optimum objective to be searched (default None).
 Returns
 Return type
See also

_repr_html_
(self)¶

class
cobra.core.
Model
(id_or_model=None, name=None)[source]¶ Bases:
cobra.core.object.Object
Class representation for a cobra model
 Parameters
id_or_model (Model, string) – Either an existing Model object in which case a new model object is instantiated with the same properties as the original model, or an identifier to associate with the model as a string.
name (string) – Human readable name for the model

reactions
¶ A DictList where the key is the reaction identifier and the value a Reaction
 Type

metabolites
¶ A DictList where the key is the metabolite identifier and the value a Metabolite
 Type

__setstate__
(self, state)¶ Make sure all cobra.Objects in the model point to the model.

__getstate__
(self)¶ Get state for serialization.
Ensures that the context stack is cleared prior to serialization, since partial functions cannot be pickled reliably.

property
solver
(self)¶ Get or set the attached solver instance.
The associated the solver object, which manages the interaction with the associated solver, e.g. glpk.
This property is useful for accessing the optimization problem directly and to define additional nonmetabolic constraints.
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> new = model.problem.Constraint(model.objective.expression, >>> lb=0.99) >>> model.solver.add(new)

property
tolerance
(self)¶

property
description
(self)¶

get_metabolite_compartments
(self)¶ Return all metabolites’ compartments.

property
compartments
(self)¶

property
medium
(self)¶

__add__
(self, other_model)¶ Add the content of another model to this model (+).
The model is copied as a new object, with a new model identifier, and copies of all the reactions in the other model are added to this model. The objective is the sum of the objective expressions for the two models.

__iadd__
(self, other_model)¶ Incrementally add the content of another model to this model (+=).
Copies of all the reactions in the other model are added to this model. The objective is the sum of the objective expressions for the two models.

copy
(self)¶ Provides a partial ‘deepcopy’ of the Model. All of the Metabolite, Gene, and Reaction objects are created anew but in a faster fashion than deepcopy

add_metabolites
(self, metabolite_list)¶ Will add a list of metabolites to the model object and add new constraints accordingly.
The change is reverted upon exit when using the model as a context.
 Parameters
metabolite_list (A list of cobra.core.Metabolite objects) –

remove_metabolites
(self, metabolite_list, destructive=False)¶ Remove a list of metabolites from the the object.
The change is reverted upon exit when using the model as a context.

add_reaction
(self, reaction)¶ Will add a cobra.Reaction object to the model, if reaction.id is not in self.reactions.
 Parameters
reaction (cobra.Reaction) – The reaction to add
(0.6) Use ~cobra.Model.add_reactions instead (Deprecated) –

add_boundary
(self, metabolite, type='exchange', reaction_id=None, lb=None, ub=None, sbo_term=None)¶ Add a boundary reaction for a given metabolite.
There are three different types of predefined boundary reactions: exchange, demand, and sink reactions. An exchange reaction is a reversible, unbalanced reaction that adds to or removes an extracellular metabolite from the extracellular compartment. A demand reaction is an irreversible reaction that consumes an intracellular metabolite. A sink is similar to an exchange but specifically for intracellular metabolites, i.e., a reversible reaction that adds or removes an intracellular metabolite.
If you set the reaction type to something else, you must specify the desired identifier of the created reaction along with its upper and lower bound. The name will be given by the metabolite name and the given type.
 Parameters
metabolite (cobra.Metabolite) – Any given metabolite. The compartment is not checked but you are encouraged to stick to the definition of exchanges and sinks.
type (str, {"exchange", "demand", "sink"}) – Using one of the predefined reaction types is easiest. If you want to create your own kind of boundary reaction choose any other string, e.g., ‘myboundary’.
reaction_id (str, optional) – The ID of the resulting reaction. This takes precedence over the autogenerated identifiers but beware that it might make boundary reactions harder to identify afterwards when using model.boundary or specifically model.exchanges etc.
lb (float, optional) – The lower bound of the resulting reaction.
ub (float, optional) – The upper bound of the resulting reaction.
sbo_term (str, optional) – A correct SBO term is set for the available types. If a custom type is chosen, a suitable SBO term should also be set.
 Returns
The created boundary reaction.
 Return type
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> demand = model.add_boundary(model.metabolites.atp_c, type="demand") >>> demand.id 'DM_atp_c' >>> demand.name 'ATP demand' >>> demand.bounds (0, 1000.0) >>> demand.build_reaction_string() 'atp_c > '

add_reactions
(self, reaction_list)¶ Add reactions to the model.
Reactions with identifiers identical to a reaction already in the model are ignored.
The change is reverted upon exit when using the model as a context.
 Parameters
reaction_list (list) – A list of cobra.Reaction objects

remove_reactions
(self, reactions, remove_orphans=False)¶ Remove reactions from the model.
The change is reverted upon exit when using the model as a context.

add_groups
(self, group_list)¶ Add groups to the model.
Groups with identifiers identical to a group already in the model are ignored.
If any group contains members that are not in the model, these members are added to the model as well. Only metabolites, reactions, and genes can have groups.
 Parameters
group_list (list) – A list of cobra.Group objects to add to the model.

remove_groups
(self, group_list)¶ Remove groups from the model.
Members of each group are not removed from the model (i.e. metabolites, reactions, and genes in the group stay in the model after any groups containing them are removed).
 Parameters
group_list (list) – A list of cobra.Group objects to remove from the model.

get_associated_groups
(self, element)¶ Returns a list of groups that an element (reaction, metabolite, gene) is associated with.
 Parameters
element (cobra.Reaction, cobra.Metabolite, or cobra.Gene) –
 Returns
All groups that the provided object is a member of
 Return type
list of cobra.Group

add_cons_vars
(self, what, **kwargs)¶ Add constraints and variables to the model’s mathematical problem.
Useful for variables and constraints that can not be expressed with reactions and simple lower and upper bounds.
Additions are reversed upon exit if the model itself is used as context.
 Parameters
what (list or tuple of optlang variables or constraints.) – The variables or constraints to add to the model. Must be of class optlang.interface.Variable or optlang.interface.Constraint.
**kwargs (keyword arguments) – Passed to solver.add()

remove_cons_vars
(self, what)¶ Remove variables and constraints from the model’s mathematical problem.
Remove variables and constraints that were added directly to the model’s underlying mathematical problem. Removals are reversed upon exit if the model itself is used as context.
 Parameters
what (list or tuple of optlang variables or constraints.) – The variables or constraints to add to the model. Must be of class optlang.interface.Variable or optlang.interface.Constraint.

property
problem
(self)¶ The interface to the model’s underlying mathematical problem.
Solutions to cobra models are obtained by formulating a mathematical problem and solving it. Cobrapy uses the optlang package to accomplish that and with this property you can get access to the problem interface directly.
 Returns
The problem interface that defines methods for interacting with the problem and associated solver directly.
 Return type
optlang.interface

property
variables
(self)¶ The mathematical variables in the cobra model.
In a cobra model, most variables are reactions. However, for specific use cases, it may also be useful to have other types of variables. This property defines all variables currently associated with the model’s problem.
 Returns
A container with all associated variables.
 Return type
optlang.container.Container

property
constraints
(self)¶ The constraints in the cobra model.
In a cobra model, most constraints are metabolites and their stoichiometries. However, for specific use cases, it may also be useful to have other types of constraints. This property defines all constraints currently associated with the model’s problem.
 Returns
A container with all associated constraints.
 Return type
optlang.container.Container

property
boundary
(self)¶ Boundary reactions in the model. Reactions that either have no substrate or product.

property
exchanges
(self)¶ Exchange reactions in model. Reactions that exchange mass with the exterior. Uses annotations and heuristics to exclude nonexchanges such as sink reactions.

property
demands
(self)¶ Demand reactions in model. Irreversible reactions that accumulate or consume a metabolite in the inside of the model.

property
sinks
(self)¶ Sink reactions in model. Reversible reactions that accumulate or consume a metabolite in the inside of the model.

_populate_solver
(self, reaction_list, metabolite_list=None)¶ Populate attached solver with constraints and variables that model the provided reactions.

slim_optimize
(self, error_value=float('nan'), message=None)¶ Optimize model without creating a solution object.
Creating a full solution object implies fetching shadow prices and flux values for all reactions and metabolites from the solver object. This necessarily takes some time and in cases where only one or two values are of interest, it is recommended to instead use this function which does not create a solution object returning only the value of the objective. Note however that the optimize() function uses efficient means to fetch values so if you need fluxes/shadow prices for more than say 4 reactions/metabolites, then the total speed increase of slim_optimize versus optimize is expected to be small or even negative depending on how you fetch the values after optimization.

optimize
(self, objective_sense=None, raise_error=False)¶ Optimize the model using flux balance analysis.
 Parameters
objective_sense ({None, 'maximize' 'minimize'}, optional) – Whether fluxes should be maximized or minimized. In case of None, the previous direction is used.
raise_error (bool) –
 If true, raise an OptimizationError if solver status is not
optimal.
Notes
Only the most commonly used parameters are presented here. Additional parameters for cobra.solvers may be available and specified with the appropriate keyword argument.

repair
(self, rebuild_index=True, rebuild_relationships=True)¶ Update all indexes and pointers in a model

property
objective
(self)¶ Get or set the solver objective
Before introduction of the optlang based problems, this function returned the objective reactions as a list. With optlang, the objective is not limited a simple linear summation of individual reaction fluxes, making that return value ambiguous. Henceforth, use cobra.util.solver.linear_reaction_coefficients to get a dictionary of reactions with their linear coefficients (empty if there are none)
The set value can be dictionary (reactions as keys, linear coefficients as values), string (reaction identifier), int (reaction index), Reaction or problem.Objective or sympy expression directly interpreted as objectives.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
objective_direction
(self)¶ Get or set the objective direction.
When using a HistoryManager context, this attribute can be set temporarily, reversed when exiting the context.

summary
(self, solution=None, fva=None)¶ Create a summary of the exchange fluxes of the model.
 Parameters
solution (cobra.Solution, optional) – A previous model solution to use for generating the summary. If
None
, the summary method will generate a parsimonious flux distribution (default None).fva (pandas.DataFrame or float, optional) – Whether or not to include flux variability analysis in the output. If given, fva should either be a previous FVA solution matching the model or a float between 0 and 1 representing the fraction of the optimum objective to be searched (default None).
 Returns
 Return type
cobra.ModelSummary
See also

__enter__
(self)¶ Record all future changes to the model, undoing them when a call to __exit__ is received

__exit__
(self, type, value, traceback)¶ Pop the top context manager and trigger the undo functions

merge
(self, right, prefix_existing=None, inplace=True, objective='left')¶ Merge two models to create a model with the reactions from both models.
Custom constraints and variables from right models are also copied to left model, however note that, constraints and variables are assumed to be the same if they have the same name.
 rightcobra.Model
The model to add reactions from
 prefix_existingstring
Prefix the reaction identifier in the right that already exist in the left model with this string.
 inplacebool
Add reactions from right directly to left model object. Otherwise, create a new model leaving the left model untouched. When done within the model as context, changes to the models are reverted upon exit.
 objectivestring
One of ‘left’, ‘right’ or ‘sum’ for setting the objective of the resulting model to that of the corresponding model or the sum of both.

_repr_html_
(self)¶

class
cobra.core.
Object
(id=None, name='')[source]¶ Bases:
object
Defines common behavior of object in cobra.core

property
id
(self)¶

_set_id_with_model
(self, value)¶

property
annotation
(self)¶

__getstate__
(self)¶ To prevent excessive replication during deepcopy.

__repr__
(self)¶ Return repr(self).

__str__
(self)¶ Return str(self).

property

class
cobra.core.
Reaction
(id=None, name='', subsystem='', lower_bound=0.0, upper_bound=None)[source]¶ Bases:
cobra.core.object.Object
Reaction is a class for holding information regarding a biochemical reaction in a cobra.Model object.
Reactions are by default irreversible with bounds (0.0, cobra.Configuration().upper_bound) if no bounds are provided on creation. To create an irreversible reaction use lower_bound=None, resulting in reaction bounds of (cobra.Configuration().lower_bound, cobra.Configuration().upper_bound).
 Parameters

__radd__
¶

_set_id_with_model
(self, value)¶

property
reverse_id
(self)¶ Generate the id of reverse_variable from the reaction’s id.

property
flux_expression
(self)¶ Forward flux expression
 Returns
The expression representing the the forward flux (if associated with model), otherwise None. Representing the net flux if model.reversible_encoding == ‘unsplit’ or None if reaction is not associated with a model
 Return type
sympy expression

property
forward_variable
(self)¶ An optlang variable representing the forward flux
 Returns
An optlang variable for the forward flux or None if reaction is not associated with a model.
 Return type
optlang.interface.Variable

property
reverse_variable
(self)¶ An optlang variable representing the reverse flux
 Returns
An optlang variable for the reverse flux or None if reaction is not associated with a model.
 Return type
optlang.interface.Variable

property
objective_coefficient
(self)¶ Get the coefficient for this reaction in a linear objective (float)
Assuming that the objective of the associated model is summation of fluxes from a set of reactions, the coefficient for each reaction can be obtained individually using this property. A more general way is to use the model.objective property directly.

__copy__
(self)¶

__deepcopy__
(self, memo)¶

static
_check_bounds
(lb, ub)¶

update_variable_bounds
(self)¶

property
lower_bound
(self)¶ Get or set the lower bound
Setting the lower bound (float) will also adjust the associated optlang variables associated with the reaction. Infeasible combinations, such as a lower bound higher than the current upper bound will update the other bound.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
upper_bound
(self)¶ Get or set the upper bound
Setting the upper bound (float) will also adjust the associated optlang variables associated with the reaction. Infeasible combinations, such as a upper bound lower than the current lower bound will update the other bound.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
bounds
(self)¶ Get or set the bounds directly from a tuple
Convenience method for setting upper and lower bounds in one line using a tuple of lower and upper bound. Invalid bounds will raise an AssertionError.
When using a HistoryManager context, this attribute can be set temporarily, reversed when the exiting the context.

property
flux
(self)¶ The flux value in the most recent solution.
Flux is the primal value of the corresponding variable in the model.
Warning
Accessing reaction fluxes through a Solution object is the safer, preferred, and only guaranteed to be correct way. You can see how to do so easily in the examples.
Reaction flux is retrieved from the currently defined self._model.solver. The solver status is checked but there are no guarantees that the current solver state is the one you are looking for.
If you modify the underlying model after an optimization, you will retrieve the old optimization values.
 Raises
RuntimeError – If the underlying model was never optimized beforehand or the reaction is not part of a model.
OptimizationError – If the solver status is anything other than ‘optimal’.
AssertionError – If the flux value is not within the bounds.
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> solution = model.optimize() >>> model.reactions.PFK.flux 7.477381962160283 >>> solution.fluxes.PFK 7.4773819621602833

property
reduced_cost
(self)¶ The reduced cost in the most recent solution.
Reduced cost is the dual value of the corresponding variable in the model.
Warning
Accessing reduced costs through a Solution object is the safer, preferred, and only guaranteed to be correct way. You can see how to do so easily in the examples.
Reduced cost is retrieved from the currently defined self._model.solver. The solver status is checked but there are no guarantees that the current solver state is the one you are looking for.
If you modify the underlying model after an optimization, you will retrieve the old optimization values.
 Raises
RuntimeError – If the underlying model was never optimized beforehand or the reaction is not part of a model.
OptimizationError – If the solver status is anything other than ‘optimal’.
Examples
>>> import cobra.test >>> model = cobra.test.create_test_model("textbook") >>> solution = model.optimize() >>> model.reactions.PFK.reduced_cost 8.673617379884035e18 >>> solution.reduced_costs.PFK 8.6736173798840355e18

property
metabolites
(self)¶

property
genes
(self)¶

property
gene_reaction_rule
(self)¶

property
gene_name_reaction_rule
(self)¶ Display gene_reaction_rule with names intead.
Do NOT use this string for computation. It is intended to give a representation of the rule using more familiar gene names instead of the often cryptic ids.

property
functional
(self)¶ All required enzymes for reaction are functional.
 Returns
True if the geneproteinreaction (GPR) rule is fulfilled for this reaction, or if reaction is not associated to a model, otherwise False.
 Return type

property
x
(self)¶ The flux through the reaction in the most recent solution.
Flux values are computed from the primal values of the variables in the solution.

property
y
(self)¶ The reduced cost of the reaction in the most recent solution.
Reduced costs are computed from the dual values of the variables in the solution.

property
reversibility
(self)¶ Whether the reaction can proceed in both directions (reversible)
This is computed from the current upper and lower bounds.

property
boundary
(self)¶ Whether or not this reaction is an exchange reaction.
Returns True if the reaction has either no products or reactants.

property
model
(self)¶ returns the model the reaction is a part of

_update_awareness
(self)¶ Make sure all metabolites and genes that are associated with this reaction are aware of it.

remove_from_model
(self, remove_orphans=False)¶ Removes the reaction from a model.
This removes all associations between a reaction the associated model, metabolites and genes.
The change is reverted upon exit when using the model as a context.
 Parameters
remove_orphans (bool) – Remove orphaned genes and metabolites from the model as well

delete
(self, remove_orphans=False)¶ Removes the reaction from a model.
This removes all associations between a reaction the associated model, metabolites and genes.
The change is reverted upon exit when using the model as a context.
Deprecated, use reaction.remove_from_model instead.
 Parameters
remove_orphans (bool) – Remove orphaned genes and metabolites from the model as well

__setstate__
(self, state)¶ Probably not necessary to set _model as the cobra.Model that contains self sets the _model attribute for all metabolites and genes in the reaction.
However, to increase performance speed we do want to let the metabolite and gene know that they are employed in this reaction

copy
(self)¶ Copy a reaction
The referenced metabolites and genes are also copied.

__add__
(self, other)¶ Add two reactions
The stoichiometry will be the combined stoichiometry of the two reactions, and the gene reaction rule will be both rules combined by an and. All other attributes (i.e. reaction bounds) will match those of the first reaction

__iadd__
(self, other)¶

__sub__
(self, other)¶

__isub__
(self, other)¶

__imul__
(self, coefficient)¶ Scale coefficients in a reaction by a given value
E.g. A > B becomes 2A > 2B.
If coefficient is less than zero, the reaction is reversed and the bounds are swapped.

__mul__
(self, coefficient)¶

property
reactants
(self)¶ Return a list of reactants for the reaction.

property
products
(self)¶ Return a list of products for the reaction

get_coefficient
(self, metabolite_id)¶ Return the stoichiometric coefficient of a metabolite.
 Parameters
metabolite_id (str or cobra.Metabolite) –

get_coefficients
(self, metabolite_ids)¶ Return the stoichiometric coefficients for a list of metabolites.
 Parameters
metabolite_ids (iterable) – Containing
str
or ``cobra.Metabolite``s.

add_metabolites
(self, metabolites_to_add, combine=True, reversibly=True)¶ Add metabolites and stoichiometric coefficients to the reaction. If the final coefficient for a metabolite is 0 then it is removed from the reaction.
The change is reverted upon exit when using the model as a context.
 Parameters
metabolites_to_add (dict) – Dictionary with metabolite objects or metabolite identifiers as keys and coefficients as values. If keys are strings (name of a metabolite) the reaction must already be part of a model and a metabolite with the given name must exist in the model.
combine (bool) – Describes behavior a metabolite already exists in the reaction. True causes the coefficients to be added. False causes the coefficient to be replaced.
reversibly (bool) – Whether to add the change to the context to make the change reversibly or not (primarily intended for internal use).

subtract_metabolites
(self, metabolites, combine=True, reversibly=True)¶ Subtract metabolites from a reaction.
That means add the metabolites with 1*coefficient. If the final coefficient for a metabolite is 0 then the metabolite is removed from the reaction.
Notes
A final coefficient < 0 implies a reactant.
The change is reverted upon exit when using the model as a context.
 Parameters
metabolites (dict) – Dictionary where the keys are of class Metabolite and the values are the coefficients. These metabolites will be added to the reaction.
combine (bool) – Describes behavior a metabolite already exists in the reaction. True causes the coefficients to be added. False causes the coefficient to be replaced.
reversibly (bool) – Whether to add the change to the context to make the change reversibly or not (primarily intended for internal use).

property
reaction
(self)¶ Human readable reaction string

build_reaction_string
(self, use_metabolite_names=False)¶ Generate a human readable reaction string

check_mass_balance
(self)¶ Compute mass and charge balance for the reaction
returns a dict of {element: amount} for unbalanced elements. “charge” is treated as an element in this dict This should be empty for balanced reactions.

property
compartments
(self)¶ lists compartments the metabolites are in

get_compartments
(self)¶ lists compartments the metabolites are in

_associate_gene
(self, cobra_gene)¶ Associates a cobra.Gene object with a cobra.Reaction.
 Parameters
cobra_gene (cobra.core.Gene.Gene) –

_dissociate_gene
(self, cobra_gene)¶ Dissociates a cobra.Gene object with a cobra.Reaction.
 Parameters
cobra_gene (cobra.core.Gene.Gene) –

knock_out
(self)¶ Knockout reaction by setting its bounds to zero.

build_reaction_from_string
(self, reaction_str, verbose=True, fwd_arrow=None, rev_arrow=None, reversible_arrow=None, term_split='+')¶ Builds reaction from reaction equation reaction_str using parser
Takes a string and using the specifications supplied in the optional arguments infers a set of metabolites, metabolite compartments and stoichiometries for the reaction. It also infers the reversibility of the reaction from the reaction arrow.
Changes to the associated model are reverted upon exit when using the model as a context.
 Parameters
reaction_str (string) – a string containing a reaction formula (equation)
verbose (bool) – setting verbosity of function
fwd_arrow (re.compile) – for forward irreversible reaction arrows
rev_arrow (re.compile) – for backward irreversible reaction arrows
reversible_arrow (re.compile) – for reversible reaction arrows
term_split (string) – dividing individual metabolite entries

summary
(self, solution=None, fva=None)¶ Create a summary of the reaction flux.
 Parameters
solution (cobra.Solution, optional) – A previous model solution to use for generating the summary. If
None
, the summary method will generate a parsimonious flux distribution (default None).fva (pandas.DataFrame or float, optional) – Whether or not to include flux variability analysis in the output. If given, fva should either be a previous FVA solution matching the model or a float between 0 and 1 representing the fraction of the optimum objective to be searched (default None).
 Returns
 Return type
cobra.summary.ReactionSummary
See also

__str__
(self)¶ Return str(self).

_repr_html_
(self)¶

class
cobra.core.
Group
(id, name='', members=None, kind=None)[source]¶ Bases:
cobra.core.object.Object
Manage groups via this implementation of the SBML group specification.
Group is a class for holding information regarding a pathways, subsystems, or other custom groupings of objects within a cobra.Model object.
 Parameters
id (str) – The identifier to associate with this group
name (str, optional) – A human readable name for the group
members (iterable, optional) – A DictList containing references to cobra.Modelassociated objects that belong to the group.
kind ({"collection", "classification", "partonomy"}, optional) – The kind of group, as specified for the Groups feature in the SBML level 3 package specification. Can be any of “classification”, “partonomy”, or “collection”. The default is “collection”. Please consult the SBML level 3 package specification to ensure you are using the proper value for kind. In short, members of a “classification” group should have an “isa” relationship to the group (e.g. member isa polar compound, or member isa transporter). Members of a “partonomy” group should have a “partof” relationship (e.g. member is partof glycolysis). Members of a “collection” group do not have an implied relationship between the members, so use this value for kind when in doubt (e.g. member is a gapfilled reaction, or member is involved in a disease phenotype).

KIND_TYPES
= ['collection', 'classification', 'partonomy']¶

__len__
(self)¶

property
members
(self)¶

property
kind
(self)¶

class
cobra.core.
Solution
(objective_value, status, fluxes, reduced_costs=None, shadow_prices=None, **kwargs)[source]¶ Bases:
object
A unified interface to a cobra.Model optimization solution.
Notes
Solution is meant to be constructed by get_solution please look at that function to fully understand the Solution class.

fluxes
¶ Contains the reaction fluxes (primal values of variables).
 Type
pandas.Series

reduced_costs
¶ Contains reaction reduced costs (dual values of variables).
 Type
pandas.Series

shadow_prices
¶ Contains metabolite shadow prices (dual values of constraints).
 Type
pandas.Series

get_primal_by_id
¶

__repr__
(self)¶ String representation of the solution instance.

_repr_html_
(self)¶

__getitem__
(self, reaction_id)¶ Return the flux of a reaction.
 Parameters
reaction (str) – A model reaction ID.

to_frame
(self)¶ Return the fluxes and reduced costs as a data frame


class
cobra.core.
LegacySolution
(f, x=None, x_dict=None, y=None, y_dict=None, solver=None, the_time=0, status='NA', **kwargs)[source]¶ Bases:
object
Legacy support for an interface to a cobra.Model optimization solution.

x
¶ List or Array of the fluxes (primal values).
 Type
iterable

y
¶ List or Array of the dual values.
 Type
iterable
Warning
The LegacySolution class and its interface is deprecated.

__repr__
(self)¶ String representation of the solution instance.

__getitem__
(self, reaction_id)¶ Return the flux of a reaction.
 Parameters
reaction_id (str) – A reaction ID.

dress_results
(self, model)¶ Method could be intended as a decorator.
Warning
deprecated


cobra.core.
get_solution
(model, reactions=None, metabolites=None, raise_error=False)[source]¶ Generate a solution representation of the current solver state.
 Parameters
model (cobra.Model) – The model whose reactions to retrieve values for.
reactions (list, optional) – An iterable of cobra.Reaction objects. Uses model.reactions by default.
metabolites (list, optional) – An iterable of cobra.Metabolite objects. Uses model.metabolites by default.
raise_error (bool) – If true, raise an OptimizationError if solver status is not optimal.
 Returns
 Return type
Note
This is only intended for the optlang solver interfaces and not the legacy solvers.

class
cobra.core.
Species
(id=None, name=None)[source]¶ Bases:
cobra.core.object.Object
Species is a class for holding information regarding a chemical Species
 Parameters
id (string) – An identifier for the chemical species
name (string) – A human readable name.

property
reactions
(self)¶

__getstate__
(self)¶ Remove the references to container reactions when serializing to avoid problems associated with recursion.

copy
(self)¶ When copying a reaction, it is necessary to deepcopy the components so the list references aren’t carried over.
Additionally, a copy of a reaction is no longer in a cobra.Model.
This should be fixed with self.__deepcopy__ if possible

property
model
(self)¶
cobra.flux_analysis
¶
Submodules¶
cobra.flux_analysis.deletion
¶Access unique combinations of reactions in deletion results. 











Provide a common interface for single or multiple knockouts. 





Knock out each reaction from a given list. 

Knock out each gene from a given list. 

Knock out each reaction pair from the combinations of two given lists. 

Knock out each gene pair from the combination of two given lists. 

cobra.flux_analysis.deletion.
_multi_deletion
(model, entity, element_lists, method='fba', solution=None, processes=None, **kwargs)[source]¶ Provide a common interface for single or multiple knockouts.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
entity ('gene' or 'reaction') – The entity to knockout (
cobra.Gene
orcobra.Reaction
).element_lists (list) – List of iterables ``cobra.Reaction``s or ``cobra.Gene``s (or their IDs) to be deleted.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Passed on to underlying simulation functions.
 Returns
A representation of all combinations of entity deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The gene or reaction identifiers that were knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.deletion.
single_reaction_deletion
(model, reaction_list=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each reaction from a given list.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
reaction_list (iterable, optional) – ``cobra.Reaction``s to be deleted. If not passed, all the reactions from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all single reaction deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The reaction identifier that was knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.deletion.
single_gene_deletion
(model, gene_list=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each gene from a given list.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
gene_list (iterable) – ``cobra.Gene``s to be deleted. If not passed, all the genes from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all single gene deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The gene identifier that was knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.deletion.
double_reaction_deletion
(model, reaction_list1=None, reaction_list2=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each reaction pair from the combinations of two given lists.
We say ‘pair’ here but the order order does not matter.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
reaction_list1 (iterable, optional) – First iterable of ``cobra.Reaction``s to be deleted. If not passed, all the reactions from the model are used.
reaction_list2 (iterable, optional) – Second iterable of ``cobra.Reaction``s to be deleted. If not passed, all the reactions from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all combinations of reaction deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The reaction identifiers that were knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.deletion.
double_gene_deletion
(model, gene_list1=None, gene_list2=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each gene pair from the combination of two given lists.
We say ‘pair’ here but the order order does not matter.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
gene_list1 (iterable, optional) – First iterable of ``cobra.Gene``s to be deleted. If not passed, all the genes from the model are used.
gene_list2 (iterable, optional) – Second iterable of ``cobra.Gene``s to be deleted. If not passed, all the genes from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all combinations of gene deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The gene identifiers that were knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

class
cobra.flux_analysis.deletion.
KnockoutAccessor
(pandas_obj: pd.DataFrame)[source]¶ Access unique combinations of reactions in deletion results.
This allows acces in the form of results.knockout[rxn1] or results.knockout[“rxn1_id”]. Each individual entry will return a deletion so results.knockout[rxn1, rxn2] will return two deletions (for individual knockouts of rxn1 and rxn2 respectively). Multideletions can be accessed by passing in sets like results.knockou[{rxn1, rxn2}] which denotes the double deletion of both reactions. Thus, the following are allowed index elements:
single reactions or genes (depending on whether it is a gene or reaction deletion)
single reaction IDs or gene IDs
lists of single single reaction IDs or gene IDs (will return one row for each element in the list)
sets of reactions or genes (for multideletions)
sets of reactions IDs or gene IDs
list of sets of objects or IDs (to get several multideletions)

__getitem__
(self, args: Union[Gene, List[Gene], Set[Gene], List[Set[Gene]], Reaction, List[Reaction], Set[Reaction], List[Set[Reaction]], str, List[str], Set[str], List[Set[str]]]) → pd.DataFrame[source]¶ Return the deletion result for a particular set of knocked entities.
 argscobra.Reactions, cobra.Gene, str, set, or list
The deletions to be returned. Accepts:  single reactions or genes  single reaction IDs or gene IDs  lists of single single reaction IDs or gene IDs  sets of reactions or genes  sets of reactions IDs or gene IDs  list of sets of objects or IDs See the docs for usage examples.
 pd.DataFrame
The deletion result where the chosen entities have been deleted. Each row denotes a deletion.
cobra.flux_analysis.fastcc
¶Provide an implementation of FASTCC.

Perform the LP required for FASTCC. 

Flip the coefficients for optimizing in reverse direction. 

Check consistency of a metabolic network using FASTCC 1. 

cobra.flux_analysis.fastcc.
_find_sparse_mode
(model, rxns, flux_threshold, zero_cutoff)[source]¶ Perform the LP required for FASTCC.
 Parameters
model (cobra.core.Model) – The cobra model to perform FASTCC on.
rxns (list of cobra.core.Reactions) – The reactions to use for LP.
flux_threshold (float) – The upper threshold an auxiliary variable can have.
zero_cutoff (float) – The cutoff below which flux is considered zero.
 Returns
result – The list of reactions to consider as consistent.
 Return type

cobra.flux_analysis.fastcc.
_flip_coefficients
(model, rxns)[source]¶ Flip the coefficients for optimizing in reverse direction.

cobra.flux_analysis.fastcc.
fastcc
(model, flux_threshold=1.0, zero_cutoff=None)[source]¶ Check consistency of a metabolic network using FASTCC 1.
FASTCC (Fast Consistency Check) is an algorithm for rapid and efficient consistency check in metabolic networks. FASTCC is a pure LP implementation and is low on computation resource demand. FASTCC also circumvents the problem associated with reversible reactions for the purpose. Given a global model, it will generate a consistent global model i.e., remove blocked reactions. For more details on FASTCC, please check 1.
 Parameters
model (cobra.Model) – The constraintbased model to operate on.
flux_threshold (float, optional (default 1.0)) – The flux threshold to consider.
zero_cutoff (float, optional) – The cutoff to consider for zero flux (default model.tolerance).
 Returns
The consistent constraintbased model.
 Return type
Notes
The LP used for FASTCC is like so: maximize: sum_{i in J} z_i s.t. : z_i in [0, varepsilon] forall i in J, z_i in mathbb{R}_+
v_i ge z_i forall i in J Sv = 0 v in B
References
cobra.flux_analysis.gapfilling
¶

Perform gapfilling on a model. 

class
cobra.flux_analysis.gapfilling.
GapFiller
(model, universal=None, lower_bound=0.05, penalties=None, exchange_reactions=False, demand_reactions=True, integer_threshold=1e06)[source]¶ Bases:
object
Class for performing gap filling.
This class implements gap filling based on a mixedinteger approach, very similar to that described in 1 and the ‘nogrowth but growth’ part of [2]_ but with minor adjustments. In short, we add indicator variables for using the reactions in the universal model, z_i and then solve problem
minimize sum_i c_i * z_i s.t. Sv = 0
v_o >= t lb_i <= v_i <= ub_i v_i = 0 if z_i = 0
where lb, ub are the upper, lower flux bounds for reaction i, c_i is a cost parameter and the objective v_o is greater than the lower bound t. The default costs are 1 for reactions from the universal model, 100 for exchange (uptake) reactions added and 1 for added demand reactions.
Note that this is a mixedinteger linear program and as such will expensive to solve for large models. Consider using alternatives [3]_ such as CORDA instead [4,5]_.
 Parameters
model (cobra.Model) – The model to perform gap filling on.
universal (cobra.Model) – A universal model with reactions that can be used to complete the model.
lower_bound (float) – The minimally accepted flux for the objective in the filled model.
penalties (dict, None) – A dictionary with keys being ‘universal’ (all reactions included in the universal model), ‘exchange’ and ‘demand’ (all additionally added exchange and demand reactions) for the three reaction types. Can also have reaction identifiers for reaction specific costs. Defaults are 1, 100 and 1 respectively.
integer_threshold (float) – The threshold at which a value is considered nonzero (aka integrality threshold). If gapfilled models fail to validate, you may want to lower this value.
exchange_reactions (bool) – Consider adding exchange (uptake) reactions for all metabolites in the model.
demand_reactions (bool) – Consider adding demand reactions for all metabolites.
References
 1
Reed, Jennifer L., Trina R. Patel, Keri H. Chen, Andrew R. Joyce, Margaret K. Applebee, Christopher D. Herring, Olivia T. Bui, Eric M. Knight, Stephen S. Fong, and Bernhard O. Palsson. “Systems Approach to Refining Genome Annotation.” Proceedings of the National Academy of Sciences 103, no. 46 (2006): 17480–17484.
[2] Kumar, Vinay Satish, and Costas D. Maranas. “GrowMatch: An Automated Method for Reconciling In Silico/In Vivo Growth Predictions.” Edited by Christos A. Ouzounis. PLoS Computational Biology 5, no. 3 (March 13, 2009): e1000308. doi:10.1371/journal.pcbi.1000308.
[3] http://opencobra.github.io/cobrapy/tags/gapfilling/
[4] Schultz, André, and Amina A. Qutub. “Reconstruction of TissueSpecific Metabolic Networks Using CORDA.” Edited by Costas D. Maranas. PLOS Computational Biology 12, no. 3 (March 4, 2016): e1004808. doi:10.1371/journal.pcbi.1004808.
[5] Diener, Christian https://github.com/cdiener/corda

extend_model
(self, exchange_reactions=False, demand_reactions=True)[source]¶ Extend gapfilling model.
Add reactions from universal model and optionally exchange and demand reactions for all metabolites in the model to perform gapfilling on.

update_costs
(self)[source]¶ Update the coefficients for the indicator variables in the objective.
Done incrementally so that second time the function is called, active indicators in the current solutions gets higher cost than the unused indicators.

add_switches_and_objective
(self)[source]¶ Update gapfilling model with switches and the indicator objective.

fill
(self, iterations=1)[source]¶ Perform the gapfilling by iteratively solving the model, updating the costs and recording the used reactions.
 Parameters
iterations (int) – The number of rounds of gapfilling to perform. For every iteration, the penalty for every used reaction increases linearly. This way, the algorithm is encouraged to search for alternative solutions which may include previously used reactions. I.e., with enough iterations pathways including 10 steps will eventually be reported even if the shortest pathway is a single reaction.
 Returns
A list of lists where each element is a list reactions that were used to gapfill the model.
 Return type
iterable
 Raises
RuntimeError – If the model fails to be validated (i.e. the original model with the proposed reactions added, still cannot get the required flux through the objective).

cobra.flux_analysis.gapfilling.
gapfill
(model, universal=None, lower_bound=0.05, penalties=None, demand_reactions=True, exchange_reactions=False, iterations=1)[source]¶ Perform gapfilling on a model.
See documentation for the class GapFiller.
 modelcobra.Model
The model to perform gap filling on.
 universalcobra.Model, None
A universal model with reactions that can be used to complete the model. Only gapfill considering demand and exchange reactions if left missing.
 lower_boundfloat
The minimally accepted flux for the objective in the filled model.
 penaltiesdict, None
A dictionary with keys being ‘universal’ (all reactions included in the universal model), ‘exchange’ and ‘demand’ (all additionally added exchange and demand reactions) for the three reaction types. Can also have reaction identifiers for reaction specific costs. Defaults are 1, 100 and 1 respectively.
 iterationsint
The number of rounds of gapfilling to perform. For every iteration, the penalty for every used reaction increases linearly. This way, the algorithm is encouraged to search for alternative solutions which may include previously used reactions. I.e., with enough iterations pathways including 10 steps will eventually be reported even if the shortest pathway is a single reaction.
 exchange_reactionsbool
Consider adding exchange (uptake) reactions for all metabolites in the model.
 demand_reactionsbool
Consider adding demand reactions for all metabolites.
 iterable
list of lists with on set of reactions that completes the model per requested iteration.
 import cobra as ct
>>> from cobra import Model >>> from cobra.flux_analysis import gapfill >>> model = ct.create_test_model("salmonella") >>> universal = Model('universal') >>> universal.add_reactions(model.reactions.GF6PTA.copy()) >>> model.remove_reactions([model.reactions.GF6PTA]) >>> gapfill(model, universal)
cobra.flux_analysis.geometric
¶Provide an implementation of geometric FBA.

Perform geometric FBA to obtain a unique, centered flux distribution. 

cobra.flux_analysis.geometric.
geometric_fba
(model, epsilon=1e06, max_tries=200, processes=None)[source]¶ Perform geometric FBA to obtain a unique, centered flux distribution.
Geometric FBA 1 formulates the problem as a polyhedron and then solves it by bounding the convex hull of the polyhedron. The bounding forms a box around the convex hull which reduces with every iteration and extracts a unique solution in this way.
 Parameters
model (cobra.Model) – The model to perform geometric FBA on.
epsilon (float, optional) – The convergence tolerance of the model (default 1E06).
max_tries (int, optional) – Maximum number of iterations (default 200).
processes (int, optional) – The number of parallel processes to run. If not explicitly passed, will be set from the global configuration singleton.
 Returns
The solution object containing all the constraints required for geometric FBA.
 Return type
References
 1
Smallbone, Kieran & Simeonidis, Vangelis. (2009). Flux balance analysis: A geometric perspective. Journal of theoretical biology.258. 3115. 10.1016/j.jtbi.2009.01.027.
cobra.flux_analysis.helpers
¶Helper functions for all flux analysis methods.
cobra.flux_analysis.loopless
¶Provides functions to remove thermodynamically infeasible loops.

Modify a model so all feasible flux distributions are loopless. 

Add constraints for CycleFreeFlux. 

Convert an existing solution to a loopless one. 

Plugin to get a loopless FVA solution from single FVA iteration. 

cobra.flux_analysis.loopless.
add_loopless
(model, zero_cutoff=None)[source]¶ Modify a model so all feasible flux distributions are loopless.
In most cases you probably want to use the much faster loopless_solution. May be used in cases where you want to add complex constraints and objecives (for instance quadratic objectives) to the model afterwards or use an approximation of Gibbs free energy directions in you model. Adds variables and constraints to a model which will disallow flux distributions with loops. The used formulation is described in [1]_. This function will modify your model.
 Parameters
model (cobra.Model) – The model to which to add the constraints.
zero_cutoff (positive float, optional) – Cutoff used for null space. Coefficients with an absolute value smaller than zero_cutoff are considered to be zero (default model.tolerance).
 Returns
 Return type
Nothing
References
 1
Elimination of thermodynamically infeasible loops in steadystate metabolic models. Schellenberger J, Lewis NE, Palsson BO. Biophys J. 2011 Feb 2;100(3):54453. doi: 10.1016/j.bpj.2010.12.3707. Erratum in: Biophys J. 2011 Mar 2;100(5):1381.

cobra.flux_analysis.loopless.
_add_cycle_free
(model, fluxes)[source]¶ Add constraints for CycleFreeFlux.

cobra.flux_analysis.loopless.
loopless_solution
(model, fluxes=None)[source]¶ Convert an existing solution to a loopless one.
Removes as many loops as possible (see Notes). Uses the method from CycleFreeFlux [1]_ and is much faster than add_loopless and should therefore be the preferred option to get loopless flux distributions.
 Parameters
model (cobra.Model) – The model to which to add the constraints.
fluxes (dict) – A dictionary {rxn_id: flux} that assigns a flux to each reaction. If not None will use the provided flux values to obtain a close loopless solution.
 Returns
A solution object containing the fluxes with the least amount of loops possible or None if the optimization failed (usually happening if the flux distribution in fluxes is infeasible).
 Return type
Notes
The returned flux solution has the following properties:
it contains the minimal number of loops possible and no loops at all if all flux bounds include zero
it has an objective value close to the original one and the same objective value id the objective expression can not form a cycle (which is usually true since it consumes metabolites)
it has the same exact exchange fluxes as the previous solution
all fluxes have the same sign (flow in the same direction) as the previous solution
References
 1
CycleFreeFlux: efficient removal of thermodynamically infeasible loops from flux distributions. Desouki AA, Jarre F, GeliusDietrich G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):215965. doi: 10.1093/bioinformatics/btv096.

cobra.flux_analysis.loopless.
loopless_fva_iter
(model, reaction, solution=False, zero_cutoff=None)[source]¶ Plugin to get a loopless FVA solution from single FVA iteration.
Assumes the following about model and reaction: 1. the model objective is set to be reaction 2. the model has been optimized and contains the minimum/maximum flux for
reaction
the model contains an auxiliary variable called “fva_old_objective” denoting the previous objective
 Parameters
model (cobra.Model) – The model to be used.
reaction (cobra.Reaction) – The reaction currently minimized/maximized.
solution (boolean, optional) – Whether to return the entire solution or only the minimum/maximum for reaction.
zero_cutoff (positive float, optional) – Cutoff used for loop removal. Fluxes with an absolute value smaller than zero_cutoff are considered to be zero (default model.tolerance).
 Returns
Returns the minimized/maximized flux through reaction if all_fluxes == False (default). Otherwise returns a loopless flux solution containing the minimum/maximum flux for reaction.
 Return type
single float or dict
cobra.flux_analysis.moma
¶Provide minimization of metabolic adjustment (MOMA).

Compute a single solution based on (linear) MOMA. 

Add constraints and objective representing for MOMA. 

cobra.flux_analysis.moma.
moma
(model, solution=None, linear=True)[source]¶ Compute a single solution based on (linear) MOMA.
Compute a new flux distribution that is at a minimal distance to a previous reference solution. Minimization of metabolic adjustment (MOMA) is generally used to assess the impact of knockouts. Thus the typical usage is to provide a wildtype flux distribution as reference and a model in knockout state.
 Parameters
model (cobra.Model) – The model state to compute a MOMAbased solution for.
solution (cobra.Solution, optional) – A (wildtype) reference solution.
linear (bool, optional) – Whether to use the linear MOMA formulation or not (default True).
 Returns
A flux distribution that is at a minimal distance compared to the reference solution.
 Return type
See also
add_moma()
add MOMA constraints and objective

cobra.flux_analysis.moma.
add_moma
(model, solution=None, linear=True)[source]¶ Add constraints and objective representing for MOMA.
This adds variables and constraints for the minimization of metabolic adjustment (MOMA) to the model.
 Parameters
model (cobra.Model) – The model to add MOMA constraints and objective to.
solution (cobra.Solution, optional) – A previous solution to use as a reference. If no solution is given, one will be computed using pFBA.
linear (bool, optional) – Whether to use the linear MOMA formulation or not (default True).
Notes
In the original MOMA 1 specification one looks for the flux distribution of the deletion (v^d) closest to the fluxes without the deletion (v). In math this means:
minimize sum_i (v^d_i  v_i)^2 s.t. Sv^d = 0
lb_i <= v^d_i <= ub_i
Here, we use a variable transformation v^t := v^d_i  v_i. Substituting and using the fact that Sv = 0 gives:
minimize sum_i (v^t_i)^2 s.t. Sv^d = 0
v^t = v^d_i  v_i lb_i <= v^d_i <= ub_i
So basically we just recenter the flux space at the old solution and then find the flux distribution closest to the new zero (center). This is the same strategy as used in cameo.
In the case of linear MOMA 2, we instead minimize sum_i abs(v^t_i). The linear MOMA is typically significantly faster. Also quadratic MOMA tends to give flux distributions in which all fluxes deviate from the reference fluxes a little bit whereas linear MOMA tends to give flux distributions where the majority of fluxes are the same reference with few fluxes deviating a lot (typical effect of L2 norm vs L1 norm).
The former objective function is saved in the optlang solver interface as
"moma_old_objective"
and this can be used to immediately extract the value of the former objective after MOMA optimization.See also
pfba()
parsimonious FBA
References
 1
Segrè, Daniel, Dennis Vitkup, and George M. Church. “Analysis of Optimality in Natural and Perturbed Metabolic Networks.” Proceedings of the National Academy of Sciences 99, no. 23 (November 12, 2002): 15112. https://doi.org/10.1073/pnas.232349399.
 2
Becker, Scott A, Adam M Feist, Monica L Mo, Gregory Hannum, Bernhard Ø Palsson, and Markus J Herrgard. “Quantitative Prediction of Cellular Metabolism with ConstraintBased Models: The COBRA Toolbox.” Nature Protocols 2 (March 29, 2007): 727.
cobra.flux_analysis.parsimonious
¶



Perform basic pFBA (parsimonious Enzyme Usage Flux Balance Analysis) 

Add pFBA objective 

cobra.flux_analysis.parsimonious.
pfba
(model, fraction_of_optimum=1.0, objective=None, reactions=None)[source]¶ Perform basic pFBA (parsimonious Enzyme Usage Flux Balance Analysis) to minimize total flux.
pFBA [1] adds the minimization of all fluxes the the objective of the model. This approach is motivated by the idea that high fluxes have a higher enzyme turnover and that since producing enzymes is costly, the cell will try to minimize overall flux while still maximizing the original objective function, e.g. the growth rate.
 Parameters
model (cobra.Model) – The model
fraction_of_optimum (float, optional) – Fraction of optimum which must be maintained. The original objective reaction is constrained to be greater than maximal_value * fraction_of_optimum.
objective (dict or model.problem.Objective) – A desired objective to use during optimization in addition to the pFBA objective. Dictionaries (reaction as key, coefficient as value) can be used for linear objectives.
reactions (iterable) – List of reactions or reaction identifiers. Implies return_frame to be true. Only return fluxes for the given reactions. Faster than fetching all fluxes if only a few are needed.
 Returns
The solution object to the optimized model with pFBA constraints added.
 Return type
References
 1
Lewis, N. E., Hixson, K. K., Conrad, T. M., Lerman, J. A., Charusanti, P., Polpitiya, A. D., Palsson, B. O. (2010). Omic data from evolved E. coli are consistent with computed optimal growth from genomescale models. Molecular Systems Biology, 6, 390. doi:10.1038/msb.2010.47

cobra.flux_analysis.parsimonious.
add_pfba
(model, objective=None, fraction_of_optimum=1.0)[source]¶ Add pFBA objective
Add objective to minimize the summed flux of all reactions to the current objective.
See also
 Parameters
model (cobra.Model) – The model to add the objective to
objective – An objective to set in combination with the pFBA objective.
fraction_of_optimum (float) – Fraction of optimum which must be maintained. The original objective reaction is constrained to be greater than maximal_value * fraction_of_optimum.
cobra.flux_analysis.phenotype_phase_plane
¶

Calculate the objective value conditioned on all combinations of 



Compute total output per input unit. 

Split metabolites into the atoms times their stoichiometric coefficients. 

Return the metabolite weight times its stoichiometric coefficient. 

Compute the total components consumption or production flux. 

Find all active carbon source reactions. 

cobra.flux_analysis.phenotype_phase_plane.
production_envelope
(model, reactions, objective=None, carbon_sources=None, points=20, threshold=None)[source]¶ Calculate the objective value conditioned on all combinations of fluxes for a set of chosen reactions
The production envelope can be used to analyze a model’s ability to produce a given compound conditional on the fluxes for another set of reactions, such as the uptake rates. The model is alternately optimized with respect to minimizing and maximizing the objective and the obtained fluxes are recorded. Ranges to compute production is set to the effective bounds, i.e., the minimum / maximum fluxes that can be obtained given current reaction bounds.
 Parameters
model (cobra.Model) – The model to compute the production envelope for.
reactions (list or string) – A list of reactions, reaction identifiers or a single reaction.
objective (string, dict, model.solver.interface.Objective, optional) – The objective (reaction) to use for the production envelope. Use the model’s current objective if left missing.
carbon_sources (list or string, optional) – One or more reactions or reaction identifiers that are the source of carbon for computing carbon (mol carbon in output over mol carbon in input) and mass yield (gram product over gram output). Only objectives with a carbon containing input and output metabolite is supported. Will identify active carbon sources in the medium if none are specified.
points (int, optional) – The number of points to calculate production for.
threshold (float, optional) – A cutoff under which flux values will be considered to be zero (default model.tolerance).
 Returns
A data frame with one row per evaluated point and
reaction id : one column per input reaction indicating the flux at each given point,
carbon_source: identifiers of carbon exchange reactions
A column for the maximum and minimum each for the following types:
flux: the objective flux
carbon_yield: if carbon source is defined and the product is a single metabolite (mol carbon product per mol carbon feeding source)
mass_yield: if carbon source is defined and the product is a single metabolite (gram product per 1 g of feeding source)
 Return type
pandas.DataFrame
Examples
>>> import cobra.test >>> from cobra.flux_analysis import production_envelope >>> model = cobra.test.create_test_model("textbook") >>> production_envelope(model, ["EX_glc__D_e", "EX_o2_e"])

cobra.flux_analysis.phenotype_phase_plane.
add_envelope
(model, reactions, grid, c_input, c_output, threshold)[source]¶

cobra.flux_analysis.phenotype_phase_plane.
total_yield
(input_fluxes, input_elements, output_flux, output_elements)[source]¶ Compute total output per input unit.
Units are typically mol carbon atoms or gram of source and product.
 Parameters
input_fluxes (list) – A list of input reaction fluxes in the same order as the
input_components
.input_elements (list) – A list of reaction components which are in turn list of numbers.
output_flux (float) – The output flux value.
output_elements (list) – A list of stoichiometrically weighted output reaction components.
 Returns
The ratio between output (mol carbon atoms or grams of product) and input (mol carbon atoms or grams of source compounds).
 Return type

cobra.flux_analysis.phenotype_phase_plane.
reaction_elements
(reaction)[source]¶ Split metabolites into the atoms times their stoichiometric coefficients.

cobra.flux_analysis.phenotype_phase_plane.
reaction_weight
(reaction)[source]¶ Return the metabolite weight times its stoichiometric coefficient.
cobra.flux_analysis.reaction
¶functions for analyzing / creating objective functions

Assesses production capacity. 

Assesses the ability of the model to provide sufficient precursors, 



Assesses the ability of the model to provide sufficient precursors for 

Assesses whether the model has the capacity to absorb the products of 

cobra.flux_analysis.reaction.
assess
(model, reaction, flux_coefficient_cutoff=0.001, solver=None)[source]¶ Assesses production capacity.
Assesses the capacity of the model to produce the precursors for the reaction and absorb the production of the reaction while the reaction is operating at, or above, the specified cutoff.
 Parameters
model (cobra.Model) – The cobra model to assess production capacity for
reaction (reaction identifier or cobra.Reaction) – The reaction to assess
flux_coefficient_cutoff (float) – The minimum flux that reaction must carry to be considered active.
solver (basestring) – Solver name. If None, the default solver will be used.
 Returns
True if the model can produce the precursors and absorb the products for the reaction operating at, or above, flux_coefficient_cutoff. Otherwise, a dictionary of {‘precursor’: Status, ‘product’: Status}. Where Status is the results from assess_precursors and assess_products, respectively.
 Return type

cobra.flux_analysis.reaction.
assess_component
(model, reaction, side, flux_coefficient_cutoff=0.001, solver=None)[source]¶ Assesses the ability of the model to provide sufficient precursors, or absorb products, for a reaction operating at, or beyond, the specified cutoff.
 Parameters
model (cobra.Model) – The cobra model to assess production capacity for
reaction (reaction identifier or cobra.Reaction) – The reaction to assess
side (basestring) – Side of the reaction, ‘products’ or ‘reactants’
flux_coefficient_cutoff (float) – The minimum flux that reaction must carry to be considered active.
solver (basestring) – Solver name. If None, the default solver will be used.
 Returns
True if the precursors can be simultaneously produced at the specified cutoff. False, if the model has the capacity to produce each individual precursor at the specified threshold but not all precursors at the required level simultaneously. Otherwise a dictionary of the required and the produced fluxes for each reactant that is not produced in sufficient quantities.
 Return type

cobra.flux_analysis.reaction.
assess_precursors
(model, reaction, flux_coefficient_cutoff=0.001, solver=None)[source]¶ Assesses the ability of the model to provide sufficient precursors for a reaction operating at, or beyond, the specified cutoff.
Deprecated: use assess_component instead
 Parameters
model (cobra.Model) – The cobra model to assess production capacity for
reaction (reaction identifier or cobra.Reaction) – The reaction to assess
flux_coefficient_cutoff (float) – The minimum flux that reaction must carry to be considered active.
solver (basestring) – Solver name. If None, the default solver will be used.
 Returns
True if the precursors can be simultaneously produced at the specified cutoff. False, if the model has the capacity to produce each individual precursor at the specified threshold but not all precursors at the required level simultaneously. Otherwise a dictionary of the required and the produced fluxes for each reactant that is not produced in sufficient quantities.
 Return type

cobra.flux_analysis.reaction.
assess_products
(model, reaction, flux_coefficient_cutoff=0.001, solver=None)[source]¶ Assesses whether the model has the capacity to absorb the products of a reaction at a given flux rate.
Useful for identifying which components might be blocking a reaction from achieving a specific flux rate.
Deprecated: use assess_component instead
 Parameters
model (cobra.Model) – The cobra model to assess production capacity for
reaction (reaction identifier or cobra.Reaction) – The reaction to assess
flux_coefficient_cutoff (float) – The minimum flux that reaction must carry to be considered active.
solver (basestring) – Solver name. If None, the default solver will be used.
 Returns
True if the model has the capacity to absorb all the reaction products being simultaneously given the specified cutoff. False, if the model has the capacity to absorb each individual product but not all products at the required level simultaneously. Otherwise a dictionary of the required and the capacity fluxes for each product that is not absorbed in sufficient quantities.
 Return type
cobra.flux_analysis.room
¶Provide regulatory on/off minimization (ROOM).

Compute a single solution based on regulatory on/off minimization (ROOM). 

Add constraints and objective for ROOM. 

cobra.flux_analysis.room.
room
(model, solution=None, linear=False, delta=0.03, epsilon=0.001)[source]¶ Compute a single solution based on regulatory on/off minimization (ROOM).
Compute a new flux distribution that minimizes the number of active reactions needed to accommodate a previous reference solution. Regulatory on/off minimization (ROOM) is generally used to assess the impact of knockouts. Thus the typical usage is to provide a wildtype flux distribution as reference and a model in knockout state.
 Parameters
model (cobra.Model) – The model state to compute a ROOMbased solution for.
solution (cobra.Solution, optional) – A (wildtype) reference solution.
linear (bool, optional) – Whether to use the linear ROOM formulation or not (default False).
delta (float, optional) – The relative tolerance range (additive) (default 0.03).
epsilon (float, optional) – The absolute tolerance range (multiplicative) (default 0.001).
 Returns
A flux distribution with minimal active reaction changes compared to the reference.
 Return type
See also
add_room()
add ROOM constraints and objective

cobra.flux_analysis.room.
add_room
(model, solution=None, linear=False, delta=0.03, epsilon=0.001)[source]¶ Add constraints and objective for ROOM.
This function adds variables and constraints for applying regulatory on/off minimization (ROOM) to the model.
 Parameters
model (cobra.Model) – The model to add ROOM constraints and objective to.
solution (cobra.Solution, optional) – A previous solution to use as a reference. If no solution is given, one will be computed using pFBA.
linear (bool, optional) – Whether to use the linear ROOM formulation or not (default False).
delta (float, optional) – The relative tolerance range which is additive in nature (default 0.03).
epsilon (float, optional) – The absolute range of tolerance which is multiplicative (default 0.001).
Notes
The formulation used here is the same as stated in the original paper 1. The mathematical expression is given below:
minimize sum_{i=1}^m y^i s.t. Sv = 0
v_min <= v <= v_max v_j = 0 j ∈ A for 1 <= i <= m v_i  y_i(v_{max,i}  w_i^u) <= w_i^u (1) v_i  y_i(v_{min,i}  w_i^l) <= w_i^l (2) y_i ∈ {0,1} (3) w_i^u = w_i + deltaw_i + epsilon w_i^l = w_i  deltaw_i  epsilon
So, for the linear version of the ROOM , constraint (3) is relaxed to 0 <= y_i <= 1.
See also
pfba()
parsimonious FBA
References
 1
Tomer Shlomi, Omer Berkman and Eytan Ruppin, “Regulatory on/off minimization of metabolic flux changes after genetic perturbations”, PNAS 2005 102 (21) 76957700; doi:10.1073/pnas.0406346102
cobra.flux_analysis.variability
¶

Initialize a global model object for multiprocessing. 



Determine the minimum and maximum possible flux value for each reaction. 

Find reactions that cannot carry any flux. 

Return a set of essential genes. 

Return a set of essential reactions. 

cobra.flux_analysis.variability.
_init_worker
(model, loopless, sense)[source]¶ Initialize a global model object for multiprocessing.

cobra.flux_analysis.variability.
flux_variability_analysis
(model, reaction_list=None, loopless=False, fraction_of_optimum=1.0, pfba_factor=None, processes=None)[source]¶ Determine the minimum and maximum possible flux value for each reaction.
 Parameters
model (cobra.Model) – The model for which to run the analysis. It will not be modified.
reaction_list (list of cobra.Reaction or str, optional) – The reactions for which to obtain min/max fluxes. If None will use all reactions in the model (default).
loopless (boolean, optional) – Whether to return only loopless solutions. This is significantly slower. Please also refer to the notes.
fraction_of_optimum (float, optional) – Must be <= 1.0. Requires that the objective value is at least the fraction times maximum objective value. A value of 0.85 for instance means that the objective has to be at least at 85% percent of its maximum.
pfba_factor (float, optional) – Add an additional constraint to the model that requires the total sum of absolute fluxes must not be larger than this value times the smallest possible sum of absolute fluxes, i.e., by setting the value to 1.1 the total sum of absolute fluxes must not be more than 10% larger than the pFBA solution. Since the pFBA solution is the one that optimally minimizes the total flux sum, the
pfba_factor
should, if set, be larger than one. Setting this value may lead to more realistic predictions of the effective flux bounds.processes (int, optional) – The number of parallel processes to run. If not explicitly passed, will be set from the global configuration singleton.
 Returns
A data frame with reaction identifiers as the index and two columns:  maximum: indicating the highest possible flux  minimum: indicating the lowest possible flux
 Return type
pandas.DataFrame
Notes
This implements the fast version as described in 1. Please note that the flux distribution containing all minimal/maximal fluxes does not have to be a feasible solution for the model. Fluxes are minimized/maximized individually and a single minimal flux might require all others to be suboptimal.
Using the loopless option will lead to a significant increase in computation time (about a factor of 100 for large models). However, the algorithm used here (see 2) is still more than 1000x faster than the “naive” version using
add_loopless(model)
. Also note that if you have included constraints that force a loop (for instance by setting all fluxes in a loop to be nonzero) this loop will be included in the solution.References
 1
Computationally efficient flux variability analysis. Gudmundsson S, Thiele I. BMC Bioinformatics. 2010 Sep 29;11:489. doi: 10.1186/1471210511489, PMID: 20920235
 2
CycleFreeFlux: efficient removal of thermodynamically infeasible loops from flux distributions. Desouki AA, Jarre F, GeliusDietrich G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):215965. doi: 10.1093/bioinformatics/btv096.

cobra.flux_analysis.variability.
find_blocked_reactions
(model, reaction_list=None, zero_cutoff=None, open_exchanges=False, processes=None)[source]¶ Find reactions that cannot carry any flux.
The question whether or not a reaction is blocked is highly dependent on the current exchange reaction settings for a COBRA model. Hence an argument is provided to open all exchange reactions.
Notes
Sink and demand reactions are left untouched. Please modify them manually.
 Parameters
model (cobra.Model) – The model to analyze.
reaction_list (list, optional) – List of reactions to consider, the default includes all model reactions.
zero_cutoff (float, optional) – Flux value which is considered to effectively be zero (default model.tolerance).
open_exchanges (bool, optional) – Whether or not to open all exchange reactions to very high flux ranges.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of reactions is large. If not explicitly passed, it will be set from the global configuration singleton.
 Returns
List with the identifiers of blocked reactions.
 Return type

cobra.flux_analysis.variability.
find_essential_genes
(model, threshold=None, processes=None)[source]¶ Return a set of essential genes.
A gene is considered essential if restricting the flux of all reactions that depend on it to zero causes the objective, e.g., the growth rate, to also be zero, below the threshold, or infeasible.
 Parameters
model (cobra.Model) – The model to find the essential genes for.
threshold (float, optional) – Minimal objective flux to be considered viable. By default this is 1% of the maximal objective.
processes (int, optional) – The number of parallel processes to run. If not passed, will be set to the number of CPUs found.
processes – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not explicitly passed, it will be set from the global configuration singleton.
 Returns
Set of essential genes
 Return type

cobra.flux_analysis.variability.
find_essential_reactions
(model, threshold=None, processes=None)[source]¶ Return a set of essential reactions.
A reaction is considered essential if restricting its flux to zero causes the objective, e.g., the growth rate, to also be zero, below the threshold, or infeasible.
 Parameters
model (cobra.Model) – The model to find the essential reactions for.
threshold (float, optional) – Minimal objective flux to be considered viable. By default this is 1% of the maximal objective.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not explicitly passed, it will be set from the global configuration singleton.
 Returns
Set of essential reactions
 Return type
Package Contents¶

Knock out each gene pair from the combination of two given lists. 

Knock out each reaction pair from the combinations of two given lists. 

Knock out each gene from a given list. 

Knock out each reaction from a given list. 

Check consistency of a metabolic network using FASTCC [1]_. 

Perform gapfilling on a model. 

Perform geometric FBA to obtain a unique, centered flux distribution. 

Convert an existing solution to a loopless one. 

Modify a model so all feasible flux distributions are loopless. 

Add constraints and objective representing for MOMA. 

Compute a single solution based on (linear) MOMA. 

Perform basic pFBA (parsimonious Enzyme Usage Flux Balance Analysis) 

Find reactions that cannot carry any flux. 

Return a set of essential genes. 

Return a set of essential reactions. 

Determine the minimum and maximum possible flux value for each reaction. 

Add constraints and objective for ROOM. 

Compute a single solution based on regulatory on/off minimization (ROOM). 

cobra.flux_analysis.
double_gene_deletion
(model, gene_list1=None, gene_list2=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each gene pair from the combination of two given lists.
We say ‘pair’ here but the order order does not matter.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
gene_list1 (iterable, optional) – First iterable of ``cobra.Gene``s to be deleted. If not passed, all the genes from the model are used.
gene_list2 (iterable, optional) – Second iterable of ``cobra.Gene``s to be deleted. If not passed, all the genes from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all combinations of gene deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The gene identifiers that were knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.
double_reaction_deletion
(model, reaction_list1=None, reaction_list2=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each reaction pair from the combinations of two given lists.
We say ‘pair’ here but the order order does not matter.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
reaction_list1 (iterable, optional) – First iterable of ``cobra.Reaction``s to be deleted. If not passed, all the reactions from the model are used.
reaction_list2 (iterable, optional) – Second iterable of ``cobra.Reaction``s to be deleted. If not passed, all the reactions from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all combinations of reaction deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The reaction identifiers that were knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.
single_gene_deletion
(model, gene_list=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each gene from a given list.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
gene_list (iterable) – ``cobra.Gene``s to be deleted. If not passed, all the genes from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all single gene deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The gene identifier that was knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.
single_reaction_deletion
(model, reaction_list=None, method='fba', solution=None, processes=None, **kwargs)[source]¶ Knock out each reaction from a given list.
 Parameters
model (cobra.Model) – The metabolic model to perform deletions in.
reaction_list (iterable, optional) – ``cobra.Reaction``s to be deleted. If not passed, all the reactions from the model are used.
method ({"fba", "moma", "linear moma", "room", "linear room"}, optional) – Method used to predict the growth rate.
solution (cobra.Solution, optional) – A previous solution to use as a reference for (linear) MOMA or ROOM.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not passed, will be set to the number of CPUs found.
kwargs – Keyword arguments are passed on to underlying simulation functions such as
add_room
.
 Returns
A representation of all single reaction deletions. The columns are ‘growth’ and ‘status’, where
 indextuple(str)
The reaction identifier that was knocked out.
 growthfloat
The growth rate of the adjusted model.
 statusstr
The solution’s status.
 Return type
pandas.DataFrame

cobra.flux_analysis.
fastcc
(model, flux_threshold=1.0, zero_cutoff=None)[source]¶ Check consistency of a metabolic network using FASTCC [1]_.
FASTCC (Fast Consistency Check) is an algorithm for rapid and efficient consistency check in metabolic networks. FASTCC is a pure LP implementation and is low on computation resource demand. FASTCC also circumvents the problem associated with reversible reactions for the purpose. Given a global model, it will generate a consistent global model i.e., remove blocked reactions. For more details on FASTCC, please check [1]_.
 Parameters
model (cobra.Model) – The constraintbased model to operate on.
flux_threshold (float, optional (default 1.0)) – The flux threshold to consider.
zero_cutoff (float, optional) – The cutoff to consider for zero flux (default model.tolerance).
 Returns
The consistent constraintbased model.
 Return type
Notes
The LP used for FASTCC is like so: maximize: sum_{i in J} z_i s.t. : z_i in [0, varepsilon] forall i in J, z_i in mathbb{R}_+
v_i ge z_i forall i in J Sv = 0 v in B
References
 1
Vlassis N, Pacheco MP, Sauter T (2014) Fast Reconstruction of Compact ContextSpecific Metabolic Network Models. PLoS Comput Biol 10(1): e1003424. doi:10.1371/journal.pcbi.1003424

cobra.flux_analysis.
gapfill
(model, universal=None, lower_bound=0.05, penalties=None, demand_reactions=True, exchange_reactions=False, iterations=1)[source]¶ Perform gapfilling on a model.
See documentation for the class GapFiller.
 modelcobra.Model
The model to perform gap filling on.
 universalcobra.Model, None
A universal model with reactions that can be used to complete the model. Only gapfill considering demand and exchange reactions if left missing.
 lower_boundfloat
The minimally accepted flux for the objective in the filled model.
 penaltiesdict, None
A dictionary with keys being ‘universal’ (all reactions included in the universal model), ‘exchange’ and ‘demand’ (all additionally added exchange and demand reactions) for the three reaction types. Can also have reaction identifiers for reaction specific costs. Defaults are 1, 100 and 1 respectively.
 iterationsint
The number of rounds of gapfilling to perform. For every iteration, the penalty for every used reaction increases linearly. This way, the algorithm is encouraged to search for alternative solutions which may include previously used reactions. I.e., with enough iterations pathways including 10 steps will eventually be reported even if the shortest pathway is a single reaction.
 exchange_reactionsbool
Consider adding exchange (uptake) reactions for all metabolites in the model.
 demand_reactionsbool
Consider adding demand reactions for all metabolites.
 iterable
list of lists with on set of reactions that completes the model per requested iteration.
 import cobra as ct
>>> from cobra import Model >>> from cobra.flux_analysis import gapfill >>> model = ct.create_test_model("salmonella") >>> universal = Model('universal') >>> universal.add_reactions(model.reactions.GF6PTA.copy()) >>> model.remove_reactions([model.reactions.GF6PTA]) >>> gapfill(model, universal)

cobra.flux_analysis.
geometric_fba
(model, epsilon=1e06, max_tries=200, processes=None)[source]¶ Perform geometric FBA to obtain a unique, centered flux distribution.
Geometric FBA [1]_ formulates the problem as a polyhedron and then solves it by bounding the convex hull of the polyhedron. The bounding forms a box around the convex hull which reduces with every iteration and extracts a unique solution in this way.
 Parameters
model (cobra.Model) – The model to perform geometric FBA on.
epsilon (float, optional) – The convergence tolerance of the model (default 1E06).
max_tries (int, optional) – Maximum number of iterations (default 200).
processes (int, optional) – The number of parallel processes to run. If not explicitly passed, will be set from the global configuration singleton.
 Returns
The solution object containing all the constraints required for geometric FBA.
 Return type
References
 1
Smallbone, Kieran & Simeonidis, Vangelis. (2009). Flux balance analysis: A geometric perspective. Journal of theoretical biology.258. 3115. 10.1016/j.jtbi.2009.01.027.

cobra.flux_analysis.
loopless_solution
(model, fluxes=None)[source]¶ Convert an existing solution to a loopless one.
Removes as many loops as possible (see Notes). Uses the method from CycleFreeFlux [1]_ and is much faster than add_loopless and should therefore be the preferred option to get loopless flux distributions.
 Parameters
model (cobra.Model) – The model to which to add the constraints.
fluxes (dict) – A dictionary {rxn_id: flux} that assigns a flux to each reaction. If not None will use the provided flux values to obtain a close loopless solution.
 Returns
A solution object containing the fluxes with the least amount of loops possible or None if the optimization failed (usually happening if the flux distribution in fluxes is infeasible).
 Return type
Notes
The returned flux solution has the following properties:
it contains the minimal number of loops possible and no loops at all if all flux bounds include zero
it has an objective value close to the original one and the same objective value id the objective expression can not form a cycle (which is usually true since it consumes metabolites)
it has the same exact exchange fluxes as the previous solution
all fluxes have the same sign (flow in the same direction) as the previous solution
References
 1
CycleFreeFlux: efficient removal of thermodynamically infeasible loops from flux distributions. Desouki AA, Jarre F, GeliusDietrich G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):215965. doi: 10.1093/bioinformatics/btv096.

cobra.flux_analysis.
add_loopless
(model, zero_cutoff=None)[source]¶ Modify a model so all feasible flux distributions are loopless.
In most cases you probably want to use the much faster loopless_solution. May be used in cases where you want to add complex constraints and objecives (for instance quadratic objectives) to the model afterwards or use an approximation of Gibbs free energy directions in you model. Adds variables and constraints to a model which will disallow flux distributions with loops. The used formulation is described in [1]_. This function will modify your model.
 Parameters
model (cobra.Model) – The model to which to add the constraints.
zero_cutoff (positive float, optional) – Cutoff used for null space. Coefficients with an absolute value smaller than zero_cutoff are considered to be zero (default model.tolerance).
 Returns
 Return type
Nothing
References
 1
Elimination of thermodynamically infeasible loops in steadystate metabolic models. Schellenberger J, Lewis NE, Palsson BO. Biophys J. 2011 Feb 2;100(3):54453. doi: 10.1016/j.bpj.2010.12.3707. Erratum in: Biophys J. 2011 Mar 2;100(5):1381.

cobra.flux_analysis.
add_moma
(model, solution=None, linear=True)[source]¶ Add constraints and objective representing for MOMA.
This adds variables and constraints for the minimization of metabolic adjustment (MOMA) to the model.
 Parameters
model (cobra.Model) – The model to add MOMA constraints and objective to.
solution (cobra.Solution, optional) – A previous solution to use as a reference. If no solution is given, one will be computed using pFBA.
linear (bool, optional) – Whether to use the linear MOMA formulation or not (default True).
Notes
In the original MOMA [1]_ specification one looks for the flux distribution of the deletion (v^d) closest to the fluxes without the deletion (v). In math this means:
minimize sum_i (v^d_i  v_i)^2 s.t. Sv^d = 0
lb_i <= v^d_i <= ub_i
Here, we use a variable transformation v^t := v^d_i  v_i. Substituting and using the fact that Sv = 0 gives:
minimize sum_i (v^t_i)^2 s.t. Sv^d = 0
v^t = v^d_i  v_i lb_i <= v^d_i <= ub_i
So basically we just recenter the flux space at the old solution and then find the flux distribution closest to the new zero (center). This is the same strategy as used in cameo.
In the case of linear MOMA [2]_, we instead minimize sum_i abs(v^t_i). The linear MOMA is typically significantly faster. Also quadratic MOMA tends to give flux distributions in which all fluxes deviate from the reference fluxes a little bit whereas linear MOMA tends to give flux distributions where the majority of fluxes are the same reference with few fluxes deviating a lot (typical effect of L2 norm vs L1 norm).
The former objective function is saved in the optlang solver interface as
"moma_old_objective"
and this can be used to immediately extract the value of the former objective after MOMA optimization.See also
pfba()
parsimonious FBA
References
 1
Segrè, Daniel, Dennis Vitkup, and George M. Church. “Analysis of Optimality in Natural and Perturbed Metabolic Networks.” Proceedings of the National Academy of Sciences 99, no. 23 (November 12, 2002): 15112. https://doi.org/10.1073/pnas.232349399.
 2
Becker, Scott A, Adam M Feist, Monica L Mo, Gregory Hannum, Bernhard Ø Palsson, and Markus J Herrgard. “Quantitative Prediction of Cellular Metabolism with ConstraintBased Models: The COBRA Toolbox.” Nature Protocols 2 (March 29, 2007): 727.

cobra.flux_analysis.
moma
(model, solution=None, linear=True)[source]¶ Compute a single solution based on (linear) MOMA.
Compute a new flux distribution that is at a minimal distance to a previous reference solution. Minimization of metabolic adjustment (MOMA) is generally used to assess the impact of knockouts. Thus the typical usage is to provide a wildtype flux distribution as reference and a model in knockout state.
 Parameters
model (cobra.Model) – The model state to compute a MOMAbased solution for.
solution (cobra.Solution, optional) – A (wildtype) reference solution.
linear (bool, optional) – Whether to use the linear MOMA formulation or not (default True).
 Returns
A flux distribution that is at a minimal distance compared to the reference solution.
 Return type
See also
add_moma()
add MOMA constraints and objective

cobra.flux_analysis.
pfba
(model, fraction_of_optimum=1.0, objective=None, reactions=None)[source]¶ Perform basic pFBA (parsimonious Enzyme Usage Flux Balance Analysis) to minimize total flux.
pFBA [1] adds the minimization of all fluxes the the objective of the model. This approach is motivated by the idea that high fluxes have a higher enzyme turnover and that since producing enzymes is costly, the cell will try to minimize overall flux while still maximizing the original objective function, e.g. the growth rate.
 Parameters
model (cobra.Model) – The model
fraction_of_optimum (float, optional) – Fraction of optimum which must be maintained. The original objective reaction is constrained to be greater than maximal_value * fraction_of_optimum.
objective (dict or model.problem.Objective) – A desired objective to use during optimization in addition to the pFBA objective. Dictionaries (reaction as key, coefficient as value) can be used for linear objectives.
reactions (iterable) – List of reactions or reaction identifiers. Implies return_frame to be true. Only return fluxes for the given reactions. Faster than fetching all fluxes if only a few are needed.
 Returns
The solution object to the optimized model with pFBA constraints added.
 Return type
References
 1
Lewis, N. E., Hixson, K. K., Conrad, T. M., Lerman, J. A., Charusanti, P., Polpitiya, A. D., Palsson, B. O. (2010). Omic data from evolved E. coli are consistent with computed optimal growth from genomescale models. Molecular Systems Biology, 6, 390. doi:10.1038/msb.2010.47

cobra.flux_analysis.
find_blocked_reactions
(model, reaction_list=None, zero_cutoff=None, open_exchanges=False, processes=None)[source]¶ Find reactions that cannot carry any flux.
The question whether or not a reaction is blocked is highly dependent on the current exchange reaction settings for a COBRA model. Hence an argument is provided to open all exchange reactions.
Notes
Sink and demand reactions are left untouched. Please modify them manually.
 Parameters
model (cobra.Model) – The model to analyze.
reaction_list (list, optional) – List of reactions to consider, the default includes all model reactions.
zero_cutoff (float, optional) – Flux value which is considered to effectively be zero (default model.tolerance).
open_exchanges (bool, optional) – Whether or not to open all exchange reactions to very high flux ranges.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of reactions is large. If not explicitly passed, it will be set from the global configuration singleton.
 Returns
List with the identifiers of blocked reactions.
 Return type

cobra.flux_analysis.
find_essential_genes
(model, threshold=None, processes=None)[source]¶ Return a set of essential genes.
A gene is considered essential if restricting the flux of all reactions that depend on it to zero causes the objective, e.g., the growth rate, to also be zero, below the threshold, or infeasible.
 Parameters
model (cobra.Model) – The model to find the essential genes for.
threshold (float, optional) – Minimal objective flux to be considered viable. By default this is 1% of the maximal objective.
processes (int, optional) – The number of parallel processes to run. If not passed, will be set to the number of CPUs found.
processes – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not explicitly passed, it will be set from the global configuration singleton.
 Returns
Set of essential genes
 Return type

cobra.flux_analysis.
find_essential_reactions
(model, threshold=None, processes=None)[source]¶ Return a set of essential reactions.
A reaction is considered essential if restricting its flux to zero causes the objective, e.g., the growth rate, to also be zero, below the threshold, or infeasible.
 Parameters
model (cobra.Model) – The model to find the essential reactions for.
threshold (float, optional) – Minimal objective flux to be considered viable. By default this is 1% of the maximal objective.
processes (int, optional) – The number of parallel processes to run. Can speed up the computations if the number of knockouts to perform is large. If not explicitly passed, it will be set from the global configuration singleton.
 Returns
Set of essential reactions
 Return type

cobra.flux_analysis.
flux_variability_analysis
(model, reaction_list=None, loopless=False, fraction_of_optimum=1.0, pfba_factor=None, processes=None)[source]¶ Determine the minimum and maximum possible flux value for each reaction.
 Parameters
model (cobra.Model) – The model for which to run the analysis. It will not be modified.
reaction_list (list of cobra.Reaction or str, optional) – The reactions for which to obtain min/max fluxes. If None will use all reactions in the model (default).
loopless (boolean, optional) – Whether to return only loopless solutions. This is significantly slower. Please also refer to the notes.
fraction_of_optimum (float, optional) – Must be <= 1.0. Requires that the objective value is at least the fraction times maximum objective value. A value of 0.85 for instance means that the objective has to be at least at 85% percent of its maximum.
pfba_factor (float, optional) – Add an additional constraint to the model that requires the total sum of absolute fluxes must not be larger than this value times the smallest possible sum of absolute fluxes, i.e., by setting the value to 1.1 the total sum of absolute fluxes must not be more than 10% larger than the pFBA solution. Since the pFBA solution is the one that optimally minimizes the total flux sum, the
pfba_factor
should, if set, be larger than one. Setting this value may lead to more realistic predictions of the effective flux bounds.processes (int, optional) – The number of parallel processes to run. If not explicitly passed, will be set from the global configuration singleton.
 Returns
A data frame with reaction identifiers as the index and two columns:  maximum: indicating the highest possible flux  minimum: indicating the lowest possible flux
 Return type
pandas.DataFrame
Notes
This implements the fast version as described in [1]_. Please note that the flux distribution containing all minimal/maximal fluxes does not have to be a feasible solution for the model. Fluxes are minimized/maximized individually and a single minimal flux might require all others to be suboptimal.
Using the loopless option will lead to a significant increase in computation time (about a factor of 100 for large models). However, the algorithm used here (see [2]_) is still more than 1000x faster than the “naive” version using
add_loopless(model)
. Also note that if you have included constraints that force a loop (for instance by setting all fluxes in a loop to be nonzero) this loop will be included in the solution.References
 1
Computationally efficient flux variability analysis. Gudmundsson S, Thiele I. BMC Bioinformatics. 2010 Sep 29;11:489. doi: 10.1186/1471210511489, PMID: 20920235
 2
CycleFreeFlux: efficient removal of thermodynamically infeasible loops from flux distributions. Desouki AA, Jarre F, GeliusDietrich G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):215965. doi: 10.1093/bioinformatics/btv096.

cobra.flux_analysis.
add_room
(model, solution=None, linear=False, delta=0.03, epsilon=0.001)[source]¶ Add constraints and objective for ROOM.
This function adds variables and constraints for applying regulatory on/off minimization (ROOM) to the model.
 Parameters
model (cobra.Model) – The model to add ROOM constraints and objective to.
solution (cobra.Solution, optional) – A previous solution to use as a reference. If no solution is given, one will be computed using pFBA.
linear (bool, optional) – Whether to use the linear ROOM formulation or not (default False).
delta (float, optional) – The relative tolerance range which is additive in nature (default 0.03).
epsilon (float, optional) – The absolute range of tolerance which is multiplicative (default 0.001).
Notes
The formulation used here is the same as stated in the original paper [1]_. The mathematical expression is given below:
minimize sum_{i=1}^m y^i s.t. Sv = 0
v_min <= v <= v_max v_j = 0 j ∈ A for 1 <= i <= m v_i  y_i(v_{max,i}  w_i^u) <= w_i^u (1) v_i  y_i(v_{min,i}  w_i^l) <= w_i^l (2) y_i ∈ {0,1} (3) w_i^u = w_i + deltaw_i + epsilon w_i^l = w_i  deltaw_i  epsilon
So, for the linear version of the ROOM , constraint (3) is relaxed to 0 <= y_i <= 1.
See also
pfba()
parsimonious FBA
References
 1
Tomer Shlomi, Omer Berkman and Eytan Ruppin, “Regulatory on/off minimization of metabolic flux changes after genetic perturbations”, PNAS 2005 102 (21) 76957700; doi:10.1073/pnas.0406346102

cobra.flux_analysis.
room
(model, solution=None, linear=False, delta=0.03, epsilon=0.001)[source]¶ Compute a single solution based on regulatory on/off minimization (ROOM).
Compute a new flux distribution that minimizes the number of active reactions needed to accommodate a previous reference solution. Regulatory on/off minimization (ROOM) is generally used to assess the impact of knockouts. Thus the typical usage is to provide a wildtype flux distribution as reference and a model in knockout state.
 Parameters
model (cobra.Model) – The model state to compute a ROOMbased solution for.
solution (cobra.Solution, optional) – A (wildtype) reference solution.
linear (bool, optional) – Whether to use the linear ROOM formulation or not (default False).
delta (float, optional) – The relative tolerance range (additive) (default 0.03).
epsilon (float, optional) – The absolute tolerance range (multiplicative) (default 0.001).
 Returns
A flux distribution with minimal active reaction changes compared to the reference.
 Return type
See also
add_room()
add ROOM constraints and objective
cobra.io
¶
Provide functions for loading and saving metabolic models.
Subpackages¶
cobra.io.web
¶Provide functionality to access remote model repositories.
cobra.io.web.abstract_model_repository
¶Provide an abstract base class that describes a remote model repository.
Define an abstract base class that describes a remote model repository. 

class
cobra.io.web.abstract_model_repository.
AbstractModelRepository
(*, url: Union[httpx.URL, str], **kwargs)[source]¶ Bases:
abc.ABC
Define an abstract base class that describes a remote model repository.

name
:str = Abstract[source]

cobra.io.web.bigg_models_repository
¶Provide a concrete implementation of the BioModels repository interface.
Define a concrete implementation of the BiGG Models repository. 

class
cobra.io.web.bigg_models_repository.
BiGGModels
(**kwargs)[source]¶ Bases:
cobra.io.web.abstract_model_repository.AbstractModelRepository
Define a concrete implementation of the BiGG Models repository.

name
:str = BiGG Models[source]

cobra.io.web.biomodels_repository
¶Provide functions for loading metabolic models over the wire.
Define a single BioModels file description. 

Define the BioModels files JSON response. 

Define a concrete implementation of the BioModels repository. 

class
cobra.io.web.biomodels_repository.
BioModelsFile
[source]¶ Bases:
pydantic.BaseModel
Define a single BioModels file description.

class
cobra.io.web.biomodels_repository.
BioModelsFilesResponse
[source]¶ Bases:
pydantic.BaseModel
Define the BioModels files JSON response.

class
cobra.io.web.biomodels_repository.
BioModels
(**kwargs)[source]¶ Bases:
cobra.io.web.abstract_model_repository.AbstractModelRepository
Define a concrete implementation of the BioModels repository.

name
:str = BioModels[source]

cobra.io.web.load
¶Provide a function load_model
to access remote model repositories.

Download an SBML model from a remote repository. 

Attempt to load a gzipcompressed SBML document from the cache. 

Attempt to load a gzipcompressed SBML document from the given repositories. 

Generate a model instance from a gzipcompressed, UTF8 encoded SBML document. 

cobra.io.web.load.
configuration