Parameter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "This note illustrates an approach to fitting the parameters of a bond\n",
    "graph model to experimental data. Insofar as the parameters are\n",
    "associated with a bond graph, they are *physically-plausible*\n",
    "<cite data-cite=\"GawCudCra20\"></cite>.\n",
    "\n",
    "    \n",
    "The approach uses a bond-graph derived from a stoichiometric model of *e.coli*\n",
    "<cite data-cite=\"OrtFlePal10\"></cite>\n",
    "(using a method described elsewhere <cite data-cite=\"Gaw20X\"></cite>)\n",
    "combined with experimental values of *reaction potential*, *reaction flux* and *species concentration* from the literature <cite data-cite=\"ParRubXu16\"></cite>.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Paths\n",
    "NeedPath=False\n",
    "if NeedPath:\n",
    "    import sys\n",
    "    sys.path += ['/usr/lib/python3/dist-packages']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Maths library\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "## BG tools\n",
    "import BondGraphTools as bgt\n",
    "\n",
    "## SVG bond graph\n",
    "import svgBondGraph as sbg\n",
    "\n",
    "## BG stoichiometric utilities\n",
    "import stoich as st\n",
    "\n",
    "## Modular bond graphs\n",
    "import modularBondGraph as mbg\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",
    "## PLotting\n",
    "import matplotlib.pyplot as plt\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",
    "quiet = True\n",
    "showMu=True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadratic programming QP.\n",
    "\\begin{align}\n",
    "\\text{minimise } & \\frac{1}{2}x^TPx + q^Tx\\\\\n",
    "\\text{subject to } & Gx \\le h\\\\\n",
    "\\text{and } & Ax = b\n",
    "\\end{align}\n",
    "\n",
    "In the case considered here, there is no equality constraint and\n",
    "\\begin{align}\n",
    "x &= \\hat{\\phi}\\\\\n",
    "P & = N N^T + \\mu I_{n_X \\times n_X}\\\\\n",
    "q &= (N \\Phi)^T \\\\\n",
    "G &= N^T\\\\\n",
    "h &= -\\Phi_{min}\n",
    "\\end{align}\n",
    "\n",
    "$\\mu > 0$ is required to give a convex QP: in essence it turns a non-unique solution for $\\phi$ into a minimum norm solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Quadratic programming stuff.\n",
    "import quadprog\n",
    "\n",
    "## Function from https://scaron.info/blog/quadratic-programming-in-python.html\n",
    "def quadprog_solve_qp(P, q, G=None, h=None, A=None, b=None):\n",
    "    qp_G = .5 * (P + P.T)   # make sure P is symmetric\n",
    "    qp_a = -q\n",
    "    if A is not None:\n",
    "        qp_C = -numpy.vstack([A, G]).T\n",
    "        qp_b = -numpy.hstack([b, h])\n",
    "        meq = A.shape[0]\n",
    "    else:  # no equality constraint\n",
    "        qp_C = -G.T\n",
    "        qp_b = -h\n",
    "        meq = 0\n",
    "    return quadprog.solve_qp(qp_G, qp_a, qp_C, qp_b, meq)[0]\n",
    "\n",
    "## Function to compute phi from Phi subject to Phi>positive number\n",
    "## NN Reduced N corresponding to known Phi\n",
    "def quadsolve_phi(N0,N1,Phi0,Phi_min=0.0,mu=1e-10):\n",
    "    \n",
    "    (n_X,n_V) = N1.shape\n",
    "    print(N1.shape)\n",
    "    P = 1.0*N0@(N0.T) + mu*np.eye(n_X)\n",
    "    q = (N0@Phi0).T\n",
    "    G = 1.0*N1.T\n",
    "    h = -Phi_min*np.ones((n_V))\n",
    "    phi = quadprog_solve_qp(P, q, G=G, h=h)\n",
    "    #Phi = -N.T@phi\n",
    "    \n",
    "    return phi\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conversion factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "To convert from kJ/mol to mV, divide by 10.4\n"
     ]
    }
   ],
   "source": [
    "Factor = st.F()/1e6\n",
    "print(f'To convert from kJ/mol to mV, divide by {1/Factor:4.3}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extract Model\n",
    "This example uses the Glycolysis and Pentose Phosphate pathways.\n",
    "\n",
    "Notes:\n",
    "\n",
    "- Reactions RPI, PGK and PGM are reversed to correspond to positive flows.\n",
    "- The resultant stoichiometric matrix $N$ relates reaction flows ($f$) to species flows ($\\dot{x}$): \n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{x} = N f\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract stoichiometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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",
      "Multiplying reaction PGK ( 54 ) by -1\n",
      "Multiplying reaction PGM ( 56 ) by -1\n",
      "Multiplying reaction RPI ( 65 ) by -1\n"
     ]
    }
   ],
   "source": [
    "sm = Extract.extract(cobraname='textbook',Remove=['_C','__' ], \n",
    "                     negReaction=['RPI','PGK','PGM'], quiet=quiet)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'GlyPPP_abg' from '/home/peterg/WEB/Github/Papers/tmp/GlyPPP_abg.py'>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "ss = Extract.choose(sm,reaction=reaction)\n",
    "\n",
    "## Create BG\n",
    "ss['name'] = name\n",
    "stbg.model(ss)\n",
    "import GlyPPP_abg\n",
    "imp.reload(GlyPPP_abg)\n",
    "s = st.stoich(GlyPPP_abg.model(),quiet=quiet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['ADP', 'ATP', 'CO2', 'G3P', 'G6P', 'H', 'H2O', 'NADP', 'NADPH', 'R5P']\n",
      "3 pathways\n",
      "0:  + PGI + PFK + FBA + TPI\n",
      "1:  + G6PDH2R + PGL + GND + RPI\n",
      "2:  - 2 PGI + 2 G6PDH2R + 2 PGL + 2 GND + TKT2 + TALA + TKT1 + 2 RPE\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\ch{ATP + G6P &<>[ PPP1 ] ADP + 2 G3P + H }\\\\\n",
       "\\ch{G6P + H2O + 2 NADP &<>[ PPP2 ] CO2 + 2 H + 2 NADPH + R5P }\\\\\n",
       "\\ch{2 H2O + 4 NADP + R5P &<>[ PPP3 ] 2 CO2 + G3P + 4 H + 4 NADPH }\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "## Set up chemostats\n",
    "chemostats = ['ADP','ATP','H','H2O','NADP','NADPH','CO2']\n",
    "chemostats += ['G6P','G3P','R5P']\n",
    "chemostats.sort()\n",
    "print(chemostats)\n",
    "sc = st.statify(s,chemostats=chemostats)\n",
    "\n",
    "sp = st.path(s,sc,pathname='PPP')\n",
    "print(st.sprintp(sc))\n",
    "disp.Latex(st.sprintrl(sp,chemformula=True))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\begin{align}\n",
      "K &=\n",
      "\\left(\\begin{matrix}1 & 0 & -2\\\\1 & 0 & 0\\\\1 & 0 & 0\\\\1 & 0 & 0\\\\0 & 1 & 2\\\\0 & 1 & 2\\\\0 & 1 & 2\\\\0 & 1 & 0\\\\0 & 0 & 1\\\\0 & 0 & 1\\\\0 & 0 & 1\\\\0 & 0 & 2\\end{matrix}\\right)\n",
      "\\end{align}\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "K &=\n",
       "\\left(\\begin{matrix}1 & 0 & -2\\\\1 & 0 & 0\\\\1 & 0 & 0\\\\1 & 0 & 0\\\\0 & 1 & 2\\\\0 & 1 & 2\\\\0 & 1 & 2\\\\0 & 1 & 0\\\\0 & 0 & 1\\\\0 & 0 & 1\\\\0 & 0 & 1\\\\0 & 0 & 2\\end{matrix}\\right)\n",
       "\\end{align}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(st.sprintl(sc,'K'))\n",
    "disp.Latex(st.sprintl(sc,'K'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract reaction potentials $\\Phi$ and deduce plausible species potentials $\\phi$.\n",
    "\n",
    "Because of the energetic constaints implied by the bond graph, the reaction potentials $\\Phi$ are \n",
    "related to the species potentials $\\phi$ by\n",
    "\\begin{equation}\n",
    "\\Phi = - N^T \\phi\n",
    "\\end{equation}\n",
    "Typically, there are more species than reactions and so $N$ has more rows than columns.\n",
    "Given the reaction potentials $\\Phi$, the species potentials can be estimated using\n",
    "the *pseudo inverse* $N^\\dagger$ of $-N^T$:\n",
    "\\begin{equation}\n",
    "\\hat{\\phi} = N^\\dagger \\Phi\n",
    "\\end{equation}\n",
    "\n",
    "Notes:\n",
    "\n",
    "- In general $\\hat{\\phi} \\ne \\phi$ but is physically plausible insofar as $-N^T \\hat{\\phi}=\\Phi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getPhi(s,Phi_hyd=0.5,phi_6PGL=None,quadprog=False):\n",
    "    \"\"\"Extract phi for given system using\n",
    "    Reaction potentials from ParRubXu16\"\"\"\n",
    "    \n",
    "    ## Reaction potentials from ParRubXu16\n",
    "    PHI = phiData.Phi_ParRubXu16_Measured()\n",
    "    \n",
    "#     Phenotype = 'Mammalian'\n",
    "#     Phenotype = 'Yeast'\n",
    "    Phenotype = 'Ecoli'\n",
    "    Phi_reac = PHI[Phenotype]\n",
    "    \n",
    "    Phi = np.zeros((len(s['reaction']),1))\n",
    "    N = copy.copy(s['N'])\n",
    "    N_0 = None\n",
    "    N_1 = None\n",
    "    Phi_0 = []\n",
    "    for j,reac in enumerate(s['reaction']):\n",
    "        if (reac in Phi_reac.keys()) and not np.isnan(Phi_reac[reac]):\n",
    "            Phi_0.append(Phi_reac[reac])\n",
    "            if N_0 is None:\n",
    "                N_0 = N[:,j]\n",
    "            else:\n",
    "                N_0 = np.vstack((N_0,N[:,j]))\n",
    "        else:\n",
    "            if N_1 is None:\n",
    "                N_1 = N[:,j]\n",
    "            else:\n",
    "                N_1 = np.vstack((N_1,N[:,j]))\n",
    "\n",
    "    Phi_0 = np.array(Phi_0)\n",
    "    #print(N_1)\n",
    "    \n",
    "    ## Compute Phi\n",
    "    N_0 = N_0.T\n",
    "    N_1 = N_1.T\n",
    "    \n",
    "    n_X,n_V = N_0.shape\n",
    "    print(f'Extracting {n_X} values of phi from {n_V} values of Phi')\n",
    "    \n",
    "    if quadprog:\n",
    "        phi = quadsolve_phi(N_0,N_1,Phi_0,Phi_min=1e-3,mu=1e-10)\n",
    "    else:\n",
    "        ## Compute Phi using pseudo inverse\n",
    "        pinvNT =  scipy.linalg.pinv(N_0.T)\n",
    "        phi = -pinvNT@Phi_0\n",
    "        \n",
    "    if phi_6PGL is not None:\n",
    "        ## Reset 6PGL\n",
    "        i_6PGL = s['species'].index('6PGL')\n",
    "        phi[i_6PGL] = phi_6PGL\n",
    "        print (f'Resetting phi_6GPL to {int(1e3*phi[i_6PGL])} mV' )\n",
    "\n",
    "    ## Sanity check\n",
    "    Phi_new = -N_0.T@phi\n",
    "    err = np.linalg.norm(Phi_new-Phi_0)\n",
    "    print(f'Phi error = {int(err*1000)}mV\\n')\n",
    "    \n",
    "    Phi = -N.T@phi\n",
    "    \n",
    "    return Phi,phi,Phi_0,Phi_reac\n",
    "     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracting 19 values of phi from 10 values of Phi\n",
      "Phi error = 0mV\n",
      "\n",
      "Minimum Phi =  -3 mV\n"
     ]
    }
   ],
   "source": [
    "Phi_,phi_est_,Phi_0_,Phi_reac_ = getPhi(s,quadprog=False)\n",
    "print('Minimum Phi = ', int(round(np.min(1e3*Phi_))), 'mV')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracting 19 values of phi from 10 values of Phi\n",
      "(19, 2)\n",
      "Phi error = 0mV\n",
      "\n",
      "Minimum Phi =  0 mV\n",
      "\n",
      "Change in phi\n",
      "1 6PGL\t 1\n",
      "12 H2O\t 1\n",
      "\n",
      "Change in Phi\n",
      "5 PGL\t 4 1 -3\n"
     ]
    }
   ],
   "source": [
    "Phi,phi_est,Phi_0,Phi_reac = getPhi(s,quadprog=True)\n",
    "print('Minimum Phi = ', int(round(np.min(1e3*Phi))), 'mV')\n",
    "\n",
    "print('\\nChange in phi')\n",
    "for i,spec in enumerate(s['species']):\n",
    "    change = int(1e3*(phi_est[i]-phi_est_[i]))\n",
    "    if not change==0:\n",
    "        print(f'{i} {spec}\\t {change}')\n",
    "    \n",
    "print('\\nChange in Phi')\n",
    "for i,reac in enumerate(s['reaction']):\n",
    "    change = int(round(1e3*(Phi[i]-Phi_[i])))\n",
    "    if not change == 0:\n",
    "        print(f'{i} {reac}\\t {change} {int(round(1e3*Phi[i]))} {int(round(1e3*Phi_[i]))}')\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracted reactions and reaction potentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{xalignat*}{2}\n",
       "\\ch{G6P &<>[ PGI ] F6P }&&(16~\\text{mV})\\;[-1.59~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{ATP + F6P &<>[ PFK ] ADP + FDP + H }&&(68~\\text{mV})\\;[-6.64~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{FDP &<>[ FBA ] DHAP + G3P }&&(20~\\text{mV})\\;[-1.93~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{DHAP &<>[ TPI ] G3P }&&(7~\\text{mV})\\;[-0.77~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{G6P + NADP &<>[ G6PDH2R ] 6PGL + H + NADPH }&&(82~\\text{mV})\\;[-7.99~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{6PGL + H2O &<>[ PGL ] 6PGC + H }&&(1~\\text{mV})\\;[-0.10~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{6PGC + NADP &<>[ GND ] CO2 + NADPH + RU5PD }&&(114~\\text{mV})\\;[-11.05~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{RU5PD &<>[ RPI ] R5P }&&(0~\\text{mV})\\;[-0.00~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{E4P + XU5PD &<>[ TKT2 ] F6P + G3P }&&(16~\\text{mV})\\;[-1.58~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{G3P + S7P &<>[ TALA ] E4P + F6P }&&(54~\\text{mV})\\;[-5.25~\\text{kJ mol}^{-1}]\n",
       "\\end{xalignat*}\n",
       "\n",
       "\\begin{xalignat*}{2}\n",
       "\\ch{R5P + XU5PD &<>[ TKT1 ] G3P + S7P }&&(4~\\text{mV})\\;[-0.39~\\text{kJ mol}^{-1}]\\\\\n",
       "\\ch{RU5PD &<>[ RPE ] XU5PD }&&(0~\\text{mV})\\;[-0.08~\\text{kJ mol}^{-1}]\n",
       "\\end{xalignat*}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "disp.Latex(st.sprintrl(s,chemformula=True,Phi=Phi,units=['mV','kJ'] ,showMu=showMu))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deduce Pathway Flows\n",
    "\n",
    "From basic stoichiometric analysis, steady-state flows can be written as:\n",
    "\\begin{align}\n",
    "f &= K_p f_p\\\\\n",
    "\\text{where }\n",
    "K_p N^{cd} &= 0\n",
    "\\end{align}\n",
    "Note that the *pathway matrix* $K_p$ is dependent on the choice of chemostats.\n",
    "\n",
    "Given a set of experimental flows $f$, an estimate $\\hat{f}_p$ of $f_p$ can be obtained from the *least-squares* formula:\n",
    "\\begin{equation}\n",
    "(K_p^T K_p) \\hat{f}_p = K_p^T f\n",
    "\\end{equation}\n",
    "\n",
    "Notes:\n",
    "\n",
    "- $v_p$ is a $n_p$ vector containg the pathways flows\n",
    "- $(K_p^T K_p)$ is a square $n_p \\times n_p$ matrix where $n_p$ is the number of pathways\n",
    "- If some flows are not measured, the corresponding rows of $K_p$ are deleted\n",
    "- the reaction flows (including the missing ones) can be estimated from $\\hat{f} = K_p \\hat{f}_p$.\n",
    "- the estimated chemostat flows are given by the non-zero elements of\n",
    "\\begin{equation}\n",
    "\\hat{\\dot{x}} = N \\hat{f}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PathwayFlux(K,reaction,Reaction,flux):\n",
    "    \n",
    "    #KK = st.singleRemove(K)\n",
    "    KK = K\n",
    "    Kp = None\n",
    "    Flux = {}\n",
    "    reac_known = []\n",
    "    #flux = phiData.ParRubXu16_flux()\n",
    "    for i,reac in enumerate(reaction):\n",
    "        if reac in flux.keys():\n",
    "            reac_known.append(reac)\n",
    "            fi = flux[reac]\n",
    "            #Ki = np.abs(KK[i,:])\n",
    "            Ki = KK[i,:]\n",
    "            #print(reac,Ki)\n",
    "            if Kp is None:\n",
    "                Kp = Ki\n",
    "                f = fi\n",
    "            else:\n",
    "                Kp = np.vstack((Kp,Ki))\n",
    "                f = np.vstack((f,fi))\n",
    "    #print(Kp)\n",
    "\n",
    "    if Kp is not None:    \n",
    "        #print(f)\n",
    "        f0 = np.linalg.solve(Kp.T@Kp,Kp.T@f)\n",
    "        for i,Reac in enumerate(Reaction):\n",
    "            Flux[Reac] = f0[i][0]\n",
    "        #print(f0)\n",
    "        f_est = Kp@f0\n",
    "        #print(Kp@f0-f)\n",
    "        \n",
    "    error = np.linalg.norm(f_est-f)/len(f)\n",
    "    print(f'Flux error = {error:.2e}')\n",
    "    \n",
    "    return Flux,f0,f_est,f,reac_known\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reaction constants (modified mass-action)\n",
    "\n",
    "The modified mass-action formula is <cite data-cite=\"GawCudCra20\"></cite>:\n",
    "\\begin{align}\n",
    "  f &= \\kappa \\left ( \\exp \\frac{\\Phi^f}{\\alpha V_N} - \\exp\n",
    "      \\frac{\\Phi^r}{\\alpha V_N} \\right )\\label{eq:MMA0}\n",
    "\\end{align}\n",
    "Thus an estimate for $\\kappa$ can be computed as:\n",
    "\\begin{align}\n",
    "\\hat{\\kappa} &= \\frac{\\hat{f}}{f_0}\\\\\n",
    "\\text{where }\n",
    "f_0 &= \\left ( \\exp \\frac{\\Phi^f}{\\alpha V_N} - \\exp\n",
    "      \\frac{\\Phi^r}{\\alpha V_N} \\right )\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reactionConstant(s,phi_est,f_est,alpha=1):\n",
    "    \n",
    "    V_N = st.V_N()\n",
    "    \n",
    "    ## Extract stoichiometry\n",
    "    N = s['N']\n",
    "    Nf = s['Nf']\n",
    "    Nr = s['Nr']\n",
    "    reaction = s['reaction']\n",
    "    \n",
    "    ## Compute Phis from estimated phi \n",
    "    Phi_ = -N.T@phi_est\n",
    "    Phi_f = Nf.T@phi_est \n",
    "    Phi_r = Nr.T@phi_est\n",
    "\n",
    "    ## Compute normalised flow rates\n",
    "    f0 = (np.exp(Phi_f/(alpha*V_N)) - np.exp(Phi_r/(alpha*V_N)))\n",
    "    #print(f0)\n",
    "    parameter = {}\n",
    "    for i,react in enumerate(reaction):\n",
    "        kap = f_est[i][0]/f0[i]\n",
    "        parameter[f'kappa_{react}'] = kap\n",
    "        #print(f'{react}: \\tPhi = {int(Phi_[i]*1000)}mV, \\tf_est = {f_est[i][0]:.2e}, \\tkappa = {kap:.2}')\n",
    "        \n",
    "    return parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show computed reaction flows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flux error = 3.10e-01\n",
      "\n",
      "\n",
      "%% LaTeX table\n",
      "\\hline\n",
      "Reaction &\t $\\Phi$ &\t $\\hat{\\Phi}$ &\t $f$ & $\\hat{f}$ & $\\hat{\\kappa}$\\\\\n",
      "\\hline\n",
      "PGI &\t 16 &\t 16 &\t 100 & 99 & 2.6e+02 \\\\\n",
      "PFK &\t 69 &\t 69 &\t 104 & 105 & 9.1e+01 \\\\\n",
      "FBA &\t 20 &\t 20 &\t 106 & 105 & 2.7e+02 \\\\\n",
      "TPI &\t 8 &\t 8 &\t 105 & 105 & 5.9e+02 \\\\\n",
      "G6PDH2R &\t -- &\t 83 &\t -- & 19 & 7.8 \\\\\n",
      "PGL &\t -- &\t 1 &\t -- & 19 & 4.9e+02 \\\\\n",
      "GND &\t 115 &\t 115 &\t 19 & 19 & 2.1 \\\\\n",
      "RPI &\t 0 &\t 0 &\t 13 & 13 & 7e+03 \\\\\n",
      "TKT2 &\t 16 &\t 16 &\t 2 & 3 & 1.5e+01 \\\\\n",
      "TALA &\t 54 &\t 54 &\t -- & 3 & 2.8 \\\\\n",
      "TKT1 &\t 4 &\t 4 &\t 5 & 3 & 1.5e+01 \\\\\n",
      "RPE &\t 1 &\t 1 &\t 6 & 6 & 1.6e+02 \\\\\n",
      "\\hline\n"
     ]
    }
   ],
   "source": [
    "Reaction = ['PPP1','PPP2','PPP3']\n",
    "#print(s['reaction'])\n",
    "## Reaction flows\n",
    "Flux = phiData.ParRubXu16_flux()\n",
    "flux = Flux['Ecoli']\n",
    "\n",
    "## Normalise flux wrt PGI = 100\n",
    "flux_PGI = flux['PGI']\n",
    "for reac in flux.keys():\n",
    "    flux[reac] *= 100/flux_PGI\n",
    "\n",
    "fluxp,f0,f_est,f,reaction = PathwayFlux(sc['K'],s['reaction'],Reaction,flux)\n",
    "\n",
    "\n",
    "## Reaction constants\n",
    "f_est = sc['K']@f0\n",
    "parameter = reactionConstant(s,phi_est,f_est)\n",
    "\n",
    "#f_est = sc['K']@f0\n",
    "j=0\n",
    "\n",
    "print('\\n\\n%% LaTeX table')\n",
    "print('\\\\hline')\n",
    "print('Reaction &\\t $\\Phi$ &\\t $\\hat{\\Phi}$ &\\t $f$ & $\\\\hat{f}$ & $\\\\hat{\\\\kappa}$\\\\\\\\')\n",
    "print('\\\\hline')\n",
    "for i,reac in enumerate(s['reaction']):\n",
    "    if reac in flux.keys():\n",
    "        ff = f'{int(round(f[j][0]))}'\n",
    "        j += 1\n",
    "    else:\n",
    "        ff = '--'\n",
    "    if reac in Phi_reac.keys():\n",
    "        PP = f'{int(round(1e3*Phi_reac[reac]))}'\n",
    "    else:\n",
    "        PP = '--'\n",
    "    kappa = 'kappa_'+reac\n",
    "    print(f'{reac} &\\t {PP} &\\t {int(round(1e3*Phi[i]))} &\\t {ff} & {int(round(f_est[i][0]))} & {parameter[kappa]:.2} \\\\\\\\')\n",
    "print('\\\\hline')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show computed chemostat flows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "%% LaTeX table\n",
      "\\hline\n",
      "Chemostat &\t flow \\\\\n",
      "\\hline\n",
      "ADP &\t 105 \\\\\n",
      "ATP &\t -105 \\\\\n",
      "CO2 &\t 19 \\\\\n",
      "G3P &\t 213 \\\\\n",
      "G6P &\t -119 \\\\\n",
      "H &\t 144 \\\\\n",
      "H2O &\t -19 \\\\\n",
      "NADP &\t -39 \\\\\n",
      "NADPH &\t 39 \\\\\n",
      "R5P &\t 10 \\\\\n",
      "\\hline\n"
     ]
    }
   ],
   "source": [
    "dx_est = s['N']@f_est\n",
    "\n",
    "print('\\n\\n%% LaTeX table')\n",
    "print('\\\\hline')\n",
    "print('Chemostat &\\t flow \\\\\\\\')\n",
    "print('\\\\hline')\n",
    "for i,spec in enumerate(s['species']):\n",
    "    if spec in chemostats:\n",
    "        print(f'{spec} &\\t {int(round(dx_est[i][0]))} \\\\\\\\')\n",
    "print('\\\\hline')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show pathway flows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "%% LaTeX table\n",
      "\\hline\n",
      "Pathway &\t $\\hat{f}_p$ \\\\\n",
      "\\hline\n",
      "PPP1 &\t 105 \\\\\n",
      "PPP2 &\t 13 \\\\\n",
      "PPP3 &\t 3 \\\\\n",
      "\\hline\n"
     ]
    }
   ],
   "source": [
    "print('\\n\\n%% LaTeX table')\n",
    "print('\\\\hline')\n",
    "print('Pathway &\\t $\\hat{f}_p$ \\\\\\\\')\n",
    "print('\\\\hline')\n",
    "for reac in fluxp.keys():\n",
    "    print(f'{reac} &\\t {int(round(fluxp[reac]))} \\\\\\\\')\n",
    "print('\\\\hline')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Species constants\n",
    "\n",
    "\\begin{equation}\n",
    "K = \\frac{\\exp \\phi}{x^\\circ} = \\frac{\\exp \\phi}{V c^\\circ}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "%% LaTeX table\n",
      "\\hline\n",
      "Species &\t $\\hat{\\phi}~mV$ & $c$ & $\\hat{K}$ \\\\\n",
      "\\hline\n",
      "6PGC & 29 &  \t0.017 & 62.4 \\\\\n",
      "ADP & -27 &  \t0.57 & 1.711 \\\\\n",
      "ATP & 27 &  \t4.7 & 0.22 \\\\\n",
      "CO2 & -30 &  \t7.6 & 0.1272 \\\\\n",
      "DHAP & -10 &  \t1.6 & 0.6075 \\\\\n",
      "E4P & -27 &  \t0.01 & 94.47 \\\\\n",
      "F6P & -21 &  \t0.097 & 10.1 \\\\\n",
      "FDP & -8 &  \t1.5 & 0.6528 \\\\\n",
      "G3P & -18 &  \t0.14 & 6.967 \\\\\n",
      "G6P & -5 &  \t0.68 & 1.474 \\\\\n",
      "NADP & 30 &  \t0.028 & 36.29 \\\\\n",
      "NADPH & -30 &  \t0.065 & 14.83 \\\\\n",
      "R5P & 5 &  \t0.028 & 35.4 \\\\\n",
      "RU5PD & 5 &  \t0.0053 & 190.8 \\\\\n",
      "S7P & 24 &  \t0.018 & 56.58 \\\\\n",
      "XU5PD & 5 &  \t0.03 & 33.6 \\\\\n",
      "\\hline\n"
     ]
    }
   ],
   "source": [
    "#imp.reload(phiData)\n",
    "\n",
    "print('\\n\\n%% LaTeX table')\n",
    "print('\\\\hline')\n",
    "print('Species &\\t $\\\\hat{\\\\phi}~mV$ & $c$ & $\\\\hat{K}$ \\\\\\\\')\n",
    "print('\\\\hline')\n",
    "\n",
    "concentration = phiData.ParRubXu16_conc()\n",
    "#concentration['H'] = 1e-7\n",
    "\n",
    "for i,spec in enumerate(s['species']):\n",
    "    if spec in concentration.keys():  \n",
    "        conc = 1e3*concentration[spec]\n",
    "        K_spec = np.exp(phi_est[i])/conc \n",
    "        print(f'{spec} & {int(round(1e3*phi_est[i]))} &  \\t{conc:.2} & {K_spec:.4} \\\\\\\\')\n",
    "#     else:\n",
    "#         print(f'{spec} &{phi_est[i]:.2} & -- & --\\\\\\\\')\n",
    "    \n",
    "print('\\\\hline')"
   ]
  }
 ],
 "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": 4
}