'MASS' SIMULATOR of GLYCOLYSIS, PPP, and Nucleotide Metabolism w/o REGULATION
This Workbook sets up, analyses, and simulates and graphs the equations that describe these networks.  It corresponds to Ch 12 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,fbp,dhap,gap,pg13,pg3,pg2,pep,pyr,lac,gl6p,go6p,ru5p,x5p,r5p,s7p,e4p,ade,ado,imp,ino,hyp,r1p,prpp,nad,nadh,amp,adp,atp,nadp,nadph,gsh,gssg,phos,h,h2o,co2,nh3};
varnames={"glu","g6p","f6p","fbp","dhap","gap","pg13","pg3","pg2","pep", "pyr","lac","gl6p","go6p","ru5p","x5p","r5p","s7p","e4p", "ade","ado","imp","ino","hyp","r1p","prpp","nad","nadh", "amp","adp","atp","nadp","nadph","gsh","gssg","phos","h", "h2o","co2","nh3"};
elemnames={"C","H","O","P","N", "S","NAD","NADP"};
MatrixForm[elem=(({{6, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3, 3, 6, 6, 5, 5, 5, 7, 4, 5, 10, 10, 10, 5, 5, 5, 0, 0, 10, 10, 10, 0, 0, 10, 20, 0, 0, 0, 1, 0}, {12, 11, 11, 10, 5, 5, 4, 4, 4, 2, 3, 5, 9, 10, 9, 9, 9, 13, 7, 5, 13, 12, 12, 4, 9, 8, 0, 1, 13, 13, 13, 0, 1, 17, 32, 1, 1, 2, 0, 3}, {6, 9, 9, 12, 6, 6, 10, 7, 7, 6, 3, 3, 9, 10, 8, 8, 8, 10, 7, 0, 4, 8, 5, 1, 8, 14, 0, 0, 7, 10, 13, 0, 0, 6, 12, 4, 0, 1, 2, 0}, {0, 1, 1, 2, 1, 1, 2, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 3, 0, 0, 1, 2, 3, 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, 5, 5, 4, 4, 4, 0, 0, 0, 0, 5, 5, 5, 0, 0, 3, 6, 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, 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, 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}}))];
MatrixForm[elem,TableHeadings->{elemnames,varnames}]
TeXForm[MatrixForm[elem,TableHeadings->{elemnames,varnames}]];

Define the reaction names

fluxes={vhk,vpgi,vpfk,vtpi,vald,vgapdh,vpgk,vpglm,veno,vpk,vldh,vg6pdh,vpglase,vgl6pdh,vr5pe,vr5pi,vtki,vtkii,vtala,vak,vampase,vada,vampda,vimpase,vpnpase,vprm,vprppsyn,vadprt,vapk,vatp,vnadh,vgssgr,vgshr,vgluin,vpyr,vlac,vade,vado,vino,vhyp,vphos,vh,vh2o,vco2, vnh3};

fluxnames={"vhk","vpgi","vpfk","vtpi","vald","vgapdh","vpgk","vpglm","veno","vpk","vldh","vg6pdh","vpglase","vgl6pdh","vr5pe","vr5pi","vtki","vtkii","vtala","vak","vampase","vada","vampda","vimpase","vpnpase","vprm","vprppsyn","vadprt","vapk","vatp","vnadh","vgssgr","vgshr","vgluin","vpyr","vlac","vade","vado","vino","vhyp","vphos", "vh","vh2o","vco2", "vnh3"};

knames={"khk","kpgi","kpfk","ktpi","kald","kgapdh","kpgk","kpglm","keno","kpk","kldh","kg6pdh","kpglase","kgl6pdh","kr5pe","kr5pi","ktki","ktkii","ktala","kak","kampase","kada","kampda","kimpase","kpnpase","kprm","kprppsyn","kadprt","kapk","katp","knadh","kgssgr","kgshr","kgluin","kpyr","klac","kade","kado","kino","khyp", "kphos","kh","kh2o", "kco2","knh3"};

Knames={"Khk","Kpgi","Kpfk","Ktpi","Kald","Kgapdh","Kpgk","Kpglm","Keno","Kpk","Kldh","Kg6pdh","Kpglase","Kgl6pdh","Kr5pe","Kr5pi","Ktki","Ktkii","Ktala","Kak","Kampase","Kada","Kampda","Kimpase","Kpnpase","Kprm","Kprppsyn","Kadprt","Kapk","Katp","Knadh","Kgssgr","Kgshr","Kgluin","Kpyr","Klac","Kade","Kado","Kino","Khyp","Kphos", "Kh","Kh2o", "Kco2","Knh3"};

Define stoichiometric matrix; there are m=41 compounds and n=47 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, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, -1, 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, 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, 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, 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, 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, -1, 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, 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, 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, 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, -1, 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, 0, 0, -1, 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, 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, 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, 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, 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, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 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, 0, 0, 0, -1, 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, 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, 0, 1, -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, 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, 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, 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, -1, 0, 0, 0, 0, 1, 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, 1, 0, 0, 0, 0, -1, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, -1, 0, 0, 0, 1, 1, 1, 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, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-1, 0, -1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 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}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -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, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 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, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -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, 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, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, -1, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0}, {1, 0, 1, 0, 0, 1, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, -1, 2, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, 0, 0, 0, -1, 0, -1, 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, 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, -1, 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, -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=38, 
	there are three time invariant pools, Dim[Left Null]] =  m-r=40-37 = 3
	there are eight pathways, Dim[Null]] =  n-r=45-37 = 8

MatrixForm[nullorg=NullSpace[stoich]];

