Understanding Regulation of Metabolism through Feasibility Analysis

Understanding cellular regulation of metabolism is a major challenge in systems biology. Thus far, the main assumption was that enzyme levels are key regulators in metabolic networks. However, regulation analysis recently showed that metabolism is rarely controlled via enzyme levels only, but through non-obvious combinations of hierarchical (gene and enzyme levels) and metabolic regulation (mass action and allosteric interaction). Quantitative analyses relating changes in metabolic fluxes to changes in transcript or protein levels have revealed a remarkable lack of understanding of the regulation of these networks. We study metabolic regulation via feasibility analysis (FA). Inspired by the constraint-based approach of Flux Balance Analysis, FA incorporates a model describing kinetic interactions between molecules. We enlarge the portfolio of objectives for the cell by defining three main physiologically relevant objectives for the cell: function, robustness and temporal responsiveness. We postulate that the cell assumes one or a combination of these objectives and search for enzyme levels necessary to achieve this. We call the subspace of feasible enzyme levels the feasible enzyme space. Once this space is constructed, we can study how different objectives may (if possible) be combined, or evaluate the conditions at which the cells are faced with a trade-off among those. We apply FA to the experimental scenario of long-term carbon limited chemostat cultivation of yeast cells, studying how metabolism evolves optimally. Cells employ a mixed strategy composed of increasing enzyme levels for glucose uptake and hexokinase and decreasing levels of the remaining enzymes. This trade-off renders the cells specialized in this low-carbon flux state to compete for the available glucose and get rid of over-overcapacity. Overall, we show that FA is a powerful tool for systems biologists to study regulation of metabolism, interpret experimental data and evaluate hypotheses.


Introduction
In their natural habitat, most microbes are exposed to constantly changing physical and chemical environments. To perform optimally in these conditions, they must finely regulate their metabolism. Understanding how microbes regulate metabolism to achieve a desired objective, or how they adapt to changing conditions, is a major challenge [1]. Quantitative analyses relating changes in metabolic fluxes to changes in transcript or protein levels have further revealed a remarkable lack of understanding of the regulation of these networks; it remains unclear how and to what extent metabolic networks are regulated through the modulation of enzyme levels [2]. Regulation analysis has shown that metabolic networks are controlled via nonobvious combinations of metabolic and hierarchical regulation [3,4].
Thus far, in systems biology two main model-based approaches are used to study metabolic regulation: top-down and bottom-up. The top-down approach employs genome-wide constraint-based modeling techniques, such as Flux Balance Analysis (FBA), to find viable intracellular flux distributions based on measured external fluxes and thermodynamical considerations. Constraint-based models have been shown useful in exploring cellular capabilities of biological systems and have enabled in silico characterization of several phenotypic features, such as growth yield under gene knockouts (see [5] for a review). However, an inherent limitation of constraint-based models is that they are based solely on stoichiometry and thus are limited to predicting steady-state flux distributions. In general, they do not contain explicit regulation terms and cannot predict the effect of gene or enzyme dosage via knock-ins or point mutations. It is however possible to constrain the solution space by incorporating series of physiological parameters [6] or additional -omics data [7], or by assuming certain objectives for the cell [8]. The list of such objectives ranges from maximization of biomass to minimization of redox potential. A systematic evaluation [9] revealed that Escherichia coli employs different objectives under different conditions.
In contrast, the bottom-up approach combines detailed kinetic models with the theorems of Metabolic Control Analysis (MCA, [10]) or Biochemical Systems Theory (BST, [11]) to study regulation. While kinetic models describe metabolic reaction rates as a function of enzyme levels and metabolite concentrations, the inverse models describing (changes in) enzyme levels required to obtain desired metabolite concentrations or reaction rates are more useful for studying regulation. However, as most kinetic models are highly nonlinear, explicit inversion is often impossible. Both within the framework of MCA and BST, a number of approximative kinetic formats [12][13][14] have therefore been proposed as a solution [15][16][17]. Although useful, these kinetic descriptions usually offer limited mechanistic insights.
In this paper, we employ a method of studying regulation, Feasibility Analysis (FA), combining elements of bottom-up and top-down approaches. FA starts from an explicit kinetic model describing the interactions between enzymes and metabolites. Inspired by the well-established constraint based approach of FBA, it then defines a number of physicochemical constraints on the cell, as well as three physiologically relevant objectives: function, robustness and temporal responsiveness, for which quantitative measures are introduced. Assuming that the cell follows one or a combination of these objectives, FA then searches for (a) set(s) of enzyme levels necessary to achieve these. Given the problem of inversion of general non-linear kinetic models, FA uses a straightforward sampling-based method, commonly used for various computational biology purposes, e.g. for ensemble modeling [18], or modeling the uncertainty in biochemical reaction networks [19,20]. For each sampled set of enzyme levels, the kinetic model is integrated to steady state and objective measures are calculated on the resulting phenotype. We call the subspace encompassing all feasible enzyme levels the feasible enzyme space. Once this space is constructed, we can study how different objectives can (if possible) be combined, or evaluate the conditions under which these objectives are traded-off.
A similar approach of using physiological constraints to find feasible sets of enzyme levels was successfully applied to identify the required changes in gene expression in yeast upon heat shock [21] and, more generally, to attain certain cellular adaptive responses [22]. This method was adapted to study general design principles of metabolic networks, employing optimization techniques to explore the space of feasible enzyme levels [21,23]. While mathematically advanced, it is derived from a specific type of approximative kinetic model (Generalized Mass Action or GMA models), which limits its general use. FA aims (1) to generalize the GMA-based analysis by defining more generic, quantitative objectives that can be evaluated for any kinetic model; and (2) to get deeper understanding of regulation by explicitly incorporating the modes of regulation (metabolic or hierarchical) under physiological constraints and objectives.
The feasible enzyme spaces found by FA can also be used to enhance currently available kinetic models. These models are usually derived starting from an ab initio selected set of kinetic interactions; subsequently, parameter values are set or estimated by fitting to a (small) number of measurements. Methods to expand/shrink the model by adding/removing interactions and inspect the feasibility of the resulting models are of great interest. Using FA, we can thus discriminate between available hypotheses on how metabolism is regulated and evaluate potential changes in model structure.
In this paper, we first describe FA in detail, listing a number of constraints and introducing quantitative measures for the proposed objectives. We then exemplify the approach using two cases: (1) an illustrative small model with tractable kinetics and (2) a larger dynamic model of yeast glycolysis [24]. For yeast glycolysis, we analyze two scenarios: the adaptation of yeast cells during long-term chemostat cultivation under carbon limitation and the regulation of hexokinase to infer robustness to the glycolytic pathway. In each case, we also perform regulation analysis to determine the modes of regulation, and inspect on the relation between the physiological objectives and hierarchical or metabolic regulation. Additionally, we employ FA to investigate putative regulatory links, by extending the corresponding metabolic model with novel interactions and studying the changes obtained in the feasible enzyme space. We end with a discussion of our results and an outlook on further applications and possible extensions of feasibility-based approaches in systems biology.

