Linear analysis near a steady-state of biochemical networks: Control analysis, correlation metrics and circuit theory

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 steady-state (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, linear-dynamic-stability analysis of the steady-state can be decomposed as A = SRT where R and S are directly related to the elasticity-coefficient matrix for the fluxes and chemical potentials in MCA, respectively; the control-coefficients for the fluxes and chemical potentials can be written in terms of RTBS and STBS respectively where matrix B is the inverse of A; the matrix S is precisely the stoichiometric matrix in FBA; and the matrix eAt plays a central role in CMC. Conclusion One key finding that emerges from this analysis is that the well-known summation theorems in MCA take different forms depending on whether metabolic steady-state is maintained by flux injection or concentration clamping. We demonstrate that if rate-limiting steps exist in a biochemical pathway, they are the steps with smallest biochemical conductances and largest flux control-coefficients. 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).

The physical concept of a NESS (also known simply as a steady-state 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 effuxes, 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 pseudo-equilibrium, 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 steady-state 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 control-coefficient 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 standard-state 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 whole-body 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]. (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 Michaelis-Menten 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]: 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.

Steady-state 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 non-dissipative and is associated with a number of conserved quantities. An essential mathematical characteristic of the closed system is that its equilibrium point is degenerate biochemical potential steady-state elasticity (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], 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 enzyme-substrate 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 first-order 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 Mth reaction or the Nth (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 elasticity-coefficient [1]. These coefficients are elements of the scaled, local, elastic- 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, elas- 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 steady-state response [24].
Hence, the scaled, steady-state elasticity-coefficient is [17] which may also be known as a response coefficient [23,25]. These coefficients are elements of the scaled, Here, (dg B) denotes the diagonal matrix containing only the diagonal components of the matrix B along its diagonal. The unscaled, steady-state, elasticity- 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. 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, control-

Control coefficients
Combining Eqs. (7) and (10), we get the scaled, flux control-coefficient [1], where δ mM is the Kronecker delta function. These coefficients are elements of the scaled, flux, control-coefficient matrix given by C = I -(dg J) -1 R T BS (dg J). The unscaled, flux, control-coefficient 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 matrices of ones and zeros, respectively. Eqs. (17) and (18)  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 elasticity-coefficients are known, the controlcoefficients can be calculated.
Using the steady-state elasticity-coefficients, we have The matrix (K ε') is also invertible (see Methods). Therefore, if the steady-state elasticity-coefficients are known, the control-coefficients 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 bio- [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 elasticity-coefficients for chemical potentials [25] , C1 0 Note the appearance of the matrices S T , S T B, and S T BS in the local and steady-state elasticity-coefficients and the control-coefficient 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 steady-state condition, Eq. (31) can also be derived using the homogeneous function theory (discussed above with flux and concentration control-coefficients). 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 elasticity-coefficients are known, the control-coefficients can be calculated. Using the steadystate elasticity-coefficients, we have The matrix (K ε') has an inverse (see Methods). Therefore, if the steady-state elasticity-coefficients are known, the control-coefficients 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.
CMC extracts information from reaction systems by using perturbations with insufficiently long correlation times.

Rate limiting step: largest flux control-coefficient and smallest biochemical conductance
Both flux control-coefficients and biochemical conductances can be used as indicators for rate-limiting steps in a biochemical network [30]. Traditionally the rate-limiting 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 first-passage 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 rate-limiting 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 control-coefficient and the smallest conductance are.
The NESS flux of reaction system (44) can be solved for analytically to yield [31,32] Therefore, the flux control-coefficient 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 rate-limiting 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].

Flux-controlled and concentration-controlled 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 con-stant-current (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 pseudo-steady-state enzyme-substrate 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 inhibi- 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 elasticity-coefficient 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 power-law expression which is essentially the integral form of the definition of local elasticity-coefficients.
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.

Solving the central, control-analysis equations
To demonstrate that the matrix, (K ε'), that appears in the central, control-analysis 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