'MASS' SIMULATOR OF GLYCOLYSIS

This Workbook sets up, analyses, and simulates and graphs the equations that describe the  glycolytic pathway.  It corresponds to chapter in Ch 10 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};
varnames={"glu","g6p","f6p","fbp","dhap","gap","pg13","pg3","pg2","pep", "pyr","lac","nad","nadh", "amp","adp","atp","phos","h", "h2o"};
elemnames={"C","H","O","P","N", "S","NAD"};
MatrixForm[elem=(({{6, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3, 3, 0, 0, 10, 10, 10, 0, 0, 0}, {12, 11, 11, 10, 5, 5, 4, 4, 4, 2, 3, 5, 0, 1, 13, 13, 13, 1, 1, 2}, {6, 9, 9, 12, 6, 6, 10, 7, 7, 6, 3, 3, 0, 0, 7, 10, 13, 4, 0, 1}, {0, 1, 1, 2, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 1, 2, 3, 1, 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, 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}}))];
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};

fluxnames={"vhk","vpgi","vpfk","vtpi","vald","vgapdh","vpgk","vpglm","veno","vpk","vldh","vamp","vapk","vpyr","vlac","vatp","vnadh","vgluin","vampin","vh","vh2o"};

knames={"khk","kpgi","kpfk","ktpi","kald","kgapdh","kpgk","kpglm","keno","kpk","kldh","kamp","kapk","kpyr","klac","katp","knadh","kgluin","kampin","kh","kh2o"};

Knames={"Khk","Kpgi","Kpfk","Ktpi","Kald","Kgapdh","Kpgk","Kpglm","Keno","Kpk","Kldh","Kamp","Kapk","Kpyr","Klac","Katp","Knadh","Kgluin","Kampin","Kh","Kh2o"};

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

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

Check elemental balancing--exchange reactions not balanced

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

The null spaces

The null spaces:  determining the pathways and the pools
	the rank of S is r=18, 
	there are two time invariant pools, Dim[Left Null]] =  m-r=20-18 = 2
	there are three time pathways, Dim[Null]] =  n-r=21-18 = 3

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

leftnull=NullSpace[Transpose[stoich]];
poolnames={"phosphate balance","nadh balance"};
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 three fluxes that are independent in the three pathways

gluin=1.12;
ststloadnadh=0.2*gluin;
ampin=.014;

MatrixForm[givenststfluxes=(({{gluin, ststloadnadh, ampin}}))];
TeXForm[MatrixForm[givenststfluxes=(({{gluin, ststloadnadh, ampin}}))]];
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;
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.;
Kvalues={Khk,Kpgi,Kpfk,Ktpi,Kald,Kgapdh,Kpgk,Kpglm,Keno,Kpk,Kldh,Kamp,Kapk,Kpyr,Klac,Katp,Knadh,Kgluin,Kampin,Kh,Kh2o};
MatrixForm[Kvalues,TableHeadings->{Knames}]
Length[Kvalues];
TeXForm[MatrixForm[Kvalues,TableHeadings->{Knames}]];

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

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;

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

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

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=gluin/(glustst*atpstst-g6pstst*adpstst/Khk);
kpgi =gluin/(g6pstst-f6pstst/Kpgi);
kpfk=gluin/(f6pstst*atpstst-fbpstst*adpstst/Kpfk);
ktpi=gluin/(dhapstst-gapstst/Ktpi);
kald=gluin/(fbpstst-gapstst*dhapstst/Kald);
kgapdh=2*gluin/(phosstst*gapstst*nadstst-nadhstst*pg13stst/Kgapdh);
kpgk=2*gluin/(pg13stst*adpstst-pg3stst*atpstst/Kpgk);
kpglm=2*gluin/(pg3stst-pg2stst/Kpglm);
keno=2*gluin/(pg2stst-pepstst/Keno);
kpk=2*gluin/(pepstst*adpstst-pyrstst*atpstst/Kpk);
kldh =(2*gluin-ststloadnadh)/(pyrstst*nadhstst-lacstst*nadstst/Kldh);
kamp=ampin/(ampstst-ampplasma/Kamp);
kapk =10^6;
klac =(2*gluin-ststloadnadh)/(lacstst-lacplasma);
kpyr =ststloadnadh/(pyrstst-pyrplasma);
katp=2*gluin/(atpstst-adpstst*phosstst/Katp);
knadh =ststloadnadh/(nadhstst-nadstst/Knadh);
kgluin=0.;
kampin=0;

