Thermodynamic constraints on the regulation of metabolic fluxes

Nutrition and metabolism are fundamental to cellular function. Metabolic activity (i.e. rates of flow, most commonly referred to as flux) is constrained by thermodynamics and regulated by the activity of enzymes. The general principles that relate biological and physical variables to metabolic control are incompletely understood. Using metabolic control analysis and computer simulations in several models of simplified metabolic pathways, we derive analytical expressions that define relationships between thermodynamics, enzyme activity, and flux control. The relationships are further analyzed in a mathematical model of glycolysis as an example of a complex biochemical pathway. We show that metabolic pathways that are very far from equilibrium are controlled by the activity of upstream enzymes. However, in general, regulation of metabolic fluxes by an enzyme has a more adaptable pattern, which relies more on distribution of free energy among reaction steps in the pathway than on the thermodynamic properties of the given enzyme. These findings show how the control of metabolic pathways is shaped by thermodynamic constraints of the given pathway.

Metabolism enables the utilization of nutritional resources to provide free energy, material, and cellular communication for functions of cells (1). Some metabolic pathways such as central carbon metabolism, which processes macronutrients or the major caloric sources in diet (proteins, fats, and carbohydrates), have been largely defined for over 50 years (2), and even genome-scale reconstructions of metabolic networks that include thousands of metabolites and reactions in different intracellular compartments are available in many unicellular organisms and metazoans (3)(4)(5). Despite the complexity of metabolic networks in living organisms, metabolic processes still follow the laws of thermodynamics. It is thus important to understand the principles that link fundamental thermodynamic quantities to the control of metabolic pathway activity (i.e. the rates of flow of materials through the network, most commonly referred to as fluxes).
Analytical frameworks have been developed to understand steady-state and dynamic behaviors of metabolic networks. The most widely applied method in computational modeling of metabolism is flux balance analysis, which assumes that the network is in steady state and that configurations of metabolic fluxes are determined by optimizing an objective function such as growth rate (6). This approach does not require that enzyme properties such as expression levels be known, but the utility of the method depends on the objective function and other assumptions such as nutrient uptake rates. Moreover, flux balance analysis does not provide information about pathway regulation and the control of metabolic flux, which requires additional knowledge such as enzyme activities or metabolite abundances, which are routinely measurable (7,8).
Another framework for understanding metabolism is metabolic control analysis (MCA), 3 originally developed in the 1970s once the biochemistry for many of the key pathways in metabolism such as glycolysis and the tricarboxylic acid cycle was established (9 -15). MCA quantitatively measures how the flow through a metabolic pathway responds to changes in parameters such as the abundance of an enzyme or availability of a nutrient. MCA defines the sensitivity of a metabolic flux to a perturbation in a given metabolic reaction (i.e. the flux control coefficients (FCCs)) and also provides a series of rigorously derived relationships between metabolic fluxes, metabolite concentrations, and enzyme activities. This framework can also be applied in computer simulations (16 -19) or in experimentation (20,21) when some but not all of the relevant variables are measured (22)(23)(24)(25).
Metabolism is subject to the laws of thermodynamics. These laws place constraints on the dynamics of metabolic reactions (26 -32) and are known to affect the control of fluxes in linear pathways (9,31,33,34). Moreover, the deviation from equilibrium at an individual reaction step has been applied as the criterion to identify rate-limiting steps in metabolic pathways (2,29). However, this rule of thumb has been challenged by MCA, at least for linear pathways (9), but a quantitative evaluation of the relationship between thermodynamics and flux control in metabolic pathways with different topologies such as the control at branching points is still lacking.
In this study, we use a set of models with representative topological structures observed in metabolism to investigate the quantitative relationships between thermodynamics and regulation of metabolic fluxes. Notably, all models are based on This work was supported by National Institutes of Health Grants R01 CA193256 and R00CA168997 and American Cancer Society Grant TBE-TBE434120 (to J. W. L.). The authors declare that they have no conflicts of interest with the contents of this article. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This article contains supporting information and Figs. S1 and S2. 1 To whom correspondence may be addressed. E-mail: ziwei.dai@duke.edu. 2 To whom correspondence may be addressed. E-mail: jason.locasale@ duke.edu. some simplifying assumptions, but some (and hopefully these) are useful. To our surprise, we find that, in both linear and branched pathways, the regulation of pathway fluxes by individual enzymes is in general loosely constrained by the deviation from thermodynamic equilibrium, or the thermodynamic driving force, of the whole pathway. Only pathways very far from equilibrium have their flux regulation strictly constrained by the thermodynamic driving force, in which all fluxes are almost fully controlled by the upstream enzymes. The relationships are further studied in a mathematical model of glycolysis as an example of a more complicated metabolic model. These results unravel simple principles of how metabolic pathways are regulated by the interaction between thermodynamics and enzyme activity.

