documentation_builder/getting_started.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting Started"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading a model and inspecting it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin with, cobrapy comes with bundled models for _Salmonella_ and _E. coli_, as well as a \"textbook\" model of _E. coli_ core metabolism. To load a test model, type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scaling...\n",
      " A: min|aij| =  1.000e+00  max|aij| =  1.000e+00  ratio =  1.000e+00\n",
      "Problem data seem to be well scaled\n"
     ]
    }
   ],
   "source": [
    "import cobra\n",
    "from cobra.io import load_model\n",
    "\n",
    "# \"iJO1366\" and \"salmonella\" are also valid arguments\n",
    "model = load_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": [
      "95\n",
      "72\n",
      "137\n"
     ]
    }
   ],
   "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://jupyter-notebook-beginner-guide.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",
       "                <td>0x07feb1af21f28</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>Number of groups</strong></td>\n",
       "                <td>0</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": [
       "<Model e_coli_core at 0x7feb1af21f28>"
      ]
     },
     "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 [0-indexing](https://en.wikipedia.org/wiki/Zero-based_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>L-Glutamate exchange</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Memory address</strong></td>\n",
       "                <td>0x7feb1b012eb8</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Stoichiometry</strong></td>\n",
       "                <td>\n",
       "                    <p style='text-align:right'>glu__L_e --></p>\n",
       "                    <p style='text-align:right'>L-Glutamate --></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": [
       "<Reaction EX_glu__L_e at 0x7feb1b012eb8>"
      ]
     },
     "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",
       "                <td>0x7feb1af41550</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",
       "                    Biomass_Ecoli_core, ATPS4r, PGK, ACKr, ADK1, GLNS, ATPM, PFK, PPS, GLNabc, PYK, PPCK, SUCOAS\n",
       "                    </td>\n",
       "            </tr>\n",
       "        </table>"
      ],
      "text/plain": [
       "<Metabolite atp_c at 0x7feb1af41550>"
      ]
     },
     "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 tab-complete 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 6-phosphate isomerase, which interconverts glucose 6-phosphate and fructose 6-phosphate. 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>glucose-6-phosphate isomerase</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Memory address</strong></td>\n",
       "                <td>0x7feb1b09c2e8</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Stoichiometry</strong></td>\n",
       "                <td>\n",
       "                    <p style='text-align:right'>g6p_c <=> f6p_c</p>\n",
       "                    <p style='text-align:right'>D-Glucose 6-phosphate <=> D-Fructose 6-phosphate</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": [
       "<Reaction PGI at 0x7feb1b09c2e8>"
      ]
     },
     "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": [
      "glucose-6-phosphate isomerase\n",
      "g6p_c <=> f6p_c\n"
     ]
    }
   ],
   "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.0\n",
      "True\n"
     ]
    }
   ],
   "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": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 < pgi < 1000.0\n",
      "Reversibility after modification: False\n",
      "Reversibility after resetting: True\n"
     ]
    }
   ],
   "source": [
    "old_bounds = pgi.bounds\n",
    "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_bounds\n",
    "print(\"Reversibility after resetting:\", pgi.reversibility)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bounds can also be modified one-at-a-time 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). If a lower bound is *higher* than upper bound (or vice versa), this will lead to an error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Upper bound prior to setting new lower bound: 1000.0\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "The lower bound must be less than or equal to the upper bound (1100 <= 1000.0).",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/v7/99h1ft7s0mzg1_83fmbc8mnm0000gn/T/ipykernel_23255/4168182522.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mold_bounds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpgi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbounds\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Upper bound prior to setting new lower bound:'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpgi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupper_bound\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mpgi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlower_bound\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1100\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Upper bound after setting new lower bound:'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpgi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupper_bound\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mpgi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbounds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mold_bounds\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/PycharmProjects/cobrapy/src/cobra/util/context.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(self, new_value)\u001b[0m\n\u001b[1;32m    109\u001b[0m             \u001b[0mcontext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpartial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mold_value\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    110\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 111\u001b[0;31m         \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_value\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    112\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    113\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/PycharmProjects/cobrapy/src/cobra/core/reaction.py\u001b[0m in \u001b[0;36mlower_bound\u001b[0;34m(self, value)\u001b[0m\n\u001b[1;32m    368\u001b[0m         \"\"\"\n\u001b[1;32m    369\u001b[0m         \u001b[0;31m# Validate bounds before setting them.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 370\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_check_bounds\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_upper_bound\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    371\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lower_bound\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    372\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate_variable_bounds\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/PycharmProjects/cobrapy/src/cobra/core/reaction.py\u001b[0m in \u001b[0;36m_check_bounds\u001b[0;34m(lb, ub)\u001b[0m\n\u001b[1;32m    292\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlb\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mub\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    293\u001b[0m             raise ValueError(\n\u001b[0;32m--> 294\u001b[0;31m                 \u001b[0;34mf\"The lower bound must be less than or equal to the upper bound \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    295\u001b[0m                 \u001b[0;34mf\"({lb} <= {ub}).\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    296\u001b[0m             )\n",
      "\u001b[0;31mValueError\u001b[0m: The lower bound must be less than or equal to the upper bound (1100 <= 1000.0)."
     ]
    }
   ],
   "source": [
    "old_bounds = pgi.bounds\n",
    "print('Upper bound prior to setting new lower bound:', pgi.upper_bound)\n",
    "pgi.lower_bound = 1100\n",
    "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_c\n",
      "{}\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' created\n",
      "unknown metabolite 'ham' created\n"
     ]
    }
   ],
   "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",
       "                <td>0x7feb1af41550</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",
       "                    Biomass_Ecoli_core, ATPS4r, PGK, ACKr, ADK1, GLNS, ATPM, PFK, PPS, GLNabc, PYK, PPCK, SUCOAS\n",
       "                    </td>\n",
       "            </tr>\n",
       "        </table>"
      ],
      "text/plain": [
       "<Metabolite atp_c at 0x7feb1af41550>"
      ]
     },
     "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": [
      "ATP\n",
      "c\n"
     ]
    }
   ],
   "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": [
      "C10H12N5O13P3\n"
     ]
    }
   ],
   "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 6-phosphate will participate in fewer reactions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({<Reaction Biomass_Ecoli_core at 0x7feb1afd4240>,\n",
       "           <Reaction G6PDH2r at 0x7feb1b04bf60>,\n",
       "           <Reaction GLCpts at 0x7feb1b050b38>,\n",
       "           <Reaction PGI at 0x7feb1b09c2e8>})"
      ]
     },
     "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):1290-307](http://dx.doi.org/doi:10.1038/nprot.2011.308).\n",
    "\n",
    "The GPR is stored as the GPR class in the gpr for a Reaction. A string representation of it is stored as the gene_reaction_rule for a Reaction object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b4025\n",
      "b4025\n"
     ]
    }
   ],
   "source": [
    "gpr = pgi.gpr\n",
    "print(gpr)\n",
    "gpr_string = pgi.gene_reaction_rule\n",
    "print(gpr_string)"
   ]
  },
  {
   "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": [
       "frozenset({<Gene b4025 at 0x7feb1afb0f98>})"
      ]
     },
     "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",
       "                <td>0x7feb1afb0f98</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\n",
       "                    </td>\n",
       "            </tr>\n",
       "        </table>"
      ],
      "text/plain": [
       "<Gene b4025 at 0x7feb1afb0f98>"
      ]
     },
     "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": [
       "frozenset({<Reaction PGI at 0x7feb1b09c2e8>})"
      ]
     },
     "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": [
       "frozenset({<Gene eggs at 0x7feb1b0b0cc0>, <Gene spam at 0x7feb1b0aacf8>})"
      ]
     },
     "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",
       "                <td>0x7feb1b0aacf8</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\n",
       "                    </td>\n",
       "            </tr>\n",
       "        </table>"
      ],
      "text/plain": [
       "<Gene spam at 0x7feb1b0aacf8>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.genes.get_by_id(\"spam\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `knock_out_model_genes` function will evaluate the GPR and set the upper and lower bounds to 0 if the reaction is knocked out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "after 1 KO: -1000 < flux_PGI < 1000\n",
      "after 2 KO:     0 < flux_PGI <    0\n"
     ]
    }
   ],
   "source": [
    "cobra.manipulation.knock_out_model_genes(\n",
    "    model, [\"spam\"])\n",
    "print(\"after 1 KO: %4d < flux_PGI < %4d\" % (pgi.lower_bound, pgi.upper_bound))\n",
    "\n",
    "cobra.manipulation.knock_out_model_genes(\n",
    "    model, [\"eggs\"])\n",
    "print(\"after 2 KO:  %4d < flux_PGI < %4d\" % (pgi.lower_bound, pgi.upper_bound))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When knocking out model genes in a context, it is reversed when leaving the context."
   ]
  },
  {
   "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 knock-out 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 knock-out 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 knock-out, 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.873922\n",
      "ACALDt blocked (bounds: (0, 0)), new growth rate 0.873922\n",
      "ACKr blocked (bounds: (0, 0)), new growth rate 0.873922\n",
      "ACONTa blocked (bounds: (0, 0)), new growth rate -0.000000\n",
      "ACONTb blocked (bounds: (0, 0)), new growth rate -0.000000\n"
     ]
    }
   ],
   "source": [
    "model = load_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_2cdba\n",
      "print objective in first context: 1.0*ATPM - 1.0*ATPM_reverse_5b752\n",
      "print objective in second context: 1.0*ACALD - 1.0*ACALD_reverse_fda2b\n",
      "objective after exiting second context: 1.0*ATPM - 1.0*ATPM_reverse_5b752\n",
      "back to original objective: 1.0*Biomass_Ecoli_core - 1.0*Biomass_Ecoli_core_reverse_2cdba\n"
     ]
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}