Abstract
Background:
Several approaches, including metabolic control analysis (MCA), flux balance analysis (FBA), correlation metric construction (CMC), and biochemical circuit theory (BCT), have been developed for the quantitative analysis of complex biochemical networks. Here, we present a comprehensive theory of linear analysis for nonequilibrium steadystate (NESS) biochemical reaction networks that unites these disparate approaches in a common mathematical framework and thermodynamic basis.
Results:
In this theory a number of relationships between key matrices are introduced: the matrix A obtained in the standard, lineardynamicstability analysis of the steadystate can be decomposed as A = SR^{T }where R and S are directly related to the elasticitycoefficient matrix for the fluxes and chemical potentials in MCA, respectively; the controlcoefficients for the fluxes and chemical potentials can be written in terms of R^{T }BS and S^{T }BS respectively where matrix B is the inverse of A; the matrix S is precisely the stoichiometric matrix in FBA; and the matrix e^{At }plays a central role in CMC.
Conclusion:
One key finding that emerges from this analysis is that the wellknown summation theorems in MCA take different forms depending on whether metabolic steadystate is maintained by flux injection or concentration clamping. We demonstrate that if ratelimiting steps exist in a biochemical pathway, they are the steps with smallest biochemical conductances and largest flux controlcoefficients. We hypothesize that biochemical networks for cellular signaling have a different strategy for minimizing energy waste and being efficient than do biochemical networks for biosynthesis. We also discuss the intimate relationship between MCA and biochemical systems analysis (BSA).
Background
Developing methodologies for quantitative analysis, both experimental and mathematical, of complex biochemical networks has become one of the central themes of postgenomic biochemistry and mathematical biology. Several disparate approaches, including metabolic control analysis (MCA) [1,2], flux balance analysis (FBA) [3,4], and correlation metric construction (CMC) [5,6], share many commonalities. The objective of this work is to provide a unifying mathematical framework and a thermodynamic basis for these approaches. The thermodynamic basis is the nonequilibrium steadystate (NESS) theory [7,8] originally developed to describe macromolecular, freeenergy transduction [913], and the mathematical methods are based on linear analysis near a NESS.
The physical concept of a NESS (also known simply as a steadystate in the MCA community) applies to that of a driven system, such as a living organism, with a constant flow of matter (transforming nutrients to waste) through the system and a corresponding dissipation of free energy [9]. Even though a NESS looks remarkably similar to a thermodynamic equilibrium, in that chemical concentrations reach stationary values, the concentrations are maintained constant in a NESS by balancing influxes and effluxes, rather than by balancing the forward and backward fluxes of each elementary reaction as in thermodynamic equilibrium. To an observer concerned only with the chemical concentrations, a NESS would seem to be a true equilibrium; but, in fact, it represents a pseudoequilibrium, where the work done to drive the system appears as heat, and deserves further consideration [9].
We distinguish between two types of linear analysis: (i) linear stability analysis [14], and (ii) control analysis (or sensitivity analysis) that addresses how the steadystate shifts in response to certain perturbations to the system. These perturbations may be deterministic – leading to MCA – or stochastic – leading to CMC. As expected, the matrix A obtained in linear stability analysis, the flux controlcoefficient matrix C in MCA, the stoichiometric matrix S of FBA, and the correlation matrix R from CMC are intimately related. We develop the relationships in the context of complex biochemical networks and their NESS thermodynamics. We show that all the key matrices in MCA and CMC can be computed if one knows the fluxes, chemical potential differences, and the stoichiometric matrix.
The theoretical concept of a NESS provides valuable mathematical and thermodynamic properties. Although many of the approaches mentioned above have provided a mathematical treatment of NESS, many have never explicitly taken the related thermodynamics into account. The key insight is to consider concentrations as measures of potentials, providing a relationship between potential differences (voltages), fluxes (currents), and resistances. It becomes immediately clear that the resistances are related to reaction effectors, such as the expression levels of enzymes, that affect both the forward and backward fluxes of a given reaction, but do not change the chemical potential difference. It is this thermodynamic perspective of a NESS that can be used to develop a comprehensive theory that unites these approaches.
In the context of the newly developed biochemical circuit theory (BCT) [15,16], the flux J of each metabolic reaction in a NESS metabolic network is decomposed into J = φ^{+ } φ^{ }and the chemical potential difference of the reaction in the NESS is Δμ = k_{B}T ln(φ_{}/φ_{+}), where k_{B }is Boltzmann's constant and T is absolute temperature. Therefore, knowing the chemical potential difference and the flux of a reaction gives
Furthermore, if the standardstate free energy of reaction Δμ^{o }is known from equilibrium thermodynamics, then determines the ratio of the metabolites' concentrations, products to reactants, in NESS.
In many theoretical approaches, metabolic networks and signaling pathways are often separated into disparate mechanisms with only peripheral interactions. This is far from reality, however, as these mechanisms are intimately intertwined and often indistinguishable. For example, in wholebody metabolism, glucose and fatty acids serve as the fuel for cellular metabolism, but also as signals for insulin secretion in pancreatic βcells and ketogenesis in hepatocytes, where cellular organelles, such as mitochondria, play a key role in the handshake between metabolism and signaling. Ultimately, our goal is to understand how these cellular systems behave in response to stress and disease [17]. It is in applying the unified mathematical framework that includes the thermodynamic perspective to these processes that we believe the most significant contributions will be made.
In application, the concepts presented here allows one who has the stoichiometric matrix, S, and measures the correlation matrix, R, to obtain the control coefficients, and vice versa (Table 1). Furthermore, by incorporating the thermodynamic aspect, one is able to take full advantage of the thermodynamic properties, such as standard free energies of formation, which are typically more complete, more consistent, and more well analyzed than the kinetic information [17]. This is advantageous in a field where the availability of kinetic parameters is a limiting factor. However, available data remains incomplete, and missing thermodynamic and kinetic data must continually be obtained [17].
Table 1. Summary and comparison table listing the relationship between the elasticity and control coefficient matrices and the stoichiometric (S) and correlation (R) matrices, where A = SR^{T }and B = A^{1}
Basic kinetic equations for biochemical networks
Let us consider a network of M biochemical reactions involving N biochemical species X_{i }(i = 1, 2, ..., N). The jth biochemical reaction (j = 1, 2, ..., M) is characterized by a set of stoichiometric coefficients and such that
Some of the integer ν's and κ's can be zero. The N × M matrix S = κ  ν is known as the stoichiometric matrix. Eq. (2) assumes that the forward and backward reaction kinetics are characterized by the constants and . In cases where the kinetic scheme involves intermediate steps (e.g., enzyme binding), each of the intermediate steps can be incorporated as reactions in the form of Eq. (2), requiring additional kinetic parameters, and the stoichiometric matrix can be constructed to include elementary reactions representing actual interaction events. Segel [18], for example, details how the MichaelisMenten mechanism is expressed in terms of Eq. (2).
For the reaction system (2), kinetic equations of the timedependent concentration changes for each biochemical species can be written according to the law of mass action as [19,20]
where we use x_{ℓ }to denote the concentration of respective X_{ℓ}, and (i = 1, 2, ..., N) are external injection fluxes to species i (also known as boundary fluxes). We use the convention that subscripts index species and superscripts index reactions. To further simplify the notation, we introduce forward and backward fluxes [20]:
Then Eq. (3) becomes dx/dt = SJ + J^{e}, in which x and J^{e }are Ndimensional column vectors, and J is an Mdimensional column vector.
In a closed reaction system, J^{e }= 0. In this case, it can be shown that thermodynamic equilibrium is the only positive stationary solution to Eq. (3): the internal fluxes J = φ_{+ } φ_{ }= 0 for each and every reaction. For a closed biochemical reaction system, since all the fluxes are necessarily zero in its unique equilibrium (for each given set of parameters; unique in the dynamic sense), all the control coefficients are necessarily zero. Therefore, in a NESS, at least one injection flux and one efflux are nonzero ( ≠ 0), or certain concentrations x_{i }are held at constant levels. The former case is referred to as "external flux injection," while the latter is referred to as "external concentration clamping." A distinction between these two "forcing" mechanisms greatly clarifies many controversial issues concerning NESS.
Steadystate concentrations
In the equilibrium state of a closed reaction system, φ_{+ }= φ_{ }and the ratio of chemical concentrations is independent of the amount of material and the initial condition:
giving rise to the concept of chemical equilibrium constants. Such a system of biochemical reactions is nondissipative and is associated with a number of conserved quantities. An essential mathematical characteristic of the closed system is that its equilibrium point is degenerate (neutrally stable on a certain manifold), i.e., there is no unique stationary solution for different initial conditions.
For an open system that approaches a NESS, the stationary solution is not usually degenerate. If we assume {} is an asymptotically stable NESS, small differences in the initial condition should lead to the same NESS. We use the notations to denote the set of concentrations of internal species and to denote the concentrations of external species which are clamped or independently varied.
In the linear regime near the NESS [17],
A = {A_{ij}} is called the linear stability matrix. We introduce the N × M matrix R = , such that A = SR^{T}. We also denote A^{1 }= B.
If the nonlinear dynamics scheme of Eq. (3) conserves certain quantities, such as total element, motif, or enzyme concentrations, S does not have full rank and A is singular [17,19,20]. In such cases, it is possible to transform A into a nonsingular matrix by replacing its linearly dependent rows with vectors in its left null space (see Ref. [18] of [20]). By doing so, one removes the redundancies from the original dynamics scheme. A more detailed mathematical analysis will be published elsewhere. Both here and below, A is understood to be the transformed nonsingular matrix.
Results and discussion
Metabolic control analysis
Eq. (3) is a general scheme for biochemical reaction networks. For metabolic reactions involving enzymes, the enzyme and the enzymesubstrate complexes are treated as additional "species." In this section, however, we study MCA, and assume that every reaction involves an enzyme which is not explicitly expressed as a species. Rather, we assume enzyme activity is absorbed into the rate constants for the forward and backward reactions. This is clearly not an accurate approximation for many enzymatic reaction mechanisms. Nevertheless, it provides a firstorder approximation for treating a metabolic network. More complex rate laws for enzymatic reactions have been discussed for biochemical systems analysis [21,22].
MCA focuses on how the NESS {} shifts in response to a perturbation to the amount of enzyme for a reaction or a perturbation to its substrate (∈ ) concentration. Without losing generality, we assume these perturbations are made to the enzyme for the M th reaction or the N th (external) species.
Elasticity coefficients
First, we consider the case where the concentration of external species N is changed: . Before reaching a new NESS, the local response is only in the flux of the reactions involving X_{N}. The immediate, local change is characterized by
which is called the scaled, local elasticitycoefficient [1]. These coefficients are elements of the scaled, local, elasticitycoefficient matrix, ϵ = (dg J)^{1 }R^{T}. Here, we use the same notation used by Heinrich and Schuster [23], where (dg v) denotes the diagonal matrix containing the components of the vector v along its diagonal. The unscaled, local, elasticitycoefficient matrix is given by ϵ' = (dg J) ϵ (dg x)^{1 }= R^{T }(dg x)^{1}. The prime symbol is used to distinguish unscaled coefficients from scaled coefficients, both here and in what is to follow.
The coefficient ϵ should be set apart from the coefficient ε, which characterizes the steadystate response [24]. When a new NESS is established, the concentrations {x_{i}} (i = 1, 2, ..., N  1) satisfy
Solving Eq. (6), we have where B_{iN }is the N th column vector of the matrix A^{1}. The new NESS established near {} is {}, and the new fluxes, J^{m }+ δ J^{m}, are given by
Hence, the scaled, steadystate elasticitycoefficient is [17]
which may also be known as a response coefficient [23,25]. These coefficients are elements of the scaled, steadystate, elasticitycoefficient matrix, ε = (dg J)^{1 }R^{T }B (dg B)^{1}. Here, (dg B) denotes the diagonal matrix containing only the diagonal components of the matrix B along its diagonal. The unscaled, steadystate, elasticitycoefficient matrix is given by ε' = (dg J) ε(dg x)^{1 }= R^{T }B (dg B)^{1 }(dg x)^{1}.
As we shall show, some connectivity theorems for the two quantities ϵ and ε are different. In principle, the ε can be experimentally measured as a system response, but obtaining ϵ can be more difficult since it is a transient response that must be measured individually.
Control coefficients
Next, we consider the case where the enzyme for the M th reaction E^{M }is changed: E^{M }→ E^{M }+ δ E^{M}.
Since we have assumed that the rate constants for reaction M, , are linearly proportional to the concentration of the enzyme catalyzing reaction M, E^{M}, we have
When E^{M }→ E^{M }+ δ E^{M}, the new NESS satisfies [17]
Solving for δx_{j}/ (j = 1, 2, ..., N), we have
From Eq. (10), we get the scaled, concentration controlcoefficient [26],
which is an element of the scaled, concentration, controlcoefficient matrix given by = BS (dg J). The unscaled, concentration, controlcoefficient matrix is = (dg x) (dg J)^{1 }=  (dg x) BS.
Combining Eqs. (7) and (10), we get the scaled, flux controlcoefficient [1],
where δ_{mM }is the Kronecker delta function. These coefficients are elements of the scaled, flux, controlcoefficient matrix given by C = I  (dg J)^{1 }R^{T }BS (dg J). The unscaled, flux, controlcoefficient matrix is C' = (dg J) C (dg J)^{1 }= I  R^{T }BS.
Summation and connectivity theorems
Among ϵ, ε, C, and , we have the following relationships. First, note the relationships
between the scaled coefficients, and similarly,
between the unscaled coefficients.
Second, from Eqs. (13)–(16), it is straightforward to show that if there is no injection flux, i.e., in NESS SJ = J^{e }= 0, we have
where 1_{i×j }and 0_{i×j }denote the i × j matricies of ones and zeros, respectively. Eqs. (17) and (18) are known as the summation theorems for the scaled controlcoefficients. For the unscaled controlcoefficients, we have the generalized summation theorems, C'K = K and , where K is the M × (M  N) loop (or nullspace) matrix, i.e., its columns form a basis of the nullspace of S so that SK = 0_{N×(M  N) }[15,16].
Euler's theorem on homogeneous functions can be used to understand the significance of the injection fluxes [17,26]. If there are no injection fluxes and a NESS is sustained by clamped concentrations, the steadystate fluxes are homogeneous functions of the enzyme concentrations with order 1. That is to say that, assuming every reaction in the system is catalyzed by an enzyme, if every enzyme concentration is simultaneously doubled, the flux in each reaction doubles. This leads to the summation theorem for the flux controlcoefficients. On the other hand, if there are no clamped concentrations and a NESS is sustained by injection fluxes, the steadystate fluxes are homogeneous functions of enzyme concentrations with order 0. A similar argument exists for steadystate concentrations as a function of enzyme concentrations, leading to the summation theorem for the concentration controlcoefficients.
Third, BSR^{T }= I ⇒ (R^{T }BS)(R^{T }B) = R^{T }B, which yields the connectivity theorems for the scaled coefficients:
and the generalized connectivity theorems for the unscaled coefficients [23]:
Finally, by combining the generalized summation and connectivity theorems, we get
which is the central equation in MCA [23]. The matrix (K ϵ') has been proved to be invertible [23,27].
Therefore, if the local elasticitycoefficients are known, the controlcoefficients can be calculated.
Using the steadystate elasticitycoefficients, we have
The matrix (K ε') is also invertible (see Methods). Therefore, if the steadystate elasticitycoefficients are known, the controlcoefficients can be calculated. This is an important result because ε can be experimentally measured as a system response, but obtaining ϵ can be more difficult since it is a transient response and must be measured individually.
Biochemical potential and its control analysis
The most notable difference between BCT and other theories of biochemical networks is its sound thermodynamics basis [15,16]. BCT articulates the fluxes J and chemical potential differences Δμ as two, complementary, key characteristics of a reaction in a NESS. Kirchhoff's flux and potential laws are manifestations of the fundamental laws of physical chemistry, namely conservation of mass and conservation of energy [15,16].
Both J and Δμ can be written in terms of the forward and backward fluxes φ_{+ }and φ_{ }as given in Eq. (1), and vice versa. In general, the flux J and the chemical potential difference Δμ are nonlinearly related [28]:
In the linear regime with small flux and chemical potential difference, J and Δμ are linearly related with the biochemical conductance being φ_{+}/k_{B}T (= φ_{}/k_{B}T) [11,16]. In the nonlinear regime, if a perturbation on a reaction is from its upstream, we can assume the φ_{ }is relatively unperturbed. Then we have δJ/δΔμ = φ_{+}/k_{B}T. Similarly, if the reaction is perturbed from its downstream, then δJ/δΔμ = φ_{}/k_{B}T.
Combining BCT and MCA, we can define the scaled, local elasticitycoefficients for chemical potentials [25]
as elements of the matrix ϵ^{Δμ }= k_{B}T (dg Δμ)^{1 }S^{T}, the scaled, steadystate elasticitycoefficients for chemical potentials
as elements of the matrix ε^{Δμ }= k_{B}T (dg Δμ)^{1 }S^{T }B (dg B)^{1}, and the scaled, potential controlcoefficients [25]
as elements of the matrix C^{Δμ }= k_{B}T (dg Δμ)^{1 }S^{T }BS (dg J), where reaction potentials and substrate concentrations in the reference state are used for normalization. The corresponding unscaled coefficients are given in matrix form as
Note the appearance of the matrices S^{T}, S^{T }B, and S^{T }BS in the local and steadystate elasticitycoefficients and the controlcoefficient for the potentials, as compared to the matricies R^{T}, R^{T }B, and R^{T }BS in the coefficients for the fluxes. Furthermore, note the relationships
Eqs. (19) and (29) in combination yield the summation [25] and connectivity theorems for the scaled coefficients, under consideration of the steadystate condition,
Eq. (31) can also be derived using the homogeneous function theory (discussed above with flux and concentration controlcoefficients). For the unscaled coefficients, we get the generalized summation and connectivity theorems by combining Eqs. (21) and (30) so that
Eq. (33) is a manifestation of Kirchhoff's loop law [15,16].
By combining the generalized summation and connectivity theorems, we get
The matrix (K ϵ') has been proved to be invertible [23,27]. Therefore, if the local elasticitycoefficients are known, the controlcoefficients can be calculated. Using the steadystate elasticitycoefficients, we have
The matrix (K ε') has an inverse (see Methods). Therefore, if the steadystate elasticitycoefficients are known, the controlcoefficients can be calculated.
Correlation metric construction
Correlation metric construction (CMC) focuses on how the NESS {} fluctuates in response to a random stochastic perturbation to the concentration level of an external species x_{N }∈ . We assume the perturbation is small. Hence, the linear analysis is appropriate.
If the random perturbation has a sufficiently long correlation time compared to the relaxations of all the reactions in the system, then we are dealing with a shift in the NESS where, according to Eq. (6),
From the time series of {δx_{i}} with an uncorrelated perturbation on δx_{N}, we have
for i, j = 1, 2, ..., N.
CMC extracts information from reaction systems by using perturbations with insufficiently long correlation times. δx_{N }is a Markov process characterized by
where k_{d }controls the relaxation time of the perturbation, a is its random amplitude with Gaussian distribution, and ξ(t) is a uncorrelated noise. δx_{N }in Eq. (38) is called an OrnsteinUhlenbeck (OU) process with autocovariance function
With this type of random perturbation, the temporal correlations between the concentration fluctuations give information on the connectivity of a network following the kinetic equation (4)
for i = 1, 2, ..., (N  1). More precisely (A. Arkin, personal communication), instead of Eq. (40) CMC follows an autoregressive (AR) process, the discrete version of an OU process, with discrete time T, 2T, ...,
in which . The discrete AR and continuous OU models give essentially the same result.
Let's rewrite the A matrix into a (N  1) × (N  1) matrix . Then the stationary solution to Eqs. (38) and (40) is
which leads to
Two limits of Eq. (42) are particularly interesting. If the correlation time of g(τ) in Eq. (39) is ≫ all the relaxation times of , then we have Eq. (37). On the other hand, if the correlation time for the random perturbation is ≪ all the relaxation times of , we have
The matrix is the fundamental solution to the linear kinetic equation (4). R_{ij }(t) is the response curve, as a function of t, for the ith species to an impulse of the jth species at time 0. The relation between the correlation matrix and the fundamental solution is expected for a linear system. For directly connected species i and j, the , whereas if i and j are not directly connected, . Analytic solutions for completely linear systems have been derived and analyzed [29].
Conclusion
Rate limiting step: largest flux controlcoefficient and smallest biochemical conductance
Both flux controlcoefficients and biochemical conductances can be used as indicators for ratelimiting steps in a biochemical network [30]. Traditionally the ratelimiting step is understood in terms of the highest "activation barrier" in the pathway. Here we use a simple example to demonstrate that these concepts are intimately related. Mathematically, the highest activation barrier is also related to the mean firstpassage time in stochastic dynamics [31,32].
Let's consider the linear pathway
under NESS with clamped concentrations x_{0 }and x_{n }for X_{0 }and X_{n}, respectively. Analogous to a continuous energy landscape, we have the discrete energy function
for j = 1, 2, ..., n. If there is a ratelimiting step, say for example that reaction j is nearly irreversible, then E_{j }is maximal and that is where the highest energy barrier is located. We now show that it is also where the largest flux controlcoefficient and the smallest conductance are.
The NESS flux of reaction system (44) can be solved for analytically to yield [31,32]
Therefore, the flux controlcoefficient is
which reaches its maximum at maximal E_{j}. The biochemical conductance, on the other hand, is
which is at its minimum when (or ) reaches its minimum. For small flux, (and ) reaches its minimum also at maximal E_{j }where concentration x_{ℓ }is the lowest.
Note that with J given, the smallest conductance also corresponds to the largest Δμ. In general, the ratelimiting step is not necessarily where the k's are smallest. The most efficient situation for reaction system (44) is when all the conductances are equal. This result is analogous to the constant torque principle suggested for molecular motors [12,13,33].
Fluxcontrolled and concentrationcontrolled biochemical networks
The analysis provided in the present paper poses the following biological question. When is a NESS of a biochemical network controlled by a constant flux injection and when is it controlled by a clamped concentration (chemical potential)? Clearly in real biological systems, the appropriate answer is a combination of both. Just as a real battery, having both finite internal resistance and conductance, in an electrical circuit is a combination of constantcurrent (zero internal conductance) and constantvoltage (zero internal resistance) ideal batteries, so too would a real biological system be a combination of constant flux injections and clamped concentrations [17]. Nevertheless, a brief discussion on this topic provides insights into the control of biochemical networks.
First, let's see what are the respective consequences of these two different types of "forcing" to a biochemical system. For a constant flux injection J to a simple unimolecular reaction in a NESS, we have the heat dissipation rate (hdr)
This indicates that for the same injection flux J a reaction with larger φ_{ }(i.e., conductance) dissipates a smaller amount of heat. For a constant clamped chemical potential Δμ, we have
We note that hdr increases with φ_{ }for the constant potential difference (concentration clamping) case. These results are analogous to an electric circuit in which the heat dissipation equals I^{2 }R and U^{2}/R, respectively.
This observation leads us to the following hypothesis. In a biochemical network for biosynthesis, the fluxes are essential for supply and demand. Hence, such networks will in general have large conductances in order to minimize the energy waste and be efficient. On the other hand, in a biochemical network for cellular signaling, the concentrations are the essential signal. Thus, in this case the conductances are generally small in order to minimize the energy waste.
Relation to the biochemical systems analysis
MCA shares an intimate relationship with the biochemical systems analysis (BSA) originally proposed by M.A. Savageau [21]. A critique of BSA can be found in [34]. BSA assumes that each metabolic reaction is catalyzed by an enzyme with rapid pseudosteadystate enzymesubstrate complexes. This assumption leads to an explicit mathematical expression for the flux J in the enzymatic reaction as a rational function of the concentration of a substrate x [18], such that
where K ≤ L due to substrate saturation or effector inhibition. The coefficients a_{i}, b_{i}, α_{0}, and β_{0 }contain the concentrations of all the other substrates. BSA provides insight to the properties of a_{i }and b_{j}. The network of metabolic reactions is then modeled [22] by a system of enzymatic reactions with a rate law given in the form of (49).
In terms of Eq. (49) the local elasticitycoefficient for the flux J can be written as
which also has the standard form of a Pade approximation [35]. Furthermore, BSA proposed [22] to approximate (49) by a powerlaw expression
which is essentially the integral form of the definition of local elasticitycoefficients.
In contrast, MCA assumes no specific rate laws for an enzymatic reaction except that the rate is linearly proportional to the enzyme concentration. This assumption is certainly valid for most biochemical reactions except when a substrate concentration is significantly less than that of an enzyme.
Methods
Solving the central, controlanalysis equations
To demonstrate that the matrix, (K ε'), that appears in the central, controlanalysis equations, i.e., Eqs. (24) and (36), is invertible, we note that since we are dealing with systems without conserved quantities, rank(S) = N. Furthermore, we note that the matrix K^{T }K is invertible because K contains linearly independent columns, and because all columns of K are in the right nullspace of S, the matrix
is also invertible [23]. Therefore, we have
Abbreviations
BCT Biochemical Circuit Theory; BSA Biochemical Systems Analysis; CMC Correlation Metric Construction; FBA Flux Balance Analysis; MCA Metabolic Control Analysis; NESS Nonequilibrium SteadyState; A N × N linear stability matrix; B N × N matrix = A^{1}; K M × (M  N) loop (or nullspace) matrix (SK = 0); R N × M correlation matrix defined such that A = SR^{T}; S N × M stoichiometric matrix; x Ndimensional vector of species' concentrations; J Mdimensional vector of reaction fluxes; Δμ Mdimensional vector of reaction potentials; ϵ M × N local, elasticitycoefficient matrix; ω M × N steadystate, elasticitycoefficient matrix; N × M concentration, controlcoefficient matrix; C M × M flux, controlcoefficient matrix; C^{Δμ }M × M biochemical potential, controlcoefficient matrix; (dg v) If v is a vector, (dg v) denotes the diagonal matrix with the components of v along its diagonal. If v is a matrix, (dg v) denotes the diagonal matrix with only the diagonal components of v along its diagonal. * Denotes steadystate; ' Distinguishes unscaled coefficients from scaled coefficients
Authors' contributions
HQ initiated this work. All authors contributed to and wrote the manuscript together. All authors read and approved the final manuscript.
Acknowledgements
We thank Professor J.B. Bassingthwaighte for continuous encouragement, and Adam Arkin, Kyung Kim, Shoudan Liang, and Brian Ingalls for helpful discussions. HQ and DAB are partially supported by NIH grant GM068610. WJH is supported by the Intramural Research Program, NIDDK, NIH.
References

Fell DA, Sauro HM: Metabolic control and its analysis.
Eur J Biochem 1985, 148:555561. PubMed Abstract  Publisher Full Text

Kholodenko BN, Westerhoff HV: The macroworld versus the microworld of biochemical regulation and control.
TIBS 1995, 20:5254. PubMed Abstract  Publisher Full Text

Edwards JS, Palsson BØ: The E. Coli MG1655 in silico metabolic geotype: its definition, characteristics, and capabilities.
Proc Natl Acad Sci USA 2000, 97:55285533. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schilling CH, Palsson BØ: The underlying pathway structure of biochemical reaction networks.
Proc Natl Acad Sci USA 1998, 95:41934198. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Arkin A, Ross J: Statistical construction of chemical reaction mechanisms from measured timeseries.
J Phys Chem 1995, 99:970979. Publisher Full Text

Samoilov M, Arkin A, Ross J: On the deduction of chemical reaction pathways from measurements of time series of concentrations.
Chaos 2001, 11:108114. PubMed Abstract  Publisher Full Text

Qian H: Phosphorylation energy hypothesis: open chemical systems and their biological functions.
Ann Rev Phys Chem 2007, 58:113142. Publisher Full Text

Qian H: Opensystem nonequilibrium steadystate: statistical thermodynamics, fluctuations and chemical oscillations.
J Phys Chem B 2006, 110:1506315074. PubMed Abstract  Publisher Full Text

Wyman J: The turning wheel: a study in steady state.
Proc Natl Acad Sci USA 1975, 72:39833987. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hill TL: Free Energy Transduction in Biology: The SteadyState Kinetic and Thermodynamic Formalism. New York: Academic Press; 1977.

Hill TL: Free Energy Transduction and Biochemical Cycle Kinetics. New York: SpringerVerlag; 1989.

Qian H: The Mathematical theory of molecular motor movement and chemomechanical energy transduction.
J Math Chem 2000, 27:219234. Publisher Full Text

Qian H: Cycle kinetics, steadystate thermodynamics and motors – a paradigm for living matter physics.
J Phys: Condens Matter 2005, 17:S3783S3794. Publisher Full Text

Murray JD: Mathematical Biology. 2nd edition. New York: Springer; 1991.

Beard DA, Liang SD, Qian H: Energy balance for analysis of complex metabolic networks.
Biophys J 2002, 83:7986. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Qian H, Beard DA, Liang SD: Stoichiometric network theory for nonequilibrium biochemical systems.
Eur J Biochem 2003, 270:415421. PubMed Abstract  Publisher Full Text

Beard DA, Qian H: Chapter 6: Biochemical Reaction Networks. In Chemical Biophysics: Quantitative Analysis of Cellular Systems. Cambridge, UK: Cambridge University Press; 2008.

Segel IH: Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steadystate Enzyme Systems. New York: John Wiley & Sons; 1975.

Alberty RA: Degrees of freedom in biochemical reaction systems at specific pH and pMg.
J Phys Chem 1992, 96:96149621. Publisher Full Text

Poland D: On the stability of mass action reactions in open systems.
J Chem Phys 1991, 94:44274439. Publisher Full Text

Savageau MA: Biochemical system analysis I.
J Theor Biol 1969, 25:365369. PubMed Abstract  Publisher Full Text

Savageau MA: Biochemical system analysis II.
J Theor Biol 1969, 25:370379. PubMed Abstract  Publisher Full Text

Heinrich R, Schuster S: The Regulation of Cellular Systems. New York: Chapman & Hall; 1996.

Westerhoff HV, Chen YD: How do enzyme activities control metabolite concentrations?
Eur J Biochem 1984, 142:425430. PubMed Abstract  Publisher Full Text

Westerhoff HV, van Dam K: Thermodynamics and control of biological freeenergy transduction. Amsterdam: Elsevier; 1987.

Giersch C: Control analysis and metabolic networks.
Eur J Biochem 1988, 174:509519. PubMed Abstract  Publisher Full Text

Reder C: Metabolic control theory: a structural approach.
J Theor Biol 1988, 135:175201. PubMed Abstract  Publisher Full Text

Beard DA, Qian H: Relationship between thermodynamic driving force and oneway fluxes in reversible processes.
PLoS One 2007, 2:e144. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Heuett WJ, Qian H: Grand canonical Markov model: a stochastic theory for open nonequilibrium biochemical networks.
J Chem Phys 2006, 124:044110. PubMed Abstract  Publisher Full Text

Fell DA: Increasing the flux in metabolic pathways: a metabolic control analysis perspective.
Biotechnol Bioeng 1998, 58:121124. PubMed Abstract  Publisher Full Text

Derrida B: Velocity and diffusion constant of a periodic onedimensional hopping model.
J Stat Phys 1983, 31:433450. Publisher Full Text

Fisher ME, Kolomeisky AB: The force exerted by a molecular motor.
Proc Natl Acad Sci USA 1999, 96:65976602. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Oster G, Wang HY: Reverse engineering a protein: the mechanochemistry of ATP synthase.
Biochim Biophys Acta 2000, 1458:482510. PubMed Abstract  Publisher Full Text

Beard DA, Qian H, Bassingthwaighte JB: Stoichiometric foundation of largescale biochemical system analysis. In Modelling in Molecular Biology. Edited by Ciobanu G, Rozenberg G. New York: Springer; 2004:119.
[Natural Computing Series].

Bender CM, Orszag SA: Advanced Mathematical Methods for Scientists and Engineers. New York: SpringerVerlag; 1999.