kvalues={khk,kpgi,kpfk,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kamp,kapk,kpyr,klac,katp,knadh,kgluin,kampin,kh,kh2o};

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

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"};
MatrixForm[pools=(({{2, 3, 3, 4, 2, 2, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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, 2, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0}, {2, 2, 2, 2, 1, 1, 0, 0, 0, 0, 0, 1, 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, 1, 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}, {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 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}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0}}))];
Dimensions[pools,2];
MatrixForm[poolsize = pools.varstst];
MatrixForm[Transpose[Append[Transpose[pools],poolsize]],TableHeadings->{totalpoolnames,Append[varnames,"size"]}]
TeXForm[MatrixForm[pools,TableHeadings->{totalpoolnames,varnames}]];
BarChart[poolsize,BarLabels->totalpoolnames]

Pool fluxes and their turnover times

Calculate the reaction that move the pools

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

Compute the fluxes in and out of pools

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

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

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

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

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

Calculate the response times of the pools

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

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

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

MatrixForm[Sqrt[poolsize],TableHeadings-> {totalpoolnames,{"sqrt"}}];
Sqrt[poolsize[[10]]-poolsize[[9]]-poolsize[[8]]];
Sqrt[poolsize[[11]]-poolsize[[7]]];

Define ratios and charges

Compute charge state of pools

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

MatrixForm[Sqrt[rdenominator.poolsize],TableHeadings->{rationames}];
Sqrt[poolsize[[10]]];

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

estimate time scales to set tfinal and any intermediate time scales

chart1=BarChart[Log[10,{khk,kpgi,kpfk,ktpi,kald,kgapdh,kpgk,kpglm,keno,kpk,kldh,kamp,kapk,kpyr,klac,katp,knadh,.001,.001,kh,kh2o}],BarLabels->knames, AspectRatio->1/5,BarGroupSpacing->.5]

Set perturbations

tfinal = 50; 
tinterm = 20.;
vgluin = 1.0*vgluin; 
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; 
vnadh = knadh*(nadh - nad/Knadh); 
vatp = katp*(atp - adp*(phos/Katp));

Compute the solution

fluxes = {vhk, vpgi, vpfk, vtpi, vald, vgapdh, vpgk, vpglm, veno, vpk, vldh, vamp, vapk, vpyr, vlac, vatp, vnadh, vgluin, vampin, vh, vh2o};
massbal = stoich.fluxes;
ic = varstst; 
solution = solveODE[massbal, ic, variables, 0, tfinal];

Section IV: 	Plots and Interpretation

Fluxes

Examine the beginning, interrim and the end flux map

tinterm = 20;
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]

Table[{i,fluxes[[i]]},{i,1,Length[fluxes]}]//MatrixForm;

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}];
{tile3,tile4}

Selected flux phase planes