leftMultMat = (({{0, 0, 0, 1/3, 0, 2/3, 2, 0}, {0, 0, 0, -1/3, 0, 1/3, 0, 0}, {0, 2/2, -1/2, 5/2-1/3, 0, 1/3, 1, 0}, {0, 1/3, 0, 2/3, 0, 0, -3, 0}, {1/6, 1/3, -1/6, 2/3, 0, 0, -3, 0}, {1/6, 1/3, -1/6, 2/3, 0, 0, -3, 1}, {1/6, 0, -1/6, 0, 1, 0, -2, 0}, {0, 0, 0, 0, 0, 0, -1, 0}}));
null=leftMultMat.nullorg;
pathwaynames={"Glycolysis - glc to 2lac/2h","Pyruvate Reduction - pyr/2h to lac","PPP cycle - glc/3h2o to 3co2/13h","Adenosine synthesis requiring 3 atp - glc/ade to ado/co2/4h/nh3","Inosine synthesis requiring 3 atp - glc/ade to ino/co2/4h/nh3 through ampase/ada","Inosine synthesis requiring 3 atp - glc/ade to ino/co2/4h/nh3 through ampda/impase","Salvage - ade/h2o to hyp/nh3","Futile Cycle - ATP consumption through ak/ampase"};
MatrixForm[nullorg,TableHeadings->{pathwaynames,fluxnames}];
MatrixForm[null,TableHeadings->{pathwaynames,fluxnames}]
TeXForm[MatrixForm[null,TableHeadings->{pathwaynames,fluxnames}]];

MatrixForm[Transpose[null],TableHeadings->{fluxnames,pathwaynames}];

(*Note; No phosphate balance with the r5p exchange added*)
leftnull=(({{0, 0, 1}, {0, 1, 0}, {1, 0, 0}})).NullSpace[Transpose[stoich]];
poolnames={"nadh balance","nadph balance","glutathione balance"};
MatrixForm[leftnull,TableHeadings->{poolnames,varnames}]
TeXForm[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}]

Set the flux steady state

fluxomic data: specify fluxes through four fluxes that are independent in the four pathways

gluin=1.12;
ststloadnadh=0.2*gluin;
ststgshr=0.42;
adaFlux = .01;
adoFlux = -.01;
akFlux = .12;
hypEXstst=0.097;
apkFlux = 0;

(*The following method is used to find independent sets of reactions that can be used to set the steady-state flux distribution with preference to reactions for which higher fidelity data exists*)
matrix = {};
remainingFluxes = {};
nullorgNonneg=pStoich[NullSpace[stoich]]-nStoich[NullSpace[stoich]];
Do[
Plus@@nullorgNonneg[[All,i]]
If[Plus@@nullorgNonneg[[All,i]] > 0,
remainingFluxes = Append[remainingFluxes,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];

MatrixForm[nullorgNonneg,TableHeadings->{{"p1","p2","p3","p4","p5","p6","p7"},fluxnames}];

(*Add vgluin*)
fluxesUsed = {};
stoichCopy1 = stoich;
stoichCopy1[[All,Position[fluxnames,"vgluin"][[1]]]] = 0;
fluxesUsed = Append[fluxesUsed,"vgluin"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[1]]][[1]]]]][[1]];
MatrixRank[matrix];

(*Find remaining fluxes*)
remainingFluxes1 = {};
nullorgNonneg1=(pStoich[NullSpace[stoichCopy1]])-nStoich[NullSpace[stoichCopy1]];
Do[
Plus@@nullorgNonneg1[[All,i]]
If[Plus@@nullorgNonneg1[[All,i]] > 0,
remainingFluxes1 = Append[remainingFluxes1,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes1];

(*Add vgshr*)
nullorgNonneg1=pStoich[NullSpace[stoichCopy1]]-nStoich[NullSpace[stoichCopy1]];
stoichCopy2 = stoichCopy1;
newPos = Position[fluxnames,"vgshr"][[1]];
If[(Plus@@nullorgNonneg1[[All,newPos]])[[1]] > 0,
stoichCopy2[[All,newPos]] = 0;
fluxesUsed = Append[fluxesUsed,"vgshr"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[2]]][[1]]]][[1]]];
MatrixRank[matrix];
];

(*Find remaining fluxes*)
remainingFluxes2 = {};
nullorgNonneg2=(pStoich[NullSpace[stoichCopy2]])-nStoich[NullSpace[stoichCopy2]];
Do[
Plus@@nullorgNonneg2[[All,i]]
If[Plus@@nullorgNonneg2[[All,i]] > 0,
remainingFluxes2 = Append[remainingFluxes2,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes2];

(*Add vnadh*)
nullorgNonneg2=pStoich[NullSpace[stoichCopy2]]-nStoich[NullSpace[stoichCopy2]];
stoichCopy3 = stoichCopy2;
newPos = Position[fluxnames,"vnadh"][[1]];
If[(Plus@@nullorgNonneg2[[All,newPos]])[[1]] > 0,
stoichCopy3[[All,newPos]] = 0;
fluxesUsed = Append[fluxesUsed,"vnadh"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[3]]][[1]]]][[1]]];
MatrixRank[matrix];
];

(*Find remaining fluxes*)
remainingFluxes3 = {};
nullorgNonneg3=(pStoich[NullSpace[stoichCopy3]])-nStoich[NullSpace[stoichCopy3]];
Do[
Plus@@nullorgNonneg3[[All,i]]
If[Plus@@nullorgNonneg3[[All,i]] > 0,
remainingFluxes3 = Append[remainingFluxes3,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes3];

(*Add vapk*)
nullorgNonneg3=pStoich[NullSpace[stoichCopy3]]-nStoich[NullSpace[stoichCopy3]];
MatrixForm[nullorgNonneg3,TableHeadings->{{"p1","p2","p3","p4","p5","p6","p7"},fluxnames}];
stoichCopy4 = stoichCopy3;
newPos = Position[fluxnames,"vapk"][[1]];
If[(Plus@@nullorgNonneg3[[All,newPos]])[[1]] > 0,
stoichCopy4[[All,newPos]] = 0;
fluxesUsed = Append[fluxesUsed,"vapk"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[4]]][[1]]]][[1]]];
MatrixRank[matrix];
];

