'MASS' SIMULATOR of GLYCOLYSIS and the PPP w/o REGULATION

This Workbook sets up, analyses, and simulates and graphs the equations that describe the glycolytic and pentose phosphate pathways.  It corresponds to chapter in Ch 11 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,gl6p,go6p,ru5p,x5p,r5p,s7p,e4p,nadp,nadph,gsh,gssg,co2};
varnames={"glu","g6p","f6p","fbp","dhap","gap","pg13","pg3","pg2","pep", "pyr","lac","nad","nadh", "amp","adp","atp","phos","h", "h2o","gl6p","go6p","ru5p","x5p","r5p","s7p","e4p","nadp","nadph","gsh","gssg","co2"};
elemnames={"C","H","O","P","N", "S","NAD","NADP"};
MatrixForm[elem=(({{6, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3, 3, 0, 0, 10, 10, 10, 0, 0, 0, 6, 6, 5, 5, 5, 7, 4, 0, 0, 10, 20, 1}, {12, 11, 11, 10, 5, 5, 4, 4, 4, 2, 3, 5, 0, 1, 13, 13, 13, 1, 1, 2, 9, 10, 9, 9, 9, 13, 7, 0, 1, 17, 32, 0}, {6, 9, 9, 12, 6, 6, 10, 7, 7, 6, 3, 3, 0, 0, 7, 10, 13, 4, 0, 1, 9, 10, 8, 8, 8, 10, 7, 0, 0, 6, 12, 2}, {0, 1, 1, 2, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 1, 2, 3, 1, 0, 0, 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, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 6, 0}, {0, 0, 0, 0, 0, 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, 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}}))];
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,vg6pdh,vpglase,vgl6pdh,vr5pe,vr5pi,vtki,vtkii,vtala,vgssgr,vgshr,vco2};

fluxnames={"vhk","vpgi","vpfk","vtpi","vald","vgapdh","vpgk","vpglm","veno","vpk","vldh","vamp","vapk","vpyr","vlac","vatp","vnadh","vgluin","vampin","vh","vh2o","vg6pdh","vpglase","vgl6pdh","vr5pe","vr5pi","vtki","vtkii","vtala","vgssgr","vgshr","vco2"};

knames={"khk","kpgi","kpfk","ktpi","kald","kgapdh","kpgk","kpglm","keno","kpk","kldh","kamp","kapk","kpyr","klac","katp","knadh","kgluin","kampin","kh","kh2o","kg6pdh","kpglase","kgl6pdh","kr5pe","kr5pi","ktki","ktkii","ktala","kgssgr","kgshr","kco2"};

Knames={"Khk","Kpgi","Kpfk","Ktpi","Kald","Kgapdh","Kpgk","Kpglm","Keno","Kpk","Kldh","Kamp","Kapk","Kpyr","Klac","Katp","Knadh","Kgluin","Kampin","Kh","Kh2o","Kg6pdh","Kpglase","Kgl6pdh","Kr5pe","Kr5pi","Ktki","Ktkii","Ktala","Kgssgr","Kgshr","Kco2"};

Define stoichiometric matrix; there are m=32 compounds and n=32 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, 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, 0, 0, 0, 0, 0, 0, 1, 1, 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, -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, 1, 1, -1, 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, 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, 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, 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, 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, -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, 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, -1, 1, 0, 0, 0, 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, -2, 0, 0, 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, 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, -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}, {1, 0, 1, 0, 0, 1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 1, 1, 0, 0, -1, 0, 1, 1, 0, 0, 0, 0, 0, 0, -1, 2, 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, 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, -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, -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, -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, 0, 0, 0, 0, 0}, {0, 0, 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, -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, -1, 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, 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, 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, -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, 1, 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=28, 
	there are three time invariant pools, Dim[Left Null]] =  m-r=32-28 = 4
	there are five pathways, Dim[Null]] =  n-r=32-28 = 4

nullorg=NullSpace[stoich];

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

