(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



Off[General::obspkg]
Off[General::shdw]


<<"BarCharts`";<<"Histograms`"
<<"PieCharts`"
<<"MultivariateStatistics`"
<<"ComputationalGeometry`"
<<"ErrorBarPlots`"
<<"Graphics`Arrow`" 
<<"PlotLegends`"


pStoich[stoich_]:=Table[Which[stoich[[i,j]]>0,stoich[[i,j]],stoich[[i,j]]<=0,0],{i,1,Length[stoich]},{j,1,Length[Transpose[stoich]]}];
pStoich::"usage"="Returns matrix with all negative values -> 0.";


nStoich[stoich_]:=Table[Which[stoich[[i,j]]<0,stoich[[i,j]],stoich[[i,j]]>=0,0],{i,1,Length[stoich]},{j,1,Length[Transpose[stoich]]}];
nStoich::"usage"="Returns matrix with all positive values -> 0.";


solveODE[eq_,ic_:Table[0,{Length[var]}],var_,starttime_:0,endtime_:100,MAXSTEPS_:10000,WORKINGPRECISION_:MachinePrecision]:=Module[{i,tmet,teq,set,counter,answer},(* If[Length[eq]!=Length[var],Print["Number of equations and iables must match"];Return[]]; *) 
tmet=tsub[var];teq=eq/. tmet;set=Table[Null,{2 Length[eq]}];Do[set[[i]]=(
\!\(\*SuperscriptBox[
RowBox[{"var", "\[LeftDoubleBracket]", "i", "\[RightDoubleBracket]"}], "\[Prime]",
MultilineFunction->None]\)[t]==teq[[i]]);,{i,1,Length[teq]}];counter=Length[teq]+1; Do[set[[counter]]=(var[[i]][0]==ic[[i]]);counter+=1;,{i,1,Length[teq]}]; 
answer=NDSolve[SetPrecision[set, WORKINGPRECISION],var,{t,starttime,endtime},MaxSteps->MAXSTEPS, WorkingPrecision-> WORKINGPRECISION];Return[answer]];
solveODE::"usage"="Solve a set of ODEs specified by eq in terms of variables, from starting time t1, to t2. Use MAXSTEPS to adjust the number of steps for the ODE solver.  \n\tresult = solveODE[eq, ic,variables, t1, t2];";


foldlist[list1_,list2_]:=Module[{l1,l2,dummy,result},l1=Length[list1];l2=Length[list2];If[l1!=l2,Return["The lists must be of equal length"]];result=Partition[Flatten[Table[{list1[[dummy]],list2[[dummy]]},{dummy,1,l1}]],2]];


plotset[eq_,mets_,n1_,n_,initial_,group_,groupname_:Null,markers_:{0, 1, 10, 100},answer1_:Null,lfont_:9,nlfont_:8,lsize_:0,psize_:3,points_:50,plrng_:0.999`]:=Module[{answer,dummy,dummy1,dummy2,element,var,all,tmet,final,power,increment1,increment2,increment3,minute,hour,rest,fullsettemp,fullset,vertices,set,begin,r,test,fiteq,slope,sval,slopemat,fred,red,gr,bl,textnum,textslope,upper,size,font,nfont,loop,t0,pone,one,ten,hundred, fullmarkers, start, end},
If[Length[markers] < 4, fullmarkers = Table[n1, {4}]; Do[fullmarkers[[dummy]] = markers[[dummy]];, {dummy, 1, Length[markers]}];, fullmarkers = markers;];
If[Min[fullmarkers] < n1, start = Min[fullmarkers]; , start = n1];
If[ Max[fullmarkers] > n, end = Max[fullmarkers], end = n];
If[answer1===Null,answer=solveODE[eq,initial,mets,start,end],answer=answer1];tmet=tsub[mets];power=Log[10,n];If[n1<1/3600,start=1/3600,start=n1];If[lsize==0,size=25 (Length[group]+3),size=lsize];If[lfont==0,font=9,font=lfont];If[nlfont==0,nfont=8,nfont=nlfont];If[lsize==0||lfont==0||nlfont==0,Print["label font = ",font,"   Number font = ",nfont,"   Image Size = ",size]];If[start>=n,Return[Break["Error! Please check the intial and final time points to make sure they are consistent with one another!"]];];fullset=Table[0,{dummy,1,points+1}];loop=1;dummy1=0;t0=start;While[loop<=points+1,dummy1=N[t0 2^(((loop-1) (Log[2,n]-Log[2,t0]))/points)];fullset[[loop]]=dummy1;loop++;];Table[var[dummy]=Flatten[Table[group[[dummy]]/. tmet/. t->fullset[[dummy1]]/. answer,{dummy1,1,Length[fullset]}]],{dummy,1,Length[group]}];

Table[pone[dummy]=Flatten[Table[group[[dummy]]/. tmet/. t->fullmarkers[[1]]/. answer,{dummy1,1,1}]],{dummy,1,Length[group]}];
Table[one[dummy]=Flatten[Table[group[[dummy]]/. tmet/. t->fullmarkers[[2]]/. answer,{dummy1,1,1}]],{dummy,1,Length[group]}];
Table[ten[dummy]=Flatten[Table[group[[dummy]]/. tmet/. t->fullmarkers[[3]]/. answer,{dummy1,1,1}]],{dummy,1,Length[group]}];
Table[hundred[dummy]=Flatten[Table[group[[dummy]]/. tmet/. t->fullmarkers[[4]]/. answer,{dummy1,1,1}]],{dummy,1,Length[group]}];

upper=Table[0,{Length[group]},{Length[group]}];r=Table[1,{dummy1,Length[group]},{dummy2,Length[group]}];slopemat=Table[1,{dummy1,Length[group]},{dummy2,Length[group]}];Do[test=Correlation[var[dummy1],var[dummy2]];r[[dummy2,dummy1]]=test;r[[dummy1,dummy2]]=test;If[Abs[test]>0.85`,fiteq=Fit[foldlist[var[dummy2],var[dummy1]],{1,x},x];slope=(fiteq/. x->1)-(fiteq/. x->0);slopemat[[dummy2,dummy1]]=1/slope;slopemat[[dummy1,dummy2]]=slope;sval=FromDigits[RealDigits[Abs[slopemat[[dummy1,dummy2]]]+0.005`]];If[slopemat[[dummy1,dummy2]]<0,textslope=StringTake[ToString[N[sval,24]],Max[4,Round[Log[10,Abs[slopemat[[dummy1,dummy2]]]]]]],textslope=StringTake[ToString[N[sval,24]],Max[4,Round[Log[10,Abs[slopemat[[dummy1,dummy2]]]]]]]];];sval=FromDigits[RealDigits[Abs[r[[dummy1,dummy2]]]+0.005`]];If[r[[dummy1,dummy2]]<0,textnum=StringTake[ToString[N[sval,24]],4];fred=1,textnum=StringTake[ToString[N[sval,24]],4];fred=0];red=1;gr=1;bl=1;Which[r[[dummy1,dummy2]]>=0.9`,red=0.8`;gr=0.8`;bl=0.8`,r[[dummy1,dummy2]]<=-0.9`,red=0.8`;gr=0.7`;bl=0.7`,r[[dummy1,dummy2]]>-0.25`&&r[[dummy1,dummy2]]<0.25`,red=1;gr=0.8`;bl=0.8`];If[Abs[r[[dummy1,dummy2]]]>0.85`,upper[[dummy1,dummy2]]=Graphics[{{RGBColor[red,gr,bl],Rectangle[{0,0},{1,1}]},{RGBColor[fred,0,0],Text[textnum<>"\n"<>textslope,{0.5`,0.5`}],BaseStyle->{FontFamily->"Arial",FontSize->nfont}}},AspectRatio->1,Frame->False],upper[[dummy1,dummy2]]=Graphics[{{RGBColor[red,gr,bl],Rectangle[{0,0},{1,1}]},{RGBColor[fred,0,0],Text[textnum,{0.5`,0.5`}],BaseStyle->{FontFamily->"Arial",FontSize->nfont}}},AspectRatio->1,Frame->False]],{dummy1,1,Length[group]},{dummy2,1,Length[group]}];
all=Table[Which[dummy1==1&&dummy2==1,Graphics[{Thickness[0.01`],RGBColor[1,1,1],Rectangle[{0,0},{1,1}]},AspectRatio->1],dummy1==dummy2,Graphics[{Thickness[0.01`],RGBColor[0,0,0],Rectangle[{0,0},{1,1}]},AspectRatio->1],dummy1==1,Graphics[Text[groupname[[dummy2-1]],{0,0},{0,0},{0,-1}],BaseStyle->{FontFamily->"Arial",FontSize->font}],dummy2==1,Graphics[Text[groupname[[dummy1-1]],{0,0}],BaseStyle->{FontFamily->"Arial",FontSize->font}],dummy1>dummy2,ListPlot[{set=foldlist[var[dummy2-1],var[dummy1-1]], foldlist[pone[dummy2-1],pone[dummy1-1]],foldlist[one[dummy2-1],one[dummy1-1]],foldlist[ten[dummy2-1],ten[dummy1-1]],foldlist[hundred[dummy2-1],hundred[dummy1-1]],{set[[1]]}, {set[[Length[set]]]}},PlotRange->{{plrng Min[var[dummy2-1]],Max[var[dummy2-1]]/plrng},{plrng Min[var[dummy1-1]],Max[var[dummy1-1]]/plrng}},DisplayFunction->Identity,Axes->False,Joined->{True, False,False,False,False,False,False},PlotRange->All,AspectRatio->1,PlotStyle->{RGBColor[0,0,0],RGBColor[0.5,0,1],RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[1,1,0],RGBColor[1,1,1],RGBColor[0,0,1]} , PlotMarkers->{{""},{\[FilledDiamond]},{\[FilledDiamond]},{\[FilledDiamond]},{\[FilledDiamond]},{\[FilledSquare]},{\[FilledSquare]}}],dummy2>dummy1,upper[[dummy1-1]][[dummy2-1]]],{dummy1,1,Length[group]+1},{dummy2,1,Length[group]+1}];Show[GraphicsGrid[all],ImageSize->{size,size}]];
(* PlotMarkers->{{\[FilledDiamond]},{\[FilledDiamond]},{\[FilledDiamond]},{\[FilledDiamond]},{\[FilledDiamond]},{\[FilledSquare]},{\[FilledSquare]}} *)
plotset::"usage"="plotset[massbaleq, mets, time_zero, time ,initial ,group ,groupname ,plotpoints, answer1 ,font ,nfont ,size ,psize]\n\n    Produces an array of phase planes, plotting the components of the \nvariable list versus one another.\nThere are three required arguments:\n\ttime_zero: first time point to be plotted\n\ttime:time interval\n\tinitial: initial conditions\n\tgroup: the list of variables to be plotted. \nFor example a user may define a list: adenosine = {ATP, ADP, AMP, ADO} so \nthat the different adenosine metabolites will be plotted against one another.\n\n\t variable names: a string list of names corresponding to the variable \nlist. The default value will refer to the variable list. When the variable \nlist consists solely of metabolites, variable names does not need to be \nspecified. It needs to be specified if the variable list consists of fluxes.\n\n\tOptional arguments include\n\tanswer1 is the solution returned from 'solution'.\n\tfont is the font size\n\n\tnfont is the number font size\n\tsize is the plot size (in*72)\n\t psize is \nthe size of the points defining the ends of the trajectory. \nAdjusted so that Transpose of plots (compared to initial vrsn are plotted) - i.e. column labels correspond to y-axes and row labels correspond to x-axes.\tThe option 'plotpoints' is a list of time points (4) which the user can specify to have marked w/ diamonds.";


tsub[x_]:=(#1->#1[t]&)/@x;
tsub::"usage"="Makes variables into time based variables";


sublist[list_,values_]:=Table[list[[i]]->values[[i]],{i,1,Length[list]}];
sublist::"usage"="Given a list and set of values, returns a substitution list.\n\tsubbed = sublist[{v1,v2,v3},{0.1,0.2,0.4}];";


createRxnBounds[rxn_,tstart_,tend_,variables_,solution_]:=
Module[{rxnMin,rxnMax,FluxValues,x},x=rxn/.tsub[variables]/.solution;
FluxValues=Table[1,{i,(tend-tstart)*100+1}];
Do[FluxValues[[i+1]]=x/.t->tstart+i/100,{i,0,(tend-tstart)*100}];rxnMin=Min[FluxValues];
rxnMax=Max[FluxValues];
rxnMin=rxnMin-Abs[.25*rxnMin] ;
rxnMax=rxnMax+Abs[.25*rxnMax];
Return[{rxnMin,rxnMax}];];
createRxnBounds::"usage"="Determine the min and max boundaries for plotting a reaction.";


createMetBounds[met_,tstart_,tend_,solution_]:=
Module[{metMin,metMax,MetValues,x},x=met/.tsub[{met}];MetValues=Table[1,{i,(tend-tstart)*100+1}];Do[MetValues[[i+1]]={x/.t->tstart+i/100}/.solution,{i,0,(tend-tstart)*100}];
metMin=Min[MetValues];
metMax=Max[MetValues];
metMin=metMin-Abs[.25*metMin];
metMax=metMax+Abs[.25*metMax];
Return[{metMin,metMax}];];
createMetBounds::"usage"=
"Determine the min and max boundaries for plotting a concentration.";


createPoolBounds[pool_,tstart_,tend_,variables_,solution_]:=
Module[{poolMin,poolMax,PoolValues,x,y},y=variables/.tsub[variables];
x=pool.y;
PoolValues=Table[1,{i,(tend-tstart)*100+1}];
Do[PoolValues[[i+1]]={x/.t->tstart+i/100}/.solution,{i,0,(tend-tstart)*100}];
poolMin=Min[PoolValues];
poolMax=Max[PoolValues];
poolMin=poolMin-Abs[.25*poolMin];
poolMax=poolMax+Abs[.25*poolMax];
Return[{poolMin,poolMax}];];
createPoolBounds::"usage"="Determine the min and max boundaries for plotting pool concentrations.";


createPoolFluxBounds[poolNum_,tstart_,tend_,fluxes_,variables_,ppstoich_,pnstoich_,solution_]:=Module[{poolFluxMinIn,poolFluxMaxIn,poolFluxMinOut,poolFluxMaxOut,PoolFluxValuesIn,PoolFluxValuesOut,xin,xout,y},
y=fluxes/.tsub[variables];
xin=ppstoich.y;
xout=-pnstoich.y;
PoolFluxValuesIn=Table[1,{i,(tend-tstart)*100+1}];
PoolFluxValuesOut=Table[1,{i,(tend-tstart)*100+1}];
Do[PoolFluxValuesIn[[i+1]]={xin[[poolNum]]/.t->tstart+i/100}/.solution,{i,0,(tend-tstart)*100}];
Do[PoolFluxValuesOut[[i+1]]={xout[[poolNum]]/.t->tstart+i/100}/.solution,{i,0,(tend-tstart)*100}];
poolFluxMinIn=Min[PoolFluxValuesIn];
poolFluxMaxIn=Max[PoolFluxValuesIn];
poolFluxMinIn=poolFluxMinIn-Abs[.25*poolFluxMinIn];
poolFluxMaxIn=poolFluxMaxIn+Abs[.25*poolFluxMaxIn];
poolFluxMinOut=Min[PoolFluxValuesOut];
poolFluxMaxOut=Max[PoolFluxValuesOut];
poolFluxMinOut=poolFluxMinOut-Abs[.25*poolFluxMinOut];
poolFluxMaxOut=poolFluxMaxOut+Abs[.25*poolFluxMaxOut];
Return[{poolFluxMinIn,poolFluxMaxIn,poolFluxMinOut,poolFluxMaxOut}];];
createPoolFluxBoudns::"usage"="Determine the min and max boundaries for plotting pool fluxes.";


createRatioBounds[rNum_,tstart_,tend_,variables_,pools_,rnumerator_,rdenominator_,solution_]:=
Module[{ratioMin,ratioMax,RatioValues,x,y,z},
z=variables/.tsub[variables];
y=pools.z;
x=rnumerator[[rNum]].y/rdenominator[[rNum]].y;
RatioValues=Table[1,{i,(tend-tstart)*100+1}];
Do[RatioValues[[i+1]]={x/.t->tstart+i/100}/.solution,{i,0,(tend-tstart)*100}];
ratioMin=Min[RatioValues];
ratioMax=Max[RatioValues];
ratioMin=ratioMin-Abs[.25*ratioMin];
ratioMax=ratioMax+Abs[.25*ratioMax];
Return[{ratioMin,ratioMax}];];
createRatioBounds::"usage"="Determine the min and max boundaries for plotting ratios.";


MASSTimeProfileMet[mets_,metnames_,solution_,startT_,endT_]:=
Module[{numMets,TPplot,metMin,metMax,max,min},
numMets=Length[mets];
metMax=Table[0,{numMets}];
metMin=Table[0,{numMets}];
For[i=1,i<=numMets,i++,
{metMin[[i]],metMax[[i]]}=createMetBounds[mets[[i]],startT,endT,solution];];
max=Max[metMax];
min=Min[metMin];
TPplot=Plot[Evaluate[Table[mets[[i]][t]/.solution,{i,1,numMets}]],{t,startT,endT},PlotLabel->"Concentrations vs time",PlotLegend->metnames,LegendPosition->{1.1,-0.4},LegendSize->0.5,PlotRange->{{startT,endT},{min,max}}];
Return[TPplot]];
MASSTimeProfileMet::"usage"="Plot time profiles of select metabolites";


MASSTimeProfileFlux[rxns_,rxnNames_,variables_,solution_,startT_,endT_]:=
Module[{TPplot,numRxns,rxnMin,rxnMax,max,min,plots},
numRxns=Length[rxns];
rxnMax=Table[0,{numRxns}];
rxnMin=Table[0,{numRxns}];
For[i=1,i<=numRxns,i++,
{rxnMin[[i]],rxnMax[[i]]}=createRxnBounds[rxns[[i]],startT,endT,variables,solution];];
max=Max[rxnMax];
min=Min[rxnMin];
TPplot=Plot[Evaluate[Table[rxns[[i]]/.tsub[variables]/.solution,{i,numRxns}]],{t,startT,endT},PlotLabel->"Fluxes vs time",PlotLegend->rxnNames,LegendPosition->{1.1,-0.4},LegendSize->0.5,PlotRange->{{startT,endT},{min,max}}];
Return[TPplot]];
MASSTimeProfileFlux::"usage"="Plot time profiles of select reactions";


MASSTimeProfilePool[pools_,poolNames_,variables_,solution_,startT_,endT_]:=
Module[{TPplot,numPools,poolMin,poolMax,max,min,x,y},
numPools=Length[pools];
poolMax=Table[0,{numPools}];
poolMin=Table[0,{numPools}];
For[i=1,i<=numPools,i++,
{poolMin[[i]],poolMax[[i]]}=createPoolBounds[pools[[i]],startT,endT,variables,solution];];
max=Max[poolMax];
min=Min[poolMin];
x=variables/.tsub[variables];
y=Table[pools[[i]].x,{i,1,numPools}];
TPplot=Plot[Evaluate[Table[y[[i]]/.solution,{i,1,numPools}]],{t,startT,endT},PlotLabel->"Pools vs time",PlotLegend->poolNames,LegendPosition->{1.1,-0.4},LegendSize->0.5,PlotRange->{{startT,endT},{min,max}}];
Return[TPplot]];
MASSTimeProfilePool::"usage"="Plot time profiles of select pools";


MASSTimeProfilePoolFlux[poolNum_,poolName_,fluxes_,variables_,ppstoich_,pnstoich_,solution_,startT_,endT_]:=
Module[{TPplot,poolMinIn,poolMinOut,poolMaxIn,poolMaxOut,max,min,xin,xout,y},
{poolMinIn,poolMaxIn,poolMinOut,poolMaxOut}=createPoolFluxBounds[poolNum,startT,endT,fluxes,variables,ppstoich,pnstoich,solution];
max=Max[{poolMaxIn,poolMaxOut}];
min=Min[{poolMinIn,poolMinOut}];
y=fluxes/.tsub[variables];
xin=ppstoich.y;
xout=-pnstoich.y;
TPplot=Plot[{xin[[poolNum]]/.solution,xout[[poolNum]]/.solution},{t,startT,endT},PlotLabel->poolName,PlotLegend->{"In","Out"},LegendPosition->{1.1,-0.4},LegendSize->0.5,PlotRange->{{startT,endT},{min,max}}];
Return[{TPplot}]];
MASSTimeProfilePoolFlux::"usage"="Plot time profiles of select fluxes of pools.";


MASSTimeProfileRatios[ratioNums_,ratioNames_,variables_,pools_,rnumerator_,rdenominator_,solution_,startT_,endT_]:=
Module[{TPplot,numRatios,ratioMin,ratioMax,max,min,x,y,z},
numRatios=Length[ratioNums];
ratioMax=Table[0,{numRatios}];
ratioMin=Table[0,{numRatios}];
For[i=1,i<=numRatios,i++,
{ratioMin[[i]],ratioMax[[i]]}=createRatioBounds[ratioNums[[i]],startT,endT,variables,pools,rnumerator,rdenominator,solution];];
max=Max[ratioMax];
min=Min[ratioMin];
z=variables/.tsub[variables];
y=pools.z;
TPplot=Plot[Table[((rnumerator[[ratioNums[[i]]]].y)/(rdenominator[[ratioNums[[i]]]].y))/.solution,{i,1,numRatios}],{t,startT,endT},PlotLabel->"Ratios vs time",PlotLegend->ratioNames,LegendPosition->{1.1,-0.4},LegendSize->0.5,PlotRange->{{startT,endT},{min,max}}];
Return[TPplot]];
MASSTimeProfileRatios::"usage"="Plot time profiles of select ratios";


MASSPhasePlaneMet[mets_,metnames_,variables_,solution_,startT_,endT_]:=
Module[{met1,met2,metname1,metname2,PPplot,min1,max1,min2,max2},
met1=mets[[1]];
met2=mets[[2]];
metname1=metnames[[1]];
metname2=metnames[[2]];
{min1,max1}=createMetBounds[met1,startT,endT,solution];
{min2,max2}=createMetBounds[met2,startT,endT,solution];
PPplot=ParametricPlot[Evaluate[{met1,met2}]/.tsub[variables]/.solution,{t,startT,endT},PlotRange->{{min1,max1},{min2,max2}},AspectRatio->1/2,AxesLabel->{metname1,metname2}];
Return[PPplot]];
MASSPhasePlaneMet::"usage"="Plot phase plane of two metabolites";


MASSPhasePlaneFlux[rxns_,rxnnames_,variables_,solution_,startT_,endT_]:=
Module[{rxn1,rxn2,rxnname1,rxnname2,PPplot,min1,max1,min2,max2},
rxn1=rxns[[1]];
rxn2=rxns[[2]];
rxnname1=rxnnames[[1]];
rxnname2=rxnnames[[2]];
{min1,max1}=createRxnBounds[rxn1,startT,endT,variables,solution];
{min2,max2}=createRxnBounds[rxn2,startT,endT,variables,solution];
PPplot=ParametricPlot[Evaluate[{rxn1,rxn2}]/.tsub[variables]/.solution,{t,startT,endT},PlotRange->{{min1,max1},{min2,max2}},AspectRatio->1/2,AxesLabel->{rxnname1,rxnname2}];
Return[PPplot]];
MASSPhasePlaneFlux::"usage"="Plot phase plane of two reactions";


MASSPhasePlanePool[pools_,poolnames_,variables_,solution_,startT_,endT_]:=
Module[{pool1,pool2,poolname1,poolname2,PPplot,min1,max1,min2,max2,x,y1,y2},
pool1=pools[[1]];
pool2=pools[[2]];
poolname1=poolnames[[1]];
poolname2=poolnames[[2]];
{min1,max1}=createPoolBounds[pool1,startT,endT,variables,solution];
{min2,max2}=createPoolBounds[pool2,startT,endT,variables,solution];
x=variables/.tsub[variables];
y1=pool1.x;
y2=pool2.x;
PPplot=ParametricPlot[Evaluate[{y1,y2}]/.solution,{t,startT,endT},PlotRange->{{min1,max1},{min2,max2}},AspectRatio->1/2,AxesLabel->{poolname1,poolname2}];
Return[PPplot]];
MASSPhasePlanePool::"usage"="Plot phase plane of two pools";


MASSPhasePlanePoolFlux[poolNum_,poolname_,fluxes_,variables_,ppstoich_,pnstoich_,solution_,startT_,endT_]:=
Module[{PPplot,min1,max1,min2,max2,y,xin,xout},
{min1,max1,min2,max2}=createPoolFluxBounds[poolNum,startT,endT,fluxes,variables,ppstoich,pnstoich,solution];
y=fluxes/.tsub[variables];
xin=ppstoich.y;
xout=-pnstoich.y;
PPplot=ParametricPlot[Evaluate[{xin[[poolNum]],xout[[poolNum]]}]/.solution,{t,startT,endT},PlotRange->{{min1,max1},{min2,max2}},AspectRatio->1/2,AxesLabel->{"In","Out"},PlotLabel->poolname];
Return[PPplot]];
MASSPhasePlanePoolFlux::"usage"="Plot phase plane of two pool fluxes";


MASSPhasePlaneRatios[ratioNums_,rationames_,variables_,pools_,rnumerator_,rdenominator_,solution_,startT_,endT_]:=
Module[{ratio1,ratio2,rationame1,rationame2,PPplot,min1,max1,min2,max2,x,y,z1,z2},
ratio1=ratioNums[[1]];
ratio2=ratioNums[[2]];
rationame1=rationames[[1]];
rationame2=rationames[[2]];
{min1,max1}=createRatioBounds[ratio1,startT,endT,variables,pools,rnumerator,rdenominator,solution];
{min2,max2}=createRatioBounds[ratio2,startT,endT,variables,pools,rnumerator,rdenominator,solution];
x=variables/.tsub[variables];
y=pools.x;
z1=rnumerator[[ratio1]].y/rdenominator[[ratio1]].y;
z2=rnumerator[[ratio2]].y/rdenominator[[ratio2]].y;
PPplot=ParametricPlot[Evaluate[{z1,z2}]/.solution,{t,startT,endT},PlotRange->{{min1,max1},{min2,max2}},AspectRatio->1/2,AxesLabel->{rationame1,rationame2}];
Return[PPplot]];
MASSPhasePlaneRatios::"usage"="Plot phase plane of two ratios";
