'MASS' SIMULATOR OF GLYCOLYSIS and PFK Enzyme

This Workbook sets up, analyses, and simulates and graphs the equations that describe the  glycolytic pathway with a PFK enzyme component.  It corresponds to chapter in Ch 14 in SB2.

Load functions

Off[General::spell]; Off[General::spell1];
SetDirectory["C:/users/aabordba.alanine0/desktop/website"];

<<functions.m

Section I: 	System Definition and Setting up the Equations

Define nodes, links and the stoich matrix

Define the the compounds and their elemental composition

variables = {glu,g6p,f6p, PFKR,PFKRAtp,PFKRAtpF6p,PFKR1,PFKR1Atp,PFKR1AtpF6p,PFKR2,PFKR2Atp,PFKR2AtpF6p,PFKR3,PFKR3Atp,PFKR3AtpF6p,PFKR4,PFKR4Atp,PFKR4AtpF6p,PFKT,PFKT1,PFKT2,PFKT3,PFKT4,fbp,dhap,gap,pg13,pg3,pg2,pep,pyr,lac,nad,nadh,amp,adp,atp,phos,h,h2o};
varnames={"glu","g6p","f6p"," PFKR","PFKRAtp","PFKRAtpF6p","PFKR1","PFKR1Atp","PFKR1AtpF6p","PFKR2","PFKR2Atp","PFKR2AtpF6p","PFKR3","PFKR3Atp","PFKR3AtpF6p","PFKR4","PFKR4Atp","PFKR4AtpF6p","PFKT","PFKT1","PFKT2","PFKT3","PFKT4","fbp","dhap","gap","pg13","pg3","pg2","pep", "pyr","lac","nad","nadh", "amp","adp","atp","phos","h", "h2o"};
elemnames={"C","H","O","P","N", "S","NAD","PFK"};
MatrixForm[elem=(({{6, 6, 6, 0, 10, 16, 10, 20, 26, 20, 30, 36, 30, 40, 46, 40, 50, 56, 0, 10, 20, 30, 40, 6, 3, 3, 3, 3, 3, 3, 3, 3, 0, 0, 10, 10, 10, 0, 0, 0}, {12, 11, 11, 0, 13, 24, 13, 26, 37, 26, 39, 50, 39, 52, 63, 52, 65, 76, 0, 13, 26, 39, 52, 10, 5, 5, 4, 4, 4, 2, 3, 5, 0, 1, 13, 13, 13, 1, 1, 2}, {6, 9, 9, 0, 13, 22, 7, 20, 29, 14, 27, 36, 21, 34, 43, 28, 41, 50, 0, 13, 26, 39, 52, 12, 6, 6, 10, 7, 7, 6, 3, 3, 0, 0, 7, 10, 13, 4, 0, 1}, {0, 1, 1, 0, 3, 4, 1, 4, 5, 2, 5, 6, 3, 6, 7, 4, 7, 8, 0, 3, 6, 9, 12, 2, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 1, 2, 3, 1, 0, 0}, {0, 0, 0, 0, 5, 5, 5, 10, 10, 10, 15, 15, 15, 20, 20, 20, 25, 25, 0, 5, 10, 15, 20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}))];
MatrixForm[elem,TableHeadings->{elemnames,varnames}]
TeXForm[MatrixForm[elem,TableHeadings->{elemnames,varnames}]];

Define the reaction names

fluxes={vhk,vpgi,vPFKR01,vPFKR02,vPFKR03,vPFKR10,vPFKR11,vPFKR12,vPFKR13,vPFKR20,vPFKR21,vPFKR22,vPFKR23,vPFKR30,vPFKR31,vPFKR32,vPFKR33,vPFKR40,vPFKR41,vPFKR42,vPFKR43,vPFKL,vPFKT1,vPFKT2,vPFKT3,vPFKT4,vtpi,vald,vgapdh,vpgk,vpglm,veno,vpk,vldh,vamp,vapk,vpyr,vlac,vatp,vnadh,vgluin,vampin,vh,vh2o};

fluxnames={"vhk","vpgi","vPFKR01","vPFKR02","vPFKR03","vPFKR10","vPFKR11","vPFKR12","vPFKR13","vPFKR20","vPFKR21","vPFKR22","vPFKR23","vPFKR30","vPFKR31","vPFKR32","vPFKR33","vPFKR40","vPFKR41","vPFKR42","vPFKR43","vPFKL","vPFKT1","vPFKT2","vPFKT3","vPFKT4","vtpi","vald","vgapdh","vpgk","vpglm","veno","vpk","vldh","vamp","vapk","vpyr","vlac","vatp","vnadh","vgluin","vampin","vh","vh2o"};

knames={"khk","kpgi","kPFKR01","kPFKR02","kPFKR03","kPFKR10","kPFKR11","kPFKR12","kPFKR13","kPFKR20","kPFKR21","kPFKR22","kPFKR23","kPFKR30","kPFKR31","kPFKR32","kPFKR33","kPFKR40","kPFKR41","kPFKR42","kPFKR43","kPFKL","kPFKT1","kPFKT2","kPFKT3","kPFKT4","ktpi","kald","kgapdh","kpgk","kpglm","keno","kpk","kldh","kamp","kapk","kpyr","klac","katp","knadh","kgluin","kampin","kh","kh2o"};

Knames={"Khk","Kpgi","KPFKR01","KPFKR02","KPFKR03","KPFKR10","KPFKR11","KPFKR12","KPFKR13","KPFKR20","KPFKR21","KPFKR22","KPFKR23","KPFKR30","KPFKR31","KPFKR32","KPFKR33","KPFKR40","KPFKR41","KPFKR42","KPFKR43","KPFKL","KPFKT1","KPFKT2","KPFKT3","KPFKT4","Ktpi","Kald","Kgapdh","Kpgk","Kpglm","Keno","Kpk","Kldh","Kamp","Kapk","Kpyr","Klac","Katp","Knadh","Kgluin","Kampin","Kh","Kh2o"};

Define stoichiometric matrix; there are m=20 compounds and n=21 reactions.

MatrixForm[stoich=(({{-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, {1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0, 1, 0, 0}, {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, -2, 0, 0, 1, 0, 0, 0, 0, 0}, {-1, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, -1, -1, -1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 1, 1, 0, 0, -1, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1}}))];
MatrixForm[stoich,TableHeadings->{varnames,fluxnames}]
TeXForm[MatrixForm[stoich,TableHeadings->{varnames,fluxnames}]];
dim=Dimensions[stoich,2];
rank=MatrixRank[stoich];

Check elemental balancing--exchange reactions not balanced

MatrixForm[elem.stoich,TableHeadings->{elemnames,fluxnames}]
TeXForm[MatrixForm[elem.stoich,TableHeadings->{elemnames,fluxnames}]];

The null spaces

The null spaces:  determining the pathways and the pools
	the rank of S is r=37, 
	there are two time invariant pools, Dim[Left Null]] =  40-r=40-37 = 3
	there are three time pathways, Dim[Null]] =  n-r=44-37 = 7