Metabolic flux and thermodynamics
We first consider unimolecular, first-order kinetics of reversible metabolic reactions for simplicity. In this framework, each reaction has one substrate (S) and one product (P) (Fig. 1). It is noteworthy that this simplified case approximates the more complicated Michaelis-Menten mechanism when the abundance of substrate is far below the K m . The forward reaction rate v ϩ ϭ k ϩ S and backward reaction rate v Ϫ ϭ k Ϫ P are linear in substrate and product concentrations. Because the rate constants of the forward and backward reactions are coupled by the equilibrium constant K:K ϭ k ϩ /k Ϫ , the net flux carried by this reaction is as follows: v ϭ v ϩ Ϫ v Ϫ ϭ k ϩ (S Ϫ P/K). Let k ϭ k ϩ for simplicity, and we have the following.
The equilibrium constant K can be further connected to the standard Gibbs free energy of this reaction.
Finally, the Gibbs free energy is determined by the standard Gibbs free energy in combination with concentrations of the substrate and product.
We note that Equation 1.1 can also be rewritten to explicitly incorporate the reaction free energy change. Here, we let g ϭ ⌬G/RT for simplicity. Thus, we have the following.
According to the connectivity theorem in MCA (13), flux control coefficients of reactions directly associated with a metabolite are coupled through local elasticity coefficients of these reactions with respect to the metabolite. For a reaction and a metabolite directly associated with this reaction (either as substrate or as product), the elasticity coefficient is defined as the partial derivative of the reaction rate with respect to the concentration of the metabolite on the logarithmic scale. The elasticity coefficients can be computed from Equations 1.1 and 1.4.

Linear pathways
Based on the definitions under "Metabolic flux and thermodynamics," we now consider a linear pathway consisting of n reactions ( Fig. 2A). Each reaction has unimolecular, linear kinetics, as described previously. The first substrate S in is converted to the end product S out by this reaction chain in n steps with n Ϫ 1 intermediary metabolites S 1 , . . . , S n Ϫ 1 . The ith reaction in this chain has rate constant k i and equilibrium constant K i . The concentrations of S in and S out are treated as fixed parameters, because these two values, together with rate constants and equilibrium constants of the enzymes, determine concentrations of all other metabolites.
Previous work has shown that in linear pathways with either first-order or zero-order (i.e. substrate-saturated) kinetics, the pattern of flux control can be completely determined if free energies of all reactions involved in the pathway are known (34). Here, we use these results as necessary background to study the relationship between thermodynamics and flux control in metabolic networks with different topological structures. Briefly, from MCA (13), there are analytical relationships between thermodynamic properties and the flux control coefficients, which quantify the relative importance of each enzyme in regulating the flux through the pathway.
To apply the summation theorem and connectivity theorems in the theory of metabolic control analysis, we assume that the A, forward and backward fluxes of a reversible reaction. S is the concentration of substrate, P is the concentration of product, v ϩ is the rate of the forward reaction, v Ϫ is the rate of the backward reaction, k ϩ is the rate constant of the forward reaction, and k Ϫ is the rate constant of the backward reaction. B, net flux and Gibbs free energy change of a reversible reaction. v is the net reaction rate, k is the rate constant, K is the equilibrium constant, ⌬G is the reaction Gibbs free energy change, R is the universal gas constant, and T is the temperature. Other variables are the same as in A.

