Advanced-SynBio-for-Cell-Factories-Course/Vanillin Production.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vanillin production\n",
    "------------------\n",
    "\n",
    "In 2010, Brochado *et al* used heuristic optimization together with flux simulations to design a vanillin producing yeast strain.\n",
    "\n",
    "Brochado, A. R., Andrejev, S., Maranas, C. D., & Patil, K. R. (2012). Impact of stoichiometry representation on simulation of genotype-phenotype relationships in metabolic networks. PLoS Computational Biology, 8(11), e1002758. doi:10.1371/journal.pcbi.1002758"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Genome-scale metabolic model\n",
    "--------------------------\n",
    "\n",
    "In their work, the authors used *iFF708* model, but recent insights in Yeast yielded newer and more complete versions. \n",
    "Becuase this algorithms should be agnostic to the model, we implement the same strategy with a newer model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "(function(global) {\n",
       "  function now() {\n",
       "    return new Date();\n",
       "  }\n",
       "\n",
       "  if (typeof (window._bokeh_onload_callbacks) === \"undefined\") {\n",
       "    window._bokeh_onload_callbacks = [];\n",
       "  }\n",
       "\n",
       "  function run_callbacks() {\n",
       "    window._bokeh_onload_callbacks.forEach(function(callback) { callback() });\n",
       "    delete window._bokeh_onload_callbacks\n",
       "    console.info(\"Bokeh: all callbacks have finished\");\n",
       "  }\n",
       "\n",
       "  function load_libs(js_urls, callback) {\n",
       "    window._bokeh_onload_callbacks.push(callback);\n",
       "    if (window._bokeh_is_loading > 0) {\n",
       "      console.log(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
       "      return null;\n",
       "    }\n",
       "    if (js_urls == null || js_urls.length === 0) {\n",
       "      run_callbacks();\n",
       "      return null;\n",
       "    }\n",
       "    console.log(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
       "    window._bokeh_is_loading = js_urls.length;\n",
       "    for (var i = 0; i < js_urls.length; i++) {\n",
       "      var url = js_urls[i];\n",
       "      var s = document.createElement('script');\n",
       "      s.src = url;\n",
       "      s.async = false;\n",
       "      s.onreadystatechange = s.onload = function() {\n",
       "        window._bokeh_is_loading--;\n",
       "        if (window._bokeh_is_loading === 0) {\n",
       "          console.log(\"Bokeh: all BokehJS libraries loaded\");\n",
       "          run_callbacks()\n",
       "        }\n",
       "      };\n",
       "      s.onerror = function() {\n",
       "        console.warn(\"failed to load library \" + url);\n",
       "      };\n",
       "      console.log(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
       "      document.getElementsByTagName(\"head\")[0].appendChild(s);\n",
       "    }\n",
       "  };\n",
       "\n",
       "  var js_urls = ['https://cdn.pydata.org/bokeh/release/bokeh-0.11.1.min.js', 'https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.11.1.min.js', 'https://cdn.pydata.org/bokeh/release/bokeh-compiler-0.11.1.min.js'];\n",
       "\n",
       "  var inline_js = [\n",
       "    function(Bokeh) {\n",
       "      Bokeh.set_log_level(\"info\");\n",
       "    },\n",
       "    \n",
       "    function(Bokeh) {\n",
       "      Bokeh.$(\"#c3220257-8cf2-4a0f-b1a6-85d36ccc8b9d\").text(\"BokehJS successfully loaded\");\n",
       "    },\n",
       "    function(Bokeh) {\n",
       "      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-0.11.1.min.css\");\n",
       "      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-0.11.1.min.css\");\n",
       "      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.11.1.min.css\");\n",
       "      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.11.1.min.css\");\n",
       "    }\n",
       "  ];\n",
       "\n",
       "  function run_inline_js() {\n",
       "    for (var i = 0; i < inline_js.length; i++) {\n",
       "      inline_js[i](window.Bokeh);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  if (window._bokeh_is_loading === 0) {\n",
       "    console.log(\"Bokeh: BokehJS loaded, going straight to plotting\");\n",
       "    run_inline_js();\n",
       "  } else {\n",
       "    load_libs(js_urls, function() {\n",
       "      console.log(\"Bokeh: BokehJS plotting callback run at\", now());\n",
       "      run_inline_js();\n",
       "    });\n",
       "  }\n",
       "}(this));"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from cameo import models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = models.bigg.iMM904"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constraints can be set in the model according to data found in the literature. The defined conditions allow the simulation of phenotypes very close to the experimental results.\n",
    "\n",
    "<img src=http://www.biomedcentral.com/content/figures/1752-0509-7-36-2.jpg/>\n",
    "**Model validation by comparing in silico prediction of the specific growth rate with experimental data**. Growth phenotypes were collected from literature and compared to simulated values for chemostat cultivations at four different conditions, nitrogen limited aerobic (green) and anaerobic (red), carbon limited aerobic (blue) and anaerobic (white).\n",
    "\n",
    "\n",
    "Österlund, T., Nookaew, I., Bordel, S., & Nielsen, J. (2013). Mapping condition-dependent regulation of metabolism in yeast through genome-scale modeling. BMC Systems Biology, 7, 36. doi:10.1186/1752-0509-7-36"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model.reactions.EX_glc__D_e.lower_bound = -13 #glucose exchange\n",
    "model.reactions.EX_o2_e.lower_bound = -3  #oxygen exchange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>reaction_id</th>\n",
       "      <th>reaction_name</th>\n",
       "      <th>lower_bound</th>\n",
       "      <th>upper_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>EX_fe2_e</td>\n",
       "      <td>Fe2+ exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>EX_glc__D_e</td>\n",
       "      <td>D-Glucose exchange</td>\n",
       "      <td>-13.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>EX_h2o_e</td>\n",
       "      <td>H2O exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>EX_h_e</td>\n",
       "      <td>H+ exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>EX_k_e</td>\n",
       "      <td>K+ exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>EX_na1_e</td>\n",
       "      <td>Sodium exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>EX_nh4_e</td>\n",
       "      <td>Ammonia exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>EX_o2_e</td>\n",
       "      <td>O2 exchange</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>EX_pi_e</td>\n",
       "      <td>Phosphate exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>EX_so4_e</td>\n",
       "      <td>Sulfate exchange</td>\n",
       "      <td>-999999.0</td>\n",
       "      <td>999999.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   reaction_id       reaction_name  lower_bound  upper_bound\n",
       "0     EX_fe2_e       Fe2+ exchange    -999999.0     999999.0\n",
       "1  EX_glc__D_e  D-Glucose exchange        -13.0     999999.0\n",
       "2     EX_h2o_e        H2O exchange    -999999.0     999999.0\n",
       "3       EX_h_e         H+ exchange    -999999.0     999999.0\n",
       "4       EX_k_e         K+ exchange    -999999.0     999999.0\n",
       "5     EX_na1_e     Sodium exchange    -999999.0     999999.0\n",
       "6     EX_nh4_e    Ammonia exchange    -999999.0     999999.0\n",
       "7      EX_o2_e         O2 exchange         -3.0     999999.0\n",
       "8      EX_pi_e  Phosphate exchange    -999999.0     999999.0\n",
       "9     EX_so4_e    Sulfate exchange    -999999.0     999999.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.medium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.39022235350798895"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.objective = model.reactions.BIOMASS_SC5_notrace #growth\n",
    "model.optimize().f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Heterologous pathway\n",
    "-------------------\n",
    "\n",
    "Vanillin is not produced by *S. cervisiae*. In their work an heterolgous pathway is inserted to allow generate a vanillin production strain. The pathway is described as:\n",
    "\n",
    "<img src=http://static-content.springer.com/image/art%3A10.1186%2F1475-2859-9-84/MediaObjects/12934_2010_Article_474_Fig1_HTML.jpg>\n",
    "**Schematic representation of the de novo VG biosynthetic pathway in S. Cerevisisae** (as designed by Hansen et al [5]). Metabolites are shown in black, enzymes are shown in black and in italic, cofactors and additional precursors are shown in red. Reactions catalyzed by heterologously introduced enzymes are shown in red. Reactions converting glucose to aromatic amino acids are represented by dashed black arrows. Metabolite secretion is represented by solid black arrows where relative thickness corresponds to relative extracellular accumulation. 3-DSH stands for 3-dedhydroshikimate, PAC stands for protocathechuic acid, PAL stands for protocatechuic aldehyde, SAM stands for S-adenosylmethionine. 3DSD stands for 3-dedhydroshikimate dehydratase, ACAR stands for aryl carboxylic acid reductase, PPTase stands for phosphopantetheine transferase, hsOMT stands for O-methyltransferase, and UGT stands for UDP-glycosyltransferase. Adapted from Hansen et al. [5].\n",
    "Brochado et al. Microbial Cell Factories 2010 9:84   doi:10.1186/1475-2859-9-84\n",
    "\n",
    "Using **cameo**, is very easy to generate a pathway and add it to a model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from cameo.strain_design.pathway_prediction import PathwayPredictor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predictor = PathwayPredictor(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<span>Pathway 1</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>equation</th>\n",
       "      <th>lower_bound</th>\n",
       "      <th>upper_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>MNXR5336</th>\n",
       "      <td>H2O + vanillin + NAD(+) &lt;=&gt; 2.0 H(+) + NADH(2-...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MNXR5340</th>\n",
       "      <td>3,4-dihydroxybenzoate + NAD(+) + H2O + formald...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MNXR14769</th>\n",
       "      <td>H2O + 3,4-dihydroxybenzoate + NAD(+) &lt;=&gt; NADH(...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                    equation  lower_bound  \\\n",
       "MNXR5336   H2O + vanillin + NAD(+) <=> 2.0 H(+) + NADH(2-...        -1000   \n",
       "MNXR5340   3,4-dihydroxybenzoate + NAD(+) + H2O + formald...        -1000   \n",
       "MNXR14769  H2O + 3,4-dihydroxybenzoate + NAD(+) <=> NADH(...        -1000   \n",
       "\n",
       "           upper_bound  \n",
       "MNXR5336          1000  \n",
       "MNXR5340          1000  \n",
       "MNXR14769         1000  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<span>Pathway 2</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>equation</th>\n",
       "      <th>lower_bound</th>\n",
       "      <th>upper_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>MNXR5336</th>\n",
       "      <td>H2O + vanillin + NAD(+) &lt;=&gt; 2.0 H(+) + NADH(2-...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MNXR5340</th>\n",
       "      <td>3,4-dihydroxybenzoate + NAD(+) + H2O + formald...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MNXR68718</th>\n",
       "      <td>H2O + 3,4-dihydroxybenzoate &lt;=&gt; 3-dehydroshiki...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                    equation  lower_bound  \\\n",
       "MNXR5336   H2O + vanillin + NAD(+) <=> 2.0 H(+) + NADH(2-...        -1000   \n",
       "MNXR5340   3,4-dihydroxybenzoate + NAD(+) + H2O + formald...        -1000   \n",
       "MNXR68718  H2O + 3,4-dihydroxybenzoate <=> 3-dehydroshiki...        -1000   \n",
       "\n",
       "           upper_bound  \n",
       "MNXR5336          1000  \n",
       "MNXR5340          1000  \n",
       "MNXR68718         1000  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<span>Pathway 3</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>equation</th>\n",
       "      <th>lower_bound</th>\n",
       "      <th>upper_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>MNXR230</th>\n",
       "      <td>3,4-dihydroxybenzoate + NADP(+) + H2O &lt;=&gt; O2 +...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MNXR5336</th>\n",
       "      <td>H2O + vanillin + NAD(+) &lt;=&gt; 2.0 H(+) + NADH(2-...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MNXR5340</th>\n",
       "      <td>3,4-dihydroxybenzoate + NAD(+) + H2O + formald...</td>\n",
       "      <td>-1000</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                   equation  lower_bound  \\\n",
       "MNXR230   3,4-dihydroxybenzoate + NADP(+) + H2O <=> O2 +...        -1000   \n",
       "MNXR5336  H2O + vanillin + NAD(+) <=> 2.0 H(+) + NADH(2-...        -1000   \n",
       "MNXR5340  3,4-dihydroxybenzoate + NAD(+) + H2O + formald...        -1000   \n",
       "\n",
       "          upper_bound  \n",
       "MNXR230          1000  \n",
       "MNXR5336         1000  \n",
       "MNXR5340         1000  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pathways = predictor.run('vanillin', max_predictions=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vanillin_pathway = pathways.pathways[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from cameo.core.pathway import Pathway"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>equation</th>\n",
       "      <th>lower_bound</th>\n",
       "      <th>upper_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>3DSD</th>\n",
       "      <td>3-dehydroshikimate --&gt; H2O + protocathechuic acid</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACAR_PPTase</th>\n",
       "      <td>ATP + protocathechuic acid + NADPH --&gt; ADP + p...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hsOMT</th>\n",
       "      <td>S-adenosyl-L-methionine + protocatechuic aldeh...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UGT</th>\n",
       "      <td>Vanillin + UDP-glucose --&gt; vanillin-B-glucoside</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                      equation  lower_bound  \\\n",
       "3DSD         3-dehydroshikimate --> H2O + protocathechuic acid          0.0   \n",
       "ACAR_PPTase  ATP + protocathechuic acid + NADPH --> ADP + p...          0.0   \n",
       "hsOMT        S-adenosyl-L-methionine + protocatechuic aldeh...          0.0   \n",
       "UGT            Vanillin + UDP-glucose --> vanillin-B-glucoside          0.0   \n",
       "\n",
       "             upper_bound  \n",
       "3DSD              1000.0  \n",
       "ACAR_PPTase       1000.0  \n",
       "hsOMT             1000.0  \n",
       "UGT               1000.0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vanillin_pathway = Pathway.from_file(\"data/vanillin_pathway.tsv\")\n",
    "vanillin_pathway.data_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can plug the pathway to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vanillin_pathway.plug_model(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from cameo import phenotypic_phase_plane"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Phenotypic phase plane can be used to analyse the theoretical yields at different growth rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "production_envelope = phenotypic_phase_plane(model, variables=[model.reactions.BIOMASS_SC5_notrace], \n",
    "                                             objective=model.reactions.EX_vnl_b_glu_c)\n",
    "production_envelope.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "DictList has no attribute or entry EX_vnl_b_glu_c",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0mTraceback (most recent call last)",
      "\u001b[1;32m/opt/conda/envs/python3.4/lib/python3.4/site-packages/cobra/core/DictList.py\u001b[0m in \u001b[0;36m__getattr__\u001b[1;34m(self, attr)\u001b[0m\n\u001b[0;32m    355\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 356\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mDictList\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_by_id\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mattr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    357\u001b[0m         \u001b[1;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/opt/conda/envs/python3.4/lib/python3.4/site-packages/cobra/core/DictList.py\u001b[0m in \u001b[0;36mget_by_id\u001b[1;34m(self, id)\u001b[0m\n\u001b[0;32m     49\u001b[0m         \u001b[1;34m\"\"\"return the element with a matching id\"\"\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 50\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__getitem__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_dict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mid\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     51\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: 'EX_vnl_b_glu_c'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[1;31mAttributeError\u001b[0mTraceback (most recent call last)",
      "\u001b[1;32m<ipython-input-16-ac7b728ba665>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m production_envelope = phenotypic_phase_plane(model, variables=[model.reactions.BIOMASS_SC5_notrace], \n\u001b[1;32m----> 2\u001b[1;33m                                              objective=model.reactions.EX_vnl_b_glu_c)\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mproduction_envelope\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/opt/conda/envs/python3.4/lib/python3.4/site-packages/cobra/core/DictList.py\u001b[0m in \u001b[0;36m__getattr__\u001b[1;34m(self, attr)\u001b[0m\n\u001b[0;32m    357\u001b[0m         \u001b[1;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    358\u001b[0m             raise AttributeError(\"DictList has no attribute or entry %s\" %\n\u001b[1;32m--> 359\u001b[1;33m                                  (attr))\n\u001b[0m\u001b[0;32m    360\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    361\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__dir__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: DictList has no attribute or entry EX_vnl_b_glu_c"
     ]
    }
   ],
   "source": [
    "production_envelope = phenotypic_phase_plane(model, variables=[model.reactions.BIOMASS_SC5_notrace], \n",
    "                                             objective=model.reactions.EX_vnl_b_glu_c)\n",
    "production_envelope.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find gene knockout targets, we use `cameo.strain_design.heuristic` package which implements the OptGene strategy.\n",
    "\n",
    "The authors used the biomass-product coupled yield (bpcy) for optimization which is the equivalent of running OptGene in non-robust mode. All simulations were computed using MOMA but because **cameo** does not implement MOMA we use it's equivalent linear version (it minimizes the absolute distance instead of the quadratic distance). The linear MOMA version is faster than the original MOMA formulation.\n",
    "\n",
    "By default, our OptGene implementation will run 20'000 evaluations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from cameo.strain_design.heuristic.evolutionary_based import OptGene\n",
    "from cameo.flux_analysis.simulation import lmoma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "optgene = OptGene(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting optimization at Thu, 26 May 2016 13:53:54\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run \"jupyter nbextension enable --py --sys-prefix widgetsnbextension\"\n",
      "/Users/joao/.virtualenvs/cameo06-py3/lib/python3.4/site-packages/ipywidgets/widgets/widget_string.py:55: UserWarning:\n",
      "\n",
      "The Latex widget is deprecated. Use Label instead\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished after 01:28:45\n"
     ]
    }
   ],
   "source": [
    "results = optgene.run(target=\"EX_vnl_b_glu_c\", \n",
    "                      biomass=\"BIOMASS_SC5_notrace\",\n",
    "                      substrate=\"EX_glc__D_e\",\n",
    "                      simulation_method=lmoma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run \"jupyter nbextension enable --py --sys-prefix widgetsnbextension\"\n",
      "/Users/joao/.virtualenvs/cameo06-py3/lib/python3.4/site-packages/ipywidgets/widgets/widget_string.py:55: UserWarning:\n",
      "\n",
      "The Latex widget is deprecated. Use Label instead\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "        <h3>OptGene Result</h3>\n",
       "        <ul>\n",
       "            <li>Simulation: lmoma<br/></li>\n",
       "            <li>Objective Function: $$bpcy = \\frac{(BIOMASS\\_SC5\\_notrace * EX\\_vnl\\_b\\_glu\\_c)}{EX\\_glc\\_\\_D\\_e}$$<br/></li>\n",
       "        </ul>\n",
       "        <div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>reactions</th>\n",
       "      <th>genes</th>\n",
       "      <th>size</th>\n",
       "      <th>fva_min</th>\n",
       "      <th>fva_max</th>\n",
       "      <th>target_flux</th>\n",
       "      <th>biomass_flux</th>\n",
       "      <th>yield</th>\n",
       "      <th>fitness</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>(OHPBAT, G6PI, PSERT, GLCGSD, OMCDCm, PGI, VAL...</td>\n",
       "      <td>((YBR196C, YBL011W, YOR184W, YHR208W, YJL121C,...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023639</td>\n",
       "      <td>0.176487</td>\n",
       "      <td>0.001896</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>(OHPBAT, GLCGSD, G6PI, PSERT, OMCDCm, PGI, VAL...</td>\n",
       "      <td>((YAR035W, YBR196C, YOR184W, YHR208W, YJL121C,...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023639</td>\n",
       "      <td>0.176487</td>\n",
       "      <td>0.001896</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>(OHPBAT, G6PI, PSERT, GLCGSD, PGI, G6PI3, RPE,...</td>\n",
       "      <td>((YBR196C, YBL011W, YOR184W, YLR432W, YJL121C,...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023565</td>\n",
       "      <td>0.177159</td>\n",
       "      <td>0.001890</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>(OHPBAT, PMI12346PH, GLCGSD, G6PI, PSERT, PGI,...</td>\n",
       "      <td>((YBR196C, YOR184W, YJL121C, YOR163W, YPR184W,...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023565</td>\n",
       "      <td>0.177159</td>\n",
       "      <td>0.001890</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>(OHPBAT, G6PI, PSERT, ALATA_L, PGI, G6PI3, RPE...</td>\n",
       "      <td>((YDR111C, YBR196C, YOR184W, YOR377W, YJL121C,...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009063</td>\n",
       "      <td>0.023205</td>\n",
       "      <td>0.180413</td>\n",
       "      <td>0.001862</td>\n",
       "      <td>0.000336</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>(OHPBAT, GLCGSD, G6PI, PSERT, PGI, AMETt2, G6P...</td>\n",
       "      <td>((YBR196C, YOR184W, YOR071C, YJL121C, YPR184W,...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023565</td>\n",
       "      <td>0.177159</td>\n",
       "      <td>0.001890</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>(OHPBAT, FAS120ACPm, PSERT, PGI, FAS80ACPm_L, ...</td>\n",
       "      <td>((YAR035W, YBR196C, YER061C, YOR184W, YJL121C,...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023565</td>\n",
       "      <td>0.177159</td>\n",
       "      <td>0.001890</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>(OHPBAT, ALLTTti, GLCGSD, G6PI, PSERT, PGI, G6...</td>\n",
       "      <td>((YAR035W, YBR196C, YOR184W, YJL121C, YPR184W,...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023565</td>\n",
       "      <td>0.177159</td>\n",
       "      <td>0.001890</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>(OHPBAT, ALLTTti, GLCGSD, G6PI, PSERT, PGI, G6...</td>\n",
       "      <td>((YBR196C, YOR184W, YBR166C, YJL121C, YPR184W,...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023494</td>\n",
       "      <td>0.177800</td>\n",
       "      <td>0.001891</td>\n",
       "      <td>0.000336</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>(GALKr, OHPBAT, G6PI, PSERT, AGTi, PGI, G6PI3,...</td>\n",
       "      <td>((YBR196C, YOR184W, YJL121C, YDL052C, YOL064C,...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.008878</td>\n",
       "      <td>0.023364</td>\n",
       "      <td>0.178975</td>\n",
       "      <td>0.001876</td>\n",
       "      <td>0.000336</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>(GALKr, OHPBAT, DKMPPD2, G6PI, PSERT, AGTi, PG...</td>\n",
       "      <td>((YLR354C, YBR196C, YOR184W, YMR009W, YER099C,...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.008878</td>\n",
       "      <td>0.023388</td>\n",
       "      <td>0.178755</td>\n",
       "      <td>0.001878</td>\n",
       "      <td>0.000336</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>(FAO141p_odd, OHPBAT, FAO141p_even, PSERT, FAO...</td>\n",
       "      <td>((YAR035W, YBR196C, YOR184W, YJL121C, YPR184W,...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023565</td>\n",
       "      <td>0.177159</td>\n",
       "      <td>0.001890</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>(OHPBAT, GLCGSD, G6PI, PSERT, ALLTTti, PGI, G6...</td>\n",
       "      <td>((YBR196C, YOR184W, YJL121C, YPR184W, YMR118C,...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.009158</td>\n",
       "      <td>0.023565</td>\n",
       "      <td>0.177159</td>\n",
       "      <td>0.001890</td>\n",
       "      <td>0.000335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>(OHPBAT, PSERT, AGTi, PGI, RPE, POLYAO3, ORNDC...</td>\n",
       "      <td>((YLR354C, YBR196C, YOR180C, YOR184W, YJL121C,...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.008878</td>\n",
       "      <td>0.023388</td>\n",
       "      <td>0.178755</td>\n",
       "      <td>0.001878</td>\n",
       "      <td>0.000336</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>37 rows × 9 columns</p>\n",
       "</div>\n",
       "        "
      ],
      "text/plain": [
       "<cameo.strain_design.heuristic.evolutionary_based.OptGeneResult at 0x119e5b9e8>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.4",
   "language": "python",
   "name": "python3.4"
  },
  "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.4.4"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}