null=(({{1, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0}, {1, 0, 1, 0, 0, 0, 0}, {1, 0, 1, 0, 0, 0, 1}, {1, 0, 1, 0, 0, 1, 0}, {1, 0, 1, 0, 1, 0, 0}, {1, 0, 1, 1, 0, 0, 0}})).NullSpace[stoich];
pathwaynames={"pyruvate reduction","amp in/out","glycolysis 1", "glycolysis 2","glycolysis 3", "glycolysis 4","glycolysis 5"};
MatrixForm[null,TableHeadings->{pathwaynames,fluxnames}]
TeXForm[MatrixForm[null,TableHeadings->{pathwaynames,fluxnames}]];

leftMultMat = (({{1, 0, 8}, {0, 1, 0}, {0, 0, 1}}));
leftnull=leftMultMat.NullSpace[Transpose[stoich]];
poolnames={"phosphate balance","nadh balance","total PFK"};
MatrixForm[leftnull,TableHeadings->{poolnames,varnames}]

Check elemental balancing of pathways and the elemental composition of the pools

MatrixForm[elem.stoich.Transpose[null],TableHeadings->{elemnames,pathwaynames}]
TeXForm[MatrixForm[elem.stoich.Transpose[null],TableHeadings->{elemnames,pathwaynames}]];
MatrixForm[elem.Transpose[leftnull],TableHeadings->{elemnames,poolnames}]

Compute the steady state concentrations and flux

Define condentration variables: pool sizes and external concentrations

infty=1000000;
NADHtotal = 0.089;
pyrplasma=0.06;
lacplasma=1.0;
ampplasma=0.0001;
hplasma = 10.^(-4.2);
h2oplasma=1;
kh=infty/10; kh2o=infty/10;

Define equilibrium constants--the irreversible are set very high,  arbitrarily at 1,000,000

Katp=infty;Kamp=infty;Kapk=1.65;Khk=850;Kpgi=0.41;Ktpi=1/17.5;Kald=0.082;Kgapdh=0.0179;Kpgk=1800;Kpglm=1/6.8;Keno=1/.59;Kpk=363000;Kldh=26300;Kpyr=1.;Klac=1.;Knadh=infty;Kgluin=infty;Kampin=infty;Kh=1.;Kh2o=1.;

Ka=.033; 
Ki=0.1;
Lallo=0.00107;

KATP=0.068;KF6P=0.1; KPFK=infty;

KPFKR01 = 1/KATP;
KPFKR02 = 1/KF6P;
KPFKR03 = KPFK;
KPFKR10 = 1/((1/4)*Ka);
KPFKR11 = 1/KATP;
KPFKR12 = 1/KF6P;
KPFKR13 = KPFK;
KPFKR20 = 1/((2/3)*Ka);
KPFKR21 = 1/KATP;
KPFKR22 = 1/KF6P;
KPFKR23 = KPFK;
KPFKR30 = 1/((3/2)*Ka);
KPFKR31 = 1/KATP;
KPFKR32 = 1/KF6P;
KPFKR33 =KPFK;
KPFKR40 = 1/(4*Ka);
KPFKR41 = 1/KATP;
KPFKR42 = 1/KF6P;
KPFKR43 = KPFK;
KPFKT1 = 1/((1/4)*Ki);
KPFKT2 = 1/((2/3)*Ki);
KPFKT3 = 1/((3/2)*Ki);
KPFKT4 = 1/(4*Ki);
KPFKL=0.0011

