BACKGROUND
metabolic network reconstructions <cit>  have been used as a basis for a number of analyses <cit>  that have provided insights into the topology  <cit> , modularity <cit> , robustness <cit> , and dynamics <cit>  of large biochemical networks.
in the constraint-based framework, the regulatory challenge for genome-scale metabolic networks has been described as a two-level process <cit> : first, regulatory mechanisms associated with transcription and translation geometrically delimit the steady-state flux space by determining which reactions can potentially carry flux; and second, regulation of gene product activity by post-translational mechanisms determines the flux state as a point location within the flux space.
the effective dimensionality of the first level of metabolic regulation has recently been shown to be small <cit> , but the effective dimensionality of the second level has yet to be assessed.
we approach this problem of cell-scale post-translational regulation in the context of presenting a method for the decomposition of the range of functional capabilities of large biochemical reaction systems.
we describe this decomposition procedure and demonstrate how it can elucidate a low number of reaction sets that account for nearly all of the range of behaviors in a cell-scale system.
RESULTS
our procedure is comprised of five main steps .
the reconstructed integrated transcriptional regulatory and metabolic network of e. coli  <cit>  was used for the analysis.
we first defined a growth environment and then used the transcriptional regulatory network to determine which reactions could be active.
this step corresponds to shrinking the flux space <cit> , and in effect reduces a 332-dimensional space to a 123-dimensional space.
the environment simulated was glucose aerobic minimal media conditions in which all media components  were allowed to vary from excess to limiting.
the analysis described herein is equally applicable to different and more complex environments–in particular, the local environment that an organism is constantly altering
monte carlo sampling is a method for generating large numbers of random allowable flux states and has been used to study the properties of metabolic flux states <cit> .
we comprehensively sampled the flux space corresponding to a growth rate of at least 90% of the maximum achievable growth rate and generated a large number  of flux vectors.
because the sampling procedure is a linear one and because we sought a basis for the sampled space, we performed principal components analysis using singular value decomposition.
the cumulative fractional eigenvalue distribution  reveals that 96% of the variation in the metabolic network flux states can be explained by seven principal components–implying that the post-translational regulatory problem is low-dimensional.
that is, by "regulating" a small number of dimensions the flux state of the entire network can be essentially set.
the implications and caveats associated with this interpretation are addressed below.
a biochemically meaningful interpretation of the eigenvectors was found through the use basis rotation methods <cit> , which are able to minimize the ambiguous association between metabolic reactions and eigenvectors in an information-preserving manner.
we rotated the top twenty eigenvectors and concentrated on the top seven .
we note three important results here.
first, following rotation the eigenvectors were comprised of distinct sets of metabolic reactions.
second, all oblique and orthogonal rotation methods tested in this study produced very similar results, indicating that a natural structure was latent in the random flux samples.
third, we computed the correlation between all pairs of rotated eigenvectors and found them to have low correlation .
we can interpret these independently-operable reaction sets, or eigenfluxes, as representing the regulatory challenge from a network perspective.
setting the seven eigenfluxes essentially positions the network flux state as a point within the flux space, and thus represents the second level of the two-level regulatory challenge <cit> .
the reaction loadings of the reactions on each eigenflux are all positive or are both positive and negative–indicating whether reactions in a set operate in a correlated  or anti-correlated  fashion.
the reactions that comprise each eigenflux are illustrated in figure 3a.
the reaction sets in figure 3a are the flux altering mechanisms utilized by the e. coli model to maintain high biomass formation in the varying environment studied.
these reaction sets can be examined in a biochemical and metabolic context:
<dig>  eigenfluxes  <dig>   <dig>   <dig>  and  <dig>  which together account for 49% of the flux variation, are metabolic "overflows."
overflow behavior in this context is due to oxygen limitation, and it allows the cell to generate atp via substrate level phosphorylation and to oxidize nadh into nad–which is in high demand.
both of these overflow-enabled mechanisms allow the cell to grow rapidly when glucose is in excess.
having a variety of overflow mechanisms allows the cell to balance its requirement to replenish nad with its need to produce energy.
the fluxes through these reactions are, in part, controlled by allosteric mechanisms.
phosphate acetyltransferase , which catalyzes the first reaction step of acetate secretion in eigenflux  <dig>  is allosterically activated by pyruvate and inhibited by nadh and nadph.
high concentrations of nadh are indicative of a redox imbalance, which eigenfluxes  <dig> and  <dig> serve to correct by oxidizing nadh.
<dig>  eigenflux  <dig> is associated with establishing the proton motive force.
by flexibly tuning the electro-chemical gradient of protons across the cellular membrane through translocation in the electron transport chain, a cell can economically provide the energy for essentially all cellular activity–such as atp synthesis, solute transport, and flagellar motility.
the cytochrome oxidases that catalyze the two opposing reactions in bimodal eigenflux  <dig> have slightly different functional behaviors.
the cytochrome bo3oxidase  is utilized under high oxygen concentrations and has a higher bioenergetic efficiency  than the other cytochrome bd oxidase .
even though cytochrome bd is energetically less efficient it is operational at low oxygen levels.
<dig>  eigenflux  <dig> is the flux variation through glycolysis, and functions with eigenflux  <dig> to determine the absolute magnitudes of the fluxes in glycolysis and the pentose phosphate pathway .
<dig>  eigenflux  <dig> is associated with the dominant tradeoff in central metabolism for providing the metabolic precursors, atp, and reducing power  to generate macromolecular building blocks against the need for reducing power  for sustaining the proton gradient.
until recently it was assumed that nearly all of the nadph needed for the biosynthesis reactions was produced in the ppp and the tca cycle, but recent experiments <cit>  have demonstrated that in glucose aerobic conditions the proton-coupled transhydrogenase reaction provides 35%–45% of the nadph needed for biosynthesis.
usage of the blue reaction in eigenflux  <dig> is consistent with this finding.
allosteric mechanisms at least in part mediate the balance between how the proton gradient is utilized and how flux is split between glycolysis and the ppp; the enzyme glucose-6-phophate-1-dehydrogenase, which catalyzes the first reaction of the ppp, is allosterically inhibited by nadh.
the flux carried through eigenfluxes in a given environment will be a result of molecular regulatory mechanisms, including mass action kinetics.
the contribution of an eigenflux to any given steady-state flux vector ν can be calculated from the equation
β = ut * ,
where u is the matrix of eigenfluxes and <ν > is the vector of mean flux values  as computed from the set of random flux samples.
as can be seen from the inverse equation
ν = <ν > + u* β,
the β values can be viewed as "tuning," or biasing, parameters–each β defining how the flux values for the reactions in its associated eigenflux are biased from their mean values.
thus one can effectively determine the flux state of the entire network with a low number of continuously varying and readily interpretable parameters.
to illustrate, we identified the two flux distributions that most oppositely utilized eigenflux  <dig> and computed their respective β vectors using equation  <dig>  the two selected flux distributions  are presented as colored β vectors in figure 3b.
the color coding allows one to quickly identify that, while fv <dig> is secreting no acetate and is utilizing the cytochrome oxidase that translocates fewer protons , fv <dig> is maximally excreting acetate and almost exclusively utilizing the other cytochrome oxidase.
furthermore, fv <dig> is predominantly utilizing the proton gradient for atp synthesis, while fv <dig> is utilizing it more for converting nadh to nadph.
the lower glycolytic flux of fv <dig> indicates a higher ppp flux and commensurately higher production of nadph, which allows the proton gradient to be used more for atp synthesis instead of nadph production.
the remaining components of β can be similarly interpreted, and altogether allow one to assess the flux state of the entire metabolic network at essentially a glance.
discussion
the faithfulness of the computational results of the presented procedure to biological reality depends critically on the completeness and integration of the molecular system reconstructions, assumptions made in the transformation of the integrated reconstructions into a model, and the manner in which the range of the system's functional capabilities are sampled.
cells are more than a combination of transcriptional regulatory and metabolic systems, so the exclusion of other cellular systems limits the scope of the decomposition that can be performed–as does the completeness of the included systems.
for instance, with other systems included, the decomposition procedure could potentially identify osmotic or movement  mechanisms in addition to purely metabolic mechanisms.
similarly, more complete regulatory information would allow a more accurate setting of reactions that can potentially carry flux.
the fixed biomass composition in the utilized model does not accurately describe a cell in all growth environments for all growth rates.
such a fixed composition limits the exploration of how the relative amounts of the biomass components can vary or how the cell can utilize different operating regimes <cit> .
the range of behaviors of large, interacting systems of molecules is defined by high-dimensional mathematical spaces that are non-trivial to fully explore.
monte carlo sampling of such spaces is not a solved problem, so the extent to which the full range of system capabilities can be sampled is directly related to the extent to which high-dimensional spaces can be computationally interrogated.
while these issues are important caveats, they are also the subject of active research and so will gradually diminish in their limiting roles.
molecular network reconstructions enable the objective, top-down assessment of regulatory challenges and functional capabilities associated with particular phenotypic states.
in the context of metabolism, the method presented herein is also applicable to kinetic and free energy parameter spaces and to concentration space  <cit> .
in general, it is applicable to any reconstructed cellular network in any environment and will aid in identifying  network regulatory challenges without the need to know detailed mechanisms or numerical parameter values.
our results also shed light on the network topology-function relationship, which is the result of a lengthy evolutionary process in varying environments.
robustness <cit> , defined as the ability to maintain specific functions in the face of varying environmental conditions, is believed to constitute a primary determinant of the topology-function relationship.
since this relationship is manifested in the range of flux states that the network can support, our investigation illuminates how the network topology confers the ability to robustly maintain a high growth rate in a dynamic environment through the use of a small number of reaction mechanisms.
similarly, the demonstrated procedure can be used to elucidate potential evolutionary mechanisms.
by modifying the sampling procedure it is possible to simulate the evolution  of an organism in a particular growth environment.
the decomposition procedure would identify the flux adjustment routes  by which evolution would be achieved.
CONCLUSIONS
with a top-down view comes an understanding of what must be controlled to attain a network flux state, and with a reductionist view comes an understanding of the mechanisms that achieve such control.
the top-down view will provide a context for the seemingly overlapping and redundant regulatory strategies that must make sense in not one environment, but in a large and varied range of environments.
the analysis method presented here will help bridge the systems and molecular biology approaches for understanding cellular regulation in the context of large biochemical networks.

 metabolic network reconstructions have been used as a basis for a number of analyses that have provided insights into the topology, modularity, robustness, and dynamics of large biochemical networks. in the constraint-based framework, the regulatory challenge for genome-scale metabolic networks has been described as a two-level process: first, regulatory mechanisms associated with transcription and translation geometrically delimit the steady-state flux space by determining which reactions can potentially carry flux; and second, regulation of gene product activity by post-translational mechanisms determines the flux state as a point location within the flux space. the effective dimensionality of the first level of metabolic regulation has recently been shown to be small, but the effective dimensionality of the second level has yet to be assessed. we approach this problem of cell-scale post-translational regulation in the context of presenting a method for the decomposition of the range of functional capabilities of large biochemical reaction systems. we describe this decomposition procedure and demonstrate how it can elucidate a low number of reaction sets that account for nearly all of the range of behaviors in a cell-scale system.  our procedure is comprised of five main steps. the reconstructed integrated transcriptional regulatory and metabolic network of e. coli was used for the analysis. we first defined a growth environment and then used the transcriptional regulatory network to determine which reactions could be active. this step corresponds to shrinking the flux space, and in effect reduces a 332-dimensional space to a 123-dimensional space. the environment simulated was glucose aerobic minimal media conditions in which all media components were allowed to vary from excess to limiting. the analysis described herein is equally applicable to different and more complex environments–in particular, the local environment that an organism is constantly altering monte carlo sampling is a method for generating large numbers of random allowable flux states and has been used to study the properties of metabolic flux states. we comprehensively sampled the flux space corresponding to a growth rate of at least 90% of the maximum achievable growth rate and generated a large number of flux vectors. because the sampling procedure is a linear one and because we sought a basis for the sampled space, we performed principal components analysis using singular value decomposition. the cumulative fractional eigenvalue distribution reveals that 96% of the variation in the metabolic network flux states can be explained by seven principal components–implying that the post-translational regulatory problem is low-dimensional. that is, by "regulating" a small number of dimensions the flux state of the entire network can be essentially set. the implications and caveats associated with this interpretation are addressed below. a biochemically meaningful interpretation of the eigenvectors was found through the use basis rotation methods, which are able to minimize the ambiguous association between metabolic reactions and eigenvectors in an information-preserving manner. we rotated the top twenty eigenvectors and concentrated on the top seven. we note three important results here. first, following rotation the eigenvectors were comprised of distinct sets of metabolic reactions. second, all oblique and orthogonal rotation methods tested in this study produced very similar results, indicating that a natural structure was latent in the random flux samples. third, we computed the correlation between all pairs of rotated eigenvectors and found them to have low correlation. we can interpret these independently-operable reaction sets, or eigenfluxes, as representing the regulatory challenge from a network perspective. setting the seven eigenfluxes essentially positions the network flux state as a point within the flux space, and thus represents the second level of the two-level regulatory challenge. the reaction loadings of the reactions on each eigenflux are all positive or are both positive and negative–indicating whether reactions in a set operate in a correlated or anti-correlated fashion. the reactions that comprise each eigenflux are illustrated in figure 3a. the reaction sets in figure 3a are the flux altering mechanisms utilized by the e. coli model to maintain high biomass formation in the varying environment studied. these reaction sets can be examined in a biochemical and metabolic context: eigenfluxes and which together account for 49% of the flux variation, are metabolic "overflows." overflow behavior in this context is due to oxygen limitation, and it allows the cell to generate atp via substrate level phosphorylation and to oxidize nadh into nad–which is in high demand. both of these overflow-enabled mechanisms allow the cell to grow rapidly when glucose is in excess. having a variety of overflow mechanisms allows the cell to balance its requirement to replenish nad with its need to produce energy. the fluxes through these reactions are, in part, controlled by allosteric mechanisms. phosphate acetyltransferase, which catalyzes the first reaction step of acetate secretion in eigenflux is allosterically activated by pyruvate and inhibited by nadh and nadph. high concentrations of nadh are indicative of a redox imbalance, which eigenfluxes and serve to correct by oxidizing nadh. eigenflux is associated with establishing the proton motive force. by flexibly tuning the electro-chemical gradient of protons across the cellular membrane through translocation in the electron transport chain, a cell can economically provide the energy for essentially all cellular activity–such as atp synthesis, solute transport, and flagellar motility. the cytochrome oxidases that catalyze the two opposing reactions in bimodal eigenflux have slightly different functional behaviors. the cytochrome bo3oxidase is utilized under high oxygen concentrations and has a higher bioenergetic efficiency than the other cytochrome bd oxidase. even though cytochrome bd is energetically less efficient it is operational at low oxygen levels. eigenflux is the flux variation through glycolysis, and functions with eigenflux to determine the absolute magnitudes of the fluxes in glycolysis and the pentose phosphate pathway. eigenflux is associated with the dominant tradeoff in central metabolism for providing the metabolic precursors, atp, and reducing power to generate macromolecular building blocks against the need for reducing power for sustaining the proton gradient. until recently it was assumed that nearly all of the nadph needed for the biosynthesis reactions was produced in the ppp and the tca cycle, but recent experiments have demonstrated that in glucose aerobic conditions the proton-coupled transhydrogenase reaction provides 35%–45% of the nadph needed for biosynthesis. usage of the blue reaction in eigenflux is consistent with this finding. allosteric mechanisms at least in part mediate the balance between how the proton gradient is utilized and how flux is split between glycolysis and the ppp; the enzyme glucose-6-phophate-1-dehydrogenase, which catalyzes the first reaction of the ppp, is allosterically inhibited by nadh. the flux carried through eigenfluxes in a given environment will be a result of molecular regulatory mechanisms, including mass action kinetics. the contribution of an eigenflux to any given steady-state flux vector ν can be calculated from the equation β = ut *, where u is the matrix of eigenfluxes and <ν > is the vector of mean flux values as computed from the set of random flux samples. as can be seen from the inverse equation ν = <ν > + u* β, the β values can be viewed as "tuning," or biasing, parameters–each β defining how the flux values for the reactions in its associated eigenflux are biased from their mean values. thus one can effectively determine the flux state of the entire network with a low number of continuously varying and readily interpretable parameters. to illustrate, we identified the two flux distributions that most oppositely utilized eigenflux and computed their respective β vectors using equation the two selected flux distributions are presented as colored β vectors in figure 3b. the color coding allows one to quickly identify that, while fv is secreting no acetate and is utilizing the cytochrome oxidase that translocates fewer protons, fv is maximally excreting acetate and almost exclusively utilizing the other cytochrome oxidase. furthermore, fv is predominantly utilizing the proton gradient for atp synthesis, while fv is utilizing it more for converting nadh to nadph. the lower glycolytic flux of fv indicates a higher ppp flux and commensurately higher production of nadph, which allows the proton gradient to be used more for atp synthesis instead of nadph production. the remaining components of β can be similarly interpreted, and altogether allow one to assess the flux state of the entire metabolic network at essentially a glance. discussion the faithfulness of the computational results of the presented procedure to biological reality depends critically on the completeness and integration of the molecular system reconstructions, assumptions made in the transformation of the integrated reconstructions into a model, and the manner in which the range of the system's functional capabilities are sampled. cells are more than a combination of transcriptional regulatory and metabolic systems, so the exclusion of other cellular systems limits the scope of the decomposition that can be performed–as does the completeness of the included systems. for instance, with other systems included, the decomposition procedure could potentially identify osmotic or movement mechanisms in addition to purely metabolic mechanisms. similarly, more complete regulatory information would allow a more accurate setting of reactions that can potentially carry flux. the fixed biomass composition in the utilized model does not accurately describe a cell in all growth environments for all growth rates. such a fixed composition limits the exploration of how the relative amounts of the biomass components can vary or how the cell can utilize different operating regimes. the range of behaviors of large, interacting systems of molecules is defined by high-dimensional mathematical spaces that are non-trivial to fully explore. monte carlo sampling of such spaces is not a solved problem, so the extent to which the full range of system capabilities can be sampled is directly related to the extent to which high-dimensional spaces can be computationally interrogated. while these issues are important caveats, they are also the subject of active research and so will gradually diminish in their limiting roles. molecular network reconstructions enable the objective, top-down assessment of regulatory challenges and functional capabilities associated with particular phenotypic states. in the context of metabolism, the method presented herein is also applicable to kinetic and free energy parameter spaces and to concentration space. in general, it is applicable to any reconstructed cellular network in any environment and will aid in identifying network regulatory challenges without the need to know detailed mechanisms or numerical parameter values. our results also shed light on the network topology-function relationship, which is the result of a lengthy evolutionary process in varying environments. robustness, defined as the ability to maintain specific functions in the face of varying environmental conditions, is believed to constitute a primary determinant of the topology-function relationship. since this relationship is manifested in the range of flux states that the network can support, our investigation illuminates how the network topology confers the ability to robustly maintain a high growth rate in a dynamic environment through the use of a small number of reaction mechanisms. similarly, the demonstrated procedure can be used to elucidate potential evolutionary mechanisms. by modifying the sampling procedure it is possible to simulate the evolution of an organism in a particular growth environment. the decomposition procedure would identify the flux adjustment routes by which evolution would be achieved.  with a top-down view comes an understanding of what must be controlled to attain a network flux state, and with a reductionist view comes an understanding of the mechanisms that achieve such control. the top-down view will provide a context for the seemingly overlapping and redundant regulatory strategies that must make sense in not one environment, but in a large and varied range of environments. the analysis method presented here will help bridge the systems and molecular biology approaches for understanding cellular regulation in the context of large biochemical networks. 