Thermodynamics in control of metabolism
pathway is in steady state, in which all reactions have identical net flux. In other words, all intermediary metabolites have balanced fluxes feeding and consuming them. Thus, we can write n Ϫ 1 equations for the steady-state constraints.
Let J denote the net flux through this pathway. The flux control coefficient, C vi J , quantifies the sensitivity of the pathway flux J to perturbation in activity (i.e. rate constant) of the ith reaction. It is defined as the ratio of relative change in the flux J to relative change in the rate of the ith reaction when an arbitrary parameter p that only affects the rate of the ith reaction and no other reactions has a small change.
If p ϭ k i , we have the following.
The flux control coefficients can be uniquely determined by solving equations derived from the summation and connectivity theorems. For linear pathways, the summation and connectivity theorems do not explicitly include the pathway flux J.

Summation theorem
Connectivity theorem According to Equation 1.5, we have the following. A, diagram of a linear pathway and related parameters. S in is the input substrate, S out is the final product, S i is the ith intermediary metabolite, k i is the rate constant of the ith reaction, K i is the equilibrium constant of the ith reaction, and J is the pathway flux. B, scatter plots comparing flux control coefficients and reaction free energy changes in a linear pathway with randomly sampled parameters. The pathway includes 10 reaction steps. C vi J is the ith flux control coefficient, and g i ϭ ⌬G i /RT quantifies the ith free energy change. C, scatter plots comparing flux control coefficients and the thermodynamic driving force (max{g i }) in a linear pathway with randomly sampled parameters. Variables are the same as in B. D, relationship between free energy of a reaction and the elasticity coefficient of its velocity toward the substrate. E, relationship between free energy of a reaction and the elasticity coefficient of its velocity toward the product. F, relationship between the free energies of two sequential reactions and the ratio of the flux control coefficients associated with the two reactions.

