'MASS' SIMULATOR OF GLYCOLYSIS

This Workbook sets up, analyses, and simulates and graphs the equations that describe the  glycolytic pathway coupled to the hemoglobin simulator.  It corresponds to chapter in Ch 13 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,nad,nadh,amp,adp,atp,phos,h,h2o,dpg23,o2,hb,hbo2,hbo22,hbo23,hbo24,dhb};
varnames={"glu","g6p","f6p","fbp","dhap","gap","pg13","pg3","pg2","pep", "pyr","lac","nad","nadh", "amp","adp","atp","phos","h", "h2o","dpg23","o2","hb","hbo2","hbo22","hbo23","hbo24","dhb"};
elemnames={"C","H","O","P","N", "S","NAD","Hb"};
MatrixForm[elem=(({{6, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3, 3, 0, 0, 10, 10, 10, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 3}, {12, 11, 11, 10, 5, 5, 4, 4, 4, 2, 3, 5, 0, 1, 13, 13, 13, 1, 1, 2, 3, 0, 0, 0, 0, 0, 0, 3}, {6, 9, 9, 12, 6, 6, 10, 7, 7, 6, 3, 3, 0, 0, 7, 10, 13, 4, 0, 1, 10, 2, 0, 2, 4, 6, 8, 10}, {0, 1, 1, 2, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 1, 2, 3, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 2}, {0, 0, 0, 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, 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, 1, 1, 1, 1}}))];
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,vamp,vapk,vpyr,vlac,vatp,vnadh,vgluin,vampin,vh,vh2o,vdpgm,vdpgase,vhbo1,vhbo2,vhbo3,vhbo4,vhbdpg,vo2};

fluxnames={"vhk","vpgi","vpfk","vtpi","vald","vgapdh","vpgk","vpglm","veno","vpk","vldh","vamp","vapk","vpyr","vlac","vatp","vnadh","vgluin","vampin","vh","vh2o","vdpgm","vdpgase","vhbo1","vhbo2","vhbo3","vhbo4","vhbdpg","vo2"};

knames={"khk","kpgi","kpfk","ktpi","kald","kgapdh","kpgk","kpglm","keno","kpk","kldh","kamp","kapk","kpyr","klac","katp","knadh","kgluin","kampin","kh","kh2o","kdpgm","kdpgase","khbo1","khbo2","khbo3","khbo4","khbdpg","ko2"};

Knames={"Khk","Kpgi","Kpfk","Ktpi","Kald","Kgapdh","Kpgk","Kpglm","Keno","Kpk","Kldh","Kamp","Kapk","Kpyr","Klac","Katp","Knadh","Kgluin","Kampin","Kh","Kh2o","Kdpgm","Kdpgase","Khbo1","Khbo2","Khbo3","Khbo4","Khbdpg","Ko2"};

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, 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, 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, -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, 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, 1, -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, -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, -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, -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, -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, 1, 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, 1, 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, -1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 0, 1, 0, 0, 0, -1, 0, 0, -1, 0, 0, -2, 0, 0, 1, 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, 1, 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, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {1, 0, 1, 0, 0, 1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 1, 1, 0, 0, -1, 0, 1, 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, 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, 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, -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, -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, 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, -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, -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, 0}, {0, 0, 0, 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}}))];
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=25, 
	there are two time invariant pools, Dim[Left Null]] =  m-r=28-25 = 3
	there are three time pathways, Dim[Null]] =  n-r=29-25 = 4

null=(({{0, 1, 0, 1}, {0, 1, 0, 0}, {0, 0, 1, 0}, {2, 1, 0, 1}})).NullSpace[stoich];
pathwaynames={"redox neural glycolysis","redox/plasma exch","amp in/out","rlsshunt"};
MatrixForm[null,TableHeadings->{pathwaynames,fluxnames}]
TeXForm[MatrixForm[null,TableHeadings->{pathwaynames,fluxnames}]];

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

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

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 three fluxes that are independent in the three pathways

gluin=1.12;
ststloadnadh=0.2*gluin;
ampin=.014;
rlsshunt= 0.196875;
MatrixForm[givenststfluxes=(({{gluin*(1-rlsshunt), ststloadnadh, ampin, gluin*rlsshunt}}))];
TeXForm[MatrixForm[givenststfluxes=(({{gluin*(1-rlsshunt), ststloadnadh, ampin, gluin*rlsshunt}}))]];
MatrixForm[null];
solution=Flatten[Transpose[givenststfluxes.null]];
MatrixForm[solution,TableHeadings->{fluxnames}]

Compute the steady state concentrations

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

