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.
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); energetic issues are not explicitly considered here although they are implicit in the bond graph models.
Gawthrop and Crampin (2018). 3. The bond graph approach is multi domain and can thus, for example, model electrochemical systems including neurodynamics (Gawthrop et al., 2017), redox reactions (Gawthrop, 2017) and transporters (Pan et al., 2019). 4. The bond graph approch is modular: subsystems can be connected using energy ports . 5. The chemostat concept Gawthrop and Crampin (2016) gives a convenient and flexible way of turning a closed system into an open system and analysing the concomitant pathway structure. 6. BondGraphTools provide a symbolic basis for describing and analysing bondgraphs within Python.
The bond graph analysis uses a number of Python modules:
## Paths
NeedPath=False
if NeedPath:
import sys
sys.path += ['/usr/lib/python3/dist-packages']
import BondGraphTools as bgt
import numpy as np
import IPython.display as disp
## Stoichiometric analysis
import stoich as st
## Export stoichiometry as bond graph
import stoichBondGraph as stbg
## Modular bond graphs
import modularBondGraph as mbg
## Extract stoichiometry from a CobraPy model
import CobraExtract as Extract
## Control outputs
quiet = True
chemformula = True
def printChem(chem):
Chem = ''
for c in chem:
Chem += f'\ch{{{c}}}, '
return Chem[:-2]
The key stoichiometric concept of pathways has already been given a bond graph interpretation <cite data-cite="Gaw17,GawCra17".
This note shows how:
Much remains to be done in exploiting the combination of the two approaches including:
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.
Chemostats provide a more flexible approach as they can be created without changing system structure.
In the following discussion, chemostats are added as appropriate to illustrate the various pathways.
In this case the ecoli core model is extracted from the CobraPy model: 'textbook'. This model corresponds to that discussed in the the textbook B.O. Palsson, Systems Biology (2015). 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.
Names of reactions and species are converted to upper case and the $c$ (cytosol) subscript removed for clarity.
sm = Extract.extract(Remove=['_C','__' ],
negReaction=['RPI','PGK','PGM','SUCOAS','FRD7'],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 FRD7 ( 23 ) by -1 Multiplying reaction PGK ( 54 ) by -1 Multiplying reaction PGM ( 56 ) by -1 Multiplying reaction RPI ( 65 ) by -1 Multiplying reaction SUCOAS ( 69 ) by -1
reaction = ['GLCPTS','PGI','PFK','FBP','FBA','TPI','GAPD','PGK','PGM','ENO','PYK']
s0 = Extract.choose(sm,reaction=reaction)
print('Reactions:', reaction)
disp.Latex(st.sprintrl(s0,chemformula=True))
Reactions: ['GLCPTS', 'PGI', 'PFK', 'FBP', 'FBA', 'TPI', 'GAPD', 'PGK', 'PGM', 'ENO', 'PYK']
s0['name'] = 'GLY_abg'
stbg.model(s0)
import GLY_abg
s = st.stoich(GLY_abg.model(),quiet=True)
## Sanity check
err = np.linalg.norm(s['N']-s0['N'])
print("Error:",err)
Error: 0.0
chemostats0 = ['ADP','ATP','H2O','PI','H']
sc0 = st.statify(s,chemostats=chemostats0)
chemostats = ['GLCD_E','PYR','NAD','NADH']
chemostats.extend(chemostats0)
sc = st.statify(s,chemostats=chemostats)
print(st.sprintp(sc))
sp = st.path(s,sc)
disp.Latex(st.sprintrl(sp,chemformula=True))
2 pathways 0: + PFK + FBP 1: + GLCPTS + PGI + PFK + FBA + TPI + 2 GAPD + 2 PGK + 2 PGM + 2 ENO + PYK
Pathway reactions:
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.
reaction = ['CS','ACONTA','ACONTB','ICDHYR','AKGDH','SUCOAS',
'FRD7',
'SUCDI','FUM','MDH',
'NADTRHD']
s0 = Extract.choose(sm,reaction=reaction)
disp.Latex(st.sprintrl(s0,chemformula=True))
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.
s0['name'] = 'TCA_abg'
stbg.model(s0)
import TCA_abg
s = st.stoich(TCA_abg.model(),quiet=True)
## Sanity check
err = np.linalg.norm(s['N']-s0['N'])
print("Error:",err)
Error: 0.0
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.
The overall reaction of the TCA cycle converts ACCOA to COA and \ch{CO2} , so these three species are also set as chemostats.
print('No chemostats')
print(st.sprintp(s))
chemostats0 = ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI','H'
,'Q8','Q8H2']
sc0 = st.statify(s,chemostats=chemostats0)
print('Basic chemostats',chemostats0)
print(st.sprintp(sc0))
chemostats = ['ACCOA','COA','CO2']
chemostats.extend(chemostats0)
sc = st.statify(s,chemostats=chemostats)
print('Chemostats',chemostats)
print(st.sprintp(sc))
No chemostats 1 pathways 0: - FRD7 + SUCDI Basic chemostats ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'H', 'Q8', 'Q8H2'] 1 pathways 0: - FRD7 + SUCDI Chemostats ['ACCOA', 'COA', 'CO2', 'ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'H', 'Q8', 'Q8H2'] 2 pathways 0: - FRD7 + SUCDI 1: + CS + ACONTA + ACONTB + ICDHYR + AKGDH + SUCOAS + FRD7 + FUM + MDH + NADTRHD
With only the basic chemostats, there is one internal cycle. This performs no conversions and so the reaction is empty.
sp0 = st.path(s,sc0)
disp.Latex(st.sprintrl(sp0,chemformula=True,split=10))
When the chemostats ['ACCOA','COA','CO2'] are added, the entire TCA cycle appears as a pathway and the overall reaction is
sp = st.path(s,sc)
disp.Latex(st.sprintrl(sp,chemformula=True,split=10))
Pathway reactions:
To take this example further, include the two reactions converting pyruvate to ACCOA.
reaction.extend(['PDH','PFL'])
s0 = Extract.choose(sm,reaction=reaction)
disp.Latex(st.sprintrl(s0,chemformula=True))
s0['name'] = 'PyrTCA_abg'
stbg.model(s0)
import PyrTCA_abg
s = st.stoich(PyrTCA_abg.model(),quiet=True)
## Sanity check
err = np.linalg.norm(s['N']-s0['N'])
print("Error:",err)
Error: 0.0
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
print('No chemostats')
print(st.sprintp(s,removeSingle=True))
chemostats0 = ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI','NADP', 'NADPH','H'
,'Q8','Q8H2','FOR']
sc0 = st.statify(s,chemostats=chemostats0)
print('Basic chemostats',chemostats0)
print(st.sprintp(sc0,removeSingle=True))
chemostats = ['PYR','CO2']
chemostats.extend(chemostats0)
sc = st.statify(s,chemostats=chemostats)
print('Chemostats:', printChem(chemostats))
print(st.sprintp(sc,removeSingle=True))
No chemostats 1 non-unit pathways 0: - FRD7 + SUCDI Basic chemostats ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI', 'NADP', 'NADPH', 'H', 'Q8', 'Q8H2', 'FOR'] 1 non-unit pathways 0: - FRD7 + SUCDI Chemostats: \ch{PYR}, \ch{CO2}, \ch{ADP}, \ch{ATP}, \ch{H2O}, \ch{NAD}, \ch{NADH}, \ch{PI}, \ch{NADP}, \ch{NADPH}, \ch{H}, \ch{Q8}, \ch{Q8H2}, \ch{FOR} 3 non-unit pathways 0: - FRD7 + SUCDI 1: + CS + ACONTA + ACONTB + ICDHYR + AKGDH + SUCOAS + FRD7 + FUM + MDH + PDH 2: + CS + ACONTA + ACONTB + ICDHYR + AKGDH + SUCOAS + FRD7 + FUM + MDH + PFL
sp = st.path(s,sc)
disp.Latex(st.sprintrl(sp,chemformula=True,split=10))
Pathway reactions:
reaction = ['NADH16','CYTBD']
s0 = Extract.choose(sm,reaction=reaction)
disp.Latex(st.sprintrl(s0,chemformula=True))
s0['name'] = 'ETC_abg'
stbg.model(s0)
import ETC_abg
##imp.reload(ETC_abg)
s = st.stoich(ETC_abg.model(),quiet=True)
err = np.linalg.norm(s['N']-s0['N'])
print("Error:",err)
s['reaction'] = s0['reaction']
Error: 0.0
print('No chemostats')
print(st.sprintp(s))
chemostats = ['NADH','NAD','O2','H2O','H','H_E']
sc = st.statify(s,chemostats=chemostats)
print('Chemostats',chemostats)
print(st.sprintp(sc))
No chemostats 0 pathways Chemostats ['NADH', 'NAD', 'O2', 'H2O', 'H', 'H_E'] 1 pathways 0: + 2 NADH16 + CYTBD
sp = st.path(s,sc)
disp.Latex(st.sprintrl(sp,chemformula=True))
Pathway reactions:
reaction = ['ATPS4R']
s0 = Extract.choose(sm,reaction=reaction)
disp.Latex(st.sprintrl(s0,chemformula=True))
s0['name'] = 'ATP_abg'
stbg.model(s0)
import ATP_abg
##imp.reload(ATP_abg)
s = st.stoich(ATP_abg.model(),quiet=True)
err = np.linalg.norm(s['N']-s0['N'])
print("Error:",err)
s['reaction'] = s0['reaction']
Error: 0.0
print('No chemostats')
print(st.sprintp(s))
chemostats = ['H','H_E','ATP','ADP','PI','H2O']
sc = st.statify(s,chemostats=chemostats)
print('Chemostats',chemostats)
print(st.sprintp(sc,removeSingle=False))
No chemostats 0 pathways Chemostats ['H', 'H_E', 'ATP', 'ADP', 'PI', 'H2O'] 1 pathways 0: + ATPS4R
sp = st.path(s,sc,removeSingle=False)
disp.Latex(st.sprintrl(sp,chemformula=True))
Pathway reactions:
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:
GlyTCA = bgt.new(name='GlyTCA') # Create system
Gly = GLY_abg.model()
TCA = PyrTCA_abg.model()
GlyTCA.add(Gly,TCA)
PYR is produced by Gly and consumed by TCA. ATP, ADP, PI, H, NAD, NADH , H2O are common to both modules.
common = ['PYR','ATP','ADP','PI','H','NAD','NADH','H2O']
print('Common species:', '{'+printChem(common)+'}')
mbg.unify(GlyTCA,common=common,quiet=quiet)
Common species: {\ch{PYR}, \ch{ATP}, \ch{ADP}, \ch{PI}, \ch{H}, \ch{NAD}, \ch{NADH}, \ch{H2O}}
## Stoichiometry of unified module
s = st.stoich(GlyTCA,quiet=True)
## Create BG
s['name'] = 'GlyTCA_abg'
stbg.model(s)
import GlyTCA_abg
chemostats0 = ['ADP', 'ATP', 'H2O', 'NAD', 'NADH', 'PI','H'
,'Q8','Q8H2','FOR']
chemostats = ['GLCD_E','CO2']
chemostats.extend(chemostats0)
sc = st.statify(s,chemostats=chemostats)
print('Chemostats:', printChem(chemostats))
print(st.sprintp(sc))
Chemostats: \ch{GLCD_E}, \ch{CO2}, \ch{ADP}, \ch{ATP}, \ch{H2O}, \ch{NAD}, \ch{NADH}, \ch{PI}, \ch{H}, \ch{Q8}, \ch{Q8H2}, \ch{FOR} 4 pathways 0: + PFK + FBP 1: - FRD7 + SUCDI 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 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
sp = st.path(s,sc)
disp.Latex(st.sprintrl(sp,chemformula=True))
Pathway reactions:
The above methods are brought together to generate the metabolic pathways by combining the modules describing: 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.
Met = bgt.new(name='Met') # Create systemGlyTCA
GlyTCA = GlyTCA_abg.model()
ETC = ETC_abg.model()
ATP = ATP_abg.model()
Met.add(GlyTCA,ETC,ATP)
## Unify species common to modules
common = ['ATP','ADP','PI','H','H_E','NAD','NADH','H2O','Q8','Q8H2']
print('Common species:', '{'+printChem(common)+'}')
mbg.unify(Met,common=common,quiet=True)
Common species: {\ch{ATP}, \ch{ADP}, \ch{PI}, \ch{H}, \ch{H_E}, \ch{NAD}, \ch{NADH}, \ch{H2O}, \ch{Q8}, \ch{Q8H2}}
s = st.stoich(Met,quiet=True)
disp.Latex(st.sprintrl(s,chemformula=True))
chemostats0 = ['ADP','ATP','H2O','PI','H']
chemostats = ['GLCD_E','CO2','O2']
chemostats.extend(chemostats0)
sc = st.statify(s,chemostats=chemostats)
print('Chemostats:', printChem(chemostats))
print(st.sprintp(sc))
Chemostats: \ch{GLCD_E}, \ch{CO2}, \ch{O2}, \ch{ADP}, \ch{ATP}, \ch{H2O}, \ch{PI}, \ch{H} 3 pathways 0: + PFK + FBP 1: - FRD7 + SUCDI 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
sp = st.path(s,sc)
disp.Latex(st.sprintrl(sp,chemformula=True))
There are three pathways:
by Palsson (2015).