Thermodynamics in control of metabolism (Eq. 2.6)
Flux control coefficients can be solved by combining Equations 2.1, 2.4, 2.5, and 2.6 (see supporting information). They can be written as either explicitly including only the reaction free energy terms or explicitly including only the rate constants and equilibrium constants.
According to Equations 2.7 and 2.8, the flux control coefficients are completely determined by kinetic and thermodynamic parameters of reactions in the pathway. Therefore, alterations in concentrations of the input substrate and the final product do not influence the flux control coefficients.
To illustrate the quantitative relationships defined in Equations 2.7 and 2.8, we consider a linear pathway consisting of 10 reactions; randomly sample 20,000 combinations of the parameters S in , S out , {k i }, and {K i }; and compute the corresponding flux control coefficients. We select the size 10 because it approximates the typical length of linear metabolic pathways. All parameters and substrate concentrations are sampled from a log-normal distribution.
We first correlate the flux control coefficient with reaction free energy change for each individual reaction (Fig. 2B). We find that the relationship between reaction free energy changes and flux control coefficients calls into question the longstanding hypothesis that reactions with the most negative free energy changes, such as the reactions catalyzed by hexokinase, phosphofructokinase, and pyruvate kinase in glycolysis, serve as rate-limiting steps of a pathway (35-38). Although previous work in metabolic control analysis has provided the theoretical framework to study the role of thermodynamics in the regulation of metabolic fluxes (34), this hypothesis is still widely accepted (39 -41). However, according to our analysis, for all reactions except the first reaction, their flux control coefficients correlate poorly with the reaction free energy changes, suggesting that a large absolute value of free energy change is dispensable for a rate-limiting step, and regulation of metabolic fluxes by an individual reaction is determined by both thermodynamics property and position along the pathway.
We next investigate how the flux control coefficients are associated with global thermodynamic properties of the entire pathway. We quantify the deviation from the thermodynamic equilibrium by the reaction free energy change that is closest to zero (i.e. the maximal value of free energy change among all reactions). When this quantity equals zero, the pathway has zero net flux (i.e. thermodynamic equilibrium). Moreover, when this quantity is close to zero, at least part of the pathway is not efficiently driven by thermodynamics, resulting in inefficient usage of enzymes with free energy changes close to zero (34). Thus, it is termed the thermodynamic driving force (TDF).
We correlate the thermodynamic driving force with flux control coefficient of each reaction step (Fig. 2C). From these simulations, a pattern emerges, in which the thermodynamic driving force determines either the upper bound or lower bound of the flux control coefficients; for the flux control coefficient of the first reaction (C vi J ), its lower bound depends on the thermodynamic driving force, whereas for the rest of the reactions, the upper bounds of their flux control coefficients depend on the thermodynamic driving force.
Analytical relationships that constrain flux control coefficients by thermodynamic driving force can also be derived from Equations 2.7 and 2.10 (see supporting information).
From Equations 2.11 and 2.12, we can derive corollary relationships for the flux control coefficients when the system is very far from equilibrium (i.e. thermodynamic driving force approaches minus infinity).
(Eq. 2.13) Equation 2.13 suggests that in a linear pathway far from thermodynamic equilibrium, the metabolic flux through the pathway is fully controlled by the enzyme catalyzing the first step. In other situations when the thermodynamic driving force has a smaller absolute value, the flux control is more evenly distributed among all enzymes in the pathway.
An intuitive interpretation of the relationship between the thermodynamic driving force and flux control coefficients can be given based on the fact that for a reaction with firstorder kinetics, the elasticity coefficients toward the sub-Thermodynamics in control of metabolism strate and product can be completely determined if the free energy of this reaction is known (Fig. 2 (D and E) and Equation 2.6). According to the connectivity theorem (Equation 2.5), in a linear pathway, the ratio of the ith flux control coefficient to the i ϩ 1th flux control coefficient can be written as a function of the free energies of the ith and i ϩ 1th reactions (Fig. 2F).
Thus, the thermodynamic driving force constrains the distribution of flux control among enzymes by constraining the elasticity coefficients that determine flux control coefficients through the connectivity theorem. When the thermodynamic driving force has a large absolute value (i.e. TDF approaches minus infinity; bottom left corner in Fig. 2F), for any reaction in the pathway, its elasticity coefficient with respect to the product is very close to zero, which means that altering the concentration of the product has negligible influence on the rate of this reaction. In this situation, perturbation of any downstream enzyme has a minimal effect on the first reaction, whose rate is equal to the overall flux, because the effect must be propagated through the metabolite S 1 , which serves as the product of the first reaction. On the other hand, when the thermodynamic driving force is close to zero (top right corner in Fig. 2F), the ratio between two sequential flux control coefficients is more flexible, resulting in what we term as a more adaptable pattern of flux control in these pathways.
To summarize, for a linear metabolic pathway with first-order kinetics, we have derived analytical relationships between the flux control coefficients and thermodynamic properties of both individual reactions and the entire pathway. We have shown that the flux control coefficients correlate poorly with free energy changes of the corresponding reactions but that the thermodynamic driving force of the whole pathway places bounds on how much any enzyme can control pathway flux. However, when the pathway is very far away from thermodynamic equilibrium, the metabolic flux through this pathway is completely controlled by the first reaction. We next investigate whether these principles are conserved in metabolic networks with more complex topologies.

Branch point with two upstream fluxes
We next consider a metabolic network with two converging fluxes, J 1 and J 2 , at a branch point (Fig. 3A). Through the two fluxes, the metabolite at the branch point, S BP , is produced by two input substrates, S in,1 and S in,2 . The final product S out is produced from S BP with the flux J 1 ϩ J 2 at the steady state. The three reactions included in this network have rate constants {k 1 , k 2 , k 3 } and equilibrium constants {K 1 , K 2 , K 3 }. At the steady state, the fluxes J 1 and J 2 and the concentration of S BP can be solved from the rate constants, equilibrium constants, and concentrations of the input substrates and end product based on the kinetic rules. We can then determine the flux control coefficients from the summation and the connectivity theorems. Here, we have three fluxes and three reactions in the network, giving nine flux control coefficients. It is worth noting that Equation 3.5 is not sufficient for determining all values of the flux control coefficients because it provides six linear equations for nine unknowns. To uniquely calculate every element of C v J , three additional equations from the connectivity theorem are needed.

