*Note: this is the EcoliCoreModelEnergy.ipynb notebook. The PDF version "The Escherichia coli Core Model: Modular Energetic Bond Graph Analysis of Glycolysis and Pentose Phosphate Pathways" is available here.*
As discussed in a companion notebook, the Network Thermodynamics/Bond Graph approach of Oster, Perlson and Katchalsky (1971,1973) extended by Gawthrop and Crampin (2014,2016,2017) to modelling biomolecular systems of interest to systems biologists developed independently from the stoichiometric approach .
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. This was pointed out by Cellier and Greifeneder (2009). This means that the two approaches are complementary and each can build on the strengths of the other.
In particular, as discussed here, the Bond Graph approach adds energy to stoichiometry.
This notebook focuses on building modular models of metabolism and consequent pathway analysis based on the Escherichia coli Core Model (Orth, Fleming and Palsson,2010); in particular, the Glycolysis and Pentose Phosphate portion is extracted and analysed. Following the discussion in the textbook of Garrett and Grisham (2017), section 22.6d, various possible pathways are examined by choosing appropriate chemostats and flowstats. (Gawthrop and Crampin, 2018)
Assuming steady-state conditions, the corresponding pathway potentials (Gawthrop 2017) are derived.
The bond graph analysis uses a number of Python modules:
## Maths library
import numpy as np
## BG tools
import BondGraphTools as bgt
## BG stoichiometric utilities
import stoich as st
## Stoichiometric conversion
import CobraExtract as Extract
import stoichBondGraph as stbg
## Potentials
import phiData
## Faraday constant
import scipy.constants as con
F = con.physical_constants['Faraday constant'][0]
## Display
import IPython.display as disp
## Allow output from within functions
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
## Units etc
factor = 1
units = ' V'
## Control output
quiet = True
computePhi = True
showMu = True
To perform energetic analysis it is necessary to have values of the chemical potential of the species involved. One way of this is to use experimentally derived value of species potentials at standard conditions and then derive potentials corresponding to the concentrations of the species. Another approach used here, is to take experimental values of reaction potentials $\Phi$ (Park et al., 2016) and derive a consistent set of species potentials $\phi$ using $\phi = -N^\dagger \Phi$ where $N$ is the stoichiometric matrix of the reaction system and $\dagger$ denotes pseudo inverse.
def getPhi(s):
"""Extract phi for given system using
Reaction potentials from ParRubXu16"""
## Reaction potentials from ParRubXu16
PHI = phiData.Phi_ParRubXu16()
Phi_reac = PHI['Ecoli']
## Reaction potential (33.9e3) from GarGri17
Phi_reac['GLCPTS'] = 33.9e3/F - Phi_reac['PYK']
print('Setting Phi for reaction GLCPTS to', int(Phi_reac['GLCPTS']*1000),'mV.')
Phi = np.zeros((len(s['reaction']),1))
for i,reac in enumerate(s['reaction']):
if reac in Phi_reac.keys():
Phi[i] = Phi_reac[reac]
else:
min = 0.01 # 10mV
print('Setting Phi for reaction','\\ch{'+reac+'}','to', min*1000, 'mV. \n')
Phi[i] = min
pinvN = np.linalg.pinv(s['N'].T)
phi = -pinvN@Phi
return Phi,phi
sm = Extract.extract(cobraname='textbook',Remove=['_C','__' ], negReaction=['RPI'], quiet=quiet)
Extracting stoichiometric matrix from: textbook Cobra Model name: e_coli_core BondGraphTools name: e_coli_core_abg Extract.Integer only handles one non-integer per reaction Multiplying reaction BIOMASS_ECOLIORE ( 12 ) by 0.6684491978609626 to avoid non-integer species 3PG ( 2 ) Multiplying reaction CYTBD ( 15 ) by 2.0 to avoid non-integer species O2 ( 55 ) Multiplying reaction RPI ( 65 ) by -1
name = 'GlyPPP_abg'
reaction = ['GLCPTS','PGI','PFK','FBA','TPI','GAPD','PGK','PGM','ENO','PYK']
reaction += ['G6PDH2R','PGL','GND','RPI','TKT2','TALA','TKT1','RPE']
sGlyPPP = Extract.choose(sm,reaction=reaction)
Phi,phi = getPhi(sGlyPPP)
sGlyPPP['name'] = name
stbg.model(sGlyPPP)
import GlyPPP_abg
Setting Phi for reaction GLCPTS to 277 mV. Setting Phi for reaction \ch{G6PDH2R} to 10.0 mV. Setting Phi for reaction \ch{PGL} to 10.0 mV.
See Gawthrop (2017) for a discussion of these two quantities.
disp.Latex(st.sprintrl(sGlyPPP,chemformula=True,Phi=Phi,showMu=showMu))
## Analyse pathways defined by chemostats and flowstats
def ch(name):
return '\\ch{'+name+'}'
def energetics(s,sp,phi):
"""Reaction energetics.
"""
## Phi for all reactions
Phi = -s['N'].T@phi
##Phi for pathway
## I is the relevant indices of phi
I = []
for spec in sp['species']:
i = s['species'].index(spec)
I.append(i)
Phip = -sp['N'].T@phi[I]
return Phi,Phip
def pathway(bg,phi,chemostats,flowstats=[],computePhi=False,verbose=False):
""" Analyse pathways
"""
print('Chemostats:',sorted(chemostats))
print('Flowstats:', sorted(flowstats))
## Stoichiometry
## Create stoichiometry from bond graph.
s = st.stoich(bg,quiet=True)
## Stoichiometry with chemostats
sc = st.statify(s,chemostats=chemostats,flowstats=flowstats)
## Pathway stoichiometry
sp = st.path(s,sc)
## Print info
if verbose:
for stat in sorted(chemostats):
print(ch(stat)+',')
## Energetics
if computePhi:
Phi,Phip = energetics(s,sp,phi)
#print('Phi units: kJ/mol')
# fac = -F/1000
# units='~\si{\kilo\joule\per\mol}'
units = '~\si{\volt}'
print(st.sprintp(sc))
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=Phip,showMu=showMu))
#return s,sc,sp,Phi*fac,Phip*fac,units
return s,sc,sp,Phip
else:
print(st.sprintrl(sp,chemformula=True))
Phip = 0
return s,sc,sp,Phip
print('Glycolysis')
import GlyPPP_abg
chemostats = ['H2O','H']
chemostats += ['ADP','ATP','PI']
chemostats += ['G6P','PYR','NAD','NADH']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
Glycolysis Chemostats: ['ADP', 'ATP', 'G6P', 'H', 'H2O', 'NAD', 'NADH', 'PI', 'PYR'] Flowstats: [] 1 pathways 0: + PGI + PFK + FBA + TPI + 2 GAPD - 2 PGK - 2 PGM + 2 ENO + 2 PYK
\citet[\S~18.2]{GarGri17}.
that the reaction proceeds in the forward direction.
print('R5P and NADPH generation')
chemostats = ['H2O','H']
chemostats += ['NADP','NADPH','CO2']
chemostats += ['G6P','R5P']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
R5P and NADPH generation Chemostats: ['CO2', 'G6P', 'H', 'H2O', 'NADP', 'NADPH', 'R5P'] Flowstats: [] 1 pathways 0: + G6PDH2R + PGL + GND + RPI
\ch{NADPH} synthesis discussed in comment 1 of Garrett and Grisham (2017), p787.
that the reaction proceeds in the forward direction.
print('R5P generation')
chemostats = ['H2O','H']
chemostats += ['G6P','R5P']
chemostats += ['ADP','ATP']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
R5P generation Chemostats: ['ADP', 'ATP', 'G6P', 'H', 'H2O', 'R5P'] Flowstats: [] 1 pathways 0: - 5 PGI - PFK - FBA - TPI - 4 RPI + 2 TKT2 + 2 TALA + 2 TKT1 + 4 RPE
Garrett and Grisham (2017), p787.
that the reaction proceeds in the reverse direction.
import GlyPPP_abg
chemostats = ['H2O','H']
chemostats += ['G6P']
chemostats += ['NADP','NADPH','CO2']
chemostats += ['ATP','ADP']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NADP', 'NADPH'] Flowstats: [] 1 pathways 0: - 5 PGI - PFK - FBA - TPI + 6 G6PDH2R + 6 PGL + 6 GND + 2 RPI + 2 TKT2 + 2 TALA + 2 TKT1 + 4 RPE
\ch{NADPH} synthesis discussed in comment 3 of Garrett and Grisham (2017), p787.
that the reaction proceeds in the forward direction.
import GlyPPP_abg
chemostats = ['H2O','H']
chemostats += ['NADP','NADPH','CO2']
chemostats += ['G6P']
chemostats += ['ADP','ATP','PI']
chemostats += ['PYR','NAD','NADH']
flowstats = ['PGI']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,flowstats=flowstats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NAD', 'NADH', 'NADP', 'NADPH', 'PI', 'PYR'] Flowstats: ['PGI'] 2 pathways 0: + PGI 1: + 2 PFK + 2 FBA + 2 TPI + 5 GAPD - 5 PGK - 5 PGM + 5 ENO + 5 PYK + 3 G6PDH2R + 3 PGL + 3 GND + RPI + TKT2 + TALA + TKT1 + 2 RPE
\ch{NADPH} and \ch{ATP} synthesis discussed in comment 4 of Garrett and Grisham (2017), p787.
that the reaction proceeds in the forward direction.
The pathways may also be isolated by using appropriate (zero-flow) flowstats. The comments for each section are the same as in the previous section.
import GlyPPP_abg
Chemostats = ['G6P','ADP','ATP','CO2','H','H2O','NAD','NADH','NADP','NADPH','PI','PYR']
by replacing the two connecting reactions (G6PDH2R and TKT2) by flowstats.
print('Glycolysis')
chemostats = Chemostats
flowstats = ['G6PDH2R','TKT2']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,flowstats=flowstats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
Glycolysis Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NAD', 'NADH', 'NADP', 'NADPH', 'PI', 'PYR'] Flowstats: ['G6PDH2R', 'TKT2'] 3 pathways 0: + PGI + PFK + FBA + TPI + 2 GAPD - 2 PGK - 2 PGM + 2 ENO + 2 PYK 1: + G6PDH2R 2: + TKT2
product \ch{R5P} is added to the chemostat list.
print('R5P and NADPH generation')
chemostats = Chemostats + ['R5P']
flowstats = ['PGI','TKT2']
#s,sc,sp,Phip,Phi,Phip,units = pathway(GlyPPP_abg.model(),phi,chemostats,flowstats=flowstats,computePhi=computePhi)
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,flowstats=flowstats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
R5P and NADPH generation Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NAD', 'NADH', 'NADP', 'NADPH', 'PI', 'PYR', 'R5P'] Flowstats: ['PGI', 'TKT2'] 3 pathways 0: + PGI 1: + G6PDH2R + PGL + GND + RPI 2: + TKT2
product \ch{R5P} is added to the chemostat list.
print('R5P generation')
chemostats = Chemostats + ['R5P']
flowstats = ['GAPD','G6PDH2R']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,flowstats=flowstats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
R5P generation Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NAD', 'NADH', 'NADP', 'NADPH', 'PI', 'PYR', 'R5P'] Flowstats: ['G6PDH2R', 'GAPD'] 3 pathways 0: + GAPD 1: + G6PDH2R 2: - 5 PGI - PFK - FBA - TPI - 4 RPI + 2 TKT2 + 2 TALA + 2 TKT1 + 4 RPE
print('NADPH generation')
chemostats = Chemostats
flowstats = ['GAPD']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,flowstats=flowstats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
NADPH generation Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NAD', 'NADH', 'NADP', 'NADPH', 'PI', 'PYR'] Flowstats: ['GAPD'] 2 pathways 0: + GAPD 1: - 5 PGI - PFK - FBA - TPI + 6 G6PDH2R + 6 PGL + 6 GND + 2 RPI + 2 TKT2 + 2 TALA + 2 TKT1 + 4 RPE
This pathway is isolated by setting PGI as flowstat.
print('NADPH and ATP generation')
chemostats = Chemostats
flowstats = ['PGI']
s,sc,sp,Phip = pathway(GlyPPP_abg.model(),phi,chemostats,flowstats=flowstats,computePhi=computePhi)
disp.Latex(st.sprintrl(sp,chemformula=True,Phi=factor*Phip,showMu=showMu))
NADPH and ATP generation Chemostats: ['ADP', 'ATP', 'CO2', 'G6P', 'H', 'H2O', 'NAD', 'NADH', 'NADP', 'NADPH', 'PI', 'PYR'] Flowstats: ['PGI'] 2 pathways 0: + PGI 1: + 2 PFK + 2 FBA + 2 TPI + 5 GAPD - 5 PGK - 5 PGM + 5 ENO + 5 PYK + 3 G6PDH2R + 3 PGL + 3 GND + RPI + TKT2 + TALA + TKT1 + 2 RPE