(*Find remaining fluxes*)
remainingFluxes4 = {};
nullorgNonneg4=(pStoich[NullSpace[stoichCopy4]])-nStoich[NullSpace[stoichCopy4]];
Do[
Plus@@nullorgNonneg4[[All,i]]
If[Plus@@nullorgNonneg4[[All,i]] > 0,
remainingFluxes4 = Append[remainingFluxes4,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes4];

(*Add vado*)
nullorgNonneg4=pStoich[NullSpace[stoichCopy4]]-nStoich[NullSpace[stoichCopy4]];
MatrixForm[nullorgNonneg4,TableHeadings->{{"p1","p2","p3","p4","p5","p6","p7"},fluxnames}];
stoichCopy5 = stoichCopy4;
newPos = Position[fluxnames,"vado"][[1]];
If[(Plus@@nullorgNonneg4[[All,newPos]])[[1]] > 0,
stoichCopy5[[All,newPos]] = 0;
fluxesUsed = Append[fluxesUsed,"vado"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[5]]][[1]]]][[1]]];
MatrixRank[matrix];
];

(*Find remaining fluxes*)
remainingFluxes5 = {};
nullorgNonneg5=(pStoich[NullSpace[stoichCopy5]])-nStoich[NullSpace[stoichCopy5]];
Do[
Plus@@nullorgNonneg5[[All,i]]
If[Plus@@nullorgNonneg5[[All,i]] > 0,
remainingFluxes5 = Append[remainingFluxes5,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes5];

(*Add ahyp*)
nullorgNonneg5=pStoich[NullSpace[stoichCopy5]]-nStoich[NullSpace[stoichCopy5]];
MatrixForm[nullorgNonneg5,TableHeadings->{{"p1","p2","p3","p4","p5","p6","p7"},fluxnames}];
stoichCopy6 = stoichCopy5;
newPos = Position[fluxnames,"vhyp"][[1]];
If[(Plus@@nullorgNonneg5[[All,newPos]])[[1]] > 0,
stoichCopy6[[All,newPos]] = 0;
fluxesUsed = Append[fluxesUsed,"vhyp"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[6]]][[1]]]][[1]]];
MatrixRank[matrix];
];

(*Find remaining fluxes*)
remainingFluxes6 = {};
nullorgNonneg6=(pStoich[NullSpace[stoichCopy6]])-nStoich[NullSpace[stoichCopy6]];
Do[
Plus@@nullorgNonneg6[[All,i]]
If[Plus@@nullorgNonneg6[[All,i]] > 0,
remainingFluxes6 = Append[remainingFluxes6,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes6];

(*Add ak*)
nullorgNonneg6=pStoich[NullSpace[stoichCopy6]]-nStoich[NullSpace[stoichCopy6]];
MatrixForm[nullorgNonneg6,TableHeadings->{{"p1","p2","p3","p4","p5","p6","p7"},fluxnames}];
stoichCopy7 = stoichCopy6;
newPos = Position[fluxnames,"vak"][[1]];
If[(Plus@@nullorgNonneg6[[All,newPos]])[[1]] > 0,
stoichCopy7[[All,newPos]] = 0;
fluxesUsed = Append[fluxesUsed,"vak"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[7]]][[1]]]][[1]]];
MatrixRank[matrix];
];