infty=1000000;
pO2conv=2.864*10^-4;
Katp=infty;Kamp=infty;Kapk=1.65;Khk=850;Kpgi=0.41;Kpfk=310;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.;
Kdpgm=infty;Kdpgase=infty;
Khbo1=0.012/pO2conv;Khbo2=0.021/pO2conv;Khbo3=0.051/pO2conv;Khbo4=0.37/pO2conv;Khbdpg=1/4;Ko2=1;
Kvalues={Khk,Kpgi,Kpfk,Ktpi,Kald,Kgapdh,Kpgk,Kpglm,Keno,Kpk,Kldh,Kamp,Kapk,Kpyr,Klac,Katp,Knadh,Kgluin,Kampin,Kh,Kh2o,Kdpgm,Kdpgase,Khbo1,Khbo2,Khbo3,Khbo4,Khbdpg,Ko2};
MatrixForm[Kvalues,TableHeadings->{Knames}]
Length[Kvalues];
TeXForm[MatrixForm[Kvalues,TableHeadings->{Knames}]];

Define condentration variables: pool sizes and external concentrations

NADHtotal = 0.089;
pyrplasma=0.06;
lacplasma=1.0;
ampplasma=0.0001;
hplasma = 10.^(-4.2);
h2oplasma=1;
kh=infty/10; kh2o=infty/10;
o2plasma=((100-40)/2+40)*pO2conv;
hbtot=7.3;

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;

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

h2ostst=h2oplasma+solution[[21]]/kh2o;
dpg23stst=3.1;
o2stst=o2plasma;

eqn1 =Khbo1*hb0*o2stst-hb1;
eqn2 =Khbo2*hb1*o2stst-hb2;
eqn3 =Khbo3*hb2*o2stst-hb3;
eqn4 =Khbo4*hb3*o2stst-hb4;
eqn5 =Khbdpg*hb0*dpg23stst-hbd; (* 23dpg binding to hemoglobin independent of the o2 binding state *)
eqn6=hb0+hb1+hb2+hb3+hb4+hbd-hbtot; (* total mass balance on hemoglobin *)

ststsolution=NSolve[{eqn1==0,eqn2==0,eqn3==0,eqn4==0,eqn5==0,eqn6==0},{hb0,hb1,hb2,hb3,hb4,hbd}];
{hbstst,hbo2stst,hbo22stst,hbo23stst,hbo24stst,dhbstst}={hb0,hb1,hb2,hb3,hb4,hbd}/.ststsolution[[1]];
varstst={glustst,g6pstst,f6pstst,fbpstst,dhapstst,gapstst,pg13stst,pg3stst,pg2stst,pepstst,pyrstst,lacstst,nadstst,nadhstst,ampstst,adpstst,atpstst,phosstst,hstst,h2ostst,dpg23stst,o2stst,hbstst,hbo2stst,hbo22stst,hbo23stst,hbo24stst,dhbstst};

MatrixForm[varstst,TableHeadings->{varnames}]
TeXForm[MatrixForm[varstst,TableHeadings->{varnames}]];
BarChart[varstst,BarLabels->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.;

Gvalues={Ghk,Gpgi,Gpfk,Gtpi,Gald,Ggapdh,Gpgk,Gpglm,Geno,Gpk,Gldh,Gamp,Gapk,Gpyr,Glac,Gatp,Gnadh,Ggluin,Gampin,Gh,Gh2o,0,0,Khbo1,Khbo2,Khbo3,Khbo4,Khbdpg,Ko2};

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);
kamp=solution[[12]]/(ampstst-ampplasma/Kamp);
kapk =10^6;
kpyr =solution[[14]]/(pyrstst-pyrplasma/Kpyr);
klac =solution[[15]]/(lacstst-lacplasma/Klac);
katp=solution[[16]]/(atpstst-adpstst*phosstst/Katp);
knadh =solution[[17]]/(nadhstst-nadstst/Knadh);
kgluin=0.;
kampin=0;
kdpgase=solution[[22]]/(dpg23stst);
kdpgm=solution[[23]]/(pg13stst);
khbo1=506935.270702259;
khbo2=511077.050923776;
khbo3=509243.459567699;
khbo4=501595.340624411;
khbdpg=519612.560391792;
ko2=509725.707725914;

kvalues={khk,kpgi,kpfk,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kamp,kapk,kpyr,klac,katp,knadh,kgluin,kampin,kh,kh2o,kdpgm,kdpgase,khbo1,khbo2,khbo3,khbo4,khbdpg,ko2};