{graph1vp,graph2vp,graph3vp,graph4vp,graph5vp,graph6vp}={
ParametricPlot[Evaluate[{vatp,vamp}]/.tsub[variables]/.solution,
{t,0,tfinal}, PlotRange->{{1.5,3.5},{0,.08}},AspectRatio->1/1],
ParametricPlot[Evaluate[{vatp,vpyr}]/.tsub[variables]/.solution,
{t,0,tfinal}, PlotRange->{{1.,3.5},{0.15,.3}},AspectRatio->1/1],
ParametricPlot[{vamp,vpyr}/.tsub[variables]/.solution,
{t,0,tfinal}, PlotRange->{{0,.1},{0.15,.3}},AspectRatio->1/1],ParametricPlot[Evaluate[{vgapdh,vldh}]/.tsub[variables]/.solution,{t,0,tfinal},PlotRange->{{1.6,2.2},{1.4,2.0}},AspectRatio->1/1],ParametricPlot[Evaluate[{vldh,vnadh}]/.tsub[variables]/.solution,{t,0,tfinal},PlotRange->{{1.4,2.2},{0.1,.3}},AspectRatio->1/1],ParametricPlot[Evaluate[{vhk,vpk}]/.tsub[variables]/.solution,{t,0,tfinal},PlotRange->{{0.75,1.25},{1.5,2.5}},AspectRatio->1/1]}

Select time responses of fluxes