(*Note; No phosphate balance with the r5p exchange added*)
leftnull=(({{0, 0, 0, 1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {1, 0, 0, 0}})).NullSpace[Transpose[stoich]];
poolnames={"nadh balance","phosphate blance","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;
ampin=.014;
MatrixForm[givenststfluxes=({{gluin-ststgshr/6, ststloadnadh, ampin, ststgshr/6}})  ];
TeXForm[MatrixForm[givenststfluxes]];

solution=Flatten[Transpose[givenststfluxes.null]];
nullstst=Append[null,solution];

MatrixForm[Transpose[nullstst],TableHeadings->{fluxnames,Append[Flatten[givenststfluxes],"stst"]}]

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;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;
Kg6pdh= 1000.; Kpglase=1000.; Kgl6pdh =1000.;Kr5pe=3.0; Kr5pi=2.57;Ktki=1.2;Ktkii=10.3;Ktala=1.05;Kgssgr=100.;Kgshr=2;Kco2=1; 
(*With new concentrations used in the integrated model Ktki needs to be changed to prevent negative PERC*)
Ktki = 3;
Kvalues={Khk,Kpgi,Kpfk,Ktpi,Kald,Kgapdh,Kpgk,Kpglm,Keno,Kpk,Kldh,Kamp,Kapk,Kpyr,Klac,Katp,Knadh,Kgluin,Kampin,Kh,Kh2o,Kg6pdh,Kpglase,Kgl6pdh,Kr5pe,Kr5pi,Ktki,Ktkii,Ktala,Kgssgr,Kgshr,Kco2};
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;

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;

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

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;

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

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

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.;
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);
Gco2 = 1.;

Gvalues={Ghk,Gpgi,Gpfk,Gtpi,Gald,Ggapdh,Gpgk,Gpglm,Geno,Gpk,Gldh,Gamp,Gapk,Gpyr,Glac,Gatp,Gnadh,Ggluin,Gampin,Gh,Gh2o,Gg6pdh,Gpglase,Ggl6pdh,Gr5pe,Gr5pi,Gtki,Gtkii,Gtala,Ggssgr,Ggshr,Gco2};
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);
klac =solution[[15]]/(lacstst-lacplasma);
katp=solution[[16]]/(atpstst-adpstst*phosstst/Katp);
knadh =ststloadnadh/(nadhstst-nadstst/Knadh);
kgluin=0.;
kampin=0.;
kg6pdh=solution[[22]]/(g6pstst*nadpstst-gl6pstst*nadphstst/Kg6pdh);
kpglase=solution[[23]]/(gl6pstst-go6pstst/Kpglase);
kgl6pdh=solution[[24]]/(go6pstst*nadpstst-ru5pstst*nadphstst/Kgl6pdh);
kr5pe=solution[[25]]/(ru5pstst-x5pstst/Kr5pe);
kr5pi=solution[[26]]/(ru5pstst-r5pstst/Kr5pi);
ktki=solution[[27]]/(x5pstst*r5pstst-s7pstst*gapstst/Ktki);
ktkii=solution[[28]]/(x5pstst*e4pstst-f6pstst*gapstst/Ktkii);
ktala=solution[[29]]/(s7pstst*gapstst-e4pstst*f6pstst/Ktala);
kgssgr=solution[[30]]/(nadphstst*gssgstst-nadpstst*gshstst^2/Kgssgr);
kgshr=solution[[31]]/(gshstst^2-gssgstst/Kgshr);

kvalues={khk,kpgi,kpfk,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kamp,kapk,kpyr,klac,katp,knadh,kgluin,kampin,kh,kh2o,kg6pdh,kpglase,kgl6pdh,kr5pe,kr5pi,ktki,ktkii,ktala,kgssgr,kgshr,kco2};

BarChart[Log[10,Sort[{khk,kpgi,kpfk,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kamp,kapk,kpyr,klac,katp,knadh,kh,kh2o,kg6pdh,kpglase,kgl6pdh,kr5pe,kr5pi,ktki,ktkii,ktala,kgssgr,kgshr,kco2},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/Kh2o);

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

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);


fluxstst= Flatten[fluxes/.sublist[variables,varstst]];

Print["The steady state fluxes are:  "]
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.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

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, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, 8/3, 8/3, 8/3, 10/3, 7/3, 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, 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, 2, 1, 0, 0, 0, 0, 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, 2, 2, 5/3, 5/3, 5/3, 7/3, 4/3, 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, 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, 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, 1, 1, 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, 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}, {6, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 4, 4, 4, 8, 2, 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, 1, 1, 2, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 1, 1, 1, 1, 1, 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, 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, 1, 1, 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}]
Dimensions[pstoich,2]

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,32}],{i,1,17}]
ppstoich//MatrixForm;

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

MatrixForm[fluxnet ={0,0,0,0,0,0,0,0,0,0,0,0,0,".",".",".","."}];
Do[fluxnet[[i]]=ppstoich[[i]].fluxstst,{i,1,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

totalpoolnames={"GP+","GP-","AP+","AP-","GR+","GR-","nadh","glutathione","nadph","P+","P-","GRP+","GRP-","Ptot","nadh-tot","glutathione-tot","nadph-tot"};
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={0,0,0,0,0,0,0,0};

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

gluin = 1.*gluin; 
ampin = 1.*ampin; 

pyrplasma = 1.*0.06; 
lacplasma = 1.*1.0; 
ampplasma = 1.*0.0001; 
vamp = kamp*(amp - ampplasma/Kamp); 
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, vamp, vapk, vpyr, vlac, vatp, vnadh, vgluin, vampin, vh, vh2o, vg6pdh, vpglase, vgl6pdh, vr5pe, vr5pi, vtki, vtkii, vtala, vgssgr, vgshr, vco2};

Compute the solution

massbal = stoich . fluxes; 
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,vamp}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{1.5,3.5},{0,.09}},AspectRatio->1/1];
graph2vp=ParametricPlot[Evaluate[{vatp,vpyr}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{1.,3.5},{0.1,.3}},AspectRatio->1/1];
graph3vp=ParametricPlot[Evaluate[{vamp,vpyr}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{0,.1},{0.1,.3}},AspectRatio->1/1];
graph4vp=ParametricPlot[Evaluate[{vtkii,vgl6pdh}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{0.,.1},{.2,.22}},AspectRatio->1/1];
graph5vp=ParametricPlot[Evaluate[{vgssgr,vgl6pdh}]/.tsub[variables]/.solution, {t,0,tfinal}, PlotRange->{{.41,.43},{.2,.22}},AspectRatio->1/1];
graph6vp=ParametricPlot[Evaluate[{vgapdh,vldh}]/.tsub[variables]/.solution,{t,0,tfinal},PlotRange->{{1.6,2.2},{1.5,2.}},AspectRatio->1/1 ];
graph7vp=ParametricPlot[Evaluate[{vgapdh,vnadh}]/.tsub[variables]/.solution,{t,0,tfinal},PlotRange->{{1.6,2.2},{0.0,.3}},AspectRatio->1/1 ];

{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]}]

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 time responses of concentrations

Table[{MASSTimeProfileMet[{variables[[i]]},{varnames[[i]]},solution,0,tfinal]},{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,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}]

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 -> {totalpool.names, {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

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

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

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