BarChart[Log[10,Sort[{khk,kpgi,kpfk,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kamp,kapk,kpyr,klac,katp,knadh,kh,kh2o,kdpgm,kdpgase,khbo1,khbo2,khbo3,khbo4,khbdpg,ko2},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);
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);

vdpgm=kdpgm*pg13;
vdpgase=kdpgase*dpg23;
vhbo1=khbo1*(hb*o2-hbo2/Khbo1);
vhbo2=khbo2*(hbo2*o2-hbo22/Khbo2);
vhbo3=khbo3*(hbo22*o2-hbo23/Khbo3);
vhbo4=khbo4*(hbo23*o2-hbo24/Khbo4);
vhbdpg=khbdpg*(hb*dpg23-dhb/Khbdpg);
vo2=ko2*(o2-o2plasma);

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

Numerical checks

Numerical check

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

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","O2 on Hb","tot Hb cap","dhb amount","tot 23dpg"};
MatrixForm[pools=(({{2, 3, 3, 4, 2, 2, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1}, {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, 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, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {2, 2, 2, 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, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 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, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1}, {0, 1, 1, 2, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 2}, {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, 1, 2, 3, 4, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 4, 4, 4, 4}, {0, 0, 0, 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, 0, 0, 0, 0, 0, 0, 1}}))];
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,29}],{i,1,15}]
ppstoich//MatrixForm;

pnstoich=pstoich;
Do[Do[pnstoich[[i,j]]=If[pstoich[[i,j]]<0,pstoich[[i,j]],0],{j,1,29}],{i,1,15}]
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,Infinity,0,Infinity}];
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,1,9}]
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,12,12}]
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,14,14}]

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,12}]
Do[fluxnet[[i]]=ppstoich[[i]].fluxstst,{i,14,14}]

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","hb_usage","dhb_ratio"};
MatrixForm[rnumerator =(({{1, 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, 1, 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, 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}, {0, 0, 1, 1, 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, 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, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}))];
ratiovars={0,0,0,0,0,0,0};
MatrixForm[Do[ratiovars[[i]]=(rnumerator[[i]].poolsize/rdenominator[[i]].poolsize),{i,1,7}]];
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 (oscillation)

Set perturbations

tfinal = 1; 

o2plasma = (30*Cos[120*Pi*t] + 70)*pO2conv;
vo2 = ko2*(o2 - o2plasma);
massbal = stoich.fluxes;

Compute the solution

ic = varstst; 
solution = solveODE[massbal, ic, variables, 0, tfinal,20000];

Section IV: 	Plots and Interpretation (oscillation)

Fluxes

Examine the beginning, interrim and the end flux map

tHigh = 0.018;
tLow = 0.027;
tflux = fluxes /. tsub[variables];
tset1 = Table[Flatten[{tflux[[i]] /. t -> 0, tflux[[i]] /. t -> tHigh, tflux[[i]] /. t -> tLow, tflux[[i]] /. t -> tfinal} /. solution], {i, 1, Length[fluxes]}]; 
MatrixForm[tset1, TableHeadings -> {fluxnames, {0, tHigh, tLow,tfinal}}]
BarChart[Transpose[tset1][[1]], Transpose[tset1][[2]], Transpose[tset1][[3]],Transpose[tset1][[4]],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}

graph8vp=MASSPhasePlaneFlux[{vdpgase,vhbo1},{"vdpgase","vhbo1"},variables,solution,0,tLow];
graph9vp=MASSPhasePlaneFlux[{vhbo4,vo2},{"vdpgase","vhbo1"},variables,solution,tLow,tfinal];

{graph8vp,graph9vp}

Select time responses of fluxes

graph1v=MASSTimeProfileFlux[{vatp,vamp,vpyr},{"vatp","vamp","vpyr"},variables,solution,0,tfinal];
graph2v= MASSTimeProfileFlux[{vgapdh,vldh,vnadh},{"vgapdh","vldh","vnadh"},variables,solution,0,tfinal];
{graph1v, graph2v}

graph4v=MASSTimeProfileFlux[{vdpgm,vdpgase},{"vdpgm","vdpgase"},variables,solution,0,0.1];
graph5v=MASSTimeProfileFlux[{vhbo2,vhbo3,vhbo4},{"vhbo2","vhbo3","vhbo4"},variables,solution,0,0.1];
graph8v=MASSTimeProfileFlux[{vhbdpg},{"vhbdpg"},variables,solution,0,0.1];
{graph4v,graph5v,graph8v}

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 -> tHigh, tvar[[i]] /. t -> tLow, tvar[[i]] /. t -> tfinal} /. solution], {i, 1, Length[variables]}]; 
MatrixForm[tset, TableHeadings -> {variables, {0, tHigh, tLow, tfinal}}]
BarChart[Transpose[tset][[1]], Transpose[tset][[2]], Transpose[tset][[3]], Transpose[tset][[4]], BarLabels -> varnames]

Selected time responses of concentrations