(*Find remaining fluxes*)
remainingFluxes7 = {};
nullorgNonneg7=(pStoich[NullSpace[stoichCopy7]])-nStoich[NullSpace[stoichCopy7]];
Do[
Plus@@nullorgNonneg7[[All,i]]
If[Plus@@nullorgNonneg7[[All,i]] > 0,
remainingFluxes7 = Append[remainingFluxes7,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes7];

(*Add ada*)
nullorgNonneg7=pStoich[NullSpace[stoichCopy7]]-nStoich[NullSpace[stoichCopy7]];
MatrixForm[nullorgNonneg7,TableHeadings->{{"p1","p2","p3","p4","p5","p6","p7"},fluxnames}];
stoichCopy8 = stoichCopy7;
newPos = Position[fluxnames,"vada"][[1]];
If[(Plus@@nullorgNonneg7[[All,newPos]])[[1]] > 0,
stoichCopy8[[All,newPos]] = 0;
fluxesUsed = Append[fluxesUsed,"vada"];
matrix = Append[matrix,Transpose[nullorg][[Position[fluxnames,fluxesUsed[[8]]][[1]]]][[1]]];
MatrixRank[matrix];
];

(*Find remaining fluxes*)
remainingFluxes8 = {};
nullorgNonneg8=(pStoich[NullSpace[stoichCopy8]])-nStoich[NullSpace[stoichCopy8]];
Do[
Plus@@nullorgNonneg8[[All,i]]
If[Plus@@nullorgNonneg8[[All,i]] > 0,
remainingFluxes8 = Append[remainingFluxes8,fluxnames[[i]]];
];,
{i,1,Dimensions[stoich,2][[2]]}
];
MatrixForm[remainingFluxes8];

MatrixForm[stoichCopy7];
MatrixForm[fluxesUsed];

MatrixForm[matrix,TableHeadings->{fluxesUsed,{"p1","p2","p3","p4","p5","p6","p7","p8"}}];
fluxesstst = (({{gluin}, {ststgshr}, {ststloadnadh}, {apkFlux}, {adoFlux}, {hypEXstst}, {akFlux}, {adaFlux}}));
MatrixForm[Transpose[givenststfluxes= Transpose[Inverse[matrix].fluxesstst]], TableHeadings->{pathwaynames}];

TeXForm[MatrixForm[givenststfluxes]];
MatrixForm[nullorg,TableHeadings->{{"p1","p2","p3","p4","p5","p6","p7","p8"},fluxnames}];
solution=Flatten[Transpose[givenststfluxes.nullorg]];
solution = Chop[solution,10^-10];
MatrixForm[solution,TableHeadings->{fluxnames}]
MatrixForm[stoich.solution];

Compute the steady state concentrations

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

infty=1000000;
Katp=infty;Kamp=infty;Kapk=1.65;Khk=850;Kpgi=0.41;Kpfk=310;Kald=0.082;Kgapdh=0.0179;Kpgk=1800;
Ktpi=1/17.5;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;
Kg6pdh= 1000.; Kpglase=1000.; Kgl6pdh =1000.;Kr5pe=3.0; Kr5pi=2.57;Ktkii=10.3;Ktala=1.05;Kgssgr=100.;Kgshr=2;Kco2=1; 

Kak=infty;
Kampase=infty; (*1 in Neema's model*)
Kampda=infty; (*.1 in Neema's model*)
Kimpase=infty; (*1000 in Neema's model*)
Kada=infty; (*1 in Neema's model*)

Katpgen=infty;
Kprppsyn=infty; (*605 in Neema's model*)
Kadprt=infty; (*100000 in Neema's model, almost same*)
Kado=1; (*.01 in Neema's model*)
Kade=1;
Kino=1;
Khyp=1;
Kamp=1;
Kh=1;
Kh2o=1;
Kphos=1;
Knh3=1;
Kapk = 1.65;

(*The following equilibrium constants were changed to make positive PERCs*)
Kpnpase=1/0.09; (*.11 in Neema's model, .09 in Joshi's model*)
Kprm=1/13.30; (*13.3 in Neema's model*)
Ktki=3; (*1.2 in Neema's model*)

Kvalues={Khk,Kpgi,Kpfk,Ktpi,Kald,Kgapdh,Kpgk,Kpglm,Keno,Kpk,Kldh,Kg6pdh,Kpglase,Kgl6pdh,Kr5pe,Kr5pi,Ktki,Ktkii,Ktala,Kak,Kampase,Kada,Kampda,Kimpase,Kpnpase,Kprm,Kprppsyn,Kadprt,Kapk,Katp,Knadh,Kgssgr,Kgshr,Kgluin,Kpyr,Klac,Kade,Kado,Kino,Khyp,Kphos,Kh,Kh2o, Kco2, Knh3};
MatrixForm[Kvalues,TableHeadings->{Knames}]
Length[Kvalues];
TeXForm[MatrixForm[Kvalues,TableHeadings->{Knames}]];

Define concentration variables: pool sizes and external concentrations

NADHtotal = 0.089;
NADPHtotal = 0.066;
Glutathionetotal = 3.32;
pyrplasma=0.06;
lacplasma=1.0;
ampplasma=0.0001;
hplasma = 10.^(-4.2);
h2oplasma=1;
co2plasma=1;
kh=infty/10; 
kh2o=infty/10;
kco2=infty/10;

Define steady state concentration values

atpstst=1.6;adpstst=0.29;

ampstst=Kapk*adpstst*adpstst/atpstst; (*This could be changed to allow apk flux without an infinite reaction rate*)

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;

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

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

gl6pstst=0.001754242723;go6pstst=0.037475258;ru5pstst=0.0049367903;x5pstst=0.014784196;r5pstst=0.012668936;s7pstst=0.023987984;e4pstst=0.0050750696;

nadpstst=0.0002;nadphstst= NADPHtotal - nadpstst;
gshstst=3.2;gssgstst= Glutathionetotal - gshstst;
co2stst=co2plasma+solution[[44]]/kco2 ;

phosstst=2.5; 
r5pstst=0.00494;
adestst = 0.001;
adostst = 0.0012;
impstst = 0.01;
inostst = 0.001;
hypstst = 0.002;
r1pstst = 0.06;
prppstst = 0.005;

(*Other literature and estimated concentrations*)
nh3stst = .091; (*91uM from Harris et. al. 1989 American Physiological Society*) 

varstst={glustst,g6pstst,f6pstst,fbpstst,dhapstst,gapstst,pg13stst,pg3stst,pg2stst,pepstst,pyrstst,lacstst,gl6pstst,go6pstst,ru5pstst,x5pstst,r5pstst,s7pstst,e4pstst,adestst,adostst,impstst,inostst,hypstst,r1pstst,prppstst,nadstst,nadhstst,ampstst,adpstst,atpstst,nadpstst,nadphstst,gshstst,gssgstst,phosstst,hstst,h2ostst,co2stst, nh3stst};

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

Compute the PERCs

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.;
Gg6pdh=gl6pstst*nadphstst/(g6pstst*nadpstst);
Gpglase=go6pstst/gl6pstst;
Ggl6pdh=ru5pstst*nadphstst/(go6pstst*nadpstst);
Gr5pe= x5pstst/ru5pstst;
Gr5pi=r5pstst/ru5pstst;
Gtki=s7pstst*gapstst/(x5pstst*r5pstst);
Gtkii=f6pstst*gapstst/(x5pstst*e4pstst);
Gtala=e4pstst*f6pstst/(s7pstst*gapstst);
Ggssgr= gshstst*gshstst*nadpstst/(nadphstst*gssgstst);
Ggshr=gssgstst/(gshstst*gshstst);
Gak=adpstst*ampstst/(adostst*atpstst);
Gampase=adostst/ampstst;
Gampda=impstst/ampstst;
Gimpase=inostst/impstst;
Gada=inostst/adostst;
Gpnpase=hypstst*r1pstst/(inostst);
Gprm=r5pstst/r1pstst;
Gprppsyn=prppstst*ampstst/(r5pstst*atpstst);
Gadprt=ampstst/(prppstst*adestst);
Gado=1.;
Gade=1.;
Gino=1.;
Ghyp=1.;
Gamp=1.;
Gphos = 1.;
Gh=1.;
Gh2o=1.;
Gco2 = 1.;
Gnh3 = 1.;

Gvalues={Ghk,Gpgi,Gpfk,Gtpi,Gald,Ggapdh,Gpgk,Gpglm,Geno,Gpk,Gldh,Gg6pdh,Gpglase,Ggl6pdh,Gr5pe,Gr5pi,Gtki,Gtkii,Gtala,Gak,Gampase,Gada,Gampda,Gimpase,Gpnpase,Gprm,Gprppsyn,Gadprt,Gapk,Gatp,Gnadh,Ggssgr,Ggshr,Ggluin,Gpyr,Glac,Gade,Gado,Gino,Ghyp,Gphos,Gh,Gh2o,Gco2, Gnh3};
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

khk=solution[[1]]/(glustst*atpstst-g6pstst*adpstst/Khk);
kpgi =solution[[2]]/(g6pstst-f6pstst/Kpgi);
kpfk=solution[[3]]/(f6pstst*atpstst-fbpstst*adpstst/Kpfk);
ktpi=solution[[4]]/(dhapstst-gapstst/Ktpi);
kald=solution[[5]]/(fbpstst-gapstst*dhapstst/Kald);
kgapdh=solution[[6]]/(phosstst*gapstst*nadstst-nadhstst*pg13stst/Kgapdh);
kpgk=solution[[7]]/(pg13stst*adpstst-pg3stst*atpstst/Kpgk);
kpglm=solution[[8]]/(pg3stst-pg2stst/Kpglm);
keno=solution[[9]]/(pg2stst-pepstst/Keno);
kpk=solution[[10]]/(pepstst*adpstst-pyrstst*atpstst/Kpk);
kldh =solution[[11]]/(pyrstst*nadhstst-lacstst*nadstst/Kldh);
kg6pdh=solution[[12]]/(g6pstst*nadpstst-gl6pstst*nadphstst/Kg6pdh);
kpglase=solution[[13]]/(gl6pstst-go6pstst/Kpglase);
kgl6pdh=solution[[14]]/(go6pstst*nadpstst-ru5pstst*nadphstst/Kgl6pdh);
kr5pe=solution[[15]]/(ru5pstst-x5pstst/Kr5pe);
kr5pi=solution[[16]]/(ru5pstst-r5pstst/Kr5pi);
ktki=solution[[17]]/(x5pstst*r5pstst-s7pstst*gapstst/Ktki);
ktkii=solution[[18]]/(x5pstst*e4pstst-f6pstst*gapstst/Ktkii);
ktala=solution[[19]]/(s7pstst*gapstst-e4pstst*f6pstst/Ktala);
kak = solution[[20]]/( adostst*atpstst - adpstst*ampstst/Kak);
kampase = solution[[21]]/(ampstst-adostst/Kampase);
kada = solution[[22]]/(adostst - inostst/Kada);
kampda = solution[[23]]/(ampstst - impstst/Kampda);
kimpase = solution[[24]]/(impstst - inostst/Kimpase);
kpnpase = solution[[25]]/(inostst - hypstst*r1pstst/Kpnpase);
kprm = solution[[26]]/( r5pstst- r1pstst/Kprm);
kprppsyn = solution[[27]]/(r5pstst*atpstst - prppstst*ampstst/Kprppsyn); 
kadprt = solution[[28]]/(prppstst*adestst - ampstst*phosstst*phosstst/Kadprt);
kapk = infty;
katp=solution[[30]]/(atpstst-adpstst*phosstst/Katp);
knadh =ststloadnadh/(nadhstst-nadstst/Knadh);
kgssgr=solution[[32]]/(nadphstst*gssgstst-nadpstst*gshstst^2/Kgssgr);
kgshr=solution[[33]]/(gshstst^2-gssgstst/Kgshr);
kgluin=0.;
kpyr =solution[[35]]/(pyrstst-pyrplasma);
klac =solution[[36]]/(lacstst-lacplasma);
kado = infty/10;
kade = infty/10;
khyp = infty/10;
kino = infty/10;
kh = solution[[42]]/(hstst - hplasma/Kh);
kh2o = solution[[43]]/(h2ostst -h2oplasma /Kh2o);
kphos = infty/10;
knh3 = infty/10;

adoplasma = Kado*(adostst - solution[[38]]/kado);
adeplasma = Kade*(adestst - solution[[37]]/kade);
inoplasma = Kino*(inostst - solution[[39]]/kino);
hypplasma = Khyp*(hypstst - solution[[40]]/khyp);
phosplasma = Kphos*(phosstst - solution[[41]]/kphos);
nh3plasma = Knh3*(nh3stst - solution[[45]]/knh3);

kvalues={khk,kpgi,kpfk,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kg6pdh,kpglase,kgl6pdh,kr5pe,kr5pi,ktki,ktkii,ktala,kak,kampase,kada,kampda,kimpase,kpnpase,kprm,kprppsyn,kadprt,kapk,katp,knadh,kgssgr,kgshr,kgluin,kpyr,klac,kade,kado,kino,khyp, kphos,kh,kh2o ,kco2,knh3};


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

The rate laws

Summary of reaction properties

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

Define the rate laws

vhk=khk*(atp*glu-adp*g6p/Khk);
vpgi=kpgi*(g6p-f6p/Kpgi);
vpfk=kpfk*(f6p*atp-fbp*adp/Kpfk);
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);
vgapdh=kgapdh*(phos*gap*nad-nadh*pg13/Kgapdh);
vldh=kldh*(pyr*nadh-lac*nad/Kldh);
vg6pdh=kg6pdh*(g6p*nadp-gl6p*nadph/Kg6pdh);
vpglase=kpglase*(gl6p-go6p/Kpglase);
vgl6pdh=kgl6pdh*(go6p*nadp-ru5p*nadph/Kgl6pdh);
vr5pe=kr5pe*(ru5p-x5p/Kr5pe);
vr5pi=kr5pi*(ru5p-r5p/Kr5pi);
vtki=ktki*(x5p*r5p-s7p*gap/Ktki);
vtkii=ktkii*(x5p*e4p-f6p*gap/Ktkii);
vtala=ktala*(s7p*gap-e4p*f6p/Ktala);
vgssgr=kgssgr*(nadph*gssg-nadp*gsh^2/Kgssgr);
vgshr=kgshr*(gsh^2-gssg/Kgshr);
vco2=kco2*(co2-co2plasma);
vak= kak*(ado*atp - adp*amp/Kak); 
vampase= kampase*(amp - ado/Kampase); 
vada= kada*(ado - ino/Kada); 
vampda= kampda*(amp - imp/Kampda); 
vimpase= kimpase*(imp-ino/Kimpase); 
vpnpase= kpnpase*(ino-hyp*r1p/Kpnpase); 
vprm= kprm*(r5p-r1p/Kprm); 
vprppsyn= kprppsyn*(r5p*atp-prpp*amp*h/Kprppsyn); 
vadprt= kadprt*(prpp*ade - amp*phos*phos/Kadprt); 
vapk = kapk*(adp*adp - amp*atp/Kapk);
vatp= katp*(atp-adp*phos/Katp);
vnadh = knadh*(nadh-nad/Knadh);
vgluin=gluin;
vh = kh*(h-hplasma/Kh);
vh2o = kh2o*(h2o-h2oplasma/Kh2o);
vpyr =kpyr*(pyr-pyrplasma/Kpyr);
vlac =klac*(lac-lacplasma/Klac);
vado= kado*(ado-adoplasma/Kado);
vade= kade*(ade - adeplasma/Kade);
vhyp= khyp*(hyp - hypplasma/Khyp);
vino = kino*(ino - inoplasma/Kino);
vphos = kphos*(phos - phosplasma/Kphos);
vnh3 = knh3*(nh3 - nh3plasma/Knh3);

fluxes={vhk,vpgi,vpfk,vtpi,vald,vgapdh,vpgk,vpglm,veno,vpk,vldh,vg6pdh,vpglase,vgl6pdh,vr5pe,vr5pi,vtki,vtkii,vtala,vak,vampase,vada,vampda,vimpase,vpnpase,vprm,vprppsyn,vadprt,vapk,vatp,vnadh,vgssgr,vgshr,vgluin,vpyr,vlac,vade,vado,vino,vhyp,vphos,vh,vh2o,vco2, vnh3};

MatrixForm[fluxes,TableHeadings->{fluxnames}];
fluxstst= Flatten[fluxes/.sublist[variables,varstst]];

Print["The steady state fluxes are:  "]
MatrixForm[0.001*Round[1000*fluxstst],TableHeadings->{fluxnames}]
BarChart[fluxstst, ChartLabels -> fluxnames]

Numerical checks

Numerical check

MatrixForm[0.01*Round[100*(solution-fluxstst)], TableHeadings->{fluxnames}]

Calculate the flux balances

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

Section II: 	Defining pools and ratios

Define the pools

Compute Pools: hard and soft pools
	use the left null space for 'hard pools'
	'soft pools' defined based on biochemical intuition
	compute pool sizes based on the steady state concentrations

(*Need to determine pools by maximizing ATP for GP and maybe NADH for GR*)

totalpoolnames={"GP+","GP-","AP+","AP-","GR+","GR-","nadh","glutathione","nadph","P+","P-","GRP+","GRP-","Ptot","nadh-tot","glutathione-tot","nadph-tot"};

MatrixForm[pools=(({{2, 3, 3, 4, 2, 2, 2, 1, 1, 1, 0, 0, 3, 3, 8/3, 8/3, 8/3, 10/3, 7/3, 0, 8/3, 8/3, 8/3, 0, 8/3, 8/3, 0, 0, 8/3, 11/3, 14/3, 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, 0, 0, 0, 0, 0, 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, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {2, 2, 2, 2, 1, 1, 0, 0, 0, 0, 0, 1, 2, 2, 5/3, 5/3, 5/3, 7/3, 4/3, 0, 5/3, 5/3, 5/3, 0, 5/3, 5/3, 0, 0, 5/3, 5/3, 5/3, 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, 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, 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, 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, 1, 1, 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, 1, 2, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {6, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 4, 4, 4, 8, 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, 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, 1, 1, 2, 1, 1, 2, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 3, 0, 0, 1, 2, 3, 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, 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, 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, 1, 1, 0, 0, 0, 0, 0, 0, 0}}))];

Dimensions[pools,2];
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

MatrixForm[pstoich.fluxnames]

ppstoich=pstoich;
Do[Do[ppstoich[[i,j]]=If[pstoich[[i,j]]>0,pstoich[[i,j]],0],{j,1,Dimensions[stoich,2][[2]]}],{i,1,Dimensions[totalpoolnames][[1]]}]
ppstoich//MatrixForm;

pnstoich=pstoich;
Do[Do[pnstoich[[i,j]]=If[pstoich[[i,j]]<0,pstoich[[i,j]],0],{j,1,Dimensions[stoich,2][[2]]}],{i,1,Dimensions[totalpoolnames][[1]]}]
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 = Table[0,{i,1,Dimensions[totalpoolnames][[1]]}]];
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,1,Dimensions[totalpoolnames][[1]]}]

MatrixForm[fluxnet =Table[0,{i,1,Dimensions[totalpoolnames][[1]]}]];
Do[fluxnet[[i]]=ppstoich[[i]].fluxstst,{i,1,Dimensions[totalpoolnames][[1]]}]

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={"charge_gly_phos","charge_ade_phos","charge_gly_redox","charge_nad_redox","ratio_p_recycle","charge_gly_redoxP","charge_gsh","charge_nadp_redox"};

MatrixForm[rnumerator = (({{1, 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, 1, 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, 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, 1, 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}}))];
MatrixForm[rdenominator = (({{1, 1, 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, 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, 1, 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, 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, 1, 1, 0, 0, 0, 0}}))];

ratiovars=Table[0,{i,1,Length[rationames]}];

MatrixForm[Do[ratiovars[[i]]=(rnumerator[[i]].poolsize/rdenominator[[i]].poolsize),{i,1,Length[rationames]}]];
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

tinterm =0.01;
tfinal = 20; 

vgluin = 1.12;
vgluin = 1.*vgluin; 

pyrplasma = 1.*0.06; 
lacplasma = 1.*1.0; 
vpyr = kpyr*(pyr - pyrplasma/Kpyr); 
vlac = klac*(lac - lacplasma/Klac); 

katp 	= 1.5*katp; 
knadh 	= 1.0*knadh; 
kgshr	= 1.*kgshr;
vnadh 	= knadh*(nadh - nad/Knadh); 
vatp = katp*(atp - adp*phos/Katp);
vgshr 	= kgshr*(gsh^2 - gssg/Kgshr);
fluxes = {vhk, vpgi, vpfk, vtpi, vald, vgapdh, vpgk, vpglm, veno, vpk, vldh, vg6pdh, vpglase, vgl6pdh, vr5pe, vr5pi, vtki, vtkii, vtala, vak, vampase, vada, vampda, vimpase, vpnpase, vprm, vprppsyn, vadprt, vapk, vatp, vnadh, vgssgr, vgshr, vgluin, vpyr, vlac, vade, vado, vino, vhyp, vphos, vh, vh2o, vco2, vnh3};

Compute the solution

MatrixForm[massbal = stoich.fluxes, TableHeadings -> {varnames}];
ic = varstst;
solution = solveODE[massbal, ic, variables, 0, tfinal];

Section IV: 	Plots and Interpretation

Fluxes

Examine the beginning, interim and the end flux map

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]

Tiled flux phase planes (careful takes a long time to compute)

tile3=plotset[massbal, variables, 0, tinterm,ic,fluxes,fluxnames,{0,1,10,20}]
tile4=plotset[massbal, variables, tinterm, tfinal,ic,fluxes,fluxnames,{20,30,40,50}]

Selected flux phase planes

graph1vp=ParametricPlot[Evaluate[{vatp,vak}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{1.5,3.7},{0.04,.15}},AspectRatio->1/1,AxesLabel-> {"vatp","vak"}];
graph2vp=ParametricPlot[Evaluate[{vatp,vampase}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{1.5,3.7},{0,.5}},AspectRatio->1/1,AxesLabel-> {"vatp","vampase"}];
graph3vp=ParametricPlot[Evaluate[{vatp,(vampase- vak)}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{1.5,3.7},{-.2,.4}},AspectRatio->1/1,AxesLabel-> {"vatp","vampase-vak"}];
graph4vp=ParametricPlot[Evaluate[{vak,vampase}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{0,.2},{0,.4}},AspectRatio->1/1,AxesLabel-> {"vak","vampase"}];
graph5vp = ParametricPlot[{Evaluate[{(vak+vhk+vpfk+vprppsyn+vatp),(vpgk+vpk+vapk)}]/.tsub[variables]/.solution,Evaluate[{(vhk+vpfk+vatp),(vpgk+vpk+vapk)}]/.tsub[variables]/.solution}, {t,0,tfinal}, PlotRange->{{3,6},{3,6}},AspectRatio->1/1,AxesLabel-> {"sum atp consumption","sum atp production"},PlotStyle->{Thickness[.004],Dashed}];
graph6vp = ParametricPlot[Evaluate[{(vhk+vpfk+vatp),(vpgk+vpk+vapk)}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{0,6},{0,6}},AspectRatio->1/1,AxesLabel-> {"sum atp consumption","sum atp production"}];
{graph1vp, graph2vp, graph3vp, graph4vp,graph5vp,graph6vp}

Select time responses of fluxes

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

(*Plot solution from this model and solution from the glycolysis.ppp model*)
labels = {"vatp - gly+ppp+salvage","vatp - gly+ppp"};
colors = {Hue[0.6],Hue[1]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph1v=Plot[{vatp/.tsub[variables]/.solution,vatp/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{1,4}, Epilog->Inset[legend, {tfinal-5,1.4}] ];

labels = {"vnadh - gly+ppp+salvage","vnadh - gly+ppp"};
colors = {Hue[0.6],Hue[1]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph2v=Plot[{vnadh/.tsub[variables]/.solution,vnadh/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{.1,.4}, Epilog->Inset[legend, {tfinal-5,.14}] ];

labels = {"vgshr - gly+ppp+salvage","vgshr - gly+ppp"};
colors = {Hue[0.6],Hue[1]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph3v=Plot[{vgshr/.tsub[variables]/.solution,vgshr/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{.3,.5}, Epilog->Inset[legend, {tfinal-5,.35}] ];

labels = {"vapk - gly+ppp+salvage","vapk - gly+ppp"};
colors = {Hue[0.6],Hue[1]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph4v=Plot[{vapk/.tsub[variables]/.solution,vapk/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{-.2,.4}, Epilog->Inset[legend, {tfinal-5,.1}] ];

labels = {"vtkii - gly+ppp+salvage","vtkii - gly+ppp"};
colors = {Hue[0.6],Hue[1]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph5v=Plot[{vtkii/.tsub[variables]/.solution,vtkii/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{0.04,.15}, Epilog->Inset[legend, {tfinal-5,.05}] ];

labels = {"vpyr - gly+ppp+salvage","vpyr - gly+ppp"};
colors = {Hue[0.6],Hue[1]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph6v=Plot[{vpyr/.tsub[variables]/.solution,vpyr/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{0.1,.4}, Epilog->Inset[legend, {tfinal-5,.15}] ];
{graph1v,graph2v,graph3v,graph4v,graph5v,graph6v}

Concentrations

Get the start and end concentration 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]

Tiled concentration phase planes (careful takes a long time to compute)

tile1=plotset[massbal,variables,0,tinterm,ic,variables,varnames,{0,1,10,20}];
tile2=plotset[massbal,variables, tinterm, tfinal, ic, variables, varnames,{20,30,40,50}];
{tile1,tile2}

Selected concentration phase planes

Selected time responses of concentrations

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

labels = {"atp - glyc + PPP + salvage", "atp - glyc + PPP"};
colors = {Hue[1`],Hue[0.7`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph1c=Plot[{atp[t]/.solution,atp[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",PlotStyle->{Thickness[.005],Dashed},PlotRange->{0.8,1.6}]

labels = {"pep - glyc + PPP + salvage", "pep - glyc + PPP"};
colors = {Hue[1`],Hue[0.7`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph2c=Plot[{pep[t]/.solution,pep[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, Automatic],PlotStyle->colors,PlotRange->{0,.04}]

labels = {"r5p - glyc + PPP + salvage", "r5p - glyc + PPP"};
colors = {Hue[1`],Hue[0.7`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph3c=Plot[{r5p[t]/.solution,r5p[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, Automatic],PlotStyle->colors,PlotRange->{0,.01}]

labels = {"nadh - glyc + PPP + salvage", "nadh - glyc + PPP"};
colors = {Hue[1`],Hue[0.7`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph4c=Plot[{nadh[t]/.solution,nadh[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, Automatic],PlotStyle->colors,PlotRange->{0.02,.04}]

labels = {"nadp - glyc + PPP + salvage", "nadp - glyc + PPP"};
colors = {Hue[1`],Hue[0.7`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph5c=Plot[{nadp[t]/.solution,nadp[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, Automatic],PlotStyle->colors,PlotRange->{0,.0005}]

labels = {"amp - glyc + PPP + salvage", "amp - glyc + PPP"};
colors = {Hue[1`],Hue[0.7`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph5c=Plot[{amp[t]/.solution,amp[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",PlotStyle->{Thickness[.005],Dashed},PlotRange->{0,.6}]

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

graph1pp=ParametricPlot[Evaluate[{tpoolfluxin[[1]],tpoolfluxout[[1]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[1]],PlotRange->{{3,5},{3,5}},AspectRatio->1/2 ];
graph2pp=ParametricPlot[Evaluate[{tpoolfluxin[[2]],tpoolfluxout[[2]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[2]],PlotRange->{{1,3},{1,3}},AspectRatio->1/2 ];
graph3pp=ParametricPlot[Evaluate[{tpoolfluxin[[3]],tpoolfluxout[[3]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[3]],PlotRange->{{3,5},{3,5}},AspectRatio->1/2 ];
graph4pp=ParametricPlot[Evaluate[{tpoolfluxin[[4]],tpoolfluxout[[4]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[4]],PlotRange->{{3,5},{3,5}},AspectRatio->1/2 ];
graph5pp=ParametricPlot[Evaluate[{tpoolfluxin[[5]],tpoolfluxout[[5]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[5]],PlotRange->{{3,5},{3,5}},AspectRatio->1/2 ];
graph6pp=ParametricPlot[Evaluate[{tpoolfluxin[[6]],tpoolfluxout[[6]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[6]],PlotRange->{{1,3},{1,3}},AspectRatio->1/2 ];
graph7pp=ParametricPlot[Evaluate[{tpoolfluxin[[7]],tpoolfluxout[[7]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[7]],PlotRange->{{1,3},{1,3}},AspectRatio->1/2 ];
graph8pp=ParametricPlot[Evaluate[{tpoolfluxin[[8]],tpoolfluxout[[8]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[8]],PlotRange->{{0.7,2},{0.7,2}},AspectRatio->1/2 ];
graph9pp=ParametricPlot[Evaluate[{tpoolfluxin[[9]],tpoolfluxout[[9]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[9]],PlotRange->{{.4,.6},{.4,.6}},AspectRatio->1/2 ];
graph10pp=ParametricPlot[Evaluate[{tpoolfluxin[[10]],tpoolfluxout[[10]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[10]],PlotRange->{{2.3,2.6},{2.3,2.6}},AspectRatio->1/2 ];
graph11pp=ParametricPlot[Evaluate[{tpoolfluxin[[11]],tpoolfluxout[[11]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[11]],PlotRange->{{2,2.2},{2,2.2}},AspectRatio->1/2 ];
graph12pp=ParametricPlot[Evaluate[{tpoolfluxin[[12]],tpoolfluxout[[12]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[12]],PlotRange->{{6.5,7},{6.5,7}},AspectRatio->1/2 ];
graph13pp=ParametricPlot[Evaluate[{tpoolfluxin[[13]],tpoolfluxout[[13]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[13]],PlotRange->{{2,2.3},{2,2.3}},AspectRatio->1/2 ];

{graph1pp,graph2pp,graph3pp,graph4pp,graph5pp,graph6pp,graph7pp,graph8pp,graph9pp,graph10pp,graph11pp,graph12pp,graph13pp}

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

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

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;

Tiled pool phase planes (careful takes a long time to compute)

tile5=plotset[massbal,variables, 0, tinterm, ic, pools.variables, totalpoolnames,{0,1,10,20}];
tile6=plotset[massbal,variables, tinterm, tfinal, ic, pools.variables, totalpoolnames,{20,30,40,50}];
{tile5,tile6}

Selected pool phase planes

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

Time responses of pools

MASSTimeProfilePool[pools,totalpoolnames,(variables),solution,0,tfinal]

Ratios

Compute the ratios at the beginning and the end

tratio = {0, 0, 0, 0, 0, 0, 0, 0};
Do[tratio[[i]] = (rnumerator[[i]].tpool/rdenominator[[i]].tpool), {i, 1, 8}]
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, tfinal}}]
BarChart[Transpose[tset5][[1]], Transpose[tset5][[2]],Transpose[tset5][[3]], BarLabels -> rationames]

Compute phase planes for the ratios

graph1rp=ParametricPlot[Evaluate[{tratio[[1]],tratio[[2]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{0.6,0.8},{0.5,.9}},AspectRatio->1/2 ];
graph2rp=ParametricPlot[Evaluate[{tratio[[3]],tratio[[4]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{0.9,1.0},{0.2,.4}},AspectRatio->1/2 ];
{graph1rp, graph2rp}

Compute the time responses of the ratios

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