Connectivity theorem
In Equation 3.7, ⑀ represents the matrix consisting of the normalized elasticity coefficients.  Thermodynamics in control of metabolism (Eq. 3.9) Furthermore, the reaction free energy changes can also be written as functions of the kinetic parameters and substrate concentrations.
(Eq. 3.10) We then combine Equations 3.9 and 3.10 to express the flux control coefficients in terms of the reaction parameters and free energy terms. To simplify the expressions, we let h i ϭ e gi . Then we have the following.

(Eq. 3.11)
To illustrate the relationships between flux control coefficients and thermodynamic properties of the network, according to Equations 3.9-3.11, we simulate the flux control coefficients and reaction free energy changes based on 20,000 combinations of parameters randomly sampled from a log-normal distribution as described previously in Equation 2.9. Combinations of parameters resulting in negative fluxes are dis-carded. The distributions of flux control coefficients (Fig. 3B) show that C v2 J1 , C v1 J2 , C v3 J1 , and C v3 J2 can have absolute values much larger than 1 under certain parameter combinations, indicating that the two upstream fluxes can be dramatically altered in response to perturbation of activities of enzymes in other branches. Based on the simulation, we also compare flux control coefficients with the reaction free energy changes (Fig. 3C) and the thermodynamic driving force (Fig. 3D). Consistent with the case of the linear pathway, most of the flux control coefficients correlate poorly with free energy changes of the individual reactions catalyzed by the corresponding enzyme (Fig. 3C) but exhibit a dependence on the thermodynamic driving force (Fig. 3D). Moreover, the quantitative relationships characterizing how flux control coefficients are constrained by the thermodynamic driving force can also be analytically derived (see the supporting information).

(Eq. 3.12)
From Equation 3.12, we can derive the limits of all the flux control coefficients except for C v1 J1 ϩJ2 and C v2 J1 ϩJ2 when the thermodynamic driving force goes to minus infinity. Thus, when the pathway is far from equilibrium, the two upstream fluxes J 1 and J 2 are fully regulated by activities of the Thermodynamics in control of metabolism enzymes that generate them. Moreover, the downstream flux J 1 ϩ J 2 is also fully controlled by the two upstream enzymes. The inequalities in Equation 3.12 also indicate that the feasible regions of the flux control coefficients C v2 J1 , C v1 J2 , C v3 J1 , and C v3 J2 become semi-infinite (i.e. only bounded in one direction) when the thermodynamic driving force approaches zero (i.e. there exists at least one near-equilibrium reaction in this pathway).
To summarize, in a metabolic pathway with a branch point involving two converging fluxes, regulation of the fluxes is largely dependent on how close to thermodynamic equilibrium the pathway in entirety is. When the pathway is far away from equilibrium, all three fluxes are completely controlled by the two upstream enzymes. On the other hand, when there exists at least one nearequilibrium reaction in the pathway, there is more flexibility in the regulation of fluxes, and, unlike in the case of linear pathways, the two upstream fluxes are able to be dramatically altered even by enzymes that are not directly involved in the reaction.