Results and Discussion
Biological systems constantly adapt to their environment and regulate their metabolism for optimal performance. In this paper, we study this regulation at a system level and use feasibility analysis (FA), considering physiological constraints and a list of potential objectives. We first describe these constraints and objectives and then apply FA to analyze two illustrative cases, a toy model and a model describing the glycolysis in yeast. Figure 1 illustrates our overall approach. FA is inspired by the constraint-based approach used in FBA where an initial flux space is delimited by thermodynamic, mass balance and capacity constraints and the model is then optimized for a certain predefined objective to find the operational point or subspace (panel fig:Feasibility-FBA). Central to our FA approach, we incorporate a detailed kinetic model, taking mechanistic interactions between the enzymes, metabolites and rates quantitatively into account. The multi-dimensional space composed by enzyme levels e, which we call enzyme space, is mapped to the physiological space (containing fluxes J and metabolites x) by the parametrized kinetic model.

Feasibility Analysis
We start by considering a large range of enzyme levels as the initial enzyme space. To construct the feasible enzyme space, this set should further be constrained. However, direct measures that can be applied as constraints are generally available for the physiological space only. In theory, since the enzyme space is mapped to the physiological space with the kinetic model, constraints in one space can be translated into the other by simply inverting the kinetic model. Yet, this inversion is generally not possible in practice due to the non-linear nature of the system. To solve this, similar to [19], we use Monte Carlo (MC) sampling. For each MC sample (a point in the enzyme space) the kinetic model is simulated until it reaches a steady state, yielding the corresponding point in the physiological space.
We first apply the ''hard constraints'' (thermodynamic, mass balance etc) on the physiological space and, via the kinetic model, on the enzyme space. These constraints yield the viable enzyme and physiological space. Next, we evaluate each feasibility criterion for each of the viable physiological states. The labels ''feasible'' or ''infeasible'' are thus assigned to each state and the feasible space is constructed (panel fig:Feasibility-Feasible). Finally, hierarchical regulation analysis can be applied to inspect where metabolism is regulated mainly hierarchically or metabolically, allowing to study the relation between physiological objectives and type of regulation (see Methods for more details).

Constraints
The first step in FA is the application of the hard constraints. We take thermodynamic, stability, kinetic, capacity and total protein constraints into account (See Methods for a more formal definition of these constraints). We start by demanding that every biochemical reaction should obey thermodynamic laws. In general, this is formulated as discrete irreversibility constraints for fluxes through reactions operating far from equilibrium. When measurements on metabolites are available, thermodynamic properties such as Gibbs free energy can be calculated [25], which can be further used as continous constraints. Next, we consider stability, requiring that for each sampled point in enzyme space, the resulting model should be stable. As an approximation, this can a priori be computed by calculating the eigenvalues of the jacobian of the system at a selected steady state and requiring that all should have negative real parts. Then, owing to the available kinetic model, we take kinetic constraints into account. The relation between an enzyme, the metabolites and the rate for any reaction is constrained by its kinetic law. When extracellular fluxes are known and the entire network is considered, if any two of either enzyme, independent metabolite or intracellular flux levels are fixed, the third can be deduced using this set of laws for each enzyme in the network. Next, the capacity constraints provide upper limits for fluxes. Lastly, we assume that the cell economizes the change in total enzyme levels, so that when adapting to a new environment, the total enzyme level is kept within limited range. We also note that, though the total enzyme level is constrained, individual enzyme levels can vary independently within the allowed range.