Kvalues={Khk,Kpgi,KPFKR01,KPFKR02,KPFKR03,KPFKR10,KPFKR11,KPFKR12,KPFKR13,KPFKR20,KPFKR21,KPFKR22,KPFKR23,KPFKR30,KPFKR31,KPFKR32,KPFKR33,KPFKR40,KPFKR41,KPFKR42,KPFKR43,KPFKL,KPFKT1,KPFKT2,KPFKT3,KPFKT4,Ktpi,Kald,Kgapdh,Kpgk,Kpglm,Keno,Kpk,Kldh,Kamp,Kapk,Kpyr,Klac,Katp,Knadh,Kgluin,Kampin,Kh,Kh2o};
MatrixForm[Kvalues,TableHeadings->{Knames}]
Length[Kvalues];
TeXForm[MatrixForm[Kvalues,TableHeadings->{Knames}]];

Define steady state concentration values

atpstst=1.6;adpstst=0.29;ampstst=Kapk*adpstst*adpstst/atpstst;

glustst=1.;g6pstst=0.0486;f6pstst=0.0198;fbpstst=0.0146;dhapstst=0.16;gapstst=0.00728;pg13stst=0.000243; pg3stst=0.0773;pg2stst=0.0113;pepstst=0.017; pyrstst = 0.060301; lacstst=1.36;

phosstst=2.5; nadhstst= 0.0301;nadstst=NADHtotal-nadhstst;

PFKtot=33.*10^-6;
(* Defined in solutionPFKstst.nb module *)
enzymestates ={3.7106552487703815`*^-8,1.1264628584498815`*^-7,2.10241288178717`*^-8,3.900826330269862`*^-7,1.1841940799454375`*^-6,2.210161541978762`*^-7,1.5377788798860717`*^-6,4.668315099534904`*^-6,8.712871203769399`*^-7,2.694316745800388`*^-6,8.179277080643446`*^-6,1.5265676421604298`*^-6,1.7702502993891606`*^-6,5.374040644391512`*^-6,1.0030026461382197`*^-6,4.081720773647421`*^-11,2.612301295134349`*^-9,6.269523108322439`*^-8,6.687491315543933`*^-7,2.6749965262175732`*^-6};

varstst=ConstantArray[0,Length[variables]];
Do[varstst[[i+3]]=enzymestates[[i]],{i,1,Length[enzymestates]}]
R0tot=Sum[enzymestates[[i]],{i,1,3}];
R1tot=Sum[enzymestates[[i]],{i,4,6}];
R2tot=Sum[enzymestates[[i]],{i,7,9}];
R3tot=Sum[enzymestates[[i]],{i,10,12}];
R4tot=Sum[enzymestates[[i]],{i,13,15}];
Csum=R0tot+R1tot+R2tot+R3tot+R4tot;
fractions = {Rfrac,Tfrac,Cfrac,R0tot/Csum,R1tot/Csum,R2tot/Csum,R3tot/Csum,R4tot/Csum};

gluin=1.12;
ststloadnadh=0.2*gluin;
ampin=.014;
MatrixForm[givenststfluxes=(({{ststloadnadh, ampin, gluin*fractions[[4]], gluin*fractions[[5]], gluin*fractions[[6]], gluin*fractions[[7]], gluin*fractions[[8]]}}))];
TeXForm[MatrixForm[givenststfluxes]];
MatrixForm[null];
solution=Flatten[Transpose[givenststfluxes.null]];
MatrixForm[solution,TableHeadings->{fluxnames}]

hstst=hplasma+solution[[43]]/kh;
pHext=-Log[10,hplasma];
pHint=-Log[10,hstst];

h2ostst=h2oplasma+solution[[44]]/kh2o;

varstst[[1;;3]]={glustst,g6pstst,f6pstst};
varstst[[24;;40]]={fbpstst,dhapstst,gapstst,pg13stst,pg3stst,pg2stst,pepstst,pyrstst,lacstst,nadstst,nadhstst,ampstst,adpstst,atpstst,phosstst,hstst,h2ostst};

Rsum=Sum[enzymestates[[i]],{i,1,15}];
Tsum=Sum[enzymestates[[i]],{i,16,20}];
PFKtot-(Rsum+Tsum)

MatrixForm[varstst,TableHeadings->{varnames}]
TeXForm[MatrixForm[varstst,TableHeadings->{varnames}]];
BarChart[varstst,BarLabels->varnames];
BarChart[Log[10,varstst]];

Fractional states of enzyme

Rfrac=Rsum/PFKtot;
Tfrac=Tsum/PFKtot;
Cfrac=(enzymestates[[3]]+enzymestates[[6]]+enzymestates[[9]]+enzymestates[[12]]+enzymestates[[15]])/PFKtot;

fracnames = {"Rfrac","Tfrac","Cfrac","RAFfrac","R1AFfrac","R2AFfrac","R3AFfrac","R4AFfrac"};
MatrixForm[fractions,TableHeadings->{fracnames}]

Calculate PERCs with new steady state concentrations and define the rate laws

