{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computer exercise: Flux balance analysis (FBA) of Escherichia coli core model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Steffen Waldherr (<steffen.waldherr@univie.ac.at>)\n",
    "\n",
    "Summer school on \"Economic principles in cell physiology\", July 2022, Learning Planet Institute, Paris"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Instructions to get started:**\n",
    "- Open <https://colab.research.google.com>.\n",
    "- Choose \"Upload\" in the dialog panel and upload this file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this exercise is to run FBA on the E. coli carbon core model from the [BiGG database](http://bigg.ucsd.edu). The individual steps are as follows:\n",
    "- Download the model as SBML and solve the FBA problem with the cobra toolbox.\n",
    "- Load the network map in escher and visualize the fluxes.\n",
    "- Determine the effects of changing bounds on the uptake reactions.\n",
    "- Compute a production envelope."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main libraries that we are using in this exercise are [cobrapy](https://github.com/opencobra/cobrapy/tree/stable) and [escher](https://escher.github.io/). These can be installed from the Python package index with `pip`. When using Google Colaboratory the installation with `pip` is carried out by the following code block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "iybXAiN0T2Lc",
    "outputId": "8b5204c3-0bc7-45fe-e6ad-25c076f99736"
   },
   "outputs": [],
   "source": [
    "!pip install cobra\n",
    "!pip install escher"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download the SBML model\n",
    "\n",
    "The constraint based model that we are working with in this exercise is the \"E. coli core\" model. It is provided at the [BiGG database](http://bigg.ucsd.edu) as [`e_coli_core`](http://bigg.ucsd.edu/models/e_coli_core). The following code snippet can be run to download the SBML file directly from the database, and save it in the local working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_id = \"e_coli_core\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "QemrsXDKT6R4",
    "outputId": "ed2b1dfe-de8c-422a-a3d0-e9a8225104ce"
   },
   "outputs": [],
   "source": [
    "import urllib.request\n",
    "urllib.request.urlretrieve(\"http://bigg.ucsd.edu/static/models/\"+model_id+\".xml\", model_id + \".xml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the model with cobra and perform FBA\n",
    "\n",
    "Cobrapy can import models in various formats, one of them is SBML. Here we load the previously obtained SBML file in cobrapy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Y_i2i5psVtFn"
   },
   "outputs": [],
   "source": [
    "import cobra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "QOAbLz6BVxCU"
   },
   "outputs": [],
   "source": [
    "ecc = cobra.io.read_sbml_model(model_id + \".xml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next step, the following **tasks** should be performed:\n",
    "\n",
    "1. Run FBA by using the `ecc.optimize()` function of cobrapy. See Section 5.1 [\"Running FBA\"](https://cobrapy.readthedocs.io/en/latest/simulating.html#Running-FBA) in the cobrapy documentation. Assign the result of the optimization to a variable called `solution1`.\n",
    "2. An overview of the optimization result (biomass growth rate, uptake, and secretion fluxes) are shown by the `ecc.summary()` command. Discuss these results with respect to utilized nutrients and byproducts. You can get more details on each metabolite by looking at the `ecc.metabolites.<id>` attribute, where `<id>`should be replaced by the short identifier of the considered metabolite in the model. There is also a `.summary()` method for each metabolite that gives further information on the involved reaction fluxes.\n",
    "3. Why is the result on the glucose uptake reaction exactly 10 mmol/g/h? Look at the reaction information for the glucose exchange reaction and give your interpretation.\n",
    "4. All reaction flux values are stored in the `fluxes` attribute of the optimization result as a [pandas series](https://pandas.pydata.org/docs/reference/api/pandas.Series.html). Use the `pandas.Series.to_csv(<filename>)` command to export all fluxes as csv file. You can open the csv file within jupyter and inspect the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "u_UYLuFJWKyl"
   },
   "outputs": [],
   "source": [
    "# Task 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 436
    },
    "id": "21pHfuPIccfR",
    "outputId": "3caf15e5-ae39-4352-e6e8-7a2935b7a342"
   },
   "outputs": [],
   "source": [
    "# Task 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Looking at basic information for cytosolic atp\n",
    "ecc.metabolites.atp_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Summary of the optimization result for cytosolic atp\n",
    "ecc.metabolites.atp_c.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code picks out the index of the reaction `EX_glc__D_e` in the reaction list and displays the reaction information.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task 3\n",
    "reaction_ids = [ecc.reactions[i].id for i in range(len(ecc.reactions))]\n",
    "print(reaction_ids.index('EX_glc__D_e'))\n",
    "ecc.reactions[51]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Interpretation for the glucose uptake flux*: **TODO**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "iiJwR9NuWb61",
    "outputId": "aa4f411a-2028-49b7-a412-370e947d1993"
   },
   "outputs": [],
   "source": [
    "# Task 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Look at the saved CSV file. Can you identify the fluxes with the largest values?\n",
    "\n",
    "**Answer:** **TODO**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reaction flux illustration with escher\n",
    "\n",
    "To get a better overview of the metabolic flux distribution, it is often helpful to map the optimization result from FBA on a metabolic map. A useful tool for this is [escher](https://escher.github.io/). The following tasks can be carried out in the escher online version at <https://escher.github.io/>. If you have a local installation, you can alternatively uncomment and run the code block below.\n",
    "\n",
    "**Tasks:**\n",
    "1. In the escher GUI, load the map of the E. coli core model. In the online version, this can be started directly (use the Viewer), for the interface on a local installation you need to get the network map as a json file from the escher website and load this within the GUI with the \"Map -> Load map JSON\" menu command.\n",
    "2. Load the FBA flux values from the previously saved csv file with the menu command \"Data -> Load reaction data\".\n",
    "3. Inspect the graphical representation of the flux distribution. Which pathways carrying flux can you identify? What is the functional role of pathways not carrying flux?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 17,
     "referenced_widgets": [
      "8994822945a044d88f37198013fd943a"
     ]
    },
    "id": "HdR9oIXfWi82",
    "outputId": "d5e0c9d7-6cb8-45d8-c5df-94b455b63ec6"
   },
   "outputs": [],
   "source": [
    "# import escher\n",
    "# escher.Builder()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Answer to task 3:*\n",
    "**TODO**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jaoR986ncJk7"
   },
   "source": [
    "## Modifying bounds on exchange reactions\n",
    "\n",
    "The purpose of this step is to determine how the metabolic flux distribution, intracellularly and concerning exchange fluxes, changes when we vary bounds on the uptake reactions and presence of extracellular substrates.\n",
    "\n",
    "Reaction upper and lower bounds can be changed by assigning the `lower_bound` or `upper_bound` attribute of the reaction object. Morover, the cobrapy `Model` class contains an attribute `exchanges` which holds a list of all exchange reactions for easier access.\n",
    "\n",
    "For example, to access the lower bound of the first exchange reaction, you can use the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecc.exchanges[0].lower_bound = 0.0\n",
    "print(ecc.exchanges[0], ecc.exchanges[0].lower_bound)\n",
    "ecc.exchanges[0].lower_bound = -1.0\n",
    "print(ecc.exchanges[0], ecc.exchanges[0].lower_bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task**:\n",
    "Determine how the uptake fluxes and intracellular fluxes change when you apply one of the following modifications:\n",
    "- Nitrogen limitation by increasing the lower bound of NH4 uptake. What is the smallest value for which you expect an effect?\n",
    "- Anaerobic growth by setting the lower bound on oxygen exchange to 0.0.\n",
    "\n",
    "Compare the exchange fluxes and intracellular fluxes (using escher) to the previous solution. Does it match your expectations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulating nitrogen limitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulating anaerobic growth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Production envelopes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code requires the `numpy` package and plotting with `pylab` which are imported here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pylab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, the aim is to compute the the maximum and minimum ethanol production at different growth rates, to get a production envelope for ethanol vs. biomass production.\n",
    "\n",
    "Most of the code is already implemented in the section below, you just have to enter the part where the ethanol flux is minimized (look for `TODO` in the comments)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecc = cobra.io.read_sbml_model(model_id + \".xml\")\n",
    "growth = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85])\n",
    "etoh_min = 0*growth\n",
    "etoh_max = 0*growth\n",
    "for i, mu in enumerate(growth):\n",
    "    with ecc: # any modifications to the model are only retained within this section\n",
    "        ecc.reactions[24].lower_bound = mu # constrain biomass reaction with same lower and upper bound\n",
    "        ecc.reactions[24].upper_bound = mu\n",
    "        ecc.objective = {ecc.reactions.EX_etoh_e: 1}\n",
    "        solution = ecc.optimize()\n",
    "        etoh_max[i] = solution.objective_value\n",
    "        # TODO: Set optimization objective to minimze EX_etoh_e, using -1 instead of 1, and call optimize\n",
    "        etoh_min[i] = solution.objective_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pylab.plot(growth, etoh_max)\n",
    "pylab.plot(growth, etoh_min)\n",
    "pylab.xlabel(\"biomass [1/h]\")\n",
    "pylab.ylabel(\"ethanol [mmol/g/h]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "ecc_fba.ipynb",
   "provenance": []
  },
  "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.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "8994822945a044d88f37198013fd943a": {
     "model_module": "escher",
     "model_module_version": "1.7.3",
     "model_name": "EscherMapModel",
     "state": {
      "_dom_classes": [],
      "_loaded_map_json": null,
      "_loaded_model_json": null,
      "_model_module": "escher",
      "_model_module_version": "1.7.3",
      "_model_name": "EscherMapModel",
      "_view_count": null,
      "_view_module": "escher",
      "_view_module_version": "1.7.3",
      "_view_name": "EscherMapView",
      "allow_building_duplicate_reactions": null,
      "and_method_in_gene_reaction_rule": null,
      "canvas_size_and_loc": null,
      "cofactors": null,
      "disabled_buttons": null,
      "embedded_css": null,
      "enable_editing": null,
      "enable_keys": false,
      "enable_keys_with_tooltip": null,
      "enable_search": null,
      "enable_tooltips": false,
      "full_screen_button": {
       "enable_editing": true,
       "enable_keys": true,
       "enable_tooltips": [
        "label"
       ],
       "menu": "all",
       "scroll_behavior": "pan"
      },
      "gene_data": null,
      "gene_font_size": null,
      "height": 500,
      "hide_all_labels": null,
      "hide_secondary_metabolites": null,
      "highlight_missing": null,
      "identifiers_on_map": null,
      "layout": "IPY_MODEL_1bb6b97bf098480eba5c8f6ab9d7949a",
      "marker_radius": null,
      "menu": null,
      "metabolite_compare_style": null,
      "metabolite_data": null,
      "metabolite_no_data_color": null,
      "metabolite_no_data_size": null,
      "metabolite_scale": null,
      "metabolite_scale_preset": null,
      "metabolite_styles": null,
      "never_ask_before_quit": null,
      "primary_metabolite_radius": null,
      "reaction_compare_style": null,
      "reaction_data": null,
      "reaction_no_data_color": null,
      "reaction_no_data_size": null,
      "reaction_scale": null,
      "reaction_scale_preset": null,
      "reaction_styles": null,
      "scroll_behavior": "none",
      "secondary_metabolite_radius": null,
      "semantic_zoom": null,
      "show_gene_reaction_rules": null,
      "starting_reaction": null,
      "use_3d_transform": null,
      "zoom_to_element": null
     }
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}