Branch point with two downstream fluxes
We next consider a pathway with one branch point in which one upstream flux, J 1 , diverges to two downstream fluxes, J 2 and J 3 , and show that the principles of metabolic flux regulation and thermodynamics also hold in this case. This network includes one input substrate, S in , and two final products, S out,1 and S out,2 , which are all connected to the branch point metabolite S BP . All kinetic rules and assumptions are as described above. Thus, at the steady state, concentration of S BP can be solved from the reaction parameters and concentrations of other metabolites in the network.  Here, we also have in total nine flux control coefficients that can be solved based on the summation theorem and connectivity theorem. To investigate the relationships between flux control coefficients and thermodynamic variables, we randomly sample 20,000 combinations of reaction parameters and substrate concentrations and calculate flux control coefficients and free energy changes corresponding to the randomly sampled parameters from Equations 4.3 and 4.6. Distributions of the sampled flux control coefficients suggest that the two downstream fluxes are able to be regulated in response to variations in the activity of enzymes not directly carrying them (i.e. C v1 J2 , C v3 J2 , C v1 J3 , and C v2 J3 Thermodynamics in control of metabolism can have absolute values larger than 1; Fig. 4B). Again, we observe a poor correlation between flux control coefficients and free energy change of the corresponding reaction steps (Fig. 4C). Among the nine flux control coefficients, only C v1 J2 ϩ J3 exhibits a clear dependence on the free energy change g 1 .
Despite the poor correlation between flux control coefficients and free energy changes of individual reaction steps, seven of the nine flux control coefficients are strictly con-strained by the thermodynamic driving force (Fig. 4D). This corroborates the findings in other types of network topologies that the regulation of metabolic fluxes is constrained by the deviation of the whole pathway from thermodynamic equilibrium. We also analytically derive the quantitative relationships between the upper and lower bounds of the flux control coefficients and the thermodynamic driving force of the pathway (see the supporting information).  Thus, based on a combination of analysis and simulation, we have demonstrated that, at branch points, the regulation of fluxes by enzyme activities is constrained by the thermodynamic driving force but has little correlation with free energy changes of the individual enzymes. These findings, together with similar results in other network structures, suggest that the influences of thermodynamics on regulation of metabolic fluxes by enzyme activities primarily occur at the pathway level rather than at the individual reaction level. These principles are applicable to simple metabolic network models regardless of the structure of the metabolic network, but they also rely on the assumption that all enzymes follow first-order kinetics. We next examine whether the conclusions are conserved in metabolic networks with Michaelis-Menten kinetics.

Thermodynamics and flux control in Michaelis-Menten models
We now consider zero-order kinetic models as an extreme case of substrate-saturable enzymes and repeat the theoretical analysis and simulation in these models to characterize the relationships between thermodynamics and flux control. Briefly, although the expressions of steady-state concentrations of the intermediary metabolites and steady-state fluxes are much more complicated with this assumption, the flux control coefficients can still be written in the form of simple functions of the reaction parameters and free energy changes. However, the thermodynamic driving force no longer constrains the flux control coefficients in this case (Figs. S1 and S2). Most of the flux control coefficients in all three network structures have the form of a zero-order homogeneous function of the thermodynamic terms h i . Consequently, moving all reactions away from or toward the thermodynamic equilibrium without influencing the ratios between the h i values does not change the resulting flux control coefficients.
Because flux control is no longer constrained by the thermodynamic driving force in zero-order kinetic models, we hypothesize that in Michaelis-Menten models, which lie between the extreme cases of first-order and zero-order kinetic models, the saturation of enzymes by metabolites determines whether the flux control coefficients are still constrained by the thermodynamic driving force as they are in first-order kinetic models. In Michaelis-Menten models (28), the rate of a reaction depends on the catalytic constant (i.e. turnover number) k and the equilibrium constant K of the enzyme, the concentrations of its substrate S and product P, and the Michaelis-Menten constants K S and K P for the substrate and product, respectively.
(Eq. 5.1) We note that for Michaelis-Menten models, the elasticity coefficients toward the substrate and product concentrations differ from those in first-order kinetic models only by one term quantifying saturation of the enzyme by the given metabolite (i.e. substrate or product).