graph3c =MASSTimeProfileMet[{atp,adp,amp,glu,g6p,f6p},{"atp","adp","amp","glu","g6p","f6p"},solution,0,tfinal];
graph6c=MASSTimeProfileMet[{gap,pg13,pg2,pep},{"gap","pg13","pg2","pep"},solution,0,tfinal];
graph7c=MASSTimeProfileMet[{h},{"h"},solution,0,tfinal];
{ graph3c, graph6c, graph7c}

graph8c=MASSTimeProfileMet[{o2},{"o2"},solution,0,0.1];
graph9c=MASSTimeProfileMet[{dpg23,dhb,hbo24},{"dpg23","dhb","hbo24"},solution,0,0.1];
{graph8c,graph9c}

Pool fluxes

calculate the fluxes in and out of the pools

tpoolfluxin = ppstoich . tflux;
tset3 = Table[Flatten[{tpoolfluxin[[i]] /. t -> 0, tpoolfluxin[[i]] /. t -> tHigh,  tpoolfluxin[[i]] /. t -> tLow,tpoolfluxin[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tpoolfluxin]}]; 
matin=MatrixForm[tset3, TableHeadings -> {totalpoolnames, {0, tHigh,tLow, 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 -> tHigh,  tpoolfluxout[[i]] /. t -> tLow,tpoolfluxout[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tpoolfluxout]}]; 
matout=MatrixForm[tset4, TableHeadings -> {totalpoolnames, {0, tHigh,tLow, tfinal}}];
barout=BarChart[Transpose[tset4][[1]], Transpose[tset4][[2]], Transpose[tset4][[3]], Transpose[tset4][[4]],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,15}]

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,15}]

Table[MASSTimeProfilePoolFlux[i,totalpoolnames[[i]],fluxes,variables,ppstoich,pnstoich,solution,0,0.1],{i,12,15}]

Pools

Calculate pools at start and end of simulations

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

Selected pool phase planes

graph1pp=MASSPhasePlanePool[{pools[[1]],pools[[2]]},{totalpoolnames[[1]],totalpoolnames[[2]]},variables,solution,0,tfinal];
graph2pp=MASSPhasePlanePool[{pools[[3]],pools[[4]]},{totalpoolnames[[3]],totalpoolnames[[4]]},variables,solution,0,tfinal];
graph3pp=MASSPhasePlanePool[{pools[[5]],pools[[6]]},{totalpoolnames[[5]],totalpoolnames[[6]]},variables,solution,0,tfinal];
graph4pp=MASSPhasePlanePool[{pools[[7]],pools[[8]]},{totalpoolnames[[7]],totalpoolnames[[8]]},variables,solution,0,tfinal];
{graph1pp,graph2pp,graph3pp,graph4pp}

graph5pp=MASSPhasePlanePool[{pools[[12]],pools[[14]]},{totalpoolnames[[12]],totalpoolnames[[14]]},variables,solution,0,0.1];
{graph5pp}

Time responses of pools

graph1p=MASSTimeProfilePool[{pools[[1]],pools[[2]]},{totalpoolnames[[1]],totalpoolnames[[2]]},variables,solution,0,tfinal];
graph2p=MASSTimeProfilePool[{pools[[3]],pools[[4]]},{totalpoolnames[[3]],totalpoolnames[[4]]},variables,solution,0,tfinal];
graph3p=MASSTimeProfilePool[{pools[[5]],pools[[6]]},{totalpoolnames[[5]],totalpoolnames[[6]]},variables,solution,0,tfinal];
graph4p=MASSTimeProfilePool[{pools[[7]],pools[[8]]},{totalpoolnames[[7]],totalpoolnames[[8]]},variables,solution,0,tfinal];
{graph1p,graph2p,graph3p,graph4p}

graph5p=MASSTimeProfilePool[{pools[[12]]},{totalpoolnames[[12]]},variables,solution,0,tfinal];
graph6p=MASSTimeProfilePool[{pools[[14]]},{totalpoolnames[[14]]},variables,solution,0,tfinal];
{graph5p,graph6p}

Ratios

Compute the ratios at the beginning and the end

tratio = {0, 0, 0, 0, 0,0,0};
Do[tratio[[i]] = (rnumerator[[i]].tpool/rdenominator[[i]].tpool), {i, 1, 7}]
tset5 = Table[Flatten[{tratio[[i]] /. t -> 0, tratio[[i]] /. t -> tHigh,tratio[[i]] /. t -> tLow, 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]], Transpose[tset5][[4]], 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[{1,2},{rationames[[3]],rationames[[4]]},variables,pools,rnumerator,rdenominator,solution,0,tfinal];
{graph1rp, graph2rp}

graph3rp=MASSPhasePlaneRatios[{1,2},{rationames[[6]],rationames[[7]]},variables,pools,rnumerator,rdenominator,solution,0,tfinal];
{graph3rp}

Compute the time responses of the ratios

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