cobrapy-05-advanced-constraints-objectives.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tailored constraints, variables and objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thanks to the use of symbolic expressions via the optlang mathematical modeling package, it is relatively straight-forward to add new variables, constraints and advanced objectives that can not easily be formulated as a combination of different reactions and their corresponding upper and lower bounds. Here we demonstrate this optlang functionality which is exposed via `model.problem`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from cobra.io import read_sbml_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = read_sbml_model('data/e_coli_core.xml.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "same_flux = model.problem.Constraint(\n",
    "    model.reactions.FBA.flux_expression - model.reactions.NH4t.flux_expression,\n",
    "    lb=0,\n",
    "    ub=0)\n",
    "model.add_cons_vars(same_flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0*FBA - 1.0*FBA_reverse_84806"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.reactions.FBA.flux_expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I can maximize growth rate whilst the fluxes of reactions 'FBA' and 'NH4t' are constrained to be (near) identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.66274904774 4.66274904774 0.8551109609261567\n"
     ]
    }
   ],
   "source": [
    "solution = model.optimize()\n",
    "print(solution.fluxes.FBA, solution.fluxes.NH4t,\n",
    "      solution.objective_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objectives\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple objective such as the maximization of the flux through one or more reactions can conveniently be done by simply \n",
    "assigning to the `model.objective` property as we have seen in previous chapters, e.g.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = read_sbml_model('data/e_coli_core.xml.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8739215069684307\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    model.objective = {model.reactions.BIOMASS_Ecoli_core_w_GAM: 1}\n",
    "    model.optimize()\n",
    "    print(model.reactions.BIOMASS_Ecoli_core_w_GAM.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objectives mathematical expression is seen by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0*BIOMASS_Ecoli_core_w_GAM_reverse_712e5 + 1.0*BIOMASS_Ecoli_core_w_GAM"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.objective.expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the example problem:\n",
    "\n",
    "> **min** $\\frac{1}{2}\\left(x^2 + y^2 \\right) - y$\n",
    "\n",
    "> *subject to*\n",
    "\n",
    "> $x + y = 2$\n",
    "\n",
    "> $x \\ge 0$\n",
    "\n",
    "> $y \\ge 0$\n",
    "\n",
    "This (admittedly very artificial) problem can be visualized graphically where the optimum is indicated by the blue dot on the line of feasible solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.solver = 'cplex'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sum_two = model.problem.Constraint(\n",
    "    model.reactions.FBA.flux_expression + model.reactions.NH4t.flux_expression,\n",
    "    lb=2,\n",
    "    ub=2)\n",
    "model.add_cons_vars(sum_two)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we add the quadratic objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "quadratic_objective = model.problem.Objective(\n",
    "    0.5 * model.reactions.NH4t.flux_expression**2 + 0.5 *\n",
    "    model.reactions.FBA.flux_expression**2 -\n",
    "    model.reactions.FBA.flux_expression,\n",
    "    direction='min')\n",
    "model.objective = quadratic_objective\n",
    "solution = model.optimize(objective_sense=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5 1.5\n"
     ]
    }
   ],
   "source": [
    "print(solution.fluxes.NH4t, solution.fluxes.FBA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = read_sbml_model('data/e_coli_core.xml.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "difference = model.problem.Variable('difference')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use constraints to define what values this variable shall take"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "constraint = model.problem.Constraint(\n",
    "    model.reactions.EX_glc__D_e.flux_expression -\n",
    "    model.reactions.EX_nh4_e.flux_expression - difference,\n",
    "    lb=0,\n",
    "    ub=0)\n",
    "model.add_cons_vars([difference, constraint])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can access that difference directly during our knock-out exploration by looking at its primal value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5.234680806802543\n",
      "-5.234680806802543\n",
      "-5.234680806802543\n",
      "-1.8644444444444437\n",
      "-1.8644444444444437\n"
     ]
    }
   ],
   "source": [
    "for reaction in model.reactions[:5]:\n",
    "    with model:\n",
    "        reaction.knock_out()\n",
    "        model.optimize()\n",
    "        print(model.solver.variables.difference.primal)"
   ]
  }
 ],
 "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"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "98px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": "3",
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}