Objectives
FA continues by further constraining the viable space to obtain the feasible space, by considering a number of quantitiative, physiologically relevant objectives. Three feasibility objectives related to function, robustness, homeostasis and temporal responsiveness are proposed (for formal definitions of each of the criteria, see Methods).
Function. Biological systems have evolved to function optimally in a given environment. Within FBA, this optimal function is considered to be a flux towards a pathway, usually the growth rate; yet alternative optimality criteria such as minimization of uptake rate or redox potential provide adequate prediction of flux distribution [9]. Generalizing this, we consider a set of enzymes as functionally feasible if that set yields optimal (or near optimal) flux for a selected pathway. We also note that in FA, the total enzyme levels are constrained when maximizing flux, considering therefore the cells as optimal strategists for the use of resources, from a costbenefit point of view [26,27].
Robustness and Homeostasis. Robustness and homeostasis are two fundamental characteristics of biological systems [28, and references therein] and has long been recognized and studied from many aspects [29][30][31][32][33]. Following the definition in [28], we consider robustness as a property of systems that maintain their function under perturbations and uncertainty, and homeostasis as maintaining the state via coordinated physiological processes. Despite detailed qualitative descriptions [34] and ad hoc defined metrics, (e.g. [32,35]), a general measure to quantify robustness in metabolism is lacking. To adress this, we first concretely define state and function for a given metabolic network as metabolite levels and the flux for a selected pathway in that network, respectively. We then consider the changes in enzyme levels as perturbations. To quantify robustness and homeostasis, we propose to use the metrics defined within the framework of MCA, namely coresponse coefficients (see methods). Where MCA's control coefficients quantify relative change in one variable (state or function) upon change in another variable (perturbation), coresponse coefficients measure the ratio of relative change in two different variables (state and function) in a network, resulting from a change in a third variable (perturbation). This coefficient is especially interesting to measure robustness and homeostasis of the network, since all three entities can be in different parts of the network.
We consider a set of enzyme levels to be feasible with respect to robustness if the function is maintained (or changes marginally) upon a change of level of any of the enzymes in this set. In that case, metabolite levels are expected to change, resulting in a small co-response coefficient for robustness (D e O J x D%1). Similarly, we consider a set of enzyme levels to be homeostatically feasible, if the state is maintained (or changes marginally) upon a change of enzyme levels of any of this set, resulting in a small overall coresponse coefficient for homeostasis ( P D e O x J D%1, note the swap of indices for x and J) for a series of metabolites located on a pathway, taking into account the global coordination in the network.
Temporal responsiveness. Temporal responsiveness reflects how quickly the network responds to perturbations or external stimuli. It is based on the dynamic characteristics of (a subpart of) the system, such as the response time. From an evolutionary perspective, it is likely that certain pathways or cell types are selected based on their fast (or slow) response to changes in their environment. The key importance of dynamic properties for the cell to adapt to external stimuli has been exemplified for metabolic [36] and signaling networks [37,38]. We consider a set of enzyme levels to be feasible with respect to temporal responsiveness, if it results in a small turn-over time for a metabolite of interest.
Illustration on a small network Initially, to get insight in the shape and properties of the feasible enzyme space, we focused on a small model illustrated in Figure 2. We sampled 2 : 10 4 enzyme level triplets (e 1 ,e 2 ,e 3 ), relative to their reference values, uniformly distributed in 0ƒe i =e 0 i ƒ2. We then simulated the model to find the physiological space (the flux J and metabolite levels x 1 ,x 2 at steady state) corresponding to each triplet of enzyme levels. We then applied all constraints and finally evaluated each feasibility objective.
Constraints. For this small problem, the kinetic expressions allow to explicitly express metabolite levels as a function of enzyme levels. Starting by assuming linlog kinetics for each reaction yields: 1z ln Considering steady state mass balance, (v 1~v2 ,v 2~v3 ) and substituting the values for J 0 1,2,3 and rearranging yields: where, x=x 0 and e=e 0 are the metabolite and enzyme levels relative to their reference state. Eq. 2 describes an explicit model (metabolite concentrations as functions of enzyme levels); fluxes can be obtained by substituting Eq. 2 into Eq. 1. To construct the feasible enzyme space, we start with the thermodynamic constraint and require the steady state flux and metabolite levels to be positive. The constraints on metabolites can analytically derived from Eq. 1, and are represented in Figure 3: where e is the base of the natural logarithm. For the toy problem, all sampled enzyme level sets yielded physiological states that obey the thermodynamic and stability constraints. Finally, we constrain the total enzyme level to change by not more than 50% with respect to the reference state, noting that individual enzyme levels are allowed to vary freely within this constraint (i.e. T enz in Eq. 5 equals 0.5). By constraining the sum of all enzyme levels, around 85% of the sampled enzyme triplets remained viable. Feasibility objectives. After applying the constraints, we analyzed the remaining viable space for each feasibility criterion. As a first step, we did not use any cut-off value (e.g. T flux ,T t ) to discriminate a selected state as feasible or not; rather we visualized feasibility by assigning a color to each state according to a specific criterion (e.g. v J max ,t, e O x J ; see Methods). Function. We first consider function feasibility, by coloring each state according to the flux criterion defined in Eq. 6a. The resulting enzyme and physiological spaces are given in Figure 4(A). An immediate observation is that the flux increases as all three enzyme levels increase simultaneously (red points fall around the line e 1 =e 0 1~e 2 =e 0 2~e 3 =e 0 3 , the main diagonal). The red colored points in the right plot show the enzyme levels that allow the network to achieve roughly the top 25% of possible fluxes. Note that FA takes the cost of the enzyme into account while evaluating the flux objective; for this problem, all enzymes are at equal cost as the optimum lies around the main diagonal. Since the total enzyme level is constrained, the flux is bounded and exhibits an optimal point (indicated by black square on the 3D plot). Furthermore, by taking metabolite levels into account, FA illustrates the effect of metabolic regulation. That is, in physiological space, metabolite level x 1 changes only over a limited range (increasing around 7-fold), while x 2 can increase up to 500 fold without affecting the flux, since x 2 has no inhibitory effect on any of the rates.
Homeostasis and temporal responsiveness. Next, we analyze the homeostasis and temporal responsiveness feasibility. For homeostasis, the resulting enzyme and physiological space is presented in Figure 4(B), where the blue points in the right plot represent the enzyme levels that are homeostatically feasible, i.e. homeostasis can only be maintained if the enzymes in the network assume levels in the blue area of the enzyme space. Notably, to maintain homeostasis, all enzymes should change in concert, i.e. the blue points lie around the main diagonal (e 1 =e 0 1~e 2 =e 0 2~e 3 =e 0 3 ) in the enzyme space. Given that metabolite levels change only marginally, while the flux levels do vary, the changes in flux are mainly attributed to changes in the enzyme levels. In order for the metabolite levels to remain unchanged while the flux is increasing, all enzyme levels should increase synchronously. Comparing Figures 4(A) and 4(B) we observe an interesting trade-off between homeostasis and function. Decreasing all enzymes simultaneously is homeostatically feasible, yet functionally not (the flux decreases). Similarly, increasing all enzymes simultaneously is homeostatically feasible, yet functionally not feasible (production of the enzymes would be too costly).
For temporal responsiveness (Figure 4(C)), we find that the effect of e 1 is small compared to that of e 2 or e 3 : when either of these latter two is low enough, x 1 increases, therefore t increases (red points). Similarly, a decrease in e 3 triggers the accumulation of x 2 , which in turn increases t 2 (not shown). This indicates that temporal responsiveness of this metabolic network is regulated by the last enzyme in this pathway, i.e. that the network has a ''brake'' at the end-point.
Combining feasibility objectives. We next investigated how the three objectives can be combined. For this, we first set a cut-off value for each criterion, as opposed to scanning the entire space as performed in the previous section. The results are given in Figure 5, showing the objective space ( Figure 5(A)), the combined feasible enzyme space ( Figure 5(B)) and a number of 2D-slices at different levels of e 3 ( Figure 5(C)). In the objective space, black points represent a very small subset of the feasible states satisfying all three objectives: high levels of e 1 , e 2 and e 3 . Interestingly, low levels of e 1 , e 2 and e 3 are homeostatically feasible, yet these enzyme levels result in a low flux, therefore functionally not feasible ( Figure 5(B)). These states are especially interesting if a cell economizes on total enzyme levels. For the trade-offs, the optimal combination of objectives depends on the experimental context (see also the section ''illustration on yeast glycolysis''). Another observation from Figure 5(C) is that only high levels of e 2 , the enzyme that consumes x 1 , are feasible in terms of temporal responsiveness.
Next, we performed regulation analysis for this system and analyzed its relation with FA. To calculate the regulation coefficients (r m and r h ) for each sampled point in the enzyme space, we considered the transition from the reference state to the perturbed state (sampled point) and made use of Eq. 7. We find that exclusively hierarchically controlled states (r h *1) are homeostatically feasible. This is expected, since from the FA point of view, the homeostasis feasibility requires that the perturbation results in minimal changes in metabolite levels, and from the regulation analysis point of view the rate of a hierarchically  regulated enzyme is exclusively affected by the level of that enzyme (therefore the metabolite levels do not change). Equivalently, exclusively metabolically controlled states (r m *1) are also feasible with respect to robustness.