Thermodynamics in control of metabolism
For Michaelis-Menten models, because the rate law is more complicated than that of zero-order or first-order kinetic models, we use sets of parameters randomly generated, as previously done for first-order kinetic models to simulate the thermodynamic variables and flux control coefficients, and then compare them with those in first-order kinetic models. To quantify the saturation of enzymes, we introduce the saturation term for a reaction (Fig. 5A). This saturation term considers the total saturation of an enzyme by the substrate and the product. For a pathway, the overall level of enzyme saturation is quantified by the maximal saturation term (Fig. 5A).
To evaluate whether the thermodynamic variables and flux control coefficients in a Michaelis-Menten model are consistent with the relationships we have derived for first-order kinetic models, we compute the deviation of the flux control coefficients in the Michaelis-Menten model from those allowed in the first-order kinetic model with the same topology and thermodynamic driving force (Fig. 5B).
We first examine how the maximal saturation term affects the fraction of Michaelis-Menten models consistent with firstorder kinetic models. In all three pathway structures, we observe a gradual increase in the fraction of Michaelis-Menten models consistent with the first-order kinetic models (Fig. 5C) as the maximal saturation term decreases, suggesting that the flux control coefficients are largely constrained by the thermodynamic driving force when the enzymes in the pathway are not highly saturated. We next compare the maximal saturation term with the deviation from the results in the first-order kinetic model. We find that the deviation gradually decreases as the maximal saturation term decreases, corroborating the findings that the relationships between thermodynamics and flux control are conserved in Michaelis-Menten models if the enzymes are not highly saturated.
To summarize, here we show that the principles for thermodynamics and flux control in metabolic pathways with firstorder kinetics are still effective in Michaelis-Menten models if the enzymes are not highly saturated. It is now reasonable to ask whether the same rules apply to more complicated metabolic networks in which most enzymes not only follow the Michaelis-Menten kinetics, but also have their activities regulated by long-range interactions such as allosteric feedback.

Thermodynamics in control of metabolism Thermodynamics and flux control in a mathematical model of glycolysis
Thus, we considered a model of glycolysis to study the relationships between thermodynamics and flux control (42). We have previously constructed a mathematical model of glycolysis and used this model to study the regulation of the lactate production flux (17,19). The topological structure of this model can be approximated by a linear pathway consisting of 13 reversible reactions (Fig. 6A). Elements and interactions such as reactions with more than one substrate, allosteric feedbacks, conserved species (e.g. NAD ϩ /NADH, ATP/ADP/AMP), and irreversible reactions are also included.
We generated 50,000 sets of random parameters to simulate the glycolysis model and computed the flux control coefficients and thermodynamic driving force at the steady state for each set of parameters. We also computed the maximal saturation term and assessed how the saturation of enzymes influences the relationships between thermodynamics and flux control. We first compared the thermodynamic driving force and flux control coefficients in the glycolysis model with the values allowed in the model for a linear pathway with the same number of reactions (i.e. 13 reactions) according to Equations 2.11 and 2.12 (Fig. 6B). As expected, we find that for parameter sets yielding a high (i.e. close to 1; red dots in Fig. 6B) maximal saturation term, the glycolysis model is inconsistent with the relationships derived under the assumption of first-order kinetics. We then estimated the deviation from first-order kinetic models (Fig.  6C) and computed the fraction of models consistent with firstorder kinetic models (Fig. 6D). Consistent with the Michaelis-Menten model, we found that as the deviation decreases, the fraction of consistent models increases, suggesting that the relationships derived in first-order kinetic models with simple topological structures are largely conserved in more complicated metabolic networks.

