documentation_builder/gapfilling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gapfillling\n",
    "\n",
    "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](http://dx.doi.org/10.1371/journal.pcbi.1000308) and [Reed et al. 2006](http://www.pnas.org/content/103/46/17480.short). Cobrapy has a gap filling implementation that is very similar to that of Reed et al. where we use a mixed-integer linear program to figure out the smallest number of reactions that need to be added for a user-defined collection of reactions, i.e. a universal model. Briefly, the problem that we try to solve is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minimize: $$\\sum_i c_i * z_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "subject to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Sv = 0$$\n",
    "$$v^\\star \\geq t$$\n",
    "$$l_i\\leq v_i \\leq u_i$$\n",
    "$$v_i = 0 \\textrm{ if } z_i = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "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 user-defined cost associated with using the *i*th 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import cobra.test\n",
    "from cobra.flux_analysis import gapfill\n",
    "model = cobra.test.create_test_model(\"salmonella\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In this model D-Fructose-6-phosphate is an essential metabolite. We will remove all the reactions using it, and at them to a separate model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "universal = cobra.Model(\"universal_reactions\")\n",
    "for i in [i.id for i in model.metabolites.f6p_c.reactions]:\n",
    "    reaction = model.reactions.get_by_id(i)\n",
    "    universal.add_reaction(reaction.copy())\n",
    "    model.remove_reactions([reaction])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now, because of these gaps, the model won't grow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.optimize().objective_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "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 'no-growth but growth (NGG)' predictions of [Kumar et al. 2009](http://dx.doi.org/10.1371/journal.pcbi.1000308)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GF6PTA\n",
      "F6PP\n",
      "TKT2\n",
      "FBP\n",
      "MAN6PI\n"
     ]
    }
   ],
   "source": [
    "solution = gapfill(model, universal, demand_reactions=False)\n",
    "for reaction in solution[0]:\n",
    "    print(reaction.id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We can obtain multiple possible reaction sets by having the algorithm go through multiple iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---- Run 1 ----\n",
      "GF6PTA\n",
      "F6PP\n",
      "TKT2\n",
      "FBP\n",
      "MAN6PI\n",
      "---- Run 2 ----\n",
      "GF6PTA\n",
      "TALA\n",
      "PGI\n",
      "F6PA\n",
      "MAN6PI\n",
      "---- Run 3 ----\n",
      "GF6PTA\n",
      "F6PP\n",
      "TKT2\n",
      "FBP\n",
      "MAN6PI\n",
      "---- Run 4 ----\n",
      "GF6PTA\n",
      "TALA\n",
      "PGI\n",
      "F6PA\n",
      "MAN6PI\n"
     ]
    }
   ],
   "source": [
    "result = gapfill(model, universal, demand_reactions=False, iterations=4)\n",
    "for i, entries in enumerate(result):\n",
    "    print(\"---- Run %d ----\" % (i + 1))\n",
    "    for e in entries:\n",
    "        print(e.id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We can also instead of using the original objective, specify a given metabolite that we want the model to be able to produce."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FBP\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    model.objective = model.add_boundary(model.metabolites.f6p_c, type='demand')\n",
    "    solution = gapfill(model, universal)\n",
    "    for reaction in solution[0]:\n",
    "        print(reaction.id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Finally, note that using mixed-integer linear programming is computationally quite expensive and for larger models you may want to consider alternative [gap filling methods](http://opencobra.github.io/cobrapy/tags/gapfilling/) and [reconstruction methods](http://opencobra.github.io/cobrapy/tags/reconstruction/)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}