Regulation
for feasibility: Feedback inhibition economically maintains homeostasis. In order to assess the effect of a given regulatory mechanism (e.g. end-product feedback inhibition), we modified the initial network and inspected the changes in the objectives and feasible enzyme space. We added a feedback inhibition of x 2 on v 1 , a regulatory mechanism ubiquitous in metabolic reaction networks (the dashed line from x 2 to v 1 in Figure 2). We explored the model by changing the value of e v1 x2 from zero to 20.5 (mild inhibition), up to 25 (strong inhibition).
The effect of this additional feedback inhibition on homeostasis feasibility is presented in Figure 6. It results in a decreased range of x 2 , yet an increased range of x 1 (Figure 6(B)). For the function feasibility, to have the same flux, higher e 3 and lower e 1 levels are needed with increasing feedback inhibition. For combining homeostasis and function feasibility, more states are feasible as inhibition strength increases (the feasible volume increases by 2.5 fold as e v1 x2 changes from 0 to 25, w.r.t. initial model). With increasing feedback strength, e 3 becomes more and more hierarchically regulated, in line with the previous result on combining regulation analysis with homeostasis feasibility. Similar observations, relating the effect of adding regulatory links in a metabolic network to the network sensitivity to perturbations, are reported in [39,40]. The authors illustrated, using a frequency domain approach, that introducing feedback inhibition reduces the effect of perturbations on the output, but additionally showed that extreme feedback inhibition makes the system more sensitive to perturbations.
We also considered a possible feedforward activation of v 3 by x 1 with various strengths (dashed line in Figure 2), and its effect on the feasible enzyme space. This activation further increases the control of e 1 on the pathway flux (the flux control coefficients at the reference state are calculated as C J ei~½ 0:7 0:2 0:1 for e 1 , e 2 and e 3 respectively at e v3 x1~2 .), as v 1 produces x 1 which in turn activates e 3 . A further increase in e v3 x1 results in fewer homeostatically feasible states, illustrating that there is an optimal level of feedforward activation that maximizes the volume of homeostatically feasible space (data not shown).