Discussion
In this study, we derive quantitative relationships between thermodynamics, enzyme activity, and regulation of metabolic fluxes. For a set of example pathways, we calculate the flux control coefficients as functions of the enzyme rate constants, reaction equilibrium constants, and Gibbs free energy. Based on numerical simulation and exact analytical results, we find that in all network topologies considered, the flux control coefficients are bounded by the thermodynamic driving force, as defined by the deviation of the entire pathway from thermodynamic equilibrium. Moreover, distinct patterns of flux regulation emerge when the thermodynamic driving force approaches two extreme values: if the thermodynamic driving force is very negative (i.e. all reactions are far away from equilibrium), enzymes catalyzing the upstream reactions (i.e. reac-

Thermodynamics in control of metabolism
tions directly consume the input substrates) exert full control of all fluxes; on the other hand, if the thermodynamic driving force is close to zero (i.e. near-equilibrium reactions exist in the pathway), there is more flexibility in the pattern of flux control, and fluxes in network topologies with branch points can be dramatically altered by enzymes not directly carrying the fluxes. These findings suggest that the coupling between thermodynamics and regulation of metabolic flux occurs at the systemic level and challenges the rule of thumb that the reaction with the most negative free energy change serves as the rate-limiting step.
It is also worth noting that all analysis that we have done to derive analytical expressions here relies on the assumption of first-order kinetics. This approximates the more generalized Michaelis-Menten mechanism when the abundance of substrate is far below the K m of the enzyme. Accordingly, in a study comparing substrate concentrations and K m values in different cell types, among all substrate-enzyme interactions investigated, around 70% exhibited higher substrate concentration than the K m value, thus requiring a Michaelis-Menten or substrate-saturated, zero-order kinetics (30). For Michaelis-Menten kinetics and a glycolysis model, we have shown, that the generalizations of our findings depend on the levels of substrate saturation. We also note that the tendency of an enzyme to be operating in the linear region or to be saturated by its substrate is highly dependent on the substrate used by this enzyme. For instance, kinases use ATP as the substrate whose physiological concentration is much higher than the typical K m value of enzymes and thus are more likely to be substrate-saturated. On the other hand, other enzymes, such as chromatinmodifying enzymes, use metabolites much less abundant as substrates (43). Nevertheless, because the relationships between thermodynamics and flux control are conserved under circumstances where the enzymes are moderately saturated, we expect that our findings are applicable to more complicated biochemical processes with Michaelis-Menten kinetics.
Finally, although the network structures that we study here are much simpler than real metabolic networks, the conclusions that we derive here can still be extended to more complicated network topologies. A metabolic network without cycles can always be simplified to a set of branch points connected by linear reaction chains with varying lengths, and the enzymes in the same linear reaction chain can be treated as an entirety to calculate the flux control coefficients with respect to simultaneous change in these enzymes (44). In this case, each linear reaction chain is simplified to a single reaction step, which enables us to apply the principles we found in the simple network structures with a branch point.
To summarize, we characterize the quantitative relationship between thermodynamics and regulation of metabolic fluxes by enzymes in a metabolic network in this study. We find that the global thermodynamic property of the network (i.e. thermodynamic driving force) constrains almost all flux control coefficients in both linear and branched pathways. Specifically, fluxes in metabolic networks far away from thermodynamic equilibrium are almost fully controlled by enzymes catalyzing the upstream reactions (i.e. reactions directly consuming the input substrates). On the other hand, near-equilibrium metabolic networks allow more flexibility in the pattern of regulation.
Semi-infinite feasible regions of flux control coefficients are only allowed in branched pathways with near-equilibrium reactions. These findings highlight the importance of global thermodynamic features in constraining the pattern of regulation of metabolism.

Experimental procedures
Analytical expressions of flux control coefficients were calculated using the function LinearSolve[] in Wolfram Mathematica version 11. For the computer simulations, 20,000 combinations of random parameters were generated using the function normrnd() in MATLAB R2016b for each model. Flux control coefficients and reaction free energies were then computed using the parameter sets that generate positive flux configurations. For the glycolysis model, a model file in SBML format was downloaded from the BioModels Database (https:// www.ebi.ac.uk/biomodels-main/) 4 (45) and translated to Cϩϩ codes using the SBML translator module in the Systems Biology Workbench (http://sbw.sourceforge.net/) 4 (46). Simulation of the model was done with the ODE solver gsl_odeiv2_ step_msbdf in the GNU Scientific Library (https://www. gnu.org/software/gsl/). 4 50,000 random sets of parameters were generated by sampling V max values of the enzymes with Latin hypercube sampling. The flux control coefficient of the glycolytic flux (defined as the rate of glucose consumption) with respect to a reaction step was estimated using the finite difference approximation, in which J glycolysis is the glycolytic flux and V max,i is the maximal velocity of the ith reaction. Free energy changes for the reactions were calculated from standard reaction free energies and steady-state concentrations of the metabolites. All source codes are available at the GitHub page of the Locasale laboratory (https://github.com/LocasaleLab/MCA_thermodynamics). 4