EcoliCoreModel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Note: this is the EcoliCoreModel.ipynb notebook. The\n",
    "PDF version \"The Escherichia coli Core Model: Modular Bond Graph Analysis\"\n",
    "is available [here](EcoliCoreModel.pdf).***"
   ]
  },
  {
   "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 of interest to systems biologists 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",
    "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>; energetic issues are not explicitly considered here although they are implicit in the bond graph models.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Strengths of the stoichiometric approach\n",
    "\n",
    "1. The stoichiometric matrix provides a precise description of the structure of chemical reaction networks.\n",
    "2. Stoichiometric matrices are available for biologically significant systems including whole-cell models.\n",
    "3. Powerful linear algebra concepts such as matrix null spaces and singular-value decompositions can be applied.\n",
    "4. It provides the basis for flux-balance analysis <cite data-cite=\"OrtThiPal10\"></cite>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Strengths of the bond graph approach\n",
    "1. Bond graph provide an energy-based approach to modelling.\n",
    "2. The bond graph model can be analysed for energy flows and efficiency as described by\n",
    "<cite data-cite=\"GawCra18\">Gawthrop and Crampin (2018)</cite>.\n",
    "3. The bond graph approach is multi domain and can thus, for example, model electrochemical systems including neurodynamics <cite data-cite=\"GawSieKam17\">(Gawthrop et al., 2017)</cite>, redox reactions <cite data-cite=\"Gaw17a\">(Gawthrop, 2017)</cite>  and transporters <cite data-cite=\"PanGawTra19\">(Pan et al., 2019)</cite>.\n",
    "4. The bond graph approch is modular: subsystems can be connected using energy ports <cite data-cite=\"GawCurCra15,GawCra16\"></cite>.\n",
    "5. The chemostat concept <cite data-cite=\"GawCra16\">Gawthrop and Crampin (2016)</cite> gives a convenient and flexible way of turning a closed system into an open system and analysing the concomitant pathway structure.\n",
    "6. [BondGraphTools](https://pypi.org/project/BondGraphTools/) provide a symbolic basis for describing and analysing bondgraphs within Python."
   ]
  },
  {
   "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": [
    "import BondGraphTools as bgt\n",
    "import numpy as np\n",
    "import IPython.display as disp\n",
    "\n",
    "## Stoichiometric analysis\n",
    "import stoich as st\n",
    "\n",
    "## Export stoichiometry as bond graph\n",
    "import stoichBondGraph as stbg\n",
    "\n",
    "## Modular bond graphs\n",
    "import modularBondGraph as mbg\n",
    "\n",
    "## Extract stoichiometry from a CobraPy model\n",
    "import CobraExtract as Extract\n",
    "\n",
    "## Control outputs\n",
    "quiet = True\n",
    "chemformula = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining the two approaches\n",
    "The key stoichiometric concept of pathways has already been given a bond graph interpretation <cite data-cite=\"Gaw17,GawCra17\"</cite>.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This note shows how:\n",
    "\n",
    "1. bond graphs can be created from stoichiometric data using the ecoli core model <cite data-cite=\"OrtFlePal10\">(Orth, Fleming and Palsson,2010)</cite> as an example.\n",
    "2. bond graphs of subsystems can be extracted from such stoichiometric data\n",
    "3. pathways can be anaysed and behavior elucidated using appropriate chice of chemostats\n",
    "4. bond graphs of the extracted subsystems can be recombined in a modular fashion\n",
    "\n",
    "Much remains to be done in exploiting the combination of the two approaches including:\n",
    "\n",
    "1. A new look at flux-balance analysis.\n",
    "2. Further work on energy and efficiency of cellular systems.\n",
    "3. Using the large stoichiometric models available either directly or extracted from models in SBML or cellML format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## External metabolites\n",
    "The standard stoichiometric approach is to create open systems from closed systems by adding \"dangling reactions\" to species which connect to the outside world as external metabolites-- for example: \\ch{ATP <>}. In contrast, the bond graph approacj would declare \\ch{ATP} to be a chemostat. Thus when extracting a bond graph from a stoichiometric model, dangling reactions are deleted and the corresponding species added to a list of chemostats.\n",
    "\n",
    "Chemostats provide a more flexible approach as they can be created without changing system structure.\n",
    "\n",
    "In the following discussion, chemostats are added as appropriate to illustrate the various pathways."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extract various modules from the ecoli core model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract full ecoli core model\n",
    "In this case the ecoli core model is extracted from the CobraPy model: 'textbook'. This model corresponds to that discussed in the the textbook <cite data-cite=\"Pal15\">B.O. Palsson, *Systems Biology* (2015)</cite>.\n",
    "An integer version of the stoichiometric matrix, together with reactions and species is produced. Note that the reaction labeled 'Biomass_Ecoli_core' is not a reaction but is associated with the optimisation procedure - it can be ignored in this notebook.\n",
    "\n",
    "Names of reactions and species are converted to upper case and the $c$ (cytosol) subscript removed for clarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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"
     ]
    }
   ],
   "source": [
    "sm = Extract.extract(Remove=['_C','__' ],quiet=quiet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract Glycolysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{GLCD_E + PEP &<>[ GLCPTS ] G6P + PYR }\\\\\n",
       "\\ch{G6P &<>[ PGI ] F6P }\\\\\n",
       "\\ch{ATP + F6P &<>[ PFK ] ADP + FDP + H }\\\\\n",
       "\\ch{FDP + H2O &<>[ FBP ] F6P + PI }\\\\\n",
       "\\ch{FDP &<>[ FBA ] DHAP + G3P }\\\\\n",
       "\\ch{DHAP &<>[ TPI ] G3P }\\\\\n",
       "\\ch{G3P + NAD + PI &<>[ GAPD ] 13DPG + H + NADH }\\\\\n",
       "\\ch{3PG + ATP &<>[ PGK ] 13DPG + ADP }\\\\\n",
       "\\ch{2PG &<>[ PGM ] 3PG }\\\\\n",
       "\\ch{2PG &<>[ ENO ] H2O + PEP }\n",
       "\\end{align}\n",
       "\n",
       "\\begin{align}\n",
       "\\ch{ADP + H + PEP &<>[ PYK ] ATP + PYR }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reaction = ['GLCPTS','PGI','PFK','FBP','FBA','TPI','GAPD','PGK','PGM','ENO','PYK']\n",
    "s0 = Extract.choose(sm,reaction=reaction)\n",
    "disp.Latex(st.sprintrl(s0,chemformula=True))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 0.0\n"
     ]
    }
   ],
   "source": [
    "s0['name'] = 'GLY_abg'\n",
    "stbg.model(s0)\n",
    "import GLY_abg\n",
    "s = st.stoich(GLY_abg.model(),quiet=True)\n",
    "## Sanity check\n",
    "err = np.linalg.norm(s['N']-s0['N'])\n",
    "print(\"Error:\",err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pathway analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 pathways\n",
      "0:  + PFK + FBP\n",
      "1:  + GLCPTS + PGI + PFK + FBA + TPI + 2 GAPD - 2 PGK - 2 PGM + 2 ENO + PYK\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ATP + H2O &<>[ pr1 ] ADP + H + PI }\\\\\n",
       "\\ch{2 ADP + GLCD_E + 2 NAD + 2 PI &<>[ pr2 ] 2 ATP + 2 H + 2 H2O + 2 NADH + 2 PYR }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chemostats0 = ['ADP','ATP','H2O','PI','H']\n",
    "sc0 = st.statify(s,chemostats=chemostats0)\n",
    "chemostats = ['GLCD_E','PYR','NAD','NADH']\n",
    "chemostats.extend(chemostats0)\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "print(st.sprintp(sc))\n",
    "sp = st.path(s,sc)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pathway reactions:\n",
    "\n",
    "-  A futile cycle consuming \\ch{ATP} \n",
    "-  A pathway from \\ch{GLCD_E} to \\ch{PYR}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract TCA part of full model \n",
    "\n",
    "A list of the reactions of the TCA cycle is provided and the corresponding species and stoichiometric matrix are extracted. Note that the NAD/NADP interconversion reaction NADTRHD is included."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ACCOA + H2O + OAA &<>[ CS ] CIT + COA + H }\\\\\n",
       "\\ch{CIT &<>[ ACONTA ] ACON + H2O }\\\\\n",
       "\\ch{ACON + H2O &<>[ ACONTB ] ICIT }\\\\\n",
       "\\ch{ICIT + NADP &<>[ ICDHYR ] AKG + CO2 + NADPH }\\\\\n",
       "\\ch{AKG + COA + NAD &<>[ AKGDH ] CO2 + NADH + SUCCOA }\\\\\n",
       "\\ch{ATP + COA + SUCC &<>[ SUCOAS ] ADP + PI + SUCCOA }\\\\\n",
       "\\ch{FUM + Q8H2 &<>[ FRD7 ] Q8 + SUCC }\\\\\n",
       "\\ch{Q8 + SUCC &<>[ SUCDI ] FUM + Q8H2 }\\\\\n",
       "\\ch{FUM + H2O &<>[ FUM ] MALL }\\\\\n",
       "\\ch{MALL + NAD &<>[ MDH ] H + NADH + OAA }\n",
       "\\end{align}\n",
       "\n",
       "\\begin{align}\n",
       "\\ch{NAD + NADPH &<>[ NADTRHD ] NADH + NADP }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reaction = ['CS','ACONTA','ACONTB','ICDHYR','AKGDH','SUCOAS',\n",
    "               'FRD7',\n",
    "                'SUCDI','FUM','MDH',\n",
    "        'NADTRHD']\n",
    "s0 = Extract.choose(sm,reaction=reaction)\n",
    "disp.Latex(st.sprintrl(s0,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create bond graph and recreate stoichiometry\n",
    "\n",
    "The bond graph TCA_abg is created and written to TCA_abg.py from whence it can me imported. The stoichiometric matrix generated from TCA_abg.model() is compared with that extracted from the ecoli core model to check that all is working correctely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 0.0\n"
     ]
    }
   ],
   "source": [
    "s0['name'] = 'TCA_abg'\n",
    "stbg.model(s0)\n",
    "import TCA_abg\n",
    "s = st.stoich(TCA_abg.model(),quiet=True)\n",
    "\n",
    "## Sanity check\n",
    "err = np.linalg.norm(s['N']-s0['N'])\n",
    "print(\"Error:\",err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pathway analysis\n",
    "Thee species corresponding to ATP hydrolysis, NAD/NADH, NADP/NADPH and Q8/Q8H2 (ubiquone, but maybe FAD, I think) are set as the basic chemostats in the list chemostats0. \n",
    "\n",
    "The  overall reaction of the TCA cycle converts ACCOA to COA and \\ch{CO2} \n",
    "<cite data-cite=\"GarGri17\"></cite>, so these three species are also set as chemostats.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No chemostats\n",
      "1 pathways\n",
      "0:  + FRD7 + SUCDI\n",
      "\n",
      "Basic chemostats ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'H', 'Q8', 'Q8H2']\n",
      "1 pathways\n",
      "0:  + FRD7 + SUCDI\n",
      "\n",
      "Chemostats ['ACCOA', 'COA', 'CO2', 'ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'H', 'Q8', 'Q8H2']\n",
      "2 pathways\n",
      "0:  + FRD7 + SUCDI\n",
      "1:  + CS + ACONTA + ACONTB + ICDHYR + AKGDH - SUCOAS - FRD7 + FUM + MDH + NADTRHD\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('No chemostats')\n",
    "print(st.sprintp(s))\n",
    "chemostats0 = ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI','H'\n",
    "              ,'Q8','Q8H2']\n",
    "sc0 = st.statify(s,chemostats=chemostats0)\n",
    "print('Basic chemostats',chemostats0)\n",
    "print(st.sprintp(sc0))\n",
    "chemostats = ['ACCOA','COA','CO2']\n",
    "chemostats.extend(chemostats0)\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "print('Chemostats',chemostats)\n",
    "print(st.sprintp(sc))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With only the basic chemostats, there is one internal cycle. This performs no conversions and so the reaction is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{&<>[ pr1 ] }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp0 = st.path(s,sc0)\n",
    "disp.Latex(st.sprintrl(sp0,chemformula=True,split=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the chemostats ['ACCOA','COA','CO2'] are added, the entire TCA cycle appears as a pathway and the overall reaction is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{&<>[ pr1 ] }\\\\\n",
       "\\ch{ACCOA + ADP + 2 H2O + 3 NAD + PI + Q8 &<>[ pr2 ] ATP + 2 CO2 + COA + 2 H + 3 NADH + Q8H2 }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp = st.path(s,sc)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True,split=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pathway reactions:\n",
    "\n",
    "- A null cycle (shared with closed system)\n",
    "- Pathway from ACCOA to CO2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Include Pyruvate (Pyr) reactions: PDH and PFL\n",
    "To take this example further, include the two reactions converting pyruvate to ACCOA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ACCOA + H2O + OAA &<>[ CS ] CIT + COA + H }\\\\\n",
       "\\ch{CIT &<>[ ACONTA ] ACON + H2O }\\\\\n",
       "\\ch{ACON + H2O &<>[ ACONTB ] ICIT }\\\\\n",
       "\\ch{ICIT + NADP &<>[ ICDHYR ] AKG + CO2 + NADPH }\\\\\n",
       "\\ch{AKG + COA + NAD &<>[ AKGDH ] CO2 + NADH + SUCCOA }\\\\\n",
       "\\ch{ATP + COA + SUCC &<>[ SUCOAS ] ADP + PI + SUCCOA }\\\\\n",
       "\\ch{FUM + Q8H2 &<>[ FRD7 ] Q8 + SUCC }\\\\\n",
       "\\ch{Q8 + SUCC &<>[ SUCDI ] FUM + Q8H2 }\\\\\n",
       "\\ch{FUM + H2O &<>[ FUM ] MALL }\\\\\n",
       "\\ch{MALL + NAD &<>[ MDH ] H + NADH + OAA }\n",
       "\\end{align}\n",
       "\n",
       "\\begin{align}\n",
       "\\ch{NAD + NADPH &<>[ NADTRHD ] NADH + NADP }\\\\\n",
       "\\ch{COA + NAD + PYR &<>[ PDH ] ACCOA + CO2 + NADH }\\\\\n",
       "\\ch{COA + PYR &<>[ PFL ] ACCOA + FOR }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reaction.extend(['PDH','PFL'])\n",
    "s0 = Extract.choose(sm,reaction=reaction)\n",
    "disp.Latex(st.sprintrl(s0,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create bond graph and recreate stoichiometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 0.0\n"
     ]
    }
   ],
   "source": [
    "s0['name'] = 'PyrTCA_abg'\n",
    "stbg.model(s0)\n",
    "import PyrTCA_abg\n",
    "s = st.stoich(PyrTCA_abg.model(),quiet=True)\n",
    "\n",
    "## Sanity check\n",
    "err = np.linalg.norm(s['N']-s0['N'])\n",
    "print(\"Error:\",err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pathway analysis\n",
    "The relevant chemostats are now the substrate \\ch{PYR} and the product \\ch{CO2} together with the reactions ATP hydrolysis, NAD/NADH, NADP/NADPH and Q8/Q8H2 andthe additional product FOR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No chemostats\n",
      "1 non-unit pathways\n",
      "0:  + FRD7 + SUCDI\n",
      "\n",
      "Basic chemostats ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'NADP', 'NADPH', 'H', 'Q8', 'Q8H2', 'FOR']\n",
      "1 non-unit pathways\n",
      "0:  + FRD7 + SUCDI\n",
      "\n",
      "Chemostats ['PYR', 'CO2', 'ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'NADP', 'NADPH', 'H', 'Q8', 'Q8H2', 'FOR']\n",
      "3 non-unit pathways\n",
      "0:  + FRD7 + SUCDI\n",
      "1:  + CS + ACONTA + ACONTB + ICDHYR + AKGDH - SUCOAS - FRD7 + FUM + MDH + PDH\n",
      "2:  + CS + ACONTA + ACONTB + ICDHYR + AKGDH - SUCOAS - FRD7 + FUM + MDH + PFL\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('No chemostats')\n",
    "print(st.sprintp(s,removeSingle=True))\n",
    "chemostats0 = ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI','NADP', 'NADPH','H'\n",
    "              ,'Q8','Q8H2','FOR']\n",
    "sc0 = st.statify(s,chemostats=chemostats0)\n",
    "print('Basic chemostats',chemostats0)\n",
    "print(st.sprintp(sc0,removeSingle=True))\n",
    "chemostats = ['PYR','CO2']\n",
    "chemostats.extend(chemostats0)\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "print('Chemostats',chemostats)\n",
    "print(st.sprintp(sc,removeSingle=True))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{&<>[ pr1 ] }\\\\\n",
       "\\ch{ADP + 2 H2O + 3 NAD + NADP + PI + PYR + Q8 &<>[ pr2 ] ATP + 3 CO2 + 2 H + 3 NADH + NADPH + Q8H2 }\\\\\n",
       "\\ch{ADP + 2 H2O + 2 NAD + NADP + PI + PYR + Q8 &<>[ pr3 ] ATP + 2 CO2 + FOR + 2 H + 2 NADH + NADPH + Q8H2 }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp = st.path(s,sc)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True,split=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pathway reactions:\n",
    "\n",
    "- A null cycle (shared with closed system)\n",
    "- Pathway from PYR to CO2 via reaction PDH\n",
    "- Pathway from PYR to CO2 via reaction PFL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract the Electron Transport Chain (ETC )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{4 H + NADH + Q8 &<>[ NADH16 ] 3 H_E + NAD + Q8H2 }\\\\\n",
       "\\ch{4 H + O2 + 2 Q8H2 &<>[ CYTBD ] 2 H2O + 4 H_E + 2 Q8 }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reaction = ['NADH16','CYTBD']\n",
    "s0 = Extract.choose(sm,reaction=reaction)\n",
    "disp.Latex(st.sprintrl(s0,chemformula=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 0.0\n"
     ]
    }
   ],
   "source": [
    "s0['name'] = 'ETC_abg'\n",
    "stbg.model(s0)\n",
    "import ETC_abg\n",
    "##imp.reload(ETC_abg)\n",
    "s = st.stoich(ETC_abg.model(),quiet=True)\n",
    "err = np.linalg.norm(s['N']-s0['N'])\n",
    "print(\"Error:\",err)\n",
    "s['reaction'] = s0['reaction']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pathway analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No chemostats\n",
      "0 pathways\n",
      "\n",
      "Chemostats ['NADH', 'NAD', 'O2', 'H2O', 'H', 'H_E']\n",
      "1 pathways\n",
      "0:  + 2 NADH16 + CYTBD\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('No chemostats')\n",
    "print(st.sprintp(s))\n",
    "chemostats = ['NADH','NAD','O2','H2O','H','H_E']\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "print('Chemostats',chemostats)\n",
    "print(st.sprintp(sc))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{12 H + 2 NADH + O2 &<>[ pr1 ] 2 H2O + 10 H_E + 2 NAD }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp = st.path(s,sc)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pathway reactions:\n",
    "\n",
    "- Pathway from NADH and O2 to NAD and H2O pumping 10 protons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract ATPase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ADP + 4 H_E + PI &<>[ ATPS4R ] ATP + 3 H + H2O }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reaction = ['ATPS4R']\n",
    "s0 = Extract.choose(sm,reaction=reaction)\n",
    "disp.Latex(st.sprintrl(s0,chemformula=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 0.0\n"
     ]
    }
   ],
   "source": [
    "s0['name'] = 'ATP_abg'\n",
    "stbg.model(s0)\n",
    "import ATP_abg\n",
    "##imp.reload(ATP_abg)\n",
    "s = st.stoich(ATP_abg.model(),quiet=True)\n",
    "err = np.linalg.norm(s['N']-s0['N'])\n",
    "print(\"Error:\",err)\n",
    "s['reaction'] = s0['reaction']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pathway analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No chemostats\n",
      "0 pathways\n",
      "\n",
      "Chemostats ['H', 'H_E', 'ATP', 'ADP', 'PI', 'H2O']\n",
      "1 pathways\n",
      "0:  + ATPS4R\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('No chemostats')\n",
    "print(st.sprintp(s))\n",
    "chemostats = ['H','H_E','ATP','ADP','PI','H2O']\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "print('Chemostats',chemostats)\n",
    "print(st.sprintp(sc,removeSingle=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ADP + 4 H_E + PI &<>[ pr1 ] ATP + 3 H + H2O }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp = st.path(s,sc,removeSingle=False)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pathway reactions:\n",
    "\n",
    "- Pathway from ADP to ATP pumped by 3 protons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modularity: the Glycolysis/TCA network\n",
    "To examine the glycolysis/TCA network, it is convenient and informative to take a modular approach:  the glycolysis network and the TCA network are extracted as separate subsystems and then combined. The example shows two approaches to combining the subsystems:\n",
    "\n",
    "1. combining the stoichiometric matrices\n",
    "2. combining the bond graphs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create composite model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "GlyTCA = bgt.new(name='GlyTCA')   # Create system\n",
    "Gly = GLY_abg.model()\n",
    "TCA = PyrTCA_abg.model()\n",
    "GlyTCA.add(Gly,TCA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unify common species\n",
    "PYR is produced by Gly and consumed by TCA.\n",
    "ATP, ADP, PI, H, NAD, NADH , H2O are common to both modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "common = ['PYR','ATP','ADP','PI','H','NAD','NADH','H2O']\n",
    "mbg.unify(GlyTCA,common=common,quiet=quiet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Stoichiometry of unified module\n",
    "s = st.stoich(GlyTCA,quiet=True)\n",
    "#disp.Latex(st.sprintrl(s,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pathway analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chemostats ['GLCD_E', 'CO2', 'ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'H', 'Q8', 'Q8H2', 'FOR']\n",
      "4 pathways\n",
      "0:  + PFK + FBP\n",
      "1:  + FRD7 + SUCDI\n",
      "2:  + GLCPTS + PGI + PFK + FBA + TPI + 2 GAPD - 2 PGK - 2 PGM + 2 ENO + PYK + 2 CS + 2 ACONTA + 2 ACONTB + 2 ICDHYR + 2 AKGDH - 2 SUCOAS - 2 FRD7 + 2 FUM + 2 MDH + 2 NADTRHD + 2 PDH\n",
      "3:  + GLCPTS + PGI + PFK + FBA + TPI + 2 GAPD - 2 PGK - 2 PGM + 2 ENO + PYK + 2 CS + 2 ACONTA + 2 ACONTB + 2 ICDHYR + 2 AKGDH - 2 SUCOAS - 2 FRD7 + 2 FUM + 2 MDH + 2 NADTRHD + 2 PFL\n",
      "\n"
     ]
    }
   ],
   "source": [
    "chemostats0 = ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI','H'\n",
    "              ,'Q8','Q8H2','FOR']\n",
    "chemostats = ['GLCD_E','CO2']\n",
    "chemostats.extend(chemostats0)\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "print('Chemostats',chemostats)\n",
    "print(st.sprintp(sc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ATP + H2O &<>[ pr1 ] ADP + PI + H }\\\\\n",
       "\\ch{&<>[ pr2 ] }\\\\\n",
       "\\ch{GLCD_E + 2 Q8 + 4 ADP + 4 PI + 10 NAD + 2 H2O &<>[ pr3 ] 6 CO2 + 2 Q8H2 + 4 ATP + 6 H + 10 NADH }\\\\\n",
       "\\ch{GLCD_E + 2 Q8 + 4 ADP + 4 PI + 8 NAD + 2 H2O &<>[ pr4 ] 4 CO2 + 2 FOR + 2 Q8H2 + 4 ATP + 6 H + 8 NADH }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp = st.path(s,sc)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pathway reactions:\n",
    "\n",
    "- A futile cycle\n",
    "- A null cycle\n",
    "- Pathway from GLCD_E to CO2 via reaction PDH\n",
    "- Pathway from GLCD_E to CO2 via reaction PFL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modularity: metabolism\n",
    "The above methods are brought together to generate the metabolic pathways by combining the modules describing:\n",
    "Glycolysis, the TCA cycle, the Electron Transport Chain and ATPase. Once again, the overall bond graph is created  by combining modules and unifying common species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "Met = bgt.new(name='Met')   # Create system\n",
    "Gly = GLY_abg.model()\n",
    "TCA = PyrTCA_abg.model()\n",
    "ETC = ETC_abg.model()\n",
    "ATP = ATP_abg.model()\n",
    "Met.add(Gly,TCA,ETC,ATP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Unify species common to modules\n",
    "common = ['PYR','ATP','ADP','PI','H','H_E','NAD','NADH','H2O','Q8','Q8H2']\n",
    "mbg.unify(Met,common=common,quiet=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pathway analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{GLCD_E + PEP &<>[ GLCPTS ] G6P + PYR }\\\\\n",
       "\\ch{G6P &<>[ PGI ] F6P }\\\\\n",
       "\\ch{F6P + ATP &<>[ PFK ] FDP + ADP + H }\\\\\n",
       "\\ch{FDP + H2O &<>[ FBP ] F6P + PI }\\\\\n",
       "\\ch{FDP &<>[ FBA ] DHAP + G3P }\\\\\n",
       "\\ch{DHAP &<>[ TPI ] G3P }\\\\\n",
       "\\ch{G3P + PI + NAD &<>[ GAPD ] 13DPG + H + NADH }\\\\\n",
       "\\ch{3PG + ATP &<>[ PGK ] 13DPG + ADP }\\\\\n",
       "\\ch{2PG &<>[ PGM ] 3PG }\\\\\n",
       "\\ch{2PG &<>[ ENO ] PEP + H2O }\n",
       "\\end{align}\n",
       "\n",
       "\\begin{align}\n",
       "\\ch{PEP + ADP + H &<>[ PYK ] PYR + ATP }\\\\\n",
       "\\ch{ACCOA + OAA + H2O &<>[ CS ] CIT + COA + H }\\\\\n",
       "\\ch{CIT &<>[ ACONTA ] ACON + H2O }\\\\\n",
       "\\ch{ACON + H2O &<>[ ACONTB ] ICIT }\\\\\n",
       "\\ch{ICIT + NADP &<>[ ICDHYR ] AKG + CO2 + NADPH }\\\\\n",
       "\\ch{AKG + COA + NAD &<>[ AKGDH ] CO2 + SUCCOA + NADH }\\\\\n",
       "\\ch{COA + SUCC + ATP &<>[ SUCOAS ] SUCCOA + ADP + PI }\\\\\n",
       "\\ch{FUM + Q8H2 &<>[ FRD7 ] SUCC + Q8 }\\\\\n",
       "\\ch{SUCC + Q8 &<>[ SUCDI ] FUM + Q8H2 }\\\\\n",
       "\\ch{FUM + H2O &<>[ FUM ] MALL }\n",
       "\\end{align}\n",
       "\n",
       "\\begin{align}\n",
       "\\ch{MALL + NAD &<>[ MDH ] OAA + H + NADH }\\\\\n",
       "\\ch{NADPH + NAD &<>[ NADTRHD ] NADP + NADH }\\\\\n",
       "\\ch{COA + PYR + NAD &<>[ PDH ] ACCOA + CO2 + NADH }\\\\\n",
       "\\ch{COA + PYR &<>[ PFL ] ACCOA + FOR }\\\\\n",
       "\\ch{4 H + NADH + Q8 &<>[ NADH16 ] 3 H_E + NAD + Q8H2 }\\\\\n",
       "\\ch{O2 + 4 H + 2 Q8H2 &<>[ CYTBD ] 4 H_E + 2 H2O + 2 Q8 }\\\\\n",
       "\\ch{ADP + PI + 4 H_E &<>[ ATPS4R ] ATP + 3 H + H2O }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = st.stoich(Met,quiet=True)\n",
    "disp.Latex(st.sprintrl(s,chemformula=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chemostats ['GLCD_E', 'CO2', 'O2', 'ADP', 'ATP', 'H2O', 'PI', 'H']\n",
      "3 pathways\n",
      "0:  + PFK + FBP\n",
      "1:  + FRD7 + SUCDI\n",
      "2:  + 2 GLCPTS + 2 PGI + 2 PFK + 2 FBA + 2 TPI + 4 GAPD - 4 PGK - 4 PGM + 4 ENO + 2 PYK + 4 CS + 4 ACONTA + 4 ACONTB + 4 ICDHYR + 4 AKGDH - 4 SUCOAS - 4 FRD7 + 4 FUM + 4 MDH + 4 NADTRHD + 4 PDH + 20 NADH16 + 12 CYTBD + 27 ATPS4R\n",
      "\n"
     ]
    }
   ],
   "source": [
    "chemostats0 = ['ADP','ATP','H2O','PI','H']\n",
    "chemostats = ['GLCD_E','CO2','O2']\n",
    "chemostats.extend(chemostats0)\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "print('Chemostats',chemostats)\n",
    "print(st.sprintp(sc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ATP + H2O &<>[ pr1 ] ADP + PI + H }\\\\\n",
       "\\ch{&<>[ pr2 ] }\\\\\n",
       "\\ch{2 GLCD_E + 12 O2 + 35 ADP + 35 PI + 35 H &<>[ pr3 ] 12 CO2 + 35 ATP + 47 H2O }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp = st.path(s,sc)\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three pathways: \n",
    "\n",
    "- A futile cycle\n",
    "- A null cycle\n",
    "- Each GLCD_E, combined with \\ch{6 O2}, reverses ATP hydrolysis  (\\ch{ATP + H2O <> ADP + PI + H+}) to give 17.5 ATP molecules with an additional \\ch{6 H2O} and \\ch{6 CO2}; 17.5 ATP is the value quoted\n",
    "by <cite data-cite={Pal15}>Palsson (2015)</cite>."
   ]
  },
  {
   "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.6.9"
  },
  "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
}