Calculate mass action ratios

Ghk=g6pstst*adpstst/(glustst*atpstst);
Gpgi=f6pstst/g6pstst;
Gpfk=fbpstst*adpstst/(f6pstst*atpstst);
Gtpi=gapstst/dhapstst;
Gald=gapstst*dhapstst/fbpstst;
Ggapdh=pg13stst*nadhstst/(phosstst*gapstst*nadstst);
Gpgk=pg3stst*atpstst/(pg13stst*adpstst);
Gpglm=pg2stst/pg3stst;
Geno=pepstst/pg2stst;
Gpk=pyrstst*atpstst/(pepstst*adpstst);
Gldh=lacstst*nadstst/(pyrstst*nadhstst);
Gamp=ampplasma/ampstst;
Gapk=atpstst*ampstst/(adpstst*adpstst);
Gpyr=pyrplasma/pyrstst;
Glac=lacplasma/lacstst;
Gatp=adpstst*phosstst/atpstst;
Gnadh=nadstst/nadhstst;
Ggluin=0;
Gampin=0;
Gh=1.;
Gh2o=1.;
GPFKR01=enzymestates[[2]]/(atpstst*enzymestates[[1]]);
GPFKR02=enzymestates[[3]]/(f6pstst*enzymestates[[2]]);
GPFKR03=(fbpstst*adpstst*enzymestates[[1]])/enzymestates[[3]];
GPFKR10=enzymestates[[4]]/(ampstst*enzymestates[[1]]);
GPFKR11=enzymestates[[5]]/(atpstst*enzymestates[[4]]);
GPFKR12=enzymestates[[6]]/(f6pstst*enzymestates[[5]]);
GPFKR13=(fbpstst*adpstst*enzymestates[[4]])/enzymestates[[6]];
GPFKR20=enzymestates[[7]]/(ampstst*enzymestates[[4]]);
GPFKR21=enzymestates[[8]]/(atpstst*enzymestates[[7]]);
GPFKR22=enzymestates[[9]]/(f6pstst*enzymestates[[8]]);
GPFKR23=(fbpstst*adpstst*enzymestates[[7]])/enzymestates[[9]];
GPFKR30=enzymestates[[10]]/(ampstst*enzymestates[[7]]);
GPFKR31=enzymestates[[11]]/(atpstst*enzymestates[[10]]);
GPFKR32=enzymestates[[12]]/(f6pstst*enzymestates[[11]]);
GPFKR33=(fbpstst*adpstst*enzymestates[[10]])/enzymestates[[12]];
GPFKR40=enzymestates[[13]]/(ampstst*enzymestates[[10]]);
GPFKR41=enzymestates[[14]]/(atpstst*enzymestates[[13]]);
GPFKR42=enzymestates[[15]]/(f6pstst*enzymestates[[14]]);
GPFKR43=(fbpstst*adpstst*enzymestates[[13]])/enzymestates[[15]];
GPFKL=enzymestates[[16]]/enzymestates[[1]];
GPFKT1=enzymestates[[17]]/(atpstst*enzymestates[[16]]);
GPFKT2=enzymestates[[18]]/(atpstst*enzymestates[[17]]);
GPFKT3=enzymestates[[19]]/(atpstst*enzymestates[[18]]);
GPFKT4=enzymestates[[20]]/(atpstst*enzymestates[[19]]);

Gvalues={Ghk,Gpgi,GPFKR01,GPFKR02,GPFKR03,GPFKR10,GPFKR11,GPFKR12,GPFKR13,GPFKR20,GPFKR21,GPFKR22,GPFKR23,GPFKR30,GPFKR31,GPFKR32,GPFKR33,GPFKR40,GPFKR41,GPFKR42,GPFKR43,GPFKL,GPFKT1,GPFKT2,GPFKT3,GPFKT4,Gtpi,Gald,Ggapdh,Gpgk,Gpglm,Geno,Gpk,Gldh,Gamp,Gapk,Gpyr,Glac,Gatp,Gnadh,Ggluin,Gampin,Gh,Gh2o};

MatrixForm[Transpose[{0.001*Round[1000*Gvalues],0.001*Round[1000*Kvalues],0.001*Round[1000*Gvalues/Kvalues]}],TableHeadings->{fluxnames,{"G","K","G/K"}}];

TeXForm[MatrixForm[Transpose[{0.001*Round[1000*Gvalues],0.001*Round[1000*Kvalues],0.001*Round[1000*Gvalues/Kvalues]}],TableHeadings->{fluxnames,{"G","K","G/K"}}]];

Determine kinetic constants from the given fluxes and concentration values, set large dummy kinetic constants for PFK

khk=solution[[1]]/(glustst*atpstst-g6pstst*adpstst/Khk);
kpgi =solution[[2]]/(g6pstst-f6pstst/Kpgi);