Illustration on glycolysis in yeast
Next, we applied FA to study yeast glycolysis under evolutionary pressure, especially focusing on trade-offs between alternative objectives. We used a model describing glycolysis in yeast [24] and analyzed two scenarios: the adaptation of the yeast cells during long-term chemostat cultivation under carbon limitation and the so-called ''danger of turbo design''.
Feasibility analysis of prolonged chemostat cultivation of yeast. We first consider the scenario, where yeast cells were grown in a carbon limited chemostat for *250 generations, resulting in a number of changes in their morphology as well as in metabolite and enzyme levels reported in [41][42][43]. Using FA, we explore the three objectives, find corresponding enzyme levels and compare these with the experimental measurements from [43].
In MC sampling the enzyme levels, we appended fermentor balances to the model in [24] and the extracellular metabolites were allowed to change freely. To construct the initial space, we randomly perturbed each enzyme in the network and monitored all resulting 17 fluxes and 13 metabolites. To obtain the viable enzyme space, each stable state was recorded and lastly all feasibility criteria were calculated for each state to construct the feasible enzyme and physiological spaces. We plotted the data from [43] on top of these spaces to inspect the actual changes in the enzyme levels.
We first evaluated the hypothesis that the cells, under constant carbon influx, would economize the enzyme levels while coping with the constant carbon flux, as proposed in [41]. This hypothesis successfully predicts the enzyme levels for PGI and ALD (Figure 7(A), blue points towards to lower left corner having decreased cost). However, it fails to predict the change in enzyme level for the glucose transporter GLT and HK (Figure 7(B)). The levels of these two enzymes increase over the course of the experiment.
To explain this increase, we consider the homeostasis objective, and check the co-response coefficient of extracellular glucose and uptake flux for both enzymes (Figure 7(C)). To take the competitive advantage into account, we drop the absolute values in Eq. 6c. Cells operating in the upper right part of this plot have a competitive advantage for extracellular glucose, since these leave decreased residual glucose levels in the fermentor. Overall, we conclude that cells, being under limited substrate carbon conditions for a long time, increase the levels of those enzymes to compete for the available glucose in the environment.
To illustrate the advantage of considering the trade-off between enzyme economy and competitive ability, we designed a synthetic competition experiment. Four organisms differing by their enzyme levels are grown in a carbon limited chemostat, and the time course for each organism during this competition is simulated. The organisms are (1) wild-type, (2) only considering enzyme economy, (3) only considering competitive ability and (4) considering the trade-off between these two. The enzyme levels, relative to wild type and the time course of each organism are given in Figure 8. We find that the organism that considers the trade-off takes over the entire population in time.
Interestingly, we also predict that the evolved state is prompter for ATP response, i.e. in the evolved strain, ATP responds quicker to perturbations (Figure 7(D)). This is further confirmed by a glucose perturbation experiment (Figure 7(E), data taken from [41]). An important observation follows for PFK: based on both the function feasibility and the ATP temporal responsiveness feasibility, the level of PFK should decrease. However, if there is any decrease in the level of this enzyme, the cell can not survive in the chemostat, meaning that the cells are already at the edge of their feasible enzyme space for PFK (Figure 7(F)). Overall, the cells evolve to a state where they balance the competition for extracellular transport and getting rid of unused overcapacity. This transition makes the cells ''specialists'' in a specific condition, at the expense of loosing the capability of buffering large changes in the environment. We expect that FA will contribute to our understanding of trade-offs and the resulting evolutionary trajectories.
Feasibility analysis of alternative metabolic redesign in yeast glycolysis. To demonstrate how FA can serve to study metabolic (re)design, we consider the so-called turbo design in glycolysis. The turbo design is a general strategy followed by many catabolic pathways, consisting of first activating a substrate in a reaction that requires ATP, after which further metabolism yields a surplus of ATP [44]. In glycolysis, 2 ATP is initially invested in reactions catalysed by HK and PFK while 4 ATP are gained from reactions catalysed by PGK and PYK. The ''danger'' of this design is that when there is excess glucose, the upper part of glycolysis may run at a very fast rate that the lower part can not cope with.
This can lead to accumulation of hexoses in the upper glycolysis (G6P, F6P, F16P), even though ATP and ADP are in steady state, resulting in substrate accelerated cell death [44].
To illustrate the case, we consider the scenario where the cells are in glucose-rich conditions and inspect the homeostasis criterion of hexoses and ethanol flux (J ADH ). Figure 9(B) shows that in this initial design (Figure 9(A)) high levels of both GLT and HK (simulating a large load of substrate), are infeasible, as metabolite levels do not reach steady state. To resolve this handicap of the turbo design, we add a metabolite T6P and two reactions (tps1 and tps2) to the trehalose producing branch. We change the kinetic expression for HK, in line with [45], such that T6P inhibits HK via a feedback inhibition (Figure 9(C), see Methods for the new rate equation). The newly added reactions towards the trehalose pathway follow linear kinetics, and parameters are chosen to keep the metabolite and flux levels the same as the reference state (see caption, Figure 9). All other parameters remain the same as in [24]. A range of enzyme states that were previously infeasible become feasible with the new design (Figure 9(D)). When there is a large push of glucose, T6P acts as a ''brake'' to the glucose uptake, so that neither of the hexoses can increase uncontrollably.
Overall, our FA illustrates how a given metabolic design can be understood within the context of cellular objectives. An interesting observation on cellular trade-offs is that to overcome the danger of the turbo design, the cells have two options: increasing the capacity of reactions consuming the substrate (e.g. storage branches), or introducing T6P inhibition of HK. The first option is costly for the cell since the capacities of all enzymes in the storage pathway have to be increased. The second option is economical and homeostatically feasible, as already illustrated with the FA on the toy model. We finally speculate that evolution pushes cells to acquire this inhibition, in order to adapt to conditions where glucose levels significantly change. Note that such a ''brake'' system is not present for less-favorable carbon sources (e.g. maltose, [46]), an excess of which still results in substrate accelerated cell death.
Two main observations on FA follow from this illustration on the yeast glycolysis. First, for the prolonged chemostat scenario, although there is no prior fitting of model parameters to the experimental data, there is a remarkable quantitative correspondence between the enzyme data from long-term chemostat experiments and the prediction from FA using a dynamic model from literature. This clearly shows that FA can be used to explore the model, to evaluate alternative metabolic strategies, and to hypothesize about cellular trade-offs. Finally, mapping the experimental data clearly reveals which of the objectives is actually selected.