the control of cellular metabolic state is mediated by numerous molecular mechanisms, which in effect position the metabolic network flux state at specific locations within a mathematically-definable steady-state flux space. post-translational regulation constitutes a large class of these mechanisms, and decades of research indicate that achieving a network flux state through post-translational metabolic regulation is both a complex and complicated regulatory problem. we show that the use of monte carlo sampling of the steady-state flux space of a cell-scale metabolic system in conjunction with principal component analysis and eigenvector rotation results in a low-dimensional and biochemically interpretable decomposition of the steady flux states of the system. this result indicates an underlying simplicity and implies that the regulation of a relatively low number of reaction sets can essentially determine the flux state of the entire network in the given growth environment. our approach complements the reductionist approach to elucidation of regulatory mechanisms and facilitates the development of our understanding of global regulatory strategies in biological networks. no analysis method for the objective, top-down assessment of such regulation problems in large biochemical networks has been presented and demonstrated.
 metabolism and its regulation constitute a large fraction of the molecular activity within cells. the control of cellular metabolic state is mediated by numerous molecular mechanisms, which in effect position the metabolic network flux state at specific locations within a mathematically-definable steady-state flux space. post-translational regulation constitutes a large class of these mechanisms, and decades of research indicate that achieving a network flux state through post-translational metabolic regulation is both a complex and complicated regulatory problem. no analysis method for the objective, top-down assessment of such regulation problems in large biochemical networks has been presented and demonstrated.  we show that the use of monte carlo sampling of the steady-state flux space of a cell-scale metabolic system in conjunction with principal component analysis and eigenvector rotation results in a low-dimensional and biochemically interpretable decomposition of the steady flux states of the system. this decomposition comes in the form of a low number of small reaction sets whose flux variability accounts for nearly all of the flux variability in the entire system. this result indicates an underlying simplicity and implies that the regulation of a relatively low number of reaction sets can essentially determine the flux state of the entire network in the given growth environment.  we demonstrate how our top-down analysis of networks can be used to determine key regulatory requirements independent of specific parameters and mechanisms. our approach complements the reductionist approach to elucidation of regulatory mechanisms and facilitates the development of our understanding of global regulatory strategies in biological networks. 