ktpi=solution[[27]]/(dhapstst-gapstst/Ktpi);
kald=solution[[28]]/(fbpstst-gapstst*dhapstst/Kald);
kgapdh=solution[[29]]/(phosstst*gapstst*nadstst-nadhstst*pg13stst/Kgapdh);
kpgk=solution[[30]]/(pg13stst*adpstst-pg3stst*atpstst/Kpgk);
kpglm=solution[[31]]/(pg3stst-pg2stst/Kpglm);
keno=solution[[32]]/(pg2stst-pepstst/Keno);
kpk=solution[[33]]/(pepstst*adpstst-pyrstst*atpstst/Kpk);
kldh =solution[[34]]/(pyrstst*nadhstst-lacstst*nadstst/Kldh);
kamp=solution[[35]]/(ampstst-ampplasma/Kamp);
kapk =10^6;
kpyr =solution[[37]]/(pyrstst-pyrplasma);
klac =solution[[38]]/(lacstst-lacplasma);
katp=solution[[39]]/(atpstst-adpstst*phosstst/Katp);
knadh =solution[[40]]/(nadhstst-nadstst/Knadh);
kgluin=0.;
kampin=0;

kPFKR01=kPFKR11=kPFKR21=kPFKR31=kPFKR41=200000/atpstst;
kPFKR02=kPFKR12=kPFKR22=kPFKR32=kPFKR42=infty/f6pstst;
kPFKR03=kPFKR13=kPFKR23=kPFKR33=kPFKR43=307447.5581772384;
kPFKR10=infty;
kPFKR20=infty;
kPFKR30=infty;
kPFKR40=infty;
kPFKL=infty^2;
kPFKT1=infty;
kPFKT2=infty;
kPFKT3=infty;
kPFKT4=infty;


kvalues={khk,kpgi,kPFKR01,kPFKR02,kPFKR03,kPFKR10,kPFKR11,kPFKR12,kPFKR13,kPFKR20,kPFKR21,kPFKR22,kPFKR23,kPFKR30,kPFKR31,kPFKR32,kPFKR33,kPFKR40,kPFKR41,kPFKR42,kPFKR43,kPFKL,kPFKT1,kPFKT2,kPFKT3,kPFKT4,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kamp,kapk,kpyr,klac,katp,knadh,kgluin,kampin,kh,kh2o};

BarChart[Log[10,Sort[kvalues,Greater]]];

Summary of reaction properties

MatrixForm[Transpose[{0.001*Round[1000*Gvalues],0.001*Round[1000*Kvalues],0.001*Round[1000*Gvalues/Kvalues],Flatten[kvalues]}],TableHeadings->{fluxnames,{"G","K","G/K","k"}}]
TeXForm[MatrixForm[kvalues,TableHeadings->{knames}]];

Determine kinetic constants from the given fluxes and concentration values

vhk=khk*(atp*glu-adp*g6p/Khk);
vpgi=kpgi*(g6p-f6p/Kpgi);
vtpi=ktpi*(dhap-gap/Ktpi);
vald=kald*(fbp-gap*dhap/Kald);
vgapdh=kgapdh*(phos*gap*nad-nadh*pg13/Kgapdh);
vpgk=kpgk*(pg13*adp-pg3*atp/Kpgk);
vpglm=kpglm*(pg3-pg2/Kpglm);
veno=keno*(pg2-pep/Keno);
vpk=kpk*(pep*adp-pyr*atp/Kpk);
vldh=kldh*(pyr*nadh-lac*nad/Kldh);
vamp=kamp*(amp-ampplasma/Kamp);
vapk=kapk*(adp*adp-atp*amp/Kapk);
vpyr =kpyr*(pyr-pyrplasma/Kpyr);
vlac =klac*(lac-lacplasma/Klac);
vatp= katp*(atp-adp*phos/Katp);
vnadh = knadh*(nadh-nad/Knadh);
vgluin=gluin;
vampin=ampin;
vh = kh*(h-hplasma/Kh);
vh2o = kh2o*(h2o-h2oplasma/Kh);

vgapdh=kgapdh*(phos*gap*nad-nadh*pg13/Kgapdh);
vldh=kldh*(pyr*nadh-lac*nad/Kldh);
vnadh = knadh*(nadh-nad/Knadh);

