2001_QuantMiBi1_CobraPy_SimVal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "dPnMzwQ_IWwZ"
   },
   "source": [
    "# Seminar Quantitativ Microbiology \n",
    "# Metabolite Production with Strain Engineering\n",
    "\n",
    "## Introduction\n",
    "\n",
    "\n",
    "For a overview on Jupyter notebooks read [this review](https://www.nature.com/articles/d41586-018-07196-1). \n",
    "\n",
    "---\n",
    "\n",
    "## Tutorial Steps\n",
    "  * Set up of Python environment\n",
    "    * Basic libraries(sys, pandas, numpy, matplotlib, zipfile, cobrapy)\n",
    "  * Analysis of Genome Scale Metabolic Model\n",
    "    * Retrieval of GSMM for *P. pastoris*\n",
    "    * Flux variability of exchange reactions\n",
    "    * Minimal medium composition\n",
    "  * Experimental growth rate reproduction\n",
    "    * Familiarizing with biomass composition reactions\n",
    "    * Defining functions for correct biomass equation switch\n",
    "    * Data retrieval\n",
    "    * Simulation loop\n",
    "    * Graphical output\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XRP3j145RXOi"
   },
   "source": [
    "## Set-up compute environment\n",
    "\n",
    "Before we can analyse GSMM, we have adjust the python environment that it integrates the cobrapy toolbox and downloading the GSMM."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ewTIEAb4Rb2G"
   },
   "source": [
    "### Basic Python libraries \n",
    "Some libraries that facilitate data manipulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6wxIUGOEReDX"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "System ready\n"
     ]
    }
   ],
   "source": [
    "import sys # loading commands to control/navigate within the system architecture\n",
    "# Loading pandas, a library for data manipulation\n",
    "from os.path import join\n",
    "import xlrd\n",
    "import pandas as pd\n",
    "# import lxml\n",
    "\n",
    "# Loading numpy, a library fo manipulation of numbers\n",
    "import numpy as np\n",
    "\n",
    "# loading matplotlib, a library for visualization\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# loading cobrapy, a library dedicated to the analysis of genome scale metabolic models\n",
    "from cobra.io import read_sbml_model, write_sbml_model, load_matlab_model\n",
    "\n",
    "# loading escher for metabolic network visualization\n",
    "import escher\n",
    "from escher import Builder\n",
    "from time import sleep\n",
    "escher.rc['never_ask_before_quit'] = True\n",
    "# list of available maps\n",
    "# print(escher.list_available_maps())\n",
    "\n",
    "# loading Memote, quality assessment of GSMM\n",
    "from memote import test_model, snapshot_report\n",
    "\n",
    "print('System ready')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-08-31 16:18:28--  http://bigg.ucsd.edu/static/models/e_coli_core.xml\n",
      "Resolving bigg.ucsd.edu (bigg.ucsd.edu)... 169.228.33.117\n",
      "Connecting to bigg.ucsd.edu (bigg.ucsd.edu)|169.228.33.117|:80... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 707166 (691K) [application/xml]\n",
      "Saving to: ‘e_coli_core.xml.1’\n",
      "\n",
      "e_coli_core.xml.1   100%[===================>] 690.59K   734KB/s    in 0.9s    \n",
      "\n",
      "2020-08-31 16:18:30 (734 KB/s) - ‘e_coli_core.xml.1’ saved [707166/707166]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget http://bigg.ucsd.edu/static/models/e_coli_core.xml\n",
    "\n",
    "# loading a visualization file of the metabolic network. \n",
    "# For frequently used models, like iJO1366, Escher automatically retrieves the visualization file.  \n",
    "#!wget http://bigg.ucsd.edu/static/models/iJO1366.json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "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>0x07fb0580f1b90</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_w_GAM - 1.0*BIOMASS_Ecoli_core_w_GAM_reverse_712e5</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Compartments</strong></td>\n",
       "                <td>extracellular space, cytosol</td>\n",
       "            </tr>\n",
       "          </table>"
      ],
      "text/plain": [
       "<Model e_coli_core at 0x7fb0580f1b90>"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# generating cobra variable from SBML/xml/mat file\n",
    "myFile = 'e_coli_core.xml'\n",
    "model = read_sbml_model(myFile)\n",
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Jn805mCObF2k"
   },
   "source": [
    "### Flux variability of exchange reactions\n",
    "\n",
    "Flux balance analysis provides a single optimal solution. Mostly, there exist a number of alternative flux distributions around the optimum, which can be physiologically relevant. To identify the variability of exchange fluxes around the optimum solution 'flux variability analysis' can be performed ([Mahadevan & Schilling, 2003](http://dx.doi.org/10.1016/j.ymben.2003.09.002)). Use the following command to identify minimum and maximum flux ranges of the model for the exchange reactions.\n",
    "\n",
    "model.summary(fva=.95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2HiWV97jWhKI"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"2\" halign=\"left\">IN_FLUXES</th>\n",
       "      <th colspan=\"2\" halign=\"left\">OUT_FLUXES</th>\n",
       "      <th colspan=\"2\" halign=\"left\">OBJECTIVES</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>ID</th>\n",
       "      <th>FLUX</th>\n",
       "      <th>ID</th>\n",
       "      <th>FLUX</th>\n",
       "      <th>ID</th>\n",
       "      <th>FLUX</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>o2_e</td>\n",
       "      <td>21.799493</td>\n",
       "      <td>h2o_e</td>\n",
       "      <td>29.175827</td>\n",
       "      <td>BIOMASS_Ecoli_core_w_GAM</td>\n",
       "      <td>0.873922</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>glc__D_e</td>\n",
       "      <td>10.000000</td>\n",
       "      <td>co2_e</td>\n",
       "      <td>22.809833</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>nh4_e</td>\n",
       "      <td>4.765319</td>\n",
       "      <td>h_e</td>\n",
       "      <td>17.530865</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>pi_e</td>\n",
       "      <td>3.214895</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "<cobra.core.summary.model_summary.ModelSummary at 0x7fb057cf3950>"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution = model.optimize()\n",
    "model.summary() # fva=.95 additional argument specifies allowed deviation from the optimum \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading Map from https://escher.github.io/1-0-0/6/maps/Escherichia%20coli/e_coli_core.Core%20metabolism.json\n"
     ]
    }
   ],
   "source": [
    "# Load an Escher map\n",
    "builder.map_name = 'e_coli_core.Core metabolism'\n",
    "builder.model = model\n",
    "builder.reaction_data = solution.fluxes\n",
    "# Add some data for metabolites\n",
    "#builder.metabolite_data = solution.shadow_prices\n",
    "# Simplify the map by hiding some labels\n",
    "builder.hide_secondary_metabolites = True\n",
    "#builder.hide_all_labels = True\n",
    "builder.reaction_scale = [\n",
    "    { 'type': 'min', 'color': '#000000', 'size': 12 },\n",
    "    { 'type': 'median', 'color': '#ffffff', 'size': 20 },\n",
    "    { 'type': 'max', 'color': '#ff0000', 'size': 25 }\n",
    "]\n",
    "builder.reaction_scale_preset = 'GaBuRd'\n",
    "\n",
    "# Make all the arrows three times as thick\n",
    "builder.reaction_scale = [\n",
    "    {k: v * 3 if k == 'size' else v for k, v in x.items()}\n",
    "    for x in builder.reaction_scale\n",
    "]\n",
    "builder.save_html('escher_map_file.html')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cvSOmwrZbNIG"
   },
   "source": [
    "### Malate Export Reaction\n",
    "coding for malate export reaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mEHqAD7nbL8n"
   },
   "outputs": [],
   "source": [
    "# 1. generating random genome with defined GC-content\n",
    "Genome_Sequence, Genome_Genes = generat_Genome(Genome_Size, Genome_GGcont, Model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_Genome(Genome_Size, Genome_GGcont, Model):\n",
    "    '''\n",
    "    Function to generate random genome to map promoter activities to reactions in a cobrapy model.\n",
    "    '''\n",
    "    return''.join(random.choices(letters, k=length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATGGATCCCCGTTGA\n"
     ]
    }
   ],
   "source": [
    "Model = model\n",
    "# coding sequence construction\n",
    "# first we determine the minimum coding gene length of nucleotides to distinguish the enzymes in the model\n",
    "Enzyme_Number = len(Model.reactions)\n",
    "Gen_Minimum = np.ceil(np.log2(Enzyme_Number))\n",
    "\n",
    "# we want to represent codon triplicates, we calculate the next highest divisor of three\n",
    "Gene_Length = int(np.ceil(Gen_Minimum/3))\n",
    "CodonFile = 'CodonTriplets.csv'\n",
    "CodonTriplets = pd.read_csv(CodonFile, delimiter=';', skipinitialspace=True)\n",
    "CodonStop = CodonTriplets[['Stop' in s for s in CodonTriplets['Name']]].reset_index()\n",
    "CodonCoding = CodonTriplets.drop(CodonStop.index).reset_index()\n",
    "Gene_Coding = [random.choices(CodonCoding['Triplet'], weights=CodonCoding['Percent']) for CodonId in range(Gene_Length)]\n",
    "Gene_Stop = random.choices(CodonStop['Triplet'], weights=CodonStop['Percent'])\n",
    "Gene_Coding.append(Gene_Stop)\n",
    "Gene_Coding.insert(0,['ATG'])\n",
    "Gene_Coding = ''.join([Letter for Nest in Gene_Coding for Letter in Nest])\n",
    "print(Gene_Coding)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'CTCAAAACTTCTCCCCGTTC'"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# background genome construction\n",
    "Genome_Size = 20\n",
    "GC_cont = .5\n",
    "\n",
    "SeqNestList = [random.choices([Letter for Nest in random.choices([['G','C'],['A','T']], weights=[GC_cont, 1-GC_cont]) for Letter in Nest]) for _x in range(Genome_Size)]\n",
    "Genome_Sequence = ''.join(Letter for Nest in SeqNestList for Letter in Nest)\n",
    "Genome_Sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Int64Index([ 0,  1,  2,  3,  4,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,\n",
       "            20, 21, 22, 23, 24, 25, 26, 27, 30, 31, 32, 35, 36, 37, 38],\n",
       "           dtype='int64')"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pickle\n",
    "import os\n",
    "# generating promoter sequence\n",
    "WeightFile = os.path.join('Models','NucleotideWeightTable.pkl')\n",
    "with open(WeightFile, 'rb') as handle:\n",
    "    Nucleotides_Weight = pickle.load(handle)\n",
    "RefFile = os.path.join('Models','RefSeq.txt')\n",
    "with open(RefFile) as f:\n",
    "    RefSeq = f.read()\n",
    "Nucleotides_Weight.index = Nucleotides_Weight.index+40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'G'"
      ]
     },
     "execution_count": 239,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myAr = np.asarray([Letter for Letter in RefSeq])\n",
    "myAr[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {},
   "outputs": [],
   "source": [
    "Bases = ['A','C','G','T']\n",
    "for idx in range(len(Nucleotides_Weight)):\n",
    "#     print(Nucleotides_Weight.index[idx], Nucleotides_Weight.iloc[idx].values)\n",
    "#     print(random.choices(Bases,Nucleotides_Weight.iloc[idx].values))\n",
    "#     print(Nucleotides_Weight.index[idx])\n",
    "    np.put(myAr,Nucleotides_Weight.index[idx],random.choices(Bases,Nucleotides_Weight.iloc[idx].values))\n",
    "#     myAr[idx] = random.choices(Bases,Nucleotides_Weight.iloc[idx].values)\n",
    "\n",
    "# type(idx)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['T', 'T', 'C', 'C', 'A', 'T', 'T', 'G', 'T', 'T', 'A', 'T', 'A',\n",
       "       'A', 'G', 'G', 'A', 'G', 'G', 'C', 'A', 'G', 'T', 'A', 'C', 'T',\n",
       "       'A', 'A', 'T', 'A', 'T', 'C', 'T', 'T', 'T', 'G', 'C', 'G', 'G',\n",
       "       'G'], dtype='<U1')"
      ]
     },
     "execution_count": 241,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myAr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xpFf1wR0QX4L"
   },
   "source": [
    "## Prediction of growth rates\n",
    "### Choose the Substrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myhost.Choose_Substrate(None)\n",
    "myhost.show_BiotechSetting()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JsDpC6gVQigc"
   },
   "source": [
    "### Familiarizing with biomass composition reactions\n",
    "\n",
    "Microorganisms adapt to their substrate. Different substrates provide different energy content and require different cellular resources to become metabolized. In GSMM these differences may be represented by different equations/reactions for the substrates. In iMR1026 for *P. pastoris* there are various biomass equations for glucose, glycerol, glucose-glycerol mixtures, and methanol. When simulating a model, we have to make sure we use the right biomass equation fitting with the substrate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DBNmaT-fRpQz"
   },
   "outputs": [],
   "source": [
    "# List of all reactions with 'BIOMASS' in their name\n",
    "model.reactions.query('?..?')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "82QMmzXYR4br"
   },
   "outputs": [],
   "source": [
    "# Looking in detail to biomass with methanol; use the reaction name given to you from the list of all biomass reactions;\n",
    "model.reactions.?..?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N_eSuR-EyQ-M"
   },
   "source": [
    "### Defining functions for correct biomass equation switch\n",
    "\n",
    "For each substrate, the boundary exchange fluxes are activated and the reactions of competing substrates are disabled.\n",
    "\n",
    "Remember to write a specific function to adapt the model for glucose utilization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XqM2cco0Qpn0"
   },
   "outputs": [],
   "source": [
    "def AdaptMethanol(model, meoh_up):\n",
    "  model.objective = 'BIOMASS_meoh'\n",
    "  # setting uptake reactions right\n",
    "  model.reactions.Ex_glc_D.lower_bound = 0\n",
    "  model.reactions.Ex_glyc.lower_bound = 0\n",
    "  model.reactions.ATPM.lower_bound = .4\n",
    "  model.reactions.Ex_meoh.lower_bound = -np.abs(meoh_up)\n",
    "  # setting additional biomass composition\n",
    "  model.reactions.LIPIDS_meoh.upper_bound = 1000\n",
    "  model.reactions.PROTEINS_meoh.upper_bound = 1000\n",
    "  model.reactions.STEROLS_meoh.upper_bound = 1000\n",
    "  model.reactions.BIOMASS_meoh.upper_bound = 1000\n",
    "  # deactivating Glyc-based biomass composition\n",
    "  model.reactions.LIPIDS_glyc.upper_bound = 0\n",
    "  model.reactions.PROTEINS_glyc.upper_bound = 0\n",
    "  model.reactions.STEROLS_glyc.upper_bound = 0\n",
    "  model.reactions.BIOMASS_glyc.upper_bound = 0  \n",
    "  # deactivating Glc-based biomass composition\n",
    "  model.reactions.LIPIDS.upper_bound = 0\n",
    "  model.reactions.PROTEINS.upper_bound = 0\n",
    "  model.reactions.STEROLS.upper_bound = 0\n",
    "  model.reactions.BIOMASS.upper_bound = 0  \n",
    "  return model\n",
    "\n",
    "def AdaptGlycerol(model, glyc_up):\n",
    "  model.objective = 'BIOMASS_glyc'\n",
    "  # setting uptake reactions right\n",
    "  model.reactions.Ex_meoh.lower_bound = 0\n",
    "  model.reactions.Ex_glc_D.lower_bound = 0\n",
    "  model.reactions.ATPM.lower_bound = 2.5\n",
    "  model.reactions.Ex_glyc.lower_bound = -np.abs(glyc_up)\n",
    "  # setting additional biomass composition\n",
    "  model.reactions.LIPIDS_glyc.upper_bound = 1000\n",
    "  model.reactions.PROTEINS_glyc.upper_bound = 1000\n",
    "  model.reactions.STEROLS_glyc.upper_bound = 1000\n",
    "  model.reactions.BIOMASS_glyc.upper_bound = 1000  \n",
    "  # deactivating MeOH-based biomass composition\n",
    "  model.reactions.LIPIDS_meoh.upper_bound = 0\n",
    "  model.reactions.PROTEINS_meoh.upper_bound = 0\n",
    "  model.reactions.STEROLS_meoh.upper_bound = 0\n",
    "  model.reactions.BIOMASS_meoh.upper_bound = 0\n",
    "  # deactivating Glc-based biomass composition\n",
    "  model.reactions.LIPIDS.upper_bound = 0\n",
    "  model.reactions.PROTEINS.upper_bound = 0\n",
    "  model.reactions.STEROLS.upper_bound = 0\n",
    "  model.reactions.BIOMASS.upper_bound = 0\n",
    "  return model\n",
    "\n",
    "def AdaptGlucose(model, glc_up):\n",
    "  # setting uptake reactions right\n",
    "  model.reactions.Ex_meoh.lower_bound = 0\n",
    "  model.reactions.Ex_glyc.lower_bound = 0\n",
    "  model.reactions.Ex_glc_D.lower_bound = -np.abs(glc_up)\n",
    "  # setting additional biomass composition\n",
    "  model.reactions.LIPIDS.upper_bound = 1000\n",
    "  model.reactions.PROTEINS.upper_bound = 1000\n",
    "  model.reactions.STEROLS.upper_bound = 1000\n",
    "  model.reactions.BIOMASS.upper_bound = 1000  \n",
    "  # deactivating Glyc-based biomass composition\n",
    "  model.reactions.LIPIDS_glyc.upper_bound = 0\n",
    "  model.reactions.PROTEINS_glyc.upper_bound = 0\n",
    "  model.reactions.STEROLS_glyc.upper_bound = 0\n",
    "  model.reactions.BIOMASS_glyc.upper_bound = 0  \n",
    " # deactivating meoh-based biomass composition\n",
    "  model.reactions.LIPIDS_meoh.upper_bound = 0\n",
    "  model.reactions.PROTEINS_meoh.upper_bound = 0\n",
    "  model.reactions.STEROLS_meoh.upper_bound = 0\n",
    "  model.reactions.BIOMASS_meoh.upper_bound = 0\n",
    "  model.objective = 'BIOMASS'\n",
    "  return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "K045st1gSllg"
   },
   "source": [
    "### Data retrieval\n",
    "\n",
    "For evaluation of the growth rate prediction of the *P. pastoris* model we use experimental data from the closely related organism *Ogataea polymorpha*. The measurements in the table are extracted from [van Dijken et al. 1976](https://dx.doi.org/10.1007/bf00446560) for methanol and from [de Koning et al., 1987](https://dx.doi.org/10.1007/BF00456710) and [Moon et al., 2003](https://dx.doi.org/10.1385/ABAB:111:2:65).\n",
    "\n",
    "Data Address: [here](https://rwth-aachen.sciebo.de/s/o72jwWQWh3ame1e/download)\n",
    "\n",
    "**When using Binder:**\n",
    "Click on \"here\" and save the file. Go back to \"Home\"-Tab in the browser afterwards and navigate back to the root directory of your cloud system (i.e. you see the Jupyter notebook files). Click on the \"Upload\" button on the right side of the window and select the file you just saved. Now you should see the file \"Opol-expt-grwth_MeOH-Glyc.csv\", click again on \"Upload\" in the line of the corresponding file. Now go back to the simulation-Tab. Replace '?..?' in the cell below by entering the name of the csv file you downloaded. Then press Ctrl+Enter to load the experimental data.\n",
    "\n",
    "**When using Colab:**\n",
    "Click on \"here\" and save the file. Then click on \"Files\" on the left side of the window. Upload the file you have just saved by clicking on \"Upload\" and selecting this file. Now you should see the file \"Opol-expt-grwth_MeOH-Glyc.csv\" in the \"Files\" section. If this is not the case, click on \"Refresh\". Then right-click on this file and select \"Copy path\". Replace '?..?' in the cell below with the file by pressing Ctrl+V. Then evaluate the cell by Ctrl+Enter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mfOlv0uiCNvI"
   },
   "outputs": [],
   "source": [
    "# Create Excel file with all uptake rates and then import it instead as follows:\n",
    "# data = pd.read_excel('nameOfFile.xlsx', sheet_name='Ecol', index_col=0)\n",
    "\n",
    "data = pd.read_csv('?..?', delimiter=',|;', engine='python')\n",
    "data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aMLKO9MNe1fq"
   },
   "source": [
    "Adding measurements for glucose from [Lehnen et al., 2017](https://doi.org/10.1016/j.meteno.2017.07.001) to our basic data table. The required information is the growth rate on glucose for *H. polymorpha*. Extract the necessary information from Table 3 in the article. The reaction name for the 'Exchange'-reaction is 'Ex_glc_D', replace this in the corresponding position of '?..?' in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4mTMJvd0ZHFv"
   },
   "outputs": [],
   "source": [
    "data = data.append({'Substrate':'Glucose', 'Exchange':'?..?', 'uptake rate (mmol/gCDW/h)':?..?, 'growth rate (/h)':?..?, 'source':'Lehnen et al.'}, ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Wxe3Csg0TtYV"
   },
   "source": [
    "### Simulation loop\n",
    "For-Loop over all experimental data points.\n",
    "\n",
    "Remember to add a decision when you include glucose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "TPyYPHrXUMs1"
   },
   "outputs": [],
   "source": [
    "# At the beginning a substrate was selected, therefore only one execution is necessary, approximately as follows:\n",
    "# Substrate = myhost.var_Substrate\n",
    "# model=AdaptSubstrate(model, data.loc[['Substrate'],['uptake rate']].values)\n",
    "# growth_simulated = model.optimize()\n",
    "\n",
    "growth_simulated = [];\n",
    "# test_model = model.copy()\n",
    "# iteration over all rows in 'data'\n",
    "for index, row in data.iterrows():\n",
    "  with model as test_model:\n",
    "    print(index) # printing the row number to get feedback that everything is working\n",
    "    # selecting the right substrate in the model based on 'Substrate' in 'data'\n",
    "    if row['Substrate'] == 'Methanol':\n",
    "      model = AdaptMethanol(test_model, row['uptake rate (mmol/gCDW/h)'])\n",
    "    elif row['Substrate'] == 'Glycerol':\n",
    "      model = AdaptGlycerol(test_model, row['uptake rate (mmol/gCDW/h)'])\n",
    "    elif row['Substrate'] == 'Glucose':\n",
    "      test_model = AdaptGlucose(test_model, row['uptake rate (mmol/gCDW/h)'])\n",
    "    else:\n",
    "      print('substrate not considered')      \n",
    "  #     model.optimize()\n",
    "    growth_simulated.append(model.slim_optimize())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aN5SDn5yWJER"
   },
   "source": [
    "### Graphical output\n",
    "\n",
    "Remember to add glucose to the visualization. Add the right axis labels, and a file name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1DG-7h3qWLRi"
   },
   "outputs": [],
   "source": [
    "plt.scatter(data['growth rate (/h)'][0:7], growth_simulated[0:7], s=50, c='k', marker='o');\n",
    "plt.scatter(data['growth rate (/h)'][7], growth_simulated[7], s=50, c='k', marker='s');\n",
    "plt.scatter(data['growth rate (/h)'][8], growth_simulated[8], s=50, c='k', marker='d');\n",
    "plt.xlabel('?..?');\n",
    "plt.ylabel('?..?');\n",
    "myline = np.linspace(0,np.max(growth_simulated),10);\n",
    "plt.plot(myline,myline,'k--');\n",
    "plt.title('Growth rate comparison');\n",
    "plt.legend(['Optimum', 'Methanol (van Dijken)', 'Glycerol (deKoning, Moon)', 'Glucose (Lehnen et al.)'], loc=2);\n",
    "plt.style.use('seaborn-paper')\n",
    "\n",
    "# Saving figure\n",
    "plt.savefig('?..?.png')"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "2001_QuantMiBi1_CobraPy_SimVal.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}