labels = {fluxnames[[1]],fluxnames[[3]],fluxnames[[7]],fluxnames[[10]]};
colors = {Hue[0.5],Hue[0.6],Hue[0.8],Hue[0.9],Hue[1.0]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph1v=Plot[{vhk/.tsub[variables]/.solution,vpfk/.tsub[variables]/.solution,vpgk/.tsub[variables]/.solution,vpk/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{0,4}, Epilog->Inset[legend, {tfinal-5,2.5}] ];

labels = {fluxnames[[16]]};
colors = {Hue[0.5],Hue[0.6],Hue[0.8],Hue[0.9],Hue[1.0]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph2v=Plot[{vatp/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{0,4}, Epilog->Inset[legend, {tfinal-5,3.}] ];

labels = {fluxnames[[6]],fluxnames[[11]],fluxnames[[17]]};
colors = {Hue[0.5],Hue[0.6],Hue[0.8],Hue[0.9],Hue[1.0]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph3v=Plot[{vgapdh/.tsub[variables]/.solution,vldh/.tsub[variables]/.solution,vnadh/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{0,3}, Epilog->Inset[legend, {tfinal-5,1.}] ];

labels = {fluxnames[[12]],fluxnames[[14]]};
colors = {Hue[0.5],Hue[0.6],Hue[0.8],Hue[0.9],Hue[1.0]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph4v=Plot[{vamp/.tsub[variables]/.solution,vpyr/.tsub[variables]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",PlotStyle->colors,PlotRange->{0,0.3}, Epilog->Inset[legend, {tfinal-5,3.}] ];

{graph1v,graph2v,graph3v,graph4v}

Concentrations

Get the start and end conentration states of the simulation

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

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

labels = {"g6p","pep"};
colors = {Hue[0.5`],Hue[1.0`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph1c=Plot[{g6p[t]/.solution,pep[t]/.solution},{t,0,tinterm},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, {tinterm-5,0.045}],PlotStyle->colors,PlotRange->{{0,tinterm},{0,0.1}}];
graph2c=Plot[{g6p[t]/.solution,pep[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, {tfinal-5,0.045}],PlotStyle->colors,PlotRange->{{0,tfinal},{0,0.1}}];
{graph1c,graph2c}

labels = {"atp","amp"};
colors = {Hue[0.5`],Hue[1.0`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph3c=Plot[{atp[t]/.solution,amp[t]/.solution},{t,0,tinterm},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, {tinterm-4,0.6}],PlotStyle->colors,PlotRange->{{0,tinterm},{0,1.6}}];
graph4c=Plot[{atp[t]/.solution,amp[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, {tfinal-5,0.6}],PlotStyle->colors,PlotRange->{{0,tfinal},{0,1.6}}];
{graph3c,graph4c}

labels = {"lac"};
colors = {Hue[0.5`],Hue[1.0`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph5c=Plot[{lac[t]/.solution},{t,0,tinterm},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, {tinterm-4,1.2}],PlotStyle->colors,PlotRange->{{0,tinterm},{1.0,1.5}}];

labels = {"nad"};
colors = {Hue[0.5`],Hue[1.0`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph6c=Plot[{nad[t]/.solution},{t,0,tinterm},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, {tinterm-4,0.07}],PlotStyle->colors,PlotRange->{{0,tinterm},{0.04,0.08}}];

labels = {"pH"};
colors = {Hue[0.5`]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph7c=Plot[{-Log[10,h[t]]+3/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->Inset[legend, {tinterm-5,7.07}],PlotStyle->colors,PlotRange->{{0,tinterm},{7,7.1}}];
{graph5c,graph6c,graph7c}

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->{{1,3},{1,3}},AspectRatio->1/2 ];
graph9pp=ParametricPlot[Evaluate[{tpoolfluxin[[9]],tpoolfluxout[[9]]}]/.solution,{t,0,tfinal},PlotLabel->totalpoolnames[[9]],PlotRange->{{1,3},{1,3}},AspectRatio->1/2 ];

{graph1pp,graph2pp,graph3pp,graph4pp,graph5pp,graph6pp,graph7pp,graph8pp,graph9pp}

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

colors = {Hue[0.7],Hue[1.0]};
labels = {"in","out"};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];

graph1vpp=Plot[{tpoolfluxin[[1]]/.solution,tpoolfluxout[[1]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[1]],PlotStyle->colors, PlotRange->{3,5}, Epilog->Inset[legend, {tfinal-5,4 }] ];
graph2vpp=Plot[{tpoolfluxin[[2]]/.solution,tpoolfluxout[[2]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[2]],PlotStyle->colors, PlotRange->{1,3}, Epilog->Inset[legend, {tfinal-5,4 }] ];
graph3vpp=Plot[{tpoolfluxin[[3]]/.solution,tpoolfluxout[[3]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[3]],PlotStyle->colors, PlotRange->{3,5}, Epilog->Inset[legend, {tfinal-5,4 }] ];graph4vpp=Plot[{tpoolfluxin[[4]]/.solution,tpoolfluxout[[4]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[4]],PlotStyle->colors, PlotRange->{3,5}, Epilog->Inset[legend, {tfinal-5,4 }] ];graph5vpp=Plot[{tpoolfluxin[[5]]/.solution,tpoolfluxout[[5]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[5]],PlotStyle->colors, PlotRange->{3,5}, Epilog->Inset[legend, {tfinal-5,4 }] ];graph6vpp=Plot[{tpoolfluxin[[6]]/.solution,tpoolfluxout[[6]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[6]],PlotStyle->colors, PlotRange->{1.5,2.5}, Epilog->Inset[legend, {tfinal-5,2}] ];graph7vpp=Plot[{tpoolfluxin[[7]]/.solution,tpoolfluxout[[7]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[7]],PlotStyle->colors, PlotRange->{1.5,2.5}, Epilog->Inset[legend, {tfinal-5,2}] ];graph8vpp=Plot[{tpoolfluxin[[8]]/.solution,tpoolfluxout[[8]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[8]],PlotStyle->colors, PlotRange->{1.5,2.5}, Epilog->Inset[legend, {tfinal-5,2}] ];graph9vpp=Plot[{tpoolfluxin[[9]]/.solution,tpoolfluxout[[9]]/.solution},{t,0,tfinal},PlotLabel->totalpoolnames[[9]],PlotStyle->colors, PlotRange->{1.5,2.5}, Epilog->Inset[legend, {tfinal-5,2}] ];
{graph1vpp,graph2vpp,graph3vpp,graph4vpp,graph5vpp,graph6vpp,graph7vpp,graph8vpp,graph9vpp}

Pools

Calculate pools at start and end of simulations

should we qc/qa at t=0 relative to the pools computed above

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

Time responses of pools

colors = {Hue[0.5],Hue[0.6]};

labels = {totalpoolnames[[1]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph1p=Plot[{tpool[[1]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{2.0,4},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[2]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph2p=Plot[{tpool[[2]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{1.0,2.0},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[3]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph3p=Plot[{tpool[[3]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{2.0,4},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[4]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph4p=Plot[{tpool[[4]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{0,1.5},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[5]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph5p=Plot[{tpool[[5]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{3.5,5},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[6]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph6p=Plot[{tpool[[6]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{0.1,0.2},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[7]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph7p=Plot[{tpool[[7]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{0.02,.04},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[8]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph8p=Plot[{tpool[[8]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{2.0,4},Epilog->Inset[legend, Automatic]];
labels = {totalpoolnames[[9]]};
legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph9p=Plot[{tpool[[9]]/.solution},{t,0,tfinal},PlotLabel->"Pools vs time",PlotStyle->colors,PlotRange->{0,.2},Epilog->Inset[legend, Automatic]];

{graph1p,graph2p,graph3p,graph4p,graph5p,graph6p,graph7p,graph8p,graph9p}

Ratios

Compute the ratios at the beginning and the end

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

MatrixForm[tset5,TableHeadings->{rationames,{0,tinterm,tfinal}}];

x=50;

tratio = {0, 0, 0, 0, 0};
Do[tratio[[i]] = (rnumerator[[i]].tpool/rdenominator[[i]].tpool), {i, 1, 5}]
Manipulate[
tset5 = Table[Flatten[{tratio[[i]] /. t->x} /. solution], {i, 1, Length[tratio]}]; 
BarChart[Transpose[tset5][[1]], BarLabels -> rationames,ImageSize->{320,240
}],{x,0,50,Appearance -> "Open"}]

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.75},{0.6,0.9}},AspectRatio->1/1];
graph2rp=ParametricPlot[Evaluate[{tratio[[3]],tratio[[4]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{.95,1.0},{0.2,0.4}},AspectRatio->1/1];
graph3rp=ParametricPlot[Evaluate[{tratio[[2]],tratio[[5]]}]/.solution,{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait",PlotRange->{{0.6,0.9},{0.94,1.0}},AspectRatio->1/1];
{graph1rp,graph2rp,graph3rp}

Compute the time responses of the ratios

colors = {Hue[0.3],Hue[0.6],Hue[0.8]};

labels = {rationames[[1]]};legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph1r=Plot[{tratio[[1]]/.solution},{t,0,tfinal},PlotLabel->"Ratios vs time",PlotStyle->colors,PlotRange->{0.6,0.9},Epilog->Inset[legend, Automatic]];
labels = {rationames[[2]]};legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph2r=Plot[{tratio[[2]]/.solution},{t,0,tfinal},PlotLabel->"Ratios vs time",PlotStyle->colors,PlotRange->{0.6,0.9},Epilog->Inset[legend, Automatic]];
labels = {rationames[[3]]};legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph3r=Plot[{tratio[[3]]/.solution},{t,0,tfinal},PlotLabel->"Ratios vs time",PlotStyle->colors,PlotRange->{0.9,1.0},Epilog->Inset[legend, Automatic]];
labels = {rationames[[4]]};legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph4r=Plot[{tratio[[4]]/.solution},{t,0,tfinal},PlotLabel->"Ratios vs time",PlotStyle->colors,PlotRange->{0.2,0.4},Epilog->Inset[legend, Automatic]];
labels = {rationames[[5]]};legend = Framed[Column[Table[Style[labels[[i]],colors[[i]],Bold,"Graphics",12],{i,Length[labels]}]]];
graph5r=Plot[{tratio[[5]]/.solution},{t,0,tfinal},PlotLabel->"Ratios vs time",PlotStyle->colors,PlotRange->{0.95,1},Epilog->Inset[legend, Automatic]];

{graph1r,graph2r,graph3r,graph4r,graph5r}