vPFKR01 = kPFKR01*(PFKR*atp-PFKRAtp/KPFKR01);
vPFKR02 = kPFKR02*(f6p*PFKRAtp-PFKRAtpF6p/KPFKR02);
vPFKR03 = kPFKR03*(PFKRAtpF6p-PFKR*fbp*adp /KPFKR03);
vPFKR10 = kPFKR10*(PFKR*amp-PFKR1/KPFKR10);
vPFKR11 = kPFKR11*(PFKR1*atp-PFKR1Atp/KPFKR11);
vPFKR12 = kPFKR12*(f6p*PFKR1Atp-PFKR1AtpF6p/KPFKR12);
vPFKR13 = kPFKR13*(PFKR1AtpF6p-PFKR1*fbp*adp/KPFKR13);
vPFKR20 = kPFKR20*(PFKR1*amp-PFKR2/KPFKR20);
vPFKR21 = kPFKR21*(PFKR2*atp-PFKR2Atp/KPFKR21);
vPFKR22 = kPFKR22*(f6p*PFKR2Atp-PFKR2AtpF6p/KPFKR22);
vPFKR23 = kPFKR23*(PFKR2AtpF6p-PFKR2*fbp*adp/KPFKR23);
vPFKR30 = kPFKR30*(PFKR2*amp-PFKR3/KPFKR30);
vPFKR31 = kPFKR31*(PFKR3*atp-PFKR3Atp/KPFKR31);
vPFKR32 = kPFKR32*(f6p*PFKR3Atp-PFKR3AtpF6p/KPFKR32);
vPFKR33 = kPFKR33*(PFKR3AtpF6p-PFKR3*fbp*adp/KPFKR33);
vPFKR40 = kPFKR40*(PFKR3*amp-PFKR4/KPFKR40);
vPFKR41 = kPFKR41*(PFKR4*atp-PFKR4Atp/KPFKR41);
vPFKR42 = kPFKR42*(f6p*PFKR4Atp-PFKR4AtpF6p/KPFKR42);
vPFKR43 = kPFKR43*(PFKR4AtpF6p-PFKR4*fbp*adp/KPFKR43);
vPFKL = kPFKL*(PFKR-PFKT/KPFKL);
vPFKT1 = kPFKT1*(PFKT*atp-PFKT1/KPFKT1);
vPFKT2 = kPFKT2*(PFKT1*atp-PFKT2/KPFKT2);
vPFKT3 = kPFKT3*(PFKT2*atp-PFKT3/KPFKT3);
vPFKT4 = kPFKT4*(PFKT3*atp-PFKT4/KPFKT4);


fluxstst= Flatten[fluxes/.sublist[variables,varstst]];
MatrixForm[0.001*Round[1000*fluxstst],TableHeadings->{fluxnames}]
BarChart[fluxstst,BarLabels->fluxnames];

Numerical checks

Calculate the flux balances

Chop[MatrixForm[0.001*Round[1000*stoich.fluxstst],TableHeadings-> {varnames}]]

Section II: 	Defining Pools and Ratios

Define the pools

Compute & characterize Pools: hard and soft pools (is a similar section needed for the pathways?)
	use the left null space for 'hard pools'
	'soft pools' defined based on biochemical intuition
	compute pool sizes based on the steady state concentrations

totalpoolnames={"GP+","GP-","AP+","AP-","GR+","GR-","nadh","P+","P-","Ptot","nadh-tot","PFK-R","PFK-L"};
MatrixForm[pools=(({{2, 3, 3, 0, 0, 3, 0, 0, 3, 0, 0, 3, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0, 4, 2, 2, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 2, 2, 0, 2, 2, 0, 2, 2, 0, 2, 2, 0, 2, 2, 0, 2, 4, 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 2, 2, 2, 4, 4, 4, 6, 6, 6, 8, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0}, {2, 2, 2, 0, 0, 2, 0, 0, 2, 0, 0, 2, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 2, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 1, 1, 0, 2, 1, 0, 2, 1, 0, 2, 1, 0, 2, 1, 0, 2, 1, 0, 2, 4, 6, 8, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 1, 0, 2, 1, 0, 2, 1, 0, 2, 1, 0, 2, 1, 0, 2, 1, 0, 2, 4, 6, 8, 2, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}))];
Dimensions[pools,2];
MatrixForm[poolsize = pools.varstst];
MatrixForm[Transpose[Append[Transpose[pools],poolsize]],TableHeadings->{totalpoolnames,Append[varnames,"size"]}]
TeXForm[MatrixForm[pools,TableHeadings->{totalpoolnames,varnames}]];
BarChart[poolsize,BarLabels->totalpoolnames]

Pool fluxes and their turnover times

Calculate the reaction that move the pools

MatrixForm[pstoich=pools.stoich,TableHeadings->{totalpoolnames,fluxnames}]

Compute the fluxes in and out of pools

ppstoich=pstoich;
Do[Do[ppstoich[[i,j]]=If[pstoich[[i,j]]>0,pstoich[[i,j]],0],{j,1,40}],{i,1,11}]
ppstoich//MatrixForm;

pnstoich=pstoich;
Do[Do[pnstoich[[i,j]]=If[pstoich[[i,j]]<0,pstoich[[i,j]],0],{j,1,40}],{i,1,11}]
pnstoich//MatrixForm;

MatrixForm[Transpose[{pstoich.fluxnames,ppstoich.fluxnames,pnstoich.fluxnames}],TableHeadings->{totalpoolnames,{"net","in","out"}}];
MatrixForm[Transpose[{.01*Round[100*pstoich.fluxstst],ppstoich.fluxstst,pnstoich.fluxstst}],TableHeadings->{totalpoolnames,{"net","in","out"}}];

MatrixForm[Transpose[{pstoich.fluxnames,.01*Round[100*pstoich.fluxstst],ppstoich.fluxnames,ppstoich.fluxstst,pnstoich.fluxnames,pnstoich.fluxstst}],TableHeadings->{totalpoolnames,{"net","value","in","value","out","value"}}]

TeXForm[MatrixForm[Transpose[{.01*Round[pstoich.fluxstst],ppstoich.fluxstst,pnstoich.fluxstst}],TableHeadings->{totalpoolnames,{"net","in","out"}}]];

Calculate the response times of the pools

MatrixForm[tau ={0,0,0,0,0,0,0,0,0,Infinity,Infinity,0,0}];
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,1,9}]
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,12,13}]

