EcoliPPPnoPhi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Note: this is the EcoliPPP.ipynb notebook. The\n",
    "PDF version \"The Escherichia coli Core Model: \n",
    "Modular Energetic Bond Graph Analysis of Pentose Phosphate Pathways\"\n",
    "is available [here](EcoliPPPnoPhi.pdf).***\n",
    "\n",
    "- This is the version without reaction potentials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "The Network Thermodynamics/Bond Graph approach of \n",
    "<cite data-cite=\"OstPerKat71,OstPerKat73\">Oster, Perlson and Katchalsky (1971,1973)</cite>, extended by <cite data-cite=\"GawCra14,GawCra16,GawCra17\">Gawthrop and Crampin (2014,2016,2017)</cite>,\n",
    "to modelling biomolecular systems developed independently from the stoichiometric approach \n",
    "<cite data-cite=\"Pal06,Pal11,Pal15\"></cite>.\n",
    "\n",
    "However, the conceptual point of intersection of the two approaches is the fact that the stoichiometric matrix is the modulus of the conceptual multiport transformer linking reactions to species.\n",
    "This was pointed out by <cite data-cite=\"CelGre09\">Cellier and Greifeneder (2009)</cite>. This means that the two approaches are complementary and each can build on the strengths of the other.\n",
    "\n",
    "In particular, as discussed here, the Bond Graph approach adds energy to stoichiometry.\n",
    "\n",
    "This notebook focuses on building modular models of metabolism and consequent pathway analysis based on the Escherichia coli Core Model <cite data-cite=\"OrtFlePal10\">(Orth, Fleming and Palsson,2010)</cite>; in particular, the Glycolysis and Pentose Phosphate portion is extracted and analysed. Following the discussion in the textbook of \n",
    "<cite data-cite=\"GarGri17\">Garrett and Grisham (2017)</cite>, section 22.6d, various possible pathways are examined by choosing appropriate chemostats and flowstats.\n",
    "<cite data-cite=\"GawCra18\">(Gawthrop and Crampin, 2018)</cite>\n",
    "\n",
    "Assuming steady-state conditions, the corresponding pathway potentials <cite data-cite=\"Gaw17a\">(Gawthrop 2017)</cite> are derived.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import some python code\n",
    "The bond graph analysis uses a number of Python modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Maths library\n",
    "import numpy as np\n",
    "import scipy\n",
    "## BG tools\n",
    "import BondGraphTools as bgt\n",
    "\n",
    "## BG stoichiometric utilities\n",
    "import stoich as st\n",
    "\n",
    "## Stoichiometric conversion\n",
    "import CobraExtract as Extract\n",
    "import stoichBondGraph as stbg\n",
    "\n",
    "## Potentials\n",
    "import phiData\n",
    "\n",
    "## Faraday constant\n",
    "import scipy.constants as con\n",
    "F = con.physical_constants['Faraday constant'][0]\n",
    "\n",
    "## Display\n",
    "import IPython.display as disp\n",
    "\n",
    "import copy\n",
    "\n",
    "## Allow output from within functions\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\"\n",
    "\n",
    "import importlib as imp\n",
    "\n",
    "## Units etc\n",
    "factor = 1\n",
    "units = ['mV','kJ']\n",
    "\n",
    "## Control output\n",
    "quiet = True\n",
    "computePhi = True\n",
    "showMu = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extract the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract full ecoli core model from the CobraPy representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracting stoichiometric matrix from: textbook\n",
      "Cobra Model name: e_coli_core BondGraphTools name: e_coli_core_abg\n",
      "Extract.Integer only handles one non-integer per reaction\n",
      "Multiplying reaction BIOMASS_ECOLIORE ( 12 ) by 0.6684491978609626 to avoid non-integer species 3PG ( 2 )\n",
      "Multiplying reaction CYTBD ( 15 ) by 2.0 to avoid non-integer species O2 ( 55 )\n",
      "Multiplying reaction PGK ( 54 ) by -1\n",
      "Multiplying reaction PGM ( 56 ) by -1\n",
      "Multiplying reaction RPI ( 65 ) by -1\n",
      "Multiplying reaction SUCOAS ( 69 ) by -1\n",
      "['ACALD', 'ACALDT', 'ACKR', 'ACONTA', 'ACONTB', 'ACT2R', 'ADK1', 'AKGDH', 'AKGT2R', 'ALCD2X', 'ATPM', 'ATPS4R', 'BIOMASS_ECOLIORE', 'CO2T', 'CS', 'CYTBD', 'D_LACT2', 'ENO', 'ETOHT2R', 'FBA', 'FBP', 'FORT2', 'FORTI', 'FRD7', 'FRUPTS2', 'FUM', 'FUMT2_2', 'G6PDH2R', 'GAPD', 'GLCPTS', 'GLNS', 'GLNABC', 'GLUDY', 'GLUN', 'GLUSY', 'GLUT2R', 'GND', 'H2OT', 'ICDHYR', 'ICL', 'LDH_D', 'MALS', 'MALT2_2', 'MDH', 'ME1', 'ME2', 'NADH16', 'NADTRHD', 'NH4T', 'O2T', 'PDH', 'PFK', 'PFL', 'PGI', 'PGK', 'PGL', 'PGM', 'PIT2R', 'PPC', 'PPCK', 'PPS', 'PTAR', 'PYK', 'PYRT2', 'RPE', 'RPI', 'SUCCT2_2', 'SUCCT3', 'SUCDI', 'SUCOAS', 'TALA', 'THD2', 'TKT1', 'TKT2', 'TPI']\n"
     ]
    }
   ],
   "source": [
    "sm = Extract.extract(cobraname='textbook',Remove=['_C','__' ], \n",
    "                     negReaction=['RPI','PGK','PGM','SUCOAS'], quiet=quiet)\n",
    "print(sm['reaction'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract Glycolysis, Pentose Phosphate Pathways and TCA (using PDH and PDH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "name = 'GlyPPP_abg'\n",
    "reaction = []\n",
    "\n",
    "## Glycolysis\n",
    "reaction += ['PGI','PFK','FBA','TPI']\n",
    "\n",
    "## Pentose Phosphate\n",
    "reaction += ['G6PDH2R','PGL','GND','RPI','TKT2','TALA','TKT1','RPE']\n",
    "\n",
    "## Create submodel\n",
    "sGlyPPP = Extract.choose(sm,reaction=reaction)\n",
    "\n",
    "sGlyPPP['name'] = name\n",
    "stbg.model(sGlyPPP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create stoichiometry\n",
    "import GlyPPP_abg\n",
    "S = st.stoich(GlyPPP_abg.model(),quiet=quiet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display the extracted reactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{G6P &<>[ PGI ] F6P }\\\\\n",
       "\\ch{ATP + F6P &<>[ PFK ] ADP + FDP + H }\\\\\n",
       "\\ch{FDP &<>[ FBA ] DHAP + G3P }\\\\\n",
       "\\ch{DHAP &<>[ TPI ] G3P }\\\\\n",
       "\\ch{G6P + NADP &<>[ G6PDH2R ] 6PGL + H + NADPH }\\\\\n",
       "\\ch{6PGL + H2O &<>[ PGL ] 6PGC + H }\\\\\n",
       "\\ch{6PGC + NADP &<>[ GND ] CO2 + NADPH + RU5PD }\\\\\n",
       "\\ch{RU5PD &<>[ RPI ] R5P }\\\\\n",
       "\\ch{E4P + XU5PD &<>[ TKT2 ] F6P + G3P }\\\\\n",
       "\\ch{G3P + S7P &<>[ TALA ] E4P + F6P }\n",
       "\\end{align}\n",
       "\n",
       "\\begin{align}\n",
       "\\ch{R5P + XU5PD &<>[ TKT1 ] G3P + S7P }\\\\\n",
       "\\ch{RU5PD &<>[ RPE ] XU5PD }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "disp.Latex(st.sprintrl(sGlyPPP,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code to analyse pathways defined by chemostats and flowstats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Analyse pathways defined by chemostats and flowstats\n",
    "def ch(name):\n",
    "    return '\\\\ch{'+name+'}'\n",
    "\n",
    "\n",
    "def pathway(bg,phi,chemostats,flowstats=[],computePhi=False,verbose=False):\n",
    "    \"\"\" Analyse pathways\n",
    "    \"\"\"\n",
    "    \n",
    "    print('Chemostats:',sorted(chemostats))\n",
    "    print('Flowstats:', sorted(flowstats))\n",
    "    ## Stoichiometry\n",
    "    ## Create stoichiometry from bond graph.\n",
    "    s = st.stoich(bg,quiet=True)\n",
    "\n",
    "    ## Stoichiometry with chemostats\n",
    "    sc = st.statify(s,chemostats=chemostats,flowstats=flowstats)\n",
    "\n",
    "    ## Pathway stoichiometry\n",
    "    sp = st.path(s,sc)\n",
    "    \n",
    "    ## Print info\n",
    "    if verbose:\n",
    "        for stat in sorted(chemostats):\n",
    "            print(ch(stat)+',')\n",
    "\n",
    "    ## Energetics\n",
    "    if computePhi:\n",
    "        Phi,Phip = energetics(s,sp,phi)\n",
    "        #print('Phi units: kJ/mol')\n",
    "#         fac = -F/1000\n",
    "#         units='~\\si{\\kilo\\joule\\per\\mol}'\n",
    "        units = '~\\si{\\volt}'\n",
    "        print(st.sprintp(sc))\n",
    "        disp.Latex(st.sprintrl(sp,chemformula=True,Phi=Phip,showMu=showMu))\n",
    "        #return s,sc,sp,Phi*fac,Phip*fac,units\n",
    "        return s,sc,sp,Phip\n",
    "    else:\n",
    "        print(st.sprintrl(sp,chemformula=True))\n",
    "        Phip = 0\n",
    "        return s,sc,sp,Phip\n",
    "  \n",
    "def Pathway(S,chemostats,flowstats=[],verbose=False):\n",
    "    \"\"\" Analyse pathways\n",
    "    \"\"\"\n",
    "    \n",
    "    print('Chemostats:',sorted(chemostats))\n",
    "    print('Flowstats:', sorted(flowstats))\n",
    "    ## Stoichiometry\n",
    "    ## Create stoichiometry from bond graph.\n",
    "    #s = st.stoich(bg,quiet=True)\n",
    "    s = copy.copy(S)\n",
    "    ## Stoichiometry with chemostats\n",
    "    sc = st.statify(s,chemostats=chemostats,flowstats=flowstats)\n",
    "\n",
    "    ## Pathway stoichiometry\n",
    "    sp = st.path(s,sc)\n",
    "    \n",
    "    ## Print info\n",
    "    if verbose:\n",
    "        for stat in sorted(chemostats):\n",
    "            print(ch(stat)+',')\n",
    "\n",
    "\n",
    "    print(st.sprintrl(sp,chemformula=True))\n",
    "    \n",
    "    return s,sc,sp\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyse Pentose Phosphate Pathway with Glycolysis \n",
    "The pathways are isolated by using appropriate (zero-flow) flowstats. For compatibility with \\citet[\\S~18.2]{GarGri17} the pathways start from G6P (Glucose 6-phosphate)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Common chemostats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Chemostats(start='G6P',end=None):\n",
    "    chemostats = ['ADP','ATP','CO2','H','H2O','NADP','NADPH']\n",
    "    chemostats += [start]\n",
    "    if end is not None:\n",
    "        chemostats += end\n",
    "    return chemostats           "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \\ch{R5P} and \\ch{NADPH} generation\n",
    "- This pathway is isolated by setting PGI and TKT2 as flowstats and the\n",
    "product \\ch{R5P} is added to the chemostat list.\n",
    "\n",
    "- It is isolated from the TCA cycle by replacing the connecting reactions (PDH and PFL) by flowstats. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'stoich' from '/home/peterg/WEB/Github/Papers/GawPanCra21/stoich.py'>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R5P and NADPH generation\n",
      "Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NADP', 'NADPH', 'R5P']\n",
      "Flowstats: ['PGI', 'TKT2']\n",
      "\\begin{align}\n",
      "\\ch{G6P + H2O + 2 NADP &<>[ pr1 ] CO2 + 2 H + 2 NADPH + R5P }\n",
      "\\end{align}\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{G6P + H2O + 2 NADP &<>[ pr1 ] CO2 + 2 H + 2 NADPH + R5P }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imp.reload(st)\n",
    "print('R5P and NADPH generation')\n",
    "chemostats = Chemostats(start='G6P',end=['R5P'])\n",
    "flowstats = ['PGI','TKT2']\n",
    "s,sc,sp = Pathway(S,chemostats,flowstats=flowstats)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The pathway reaction \\ch{P_1} corresponds to the \\ch{R5P} and\n",
    "\\ch{NADPH} synthesis discussed in comment 1 of \n",
    "<cite data-cite=\"GarGri17\">Garrett and Grisham (2017)</cite>, p787.\n",
    "\n",
    "- It is isolated from the TCA cycle by replacing the connecting reactions (PDH and PFL) by flowstats. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \\ch{R5P} generation\n",
    "- This pathway is isolated by setting GAPD and G6PDH2R as flowstats and the\n",
    "product \\ch{R5P} is added to the chemostat list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R5P generation\n",
      "Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NADP', 'NADPH', 'R5P']\n",
      "Flowstats: ['G6PDH2R']\n",
      "\\begin{align}\n",
      "\\ch{ADP + H + 6 R5P &<>[ pr1 ] ATP + 5 G6P }\n",
      "\\end{align}\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ADP + H + 6 R5P &<>[ pr1 ] ATP + 5 G6P }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('R5P generation')\n",
    "chemostats = Chemostats(start='G6P',end=['R5P'])\n",
    "flowstats = ['G6PDH2R']\n",
    "s,sc,sp = Pathway(S,chemostats,flowstats=flowstats)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \\ch{NADPH} generation\n",
    "\n",
    "- This pathway is isolated by setting GAPD as a flowstat.\n",
    "\n",
    "- It is isolated from the TCA cycle by replacing the connecting reactions (PDH and PFL) by flowstats. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NADPH generation\n",
      "Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NADP', 'NADPH']\n",
      "Flowstats: []\n",
      "\\begin{align}\n",
      "\\ch{ADP + G6P + 6 H2O + 12 NADP &<>[ pr1 ] ATP + 6 CO2 + 11 H + 12 NADPH }\n",
      "\\end{align}\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ADP + G6P + 6 H2O + 12 NADP &<>[ pr1 ] ATP + 6 CO2 + 11 H + 12 NADPH }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('NADPH generation')\n",
    "chemostats = Chemostats(start='G6P')\n",
    "flowstats = []\n",
    "s,sc,sp = Pathway(S,chemostats,flowstats=flowstats)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The pathway reaction \\ch{pr1} corresponds to the \n",
    "\\ch{NADPH} synthesis discussed in comment 3 of\n",
    "<cite data-cite=\"GarGri17\">Garrett and Grisham (2017)</cite>, p787.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}