Sampling based methods
We use a Monte Carlo sampling based method to construct the feasible enzyme space. Such a sampling based approach is, intuitive, unbiased and useful as illustrated by the examples. Similar methods are frequently used for exploring biological features, to build model families [18,47], modeling of uncertainties in biochemical networks [19], robustness analysis [48], or designing synthetic networks [49]. Exploring the feasible space by sampling allows to study the trade-offs, and suboptimal behavior, frequently observed feature in biological systems [50][51][52]. Despite its advantages, sampling methods generally suffers from a number of limitations, mainly that they require lots of samples to cover the entire space, tend to waste too much effort and time on regions which are of no real interest and are therefore not scalable to large systems. In the future, our approach and the quantitative measures for feasibility can be combined with a smarter sequential sampling scheme, e.g. [53][54][55] to efficiently explore the initial space for a feasible sub-space.

Feasibility analysis to study regulation at system level
Traditionally, regulation of metabolic networks is studied either by choosing a cellular objective for a genome scale model and optimizing the flux for that objective (top-down, FBA approach), or by constructing a kinetic model with detailed molecular interactions (bottom-up) and applying theorems of MCA or BST. Our method aims to combine elements of the two approaches and allows to study objectives other than flux. Furthermore, as we propose quantitative objectives for feasibility, Figure 8. The competition experiment, to illustrate the optimal enzyme distribution considering the trade-off between enzyme economy and competitive ability for extracellular glucose. The radar plot on the left represents the enzyme levels, relative to wild-type, and the plot on the left represents the competition of each subpopulation with a specific enzyme setting as described in the radar plot. The color for each subpopulation is the same in both plots and is described in the legend. doi:10.1371/journal.pone.0039396.g008 rather than studying ''viable'' or ''lethal'' changes, we can study sub-optimality and trade-offs. By exploring the feasible enzyme space, FA allows evaluating alternative hypotheses and interpreting experimental data. FA assumes the availability of a kinetic model. Despite a long list of challenges, e.g. a high degree of non-linearity, lack of sufficient experimental data, coexistence of multiple time scales etc. [16,56], currently available information on the kinetics of individual enzymes [57] as well as the list of available reliable kinetic models is increasing [58][59][60]. Our example with the detailed kinetic model of glycolysis illustrates how FA can be applied to realistic problems. This is urgently needed, given the growing accumulation of experimental data obtained from different omics layers of cells. Tools to analyze such data are of high value.
Previous computational efforts to understand regulation of metabolism include searching for design principles using optimization principles [23], exhaustively searching and identifying enzyme-based motifs while seeking adaptive properties in a library of network topologies [61], designing synthetic networks for specific tasks [49] or the use of constraints in kinetic parameters to constrain the solution space in steady state models [62]. In particular, the approach taken by Sorribas and co-workers is similar to our FA approach, in that they also investigate feasible enzyme activity patterns leading to cellular (adaptive) responses [63]. Their analysis efficiently finds a global optimum for a given objective, under a given list of physiological constraints. Their mathematically involved approach is tightly coupled to the GMA formulation, elegantly exploits the mathematical structure of the non-convexities of the model. This coupling, in turn, limits its general use. Our proposed FA differs from [63] and [23] in two points. First, it is more general and can be used with any system that can be simulated with a model. Second, central to FA, we propose and use generic quantitative measures for cellular objectives, aiming to eliminate ad hoc definitions. This allows us to consider objectives other than thresholds on fluxes or concentrations, such as robustness/homeostasis and temporal responsiveness.

Conclusions
In this paper, we addressed the following question: being under constraints and evolutionary pressure, why and how is metabolism regulated? To answer this question, we took a top-down approach and speculated that cells, under physiological constraints, are regulated to optimize (one or a combination of) a number of objectives and can hence only assume enzyme levels falling in a socalled feasible enzyme space. We further analyzed how metabolism should be designed from a feasibility perspective, i.e. we addressed the question what are the necessary kinetic interactions in order for cells to attain an objective. Unique to our approach, we proposed quantitative metrics to measure proposed cellular objectives.
One of the fundamental characteristics of biological systems, homeostasis, requires globally coordinated regulation of enzyme levels. An interesting observation for homeostatically feasible states is that these fall in two distinct sub-regimes: a low-flux regime, where all enzymes are downregulated, and are less costly for the cell; and a high-flux regime, where all enzymes are upregulated, therefore costly for the cell. The actual regime chosen by the cell is defined with respect to the environment. In the prolonged chemostat scenario, the cell optimizes the enzyme levels for function, since the carbon influx is externally kept constant. From an enzyme budget point of view, the ubiquitously present feedback inhibition is an economical way to ensure homeostasis. This is especially important for keeping metabolite levels within limits upon a wide range of fluctuations in the environment.
In contrast to homeostasis, maintaining robustness requires a local metabolic effect, meaning that the function can still be maintained by locally adjusting the metabolite levels around specific perturbed enzymes. In line with our findings, Sauer and co-workers recently showed in yeast that alterations in enzyme capacity are buffered by converse changes in substrate metabolite concentration, thereby minimizing the difference in metabolic flux caused by the alteration [52]. In this work, we took homeostasis or robustness as objectives so that we could also study sub-optimal states and the trade-offs between various objectives. This is in contrast to previous attempts where homeostasis has been considered as constraint for the metabolic design problem [15].
Temporal responsiveness reflects a dynamic property of the system. We speculate that this objective is especially applicable to networks whose dynamic properties are of evolutionary importance, e.g. ultrasensitivity, response time etc. As an example, for signaling pathways the effect of network structure on dynamic properties has already been discussed [37,38,64]. Note, that our approach can equally well be used for any other kinetic model, although the physiological objectives may need to be customized. The objective functions we have formulated in this study are illustrations of a more general approach: it may as well be that other objectives turn out to be more relevant under different conditions. It should also be noted that, here, we proposed three ''container'' objectives that are physiologically relevant, which need to be further specified depending on the case evaluated. Additional quantifiable objectives such as overcapacity (which may be defined as the ratio of actual flux to the maximum possible flux) can easily be considered as well.
Taken together, we see that FA quantitatively evaluates alternative hypotheses, shows trade-offs between the available objectives and provides an intuitive platform to integrate the proteome information (enzyme space) with information on metabolome and fluxome (physiological space). Such an integrative approach is indispensable to analyse and interpret the increasingly available multi-omics data on regulation of metabolic networks especially when considering optimal performance or adaptation in response to external stimuli. We illustrated quantitatively via FA that there is a very limited set of enzyme set that are feasible for all the considered objectives. Similar to [51], we argued that the cells are often faced with trade-offs between alternative strategies. Furthermore, by fully exploring the initial viable space and quantitatively evaluating physiological objectives, we got insights on how the metabolic systems are designed (e.g. the ''brake'' for temporal responsiveness objective). This aspect is similar to the ''design space for biochemical systems'' concept in [65,66], but has the additional benefit of direct use of the physiological objectives, making the link from genotype to phenotype more intuitive.

Feasible enzyme space
To construct the feasible enzyme space, we first quantitatively formulate the constraints and the objectives for physiological states. Then we calculate the range of theoretically possible physiological states, and call this the viable enzyme space. We then select a feasible subspace based on the pre-defined criteria and analyze the properties of this subspace. Overall, we construct the feasible enzyme space E f for enzyme levels e as where C is the set of viable states considering the constraints, W represents the set of feasible physiological states (metabolites and fluxes (x,v)) considering the list of physiologically relevant cellular objectives, and v~v(e,x,p) is the rate of the reaction catalyzed by the enzyme e as a function of the enzyme level e, metabolite level x and set of parameters p. The constraints and objectives are detailed below. Some of the metrics for these constraints and objectives are defined with respect to a so-called ''reference state'', denoted with superscript '' 0 ''. This way, all defined entities can be measured with respect to this reference state, much like the elasticity parameters or control coefficients in MCA literature. Conventionally, the reference state can be chosen as the steady state that cells achieve when they are grown under constant, substrate limited conditions and is usually characterized by intracellular fluxes, metabolites and enzyme levels.
The feasible enzyme space constructed is, in fact, a sampling of a multidimensional space containing enzymes, metabolites and fluxes. We visualize this space by slices, i.e. 2D cross-sections. A GUI written in Matlab, and supporting functions as well as the datasets mentioned in this paper can be downloaded at: http:// bioinformatics.tudelft.nl/. The interface takes as input a dataset, calculates the feasibility criteria for selected objectives of the cell and visualizes by (selected) slices.

Constraints: Thermodynamic
constraints. Thermodynamic constraints are formulated as irreversibility constraints for fluxes through reactions operating away from equilibrium: C td~f x,vDv irr e,x,p ð Þw0g ð 4aÞ When more quantitative information is available, for example when Gibbs free energy of a reaction is known, this can also be taken into account as Note, that in a kinetic model in which the equilibrium constant is incorporated in the rate law (e.g. implicitly through the Haldanerelationship), this constraint would already be taken into account. Constraint on total enzyme level. We constrain the total enzyme level to change in a limited range, while individual enzymes can freely be interconverted: where T enz is a precision parameter (e.g. 0.1), and e 0 i is the enzyme level for reaction i at a reference state (denoted by '' 0 ''). This criterion demands that the total enzyme level stays nearly constant (e.g. can change only within 10%, when T enz = 0.1).
Cellular objectives: Function. We define states in which near-optimal flux under constrained enzyme levels is obtained as feasible: where T flux is a cut-off value for feasibility in terms of optimal flux. This criterion demands that a flux can be at most 10% (when T flux~0 :9) away of it's possible maximal flux (denoted as J max ). Robustness and Homeostasis. We consider robustness as a property that allows a system to maintain its function under perturbations and homeostasis as the coordinated physiological processes which maintain the current state [28]. In this work, the perturbations are changes in enzyme levels, states are metabolite levels and function is the flux towards a selected pathway or reaction. In order to quantify both robustness or homeostasis, we need a measure between state (metabolite levels), function (flux towards the selected enzyme/pathway) and perturbation (changes in enzyme levels). We use the co-response coefficients ( ei O yj y k ) as a measure, defined within the context of Metabolic Control Analysis (MCA) [67] as: resulting co-response coefficient should therefore be large: Second, homeostasis is considered. A state is called feasible if upon enzyme perturbation, metabolite levels do not change significantly , whereas the flux does. We formulate the feasibility related to homeostasis for a set of M metabolites as: The summation over the metabolites ensures that homeostasis is not only local in one metabolite but over a number of relevant metabolites, e.g. belonging to a pathway. Temporal responsiveness. Temporal responsiveness of metabolite levels in a metabolic network in response to perturbations is defined using the turn-over time of metabolites W temporalresponsiveness~x ,vDtvT t f g ð6dÞ with t~x 0 J 0 , x 0 and J 0 the physiological parameters at the steady state reached after a perturbation and T t is a treshold. This criterion demands that the turnover time of a metabolite should be smaller than e.g. 0.5 per unit time, when T t = 0.5.

Regulation analysis
An important question is to what extent metabolic fluxes are regulated by gene expression or by metabolic regulation. In line with the convention in regulation analysis [4], ''metabolic'' regulation is defined as change caused by concentrations of substrate(s), product(s) and modifier(s). ''Hierarchical'' changes are those caused by change in enzyme concentration, via alterations in mRNA sequestration and intracellular localization and/or rates of transcription, translation or degradation. Both types of regulation are quantified by the hierarchical and metabolic regulation coefficients (r h and r m ) defined as We consider cases where Dr h Dv0:1 as exclusively metabolically regulated, and cases where Dr m Dv0:1 as exclusively hierarchically regulated. One immediate application of information from regulation analysis is in metabolic engineering. In the first case, increasing the flux could be achieved by simply increasing the enzyme level whereas for the second case, alternative engineering strategies such as protein engineering to change the kinetic properties of the enzyme need to be considered.

Illustrative cases
Toy model. To illustrate feasibility analysis, we use a small example model with tractable kinetics, where a substrate S is converted into a product P via a linear pathway of 3 reactions and 2 intracellular metabolites ( Figure 2). The model assumes a steady state and all three rates follow linlog kinetics, allowing to calculate an explicit steady state solution for metabolites and rates in terms of enzyme levels and kinetic parameters [12]. In linlog kinetics, the rate of reaction i (v i ) is described relative to the steady state flux J 0 i as a function of the enzyme levels e i and intracellular and extracellular metabolites (x j and c k ) all relative to their steady state The reference steady state conditions (X 0 ,J 0 ) and the elasticity matrix E x composed of kinetic parameters (e) are given in Figure 2.
Glycolysis model in yeast. To study feasibility analysis applied on a real problem, we used a previously published model of glycolysis in Saccharomyces cerevisiae [24]. The kinetic expressions for each reaction and the parameters are the same as [24], and the reference state used in our work is given in Table 1. For feasibility analysis, we need to sample the enzyme levels, relative to their reference state and in the yeast model, this is performed by sampling the relative V max 'es since: For the feasibility analysis of alternative metabolic redesign in yeast glycolysis, the new kinetic expression for the HK reaction is: Lastly, in order to illustrate the competition in the fermentor, we added a growth equation for this and expressed the growth rate with simple monod-growth kinetics as: with m 0 is the growth rate at the reference conditions (and is equal to the dilution rate in the chemostat) and K g~0 :098mM, the extracellular glucose level at the reference conditions. The term between the parantheses represent the effect of total enzyme cost on growth.

Data pre-processing
In using experimental data, there were multiple measurements for a specific time point. Since the data available was insufficient to assume and fit a parametric model, we used non-parametric Gaussian kernel regression (s~20hr) to estimate the average at a specific point, taking all data into account.