MatrixForm[fluxnet ={0,0,0,0,0,0,0,0,0,".",".",0,0}];
Do[fluxnet[[i]]=ppstoich[[i]].fluxstst,{i,1,9}]
Do[fluxnet[[i]]=ppstoich[[i]].fluxstst,{i,12,13}]

MatrixForm[Transpose[{poolsize,fluxnet,tau}],TableHeadings->{totalpoolnames,{"size","net flux","tau"}}]
TeXForm[MatrixForm[Transpose[{poolsize,fluxnet,tau}],TableHeadings->{totalpoolnames,{"size","net flux","tau"}}]];

Define ratios and charges

Compute charge state of pools

rationames={"gly_phos","ade_phos","gly_redox","nad_redox","p_incorp","PFKactive"};
MatrixForm[rnumerator =(({{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}}))];MatrixForm[rdenominator =(({{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}}))];
ratiovars={0,0,0,0,0,0};
MatrixForm[Do[ratiovars[[i]]=(rnumerator[[i]].poolsize/rdenominator[[i]].poolsize),{i,1,6}]];
MatrixForm[Transpose[{rnumerator.totalpoolnames,rnumerator.poolsize,rdenominator.totalpoolnames,rdenominator.poolsize,ratiovars}],TableHeadings->{rationames,{"numerator","size(mM)","denominator","size(mM)","charge"}}]
BarChart[ratiovars];

Section III:	Simulation of the Non-Linear Mass Balance Equations

Set perturbations

tfinal = 100; 
tinterm = 15;
katp=1.5*katp;
vatp=katp*(atp - adp*phos/Katp);

Compute the solution

fluxes = {vhk, vpgi, vPFKR01, vPFKR02, vPFKR03, vPFKR10, vPFKR11, vPFKR12, vPFKR13, vPFKR20, vPFKR21, vPFKR22, vPFKR23, vPFKR30, vPFKR31, vPFKR32, vPFKR33, vPFKR40, vPFKR41, vPFKR42, vPFKR43, vPFKL, vPFKT1, vPFKT2, vPFKT3, vPFKT4, vtpi, vald, vgapdh, vpgk, vpglm, veno, vpk, vldh, vamp, vapk, vpyr, vlac, vatp, vnadh, vgluin, vampin, vh, vh2o};
massbal = stoich.Flatten[fluxes];
ic = varstst; 
solution = solveODE[massbal, ic, variables, 0, tfinal];

Section IV: 	Plots and Interpretation

Fluxes

Examine the beginning, interrim and the end flux map

tinterm = 10;
tflux = fluxes /. tsub[variables];
tset1 = Table[Flatten[{tflux[[i]] /. t -> 0, tflux[[i]] /. t -> tinterm, tflux[[i]] /. t -> tfinal} /. solution], {i, 1, Length[fluxes]}]; 
MatrixForm[tset1, TableHeadings -> {fluxnames, {0, tinterm, tfinal}}]
BarChart[Transpose[tset1][[1]], Transpose[tset1][[2]], Transpose[tset1][[3]],BarLabels -> fluxnames]

Selected flux phase planes

graph1vp=MASSPhasePlaneFlux[{vatp,vamp},{"vatp","vamp"},variables,solution,0,tfinal];
graph2vp=MASSPhasePlaneFlux[{vatp,vpyr},{"vatp","vpyr"},variables,solution,0,tfinal];
graph3vp=MASSPhasePlaneFlux[{vamp,vpyr},{"vamp","vpyr"},variables,solution,0,tfinal];
graph6vp=MASSPhasePlaneFlux[{vgapdh,vldh},{"vgapdh","vldh"},variables,solution,0,tfinal];
graph7vp=MASSPhasePlaneFlux[{vgapdh,vnadh},{"vgapdh","vnadh"},variables,solution,0,tfinal];

{graph1vp, graph2vp, graph3vp, graph6vp, graph7vp}

Select time responses of fluxes

Table[{MASSTimeProfileFlux[{fluxes[[i]]},{fluxnames[[i]]},(variables),solution,0,tfinal]},{i,1,Length[fluxes]}]

Concentrations

Get the start and end conentration states of the simulation

tvar = variables /. tsub[variables]; 
tset = Table[Flatten[{tvar[[i]] /. t -> 0, tvar[[i]] /. t -> tinterm, tvar[[i]] /. t -> tfinal} /. solution], {i, 1, Length[variables]}]; 
MatrixForm[tset, TableHeadings -> {variables, {0, tinterm, tfinal}}]
BarChart[Transpose[tset][[1]], Transpose[tset][[2]], Transpose[tset][[3]], BarLabels -> varnames]

Selected time responses of concentrations

Table[{MASSTimeProfileMet[{variables[[i]]},{varnames[[i]]},solution,0,tfinal]},{i,1,Length[variables]}]

Table[{MASSTimeProfileMet[{variables[[i]]},{varnames[[i]]},solution,0,40]},{i,1,Length[variables]}]

Pool fluxes

calculate the fluxes in and out of the pools

tpoolfluxin = ppstoich . tflux;
tset3 = Table[Flatten[{tpoolfluxin[[i]] /. t -> 0, tpoolfluxin[[i]] /. t -> tinterm, tpoolfluxin[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tpoolfluxin]}]; 
matin=MatrixForm[tset3, TableHeadings -> {totalpoolnames, {0, tinterm, tfinal}}];barin=BarChart[Transpose[tset3][[1]], Transpose[tset3][[2]], Transpose[tset3][[3]], BarLabels -> totalpoolnames];

tpoolfluxout = -pnstoich . tflux; 
tset4 = Table[Flatten[{tpoolfluxout[[i]] /. t -> 0, tpoolfluxout[[i]] /. t -> tinterm, tpoolfluxout[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tpoolfluxout]}]; 
matout=MatrixForm[tset4, TableHeadings -> {totalpoolnames, {0, tinterm, tfinal}}];
barout=BarChart[Transpose[tset4][[1]], Transpose[tset4][[2]], Transpose[tset4][[3]], BarLabels -> totalpoolnames];

{matin,matout}
{barin,barout}

flux phase planes for pools

Table[MASSPhasePlanePoolFlux[i,totalpoolnames[[i]],fluxes,variables,ppstoich,pnstoich,solution,0,tfinal],{i,1,10}]
Table[MASSPhasePlanePoolFlux[i,totalpoolnames[[i]],fluxes,variables,ppstoich,pnstoich,solution,0,tfinal],{i,12,13}]

Time responses of total pool fluxes (can break up into in/out if of interest)

Table[MASSTimeProfilePoolFlux[i,totalpoolnames[[i]],fluxes,variables,ppstoich,pnstoich,solution,0,tfinal],{i,1,13}]

Table[MASSTimeProfilePoolFlux[i,totalpoolnames[[i]],fluxes,variables,ppstoich,pnstoich,solution,0,40],{i,1,13}]

Pools

Calculate pools at start and end of simulations

tpool = pools . tvar; 
tset2 = Table[Flatten[{tpool[[i]] /. t -> 0, tpool[[i]] /. t -> tinterm, tpool[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tpool]}]; 
MatrixForm[tset2, TableHeadings -> {totalpoolnames, {0, tinterm, tfinal}}]
BarChart[Transpose[tset2][[1]], Transpose[tset2][[2]], Transpose[tset2][[3]], BarLabels -> totalpoolnames]
tpool[[1]] /. t -> 1 /. solution;

Selected pool phase planes

graph1pp=ParametricPlot[Evaluate[{tpool[[1]],tpool[[2]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{2.5,4},{1.3,1.5}},AspectRatio->1/2 ];
graph2pp=ParametricPlot[Evaluate[{tpool[[3]],tpool[[4]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{1,5},{0,1.5}},AspectRatio->1/2 ];
graph3pp=ParametricPlot[Evaluate[{tpool[[5]],tpool[[6]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{3.5,5},{0,.4}},AspectRatio->1/2 ];
graph4pp=ParametricPlot[Evaluate[{tpool[[8]],tpool[[9]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{2,4},{0,.15}},AspectRatio->1/2 ];
{graph1pp,graph2pp,graph3pp,graph4pp}

Time responses of pools

Table[MASSTimeProfilePool[{pools[[i]]},{totalpoolnames[[i]]},variables,solution,0,tfinal],{i,1,13}]

Table[MASSTimeProfilePool[{pools[[i]]},{totalpoolnames[[i]]},variables,solution,0,50],{i,1,13}]

Ratios

Compute the ratios at the beginning and the end

tratio = {0, 0, 0, 0, 0,0};
Do[tratio[[i]] = (rnumerator[[i]].tpool/rdenominator[[i]].tpool), {i, 1, 6}]
tset5 = Table[Flatten[{tratio[[i]] /. t -> 0, tratio[[i]] /. t -> tinterm, tratio[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tratio]}]; 
MatrixForm[tset5, TableHeadings -> {rationames, {0, tinterm, tfinal}}]
BarChart[Transpose[tset5][[1]], Transpose[tset5][[2]],Transpose[tset5][[3]], BarLabels -> rationames]

Compute phase planes for the ratios

graph1rp=MASSPhasePlaneRatios[{1,2},{rationames[[1]],rationames[[2]]},variables,pools,rnumerator,rdenominator,solution,0,tfinal];
graph2rp=MASSPhasePlaneRatios[{3,4},{rationames[[3]],rationames[[4]]},variables,pools,rnumerator,rdenominator,solution,0,tfinal];
graph3rp=MASSPhasePlaneRatios[{2,5},{rationames[[2]],rationames[[5]]},variables,pools,rnumerator,rdenominator,solution,0,tfinal];
graph4rp=MASSPhasePlaneRatios[{1,6},{rationames[[1]],rationames[[6]]},variables,pools,rnumerator,rdenominator,solution,0,tfinal];
{graph1rp,graph2rp,graph3rp,graph4rp}

Compute the time responses of the ratios

Table[MASSTimeProfileRatios[{i},rationames[[i]],variables,pools,rnumerator,rdenominator,solution,0,tfinal],{i,1,6}]

Table[MASSTimeProfileRatios[{i},rationames[[i]],variables,pools,rnumerator,rdenominator,solution,0,50],{i,1,6}]