Skip to main content

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

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 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 = SRTwhere 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 eAtplays 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).

Background

Developing methodologies for quantitative analysis, both experimental and mathematical, of complex biochemical networks has become one of the central themes of post-genomic 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 steady-state (NESS) theory [7, 8] originally developed to describe macromolecular, free-energy transduction [9–13], and the mathematical methods are based on linear analysis near a NESS.

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 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 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

ϕ + = J 1 − e Δ μ / k B T and ϕ − = J e Δ μ / k B T 1 − e Δ μ / k B T . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeWaaaqaaiabew9aMnaaBaaaleaacqGHRaWkaeqaaOGaeyypa0tcfa4aaSaaaeaacqWGkbGsaeaacqaIXaqmcqGHsislcqWGLbqzdaahaaqabeaacqqHuoarcqaH8oqBcqGGVaWlcqWGRbWAdaWgaaqaaiabdkeacbqabaGaemivaqfaaaaaaOqaaiabbggaHjabb6gaUjabbsgaKbqaaiabew9aMnaaBaaaleaacqGHsislaeqaaOGaeyypa0tcfa4aaSaaaeaacqWGkbGscqWGLbqzdaahaaqabeaacqqHuoarcqaH8oqBcqGGVaWlcqWGRbWAdaWgaaqaaiabdkeacbqabaGaemivaqfaaaqaaiabigdaXiabgkHiTiabdwgaLnaaCaaabeqaaiabfs5aejabeY7aTjabc+caViabdUgaRnaaBaaabaGaemOqaieabeaacqWGubavaaaaaOGaeiOla4caaaaa@5BCB@
(1)

Furthermore, if the standard-state free energy of reaction Δμois known from equilibrium thermodynamics, then e ( Δ μ − Δ μ o ) / k B T MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyzau2aaWbaaSqabeaacqGGOaakcqqHuoarcqaH8oqBcqGHsislcqqHuoarcqaH8oqBdaahaaadbeqaaiabd+gaVbaaliabcMcaPiabc+caViabdUgaRnaaBaaameaacqWGcbGqaeqaaSGaemivaqfaaaaa@3C87@ 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].

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 = SRTand 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 j th biochemical reaction (j = 1, 2, ..., M) is characterized by a set of stoichiometric coefficients ν = { ν i j } MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGae8xVd4Maeyypa0Jaei4EaSNaeqyVd42aa0baaSqaaiabdMgaPbqaaiabdQgaQbaakiabc2ha9baa@3642@ and κ = { κ i j } MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGae8NUdSMaeyypa0Jaei4EaSNaeqOUdS2aa0baaSqaaiabdMgaPbqaaiabdQgaQbaakiabc2ha9baa@3636@ such that

ν 1 j X 1 + ν 2 j X 2 + ⋯ + ν N j X N ⇌ k − j k + j κ 1 j X 1 + κ 2 j X 2 + ⋯ + κ N j X N . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyVd42aa0baaSqaaiabigdaXaqaaiabdQgaQbaakiabdIfaynaaBaaaleaacqaIXaqmaeqaaOGaey4kaSIaeqyVd42aa0baaSqaaiabikdaYaqaaiabdQgaQbaakiabdIfaynaaBaaaleaacqaIYaGmaeqaaOGaey4kaSIaeS47IWKaey4kaSIaeqyVd42aa0baaSqaaiabd6eaobqaaiabdQgaQbaakiabdIfaynaaBaaaleaacqWGobGtaeqaaOWaaCbmaeaacqWImhYGaSqaaiabdUgaRnaaDaaameaacqGHsislaeaacqWGQbGAaaaaleaacqWGRbWAdaqhaaadbaGaey4kaScabaGaemOAaOgaaaaakiabeQ7aRnaaDaaaleaacqaIXaqmaeaacqWGQbGAaaGccqWGybawdaWgaaWcbaGaeGymaedabeaakiabgUcaRiabeQ7aRnaaDaaaleaacqaIYaGmaeaacqWGQbGAaaGccqWGybawdaWgaaWcbaGaeGOmaidabeaakiabgUcaRiabl+UimjabgUcaRiabeQ7aRnaaDaaaleaacqWGobGtaeaacqWGQbGAaaGccqWGybawdaWgaaWcbaGaemOta4eabeaakiabc6caUaaa@68AC@
(2)

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 k + j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4AaS2aa0baaSqaaiabgUcaRaqaaiabdQgaQbaaaaa@2FA0@ and k − j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4AaS2aa0baaSqaaiabgkHiTaqaaiabdQgaQbaaaaa@2FAB@ . 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 Michaelis-Menten mechanism is expressed in terms of Eq. (2).

For the reaction system (2), kinetic equations of the time-dependent concentration changes for each biochemical species can be written according to the law of mass action as [19, 20]

d x i ( t ) d t = ∑ j = 1 M ( κ i j − ν i j ) ( k + j x 1 ν 1 i x 2 ν 2 i ⋯ x N ν N j − k − j x 1 κ 1 j x 2 κ 2 j ⋯ x N κ N j ) + J i e , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGKbazcqWG4baEdaWgaaqaaiabdMgaPbqabaGaeiikaGIaemiDaqNaeiykaKcabaGaemizaqMaemiDaqhaaOGaeyypa0ZaaabCaeaadaqadaqaaiabeQ7aRnaaDaaaleaacqWGPbqAaeaacqWGQbGAaaGccqGHsislcqaH9oGBdaqhaaWcbaGaemyAaKgabaGaemOAaOgaaaGccaGLOaGaayzkaaWaaeWaaeaacqWGRbWAdaqhaaWcbaGaey4kaScabaGaemOAaOgaaOGaemiEaG3aa0baaSqaaiabigdaXaqaaiabe27aUnaaDaaameaacqaIXaqmaeaacqWGPbqAaaaaaOGaemiEaG3aa0baaSqaaiabikdaYaqaaiabe27aUnaaDaaameaacqaIYaGmaeaacqWGPbqAaaaaaOGaeS47IWKaemiEaG3aa0baaSqaaiabd6eaobqaaiabe27aUnaaDaaameaacqWGobGtaeaacqWGQbGAaaaaaOGaeyOeI0Iaem4AaS2aa0baaSqaaiabgkHiTaqaaiabdQgaQbaakiabdIha4naaDaaaleaacqaIXaqmaeaacqaH6oWAdaqhaaadbaGaeGymaedabaGaemOAaOgaaaaakiabdIha4naaDaaaleaacqaIYaGmaeaacqaH6oWAdaqhaaadbaGaeGOmaidabaGaemOAaOgaaaaakiabl+UimjabdIha4naaDaaaleaacqWGobGtaeaacqaH6oWAdaqhaaadbaGaemOta4eabaGaemOAaOgaaaaaaOGaayjkaiaawMcaaiabgUcaRiabdQeaknaaDaaaleaacqWGPbqAaeaacqWGLbqzaaaabaGaemOAaOMaeyypa0JaeGymaedabaGaemyta0eaniabggHiLdGccqGGSaalaaa@87EA@
(3)

where we use xℓ to denote the concentration of respective Xℓ, and J i e MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaO0aa0baaSqaaiabdMgaPbqaaiabdwgaLbaaaaa@2FCD@ (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 J j = ϕ + j − ϕ − j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaO0aaWbaaSqabeaacqWGQbGAaaGccqGH9aqpcqaHvpGzdaqhaaWcbaGaey4kaScabaGaemOAaOgaaOGaeyOeI0Iaeqy1dy2aa0baaSqaaiabgkHiTaqaaiabdQgaQbaaaaa@38F6@ [20]:

ϕ + j = k + j ( x 1 ) ν 1 j ( x 2 ) ν 2 j ⋯ ( x N ) ν N j , ϕ − j = k − j ( x 1 ) κ 1 j ( x 2 ) κ 2 j ⋯ ( x N ) κ N j . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiqaaaqaaiabew9aMnaaDaaaleaacqGHRaWkaeaacqWGQbGAaaGccqGH9aqpcqWGRbWAdaqhaaWcbaGaey4kaScabaGaemOAaOgaaOGaeiikaGIaemiEaG3aaSbaaSqaaiabigdaXaqabaGccqGGPaqkdaahaaWcbeqaaiabe27aUnaaDaaameaacqaIXaqmaeaacqWGQbGAaaaaaOGaeiikaGIaemiEaG3aaSbaaSqaaiabikdaYaqabaGccqGGPaqkdaahaaWcbeqaaiabe27aUnaaDaaameaacqaIYaGmaeaacqWGQbGAaaaaaOGaeS47IWKaeiikaGIaemiEaG3aaSbaaSqaaiabd6eaobqabaGccqGGPaqkdaahaaWcbeqaaiabe27aUnaaDaaameaacqWGobGtaeaacqWGQbGAaaaaaOGaeiilaWcabaGaeqy1dy2aa0baaSqaaiabgkHiTaqaaiabdQgaQbaakiabg2da9iabdUgaRnaaDaaaleaacqGHsislaeaacqWGQbGAaaGccqGGOaakcqWG4baEdaWgaaWcbaGaeGymaedabeaakiabcMcaPmaaCaaaleqabaGaeqOUdS2aa0baaWqaaiabigdaXaqaaiabdQgaQbaaaaGccqGGOaakcqWG4baEdaWgaaWcbaGaeGOmaidabeaakiabcMcaPmaaCaaaleqabaGaeqOUdS2aa0baaWqaaiabikdaYaqaaiabdQgaQbaaaaGccqWIVlctcqGGOaakcqWG4baEdaWgaaWcbaGaemOta4eabeaakiabcMcaPmaaCaaaleqabaGaeqOUdS2aa0baaWqaaiabd6eaobqaaiabdQgaQbaaaaGccqGGUaGlaaaaaa@7946@

Then Eq. (3) becomes d x/dt = SJ + Je, in which x and Jeare N-dimensional column vectors, and J is an M-dimensional column vector.

In a closed reaction system, Je= 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 ( J i e MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaO0aa0baaSqaaiabdMgaPbqaaiabdwgaLbaaaaa@2FCD@ ≠ 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.

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:

( x 1 ) κ − j ( x 2 ) κ 2 j ⋯ ( x N ) κ N j ( x 1 ) ν 1 j ( x 2 ) ν 2 j ⋯ ( x N ) ν N j = k + j k − j , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqGGOaakcqWG4baEdaWgaaqaaiabigdaXaqabaGaeiykaKYaaWbaaeqabaGaeqOUdS2aa0baaeaacqGHsislaeaacqWGQbGAaaaaaiabcIcaOiabdIha4naaBaaabaGaeGOmaidabeaacqGGPaqkdaahaaqabeaacqaH6oWAdaqhaaqaaiabikdaYaqaaiabdQgaQbaaaaGaeS47IWKaeiikaGIaemiEaG3aaSbaaeaacqWGobGtaeqaaiabcMcaPmaaCaaabeqaaiabeQ7aRnaaDaaabaGaemOta4eabaGaemOAaOgaaaaaaeaacqGGOaakcqWG4baEdaWgaaqaaiabigdaXaqabaGaeiykaKYaaWbaaeqabaGaeqyVd42aa0baaeaacqaIXaqmaeaacqWGQbGAaaaaaiabcIcaOiabdIha4naaBaaabaGaeGOmaidabeaacqGGPaqkdaahaaqabeaacqaH9oGBdaqhaaqaaiabikdaYaqaaiabdQgaQbaaaaGaeS47IWKaeiikaGIaemiEaG3aaSbaaeaacqWGobGtaeqaaiabcMcaPmaaCaaabeqaaiabe27aUnaaDaaabaGaemOta4eabaGaemOAaOgaaaaaaaGccqGH9aqpjuaGdaWcaaqaaiabdUgaRnaaDaaabaGaey4kaScabaGaemOAaOgaaaqaaiabdUgaRnaaDaaabaGaeyOeI0cabaGaemOAaOgaaaaakiabcYcaSaaa@6EA5@

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 (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 { x i ∗ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdMgaPbqaaiabgEHiQaaaaaa@2FC5@ } is an asymptotically stable NESS, small differences in the initial condition should lead to the same NESS. We use the notations S MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbiGae83uamfaaa@2D0B@ to denote the set of concentrations of internal species and I MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbiGae8xsaKeaaa@2CF7@ to denote the concentrations of external species which are clamped or independently varied.

In the linear regime near the NESS [17],

d d t ( δ x i ) = ∑ j = 1 N A i j ( δ x j ) x j ∗ , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGKbazaeaacqWGKbazcqWG0baDaaGccqGGOaakcqaH0oazcqWG4baEdaWgaaWcbaGaemyAaKgabeaakiabcMcaPiabg2da9maaqahabaGaemyqae0aaSbaaSqaaiabdMgaPjabdQgaQbqabaqcfa4aaSaaaeaacqGGOaakcqaH0oazcqWG4baEdaWgaaqaaiabdQgaQbqabaGaeiykaKcabaGaemiEaG3aa0baaeaacqWGQbGAaeaacqGHxiIkaaaaaaWcbaGaemOAaOMaeyypa0JaeGymaedabaGaemOta4eaniabggHiLdGccqGGSaalaaa@4EAB@
(4)

in which δ x j = x j - x j ∗ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@2FC7@ and

A i j = ∑ ℓ = 1 M ( κ i ℓ − ν i ℓ ) ( ν j ℓ ϕ + ℓ − κ j ℓ ϕ − ℓ ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyqae0aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGH9aqpdaaeWbqaaiabcIcaOiabeQ7aRnaaDaaaleaacqWGPbqAaeaacqWItecBaaGccqGHsislcqaH9oGBdaqhaaWcbaGaemyAaKgabaGaeS4eHWgaaOGaeiykaKIaeiikaGIaeqyVd42aa0baaSqaaiabdQgaQbqaaiabloriSbaakiabew9aMnaaDaaaleaacqGHRaWkaeaacqWItecBaaGccqGHsislcqaH6oWAdaqhaaWcbaGaemOAaOgabaGaeS4eHWgaaOGaeqy1dy2aa0baaSqaaiabgkHiTaqaaiabloriSbaakiabcMcaPaWcbaGaeS4eHWMaeyypa0JaeGymaedabaGaemyta0eaniabggHiLdGccqGGUaGlaaa@57ED@

A = {A ij } is called the linear stability matrix. We introduce the N × M matrix R = { R i j = ν i j ϕ + j − κ i j ϕ − j } MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaSNaemOuai1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGH9aqpcqaH9oGBdaqhaaWcbaGaemyAaKgabaGaemOAaOgaaOGaeqy1dy2aa0baaSqaaiabgUcaRaqaaiabdQgaQbaakiabgkHiTiabeQ7aRnaaDaaaleaacqWGPbqAaeaacqWGQbGAaaGccqaHvpGzdaqhaaWcbaGaeyOeI0cabaGaemOAaOgaaOGaeiyFa0haaa@46B2@ , such that A = SRT. 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 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 { x i ∗ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdMgaPbqaaiabgEHiQaaaaaa@2FC5@ } shifts in response to a perturbation to the amount of enzyme for a reaction or a perturbation to its substrate (∈ I MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbiGae8xsaKeaaa@2CF7@ ) 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: x N ∗ → x N ∗ + δ x N MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabd6eaobqaaiabgEHiQaaakiabgkziUkabdIha4naaDaaaleaacqWGobGtaeaacqGHxiIkaaGccqGHRaWkcqaH0oazcqWG4baEdaWgaaWcbaGaemOta4eabeaaaaa@3A9B@ . 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

ϵ N m ≜ x N ∗ J m ( δ J m δ x N ) = ν N m ϕ + m − κ N m ϕ − m ϕ + m − ϕ − m = R N m J m , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWeuvgwd1utHrhAjrxySL2yaeHbJ1wBPfdmaGabciab=v=aYoaaDaaaleaacqWGobGtaeaacqWGTbqBaaGccqWICjcqjuaGdaWcaaqaaiabdIha4naaDaaabaGaemOta4eabaGaey4fIOcaaaqaaiabdQeaknaaCaaabeqaaiabd2gaTbaaaaGcdaqadaqcfayaamaalaaabaGaeqiTdqMaemOsaO0aaWbaaeqabaGaemyBa0gaaaqaaiabes7aKjabdIha4naaBaaabaGaemOta4eabeaaaaaakiaawIcacaGLPaaacqGH9aqpjuaGdaWcaaqaaiabe27aUnaaDaaabaGaemOta4eabaGaemyBa0gaaiabew9aMnaaDaaabaGaey4kaScabaGaemyBa0gaaiabgkHiTiabeQ7aRnaaDaaabaGaemOta4eabaGaemyBa0gaaiabew9aMnaaDaaabaGaeyOeI0cabaGaemyBa0gaaaqaaiabew9aMnaaDaaabaGaey4kaScabaGaemyBa0gaaiabgkHiTiabew9aMnaaDaaabaGaeyOeI0cabaGaemyBa0gaaaaakiabg2da9KqbaoaalaaabaGaemOuai1aaSbaaeaacqWGobGtcqWGTbqBaeqaaaqaaiabdQeaknaaCaaabeqaaiabd2gaTbaaaaGaeiilaWcaaa@7209@
(5)

which is called the scaled, local elasticity-coefficient [1]. These coefficients are elements of the scaled, local, elasticity-coefficient matrix, ϵ = (dg J)-1 RT. 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, elasticity-coefficient matrix is given by ϵ '= (dg J) ϵ (dg x)-1 = RT(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 steady-state response [24]. When a new NESS is established, the concentrations {x i } (i = 1, 2, ..., N - 1) satisfy

∑ j = 1 N − 1 A i j δ x j x j ∗ = − A i N δ x N x N ∗ . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaabCaeaacqWGbbqqdaWgaaWcbaGaemyAaKMaemOAaOgabeaajuaGdaWcaaqaaiabes7aKjabdIha4naaBaaabaGaemOAaOgabeaaaeaacqWG4baEdaqhaaqaaiabdQgaQbqaaiabgEHiQaaaaaGccqGH9aqpcqGHsislcqWGbbqqdaWgaaWcbaGaemyAaKMaemOta4eabeaajuaGdaWcaaqaaiabes7aKjabdIha4naaBaaabaGaemOta4eabeaaaeaacqWG4baEdaqhaaqaaiabd6eaobqaaiabgEHiQaaaaaGaeiOla4caleaacqWGQbGAcqGH9aqpcqaIXaqmaeaacqWGobGtcqGHsislcqaIXaqma0GaeyyeIuoaaaa@5113@
(6)

Solving Eq. (6), we have δ x i / x i ∗ = B i N δ x N / B N N x N ∗ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiTdqMaemiEaG3aaSbaaSqaaiabdMgaPbqabaGccqGGVaWlcqWG4baEdaqhaaWcbaGaemyAaKgabaGaey4fIOcaaOGaeyypa0JaemOqai0aaSbaaSqaaiabdMgaPjabd6eaobqabaGccqaH0oazcqWG4baEdaWgaaWcbaGaemOta4eabeaakiabc+caViabdkeacnaaBaaaleaacqWGobGtcqWGobGtaeqaaOGaemiEaG3aa0baaSqaaiabd6eaobqaaiabgEHiQaaaaaa@46D3@ where B iN is the N th column vector of the matrix A-1. The new NESS established near { x i ∗ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdMgaPbqaaiabgEHiQaaaaaa@2FC5@ } is { x i ∗ + δ x i MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdMgaPbqaaiabgEHiQaaakiabgUcaRiabes7aKjabdIha4naaBaaaleaacqWGPbqAaeqaaaaa@3556@ }, and the new fluxes, Jm+ δ Jm, are given by

δ J m = δ ϕ + m − δ ϕ − m = ∑ ℓ = 1 N ( ν ℓ m ϕ + m − κ ℓ m ϕ − m ) δ x ℓ x ℓ ∗ = ∑ ℓ = 1 N ( ν ℓ m ϕ + m − κ ℓ m ϕ − m ) B ℓ N δ x N B N N x N ∗ . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiWaaaqaaiabes7aKjabdQeaknaaCaaaleqabaGaemyBa0gaaOGaeyypa0JaeqiTdqMaeqy1dy2aa0baaSqaaiabgUcaRaqaaiabd2gaTbaakiabgkHiTiabes7aKjabew9aMnaaDaaaleaacqGHsislaeaacqWGTbqBaaaakeaacqGH9aqpaeaadaaeWbqaaiabcIcaOiabe27aUnaaDaaaleaacqWItecBaeaacqWGTbqBaaGccqaHvpGzdaqhaaWcbaGaey4kaScabaGaemyBa0gaaOGaeyOeI0IaeqOUdS2aa0baaSqaaiabloriSbqaaiabd2gaTbaakiabew9aMnaaDaaaleaacqGHsislaeaacqWGTbqBaaGccqGGPaqkjuaGdaWcaaqaaiabes7aKjabdIha4naaBaaabaGaeS4eHWgabeaaaeaacqWG4baEdaqhaaqaaiabloriSbqaaiabgEHiQaaaaaaaleaacqWItecBcqGH9aqpcqaIXaqmaeaacqWGobGta0GaeyyeIuoaaOqaaaqaaiabg2da9aqaamaaqahabaGaeiikaGIaeqyVd42aa0baaSqaaiabloriSbqaaiabd2gaTbaakiabew9aMnaaDaaaleaacqGHRaWkaeaacqWGTbqBaaGccqGHsislcqaH6oWAdaqhaaWcbaGaeS4eHWgabaGaemyBa0gaaOGaeqy1dy2aa0baaSqaaiabgkHiTaqaaiabd2gaTbaakiabcMcaPKqbaoaalaaabaGaemOqai0aaSbaaeaacqWItecBcqWGobGtaeqaaiabes7aKjabdIha4naaBaaabaGaemOta4eabeaaaeaacqWGcbGqdaWgaaqaaiabd6eaojabd6eaobqabaGaemiEaG3aa0baaeaacqWGobGtaeaacqGHxiIkaaaaaiabc6caUaWcbaGaeS4eHWMaeyypa0JaeGymaedabaGaemOta4eaniabggHiLdaaaaaa@8F54@
(7)

Hence, the scaled, steady-state elasticity-coefficient is [17]

ε N m ≜ x N ∗ J m ( δ J m δ x N ) = ∑ ℓ = 1 N ( ν ℓ m ϕ + m − κ ℓ m ϕ − m ) B ℓ N ( ϕ + m − ϕ − m ) B N N , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyTdu2aa0baaSqaaiabd6eaobqaaiabd2gaTbaakiablYLiaLqbaoaalaaabaGaemiEaG3aa0baaeaacqWGobGtaeaacqGHxiIkaaaabaGaemOsaO0aaWbaaeqabaGaemyBa0gaaaaakmaabmaajuaGbaWaaSaaaeaacqaH0oazcqWGkbGsdaahaaqabeaacqWGTbqBaaaabaGaeqiTdqMaemiEaG3aaSbaaeaacqWGobGtaeqaaaaaaOGaayjkaiaawMcaaiabg2da9KqbaoaalaaabaWaaabmaeaacqGGOaakcqaH9oGBdaqhaaqaaiabloriSbqaaiabd2gaTbaacqaHvpGzdaqhaaqaaiabgUcaRaqaaiabd2gaTbaacqGHsislcqaH6oWAdaqhaaqaaiabloriSbqaaiabd2gaTbaacqaHvpGzdaqhaaqaaiabgkHiTaqaaiabd2gaTbaacqGGPaqkcqWGcbGqdaWgaaqaaiabloriSjabd6eaobqabaaabaGaeS4eHWMaeyypa0JaeGymaedabaGaemOta4eacqGHris5aaqaaiabcIcaOiabew9aMnaaDaaabaGaey4kaScabaGaemyBa0gaaiabgkHiTiabew9aMnaaDaaabaGaeyOeI0cabaGaemyBa0gaaiabcMcaPiabdkeacnaaBaaabaGaemOta4KaemOta4eabeaaaaGaeiilaWcaaa@7249@
(8)

which may also be known as a response coefficient [23, 25]. These coefficients are elements of the scaled, steady-state, elasticity-coefficient matrix, ε= (dg J)-1 RTB (dg B)-1. Here, (dg B) denotes the diagonal matrix containing only the diagonal components of the matrix B along its diagonal. The unscaled, steady-state, elasticity-coefficient matrix is given by ε ' = (dg J) ε(dg x)-1 = RTB (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 EMis changed: EM→ EM+ δ EM.

Since we have assumed that the rate constants for reaction M, k ± M MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4AaS2aa0baaSqaaiabgglaXcqaaiabd2eanbaaaaa@3072@ , are linearly proportional to the concentration of the enzyme catalyzing reaction M, EM, we have

E M k + M ( δ k + M δ E M ) = E M k − M ( δ k − M δ E M ) = 1. MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGfbqrdaahaaqabeaacqWGnbqtaaaabaGaem4AaS2aa0baaeaacqGHRaWkaeaacqWGnbqtaaaaaOWaaeWaaKqbagaadaWcaaqaaiabes7aKjabdUgaRnaaDaaabaGaey4kaScabaGaemyta0eaaaqaaiabes7aKjabdweafnaaCaaabeqaaiabd2eanbaaaaaakiaawIcacaGLPaaacqGH9aqpjuaGdaWcaaqaaiabdweafnaaCaaabeqaaiabd2eanbaaaeaacqWGRbWAdaqhaaqaaiabgkHiTaqaaiabd2eanbaaaaGcdaqadaqcfayaamaalaaabaGaeqiTdqMaem4AaS2aa0baaeaacqGHsislaeaacqWGnbqtaaaabaGaeqiTdqMaemyrau0aaWbaaeqabaGaemyta0eaaaaaaOGaayjkaiaawMcaaiabg2da9iabigdaXiabc6caUaaa@53D7@

When EM→ EM+ δ EM, the new NESS satisfies [17]

∑ j = 1 N A i j δ x j x j ∗ = ( ν i M − κ i M ) ( ϕ + M − ϕ − M ) δ E M E M . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaabCaeaacqWGbbqqdaWgaaWcbaGaemyAaKMaemOAaOgabeaajuaGdaWcaaqaaiabes7aKjabdIha4naaBaaabaGaemOAaOgabeaaaeaacqWG4baEdaqhaaqaaiabdQgaQbqaaiabgEHiQaaaaaaaleaacqWGQbGAcqGH9aqpcqaIXaqmaeaacqWGobGta0GaeyyeIuoakiabg2da9iabcIcaOiabe27aUnaaDaaaleaacqWGPbqAaeaacqWGnbqtaaGccqGHsislcqaH6oWAdaqhaaWcbaGaemyAaKgabaGaemyta0eaaOGaeiykaKIaeiikaGIaeqy1dy2aa0baaSqaaiabgUcaRaqaaiabd2eanbaakiabgkHiTiabew9aMnaaDaaaleaacqGHsislaeaacqWGnbqtaaGccqGGPaqkjuaGdaWcaaqaaiabes7aKjabdweafnaaCaaabeqaaiabd2eanbaaaeaacqWGfbqrdaahaaqabeaacqWGnbqtaaaaaOGaeiOla4caaa@5F03@
(9)

Solving for δ x j / x j ∗ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@2FC7@ (j = 1, 2, ..., N), we have

δ x j x j ∗ = ∑ i = 1 N B j i ( ν i M − κ i M ) ( ϕ + M − ϕ − M ) δ E M E M . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqaH0oazcqWG4baEdaWgaaqaaiabdQgaQbqabaaabaGaemiEaG3aa0baaeaacqWGQbGAaeaacqGHxiIkaaaaaiabg2da9OWaaabCaeaacqWGcbGqdaWgaaWcbaGaemOAaOMaemyAaKgabeaaaeaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGobGta0GaeyyeIuoakiabcIcaOiabe27aUnaaDaaaleaacqWGPbqAaeaacqWGnbqtaaGccqGHsislcqaH6oWAdaqhaaWcbaGaemyAaKgabaGaemyta0eaaOGaeiykaKIaeiikaGIaeqy1dy2aa0baaSqaaiabgUcaRaqaaiabd2eanbaakiabgkHiTiabew9aMnaaDaaaleaacqGHsislaeaacqWGnbqtaaGccqGGPaqkjuaGdaWcaaqaaiabes7aKjabdweafnaaCaaabeqaaiabd2eanbaaaeaacqWGfbqrdaahaaqabeaacqWGnbqtaaaaaOGaeiOla4caaa@5F02@
(10)

From Eq. (10), we get the scaled, concentration control-coefficient [26],

C ^ j M ≜ E M x j ∗ ( δ x j δ E M ) = − ∑ j = 1 N B j i ( κ i M − ν i M ) ( ϕ + M − ϕ − M ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfaOafm4qamKbaKaadaqhaaqaaiabdQgaQbqaaiabd2eanbaacqWICjcqdaWcaaqaaiabdweafnaaCaaabeqaaiabd2eanbaaaeaacqWG4baEdaqhaaqaaiabdQgaQbqaaiabgEHiQaaaaaWaaeWaaeaadaWcaaqaaiabes7aKjabdIha4naaBaaabaGaemOAaOgabeaaaeaacqaH0oazcqWGfbqrdaahaaqabeaacqWGnbqtaaaaaaGaayjkaiaawMcaaiabg2da9iabgkHiTOWaaabCaeaacqWGcbGqdaWgaaWcbaGaemOAaOMaemyAaKgabeaaaeaacqWGQbGAcqGH9aqpcqaIXaqmaeaacqWGobGta0GaeyyeIuoakiabcIcaOiabeQ7aRnaaDaaaleaacqWGPbqAaeaacqWGnbqtaaGccqGHsislcqaH9oGBdaqhaaWcbaGaemyAaKgabaGaemyta0eaaOGaeiykaKIaeiikaGIaeqy1dy2aa0baaSqaaiabgUcaRaqaaiabd2eanbaakiabgkHiTiabew9aMnaaDaaaleaacqGHsislaeaacqWGnbqtaaGccqGGPaqkcqGGSaalaaa@65EE@
(11)

which is an element of the scaled, concentration, control-coefficient matrix given by C ^ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbeGaf83qamKbaKaaaaa@2CFA@ = -BS (dg J). The unscaled, concentration, control-coefficient matrix is C ^ ′ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbeGaf83qamKbaKGbauaaaaa@2D05@ = (dg x) C ^ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbeGaf83qamKbaKaaaaa@2CFA@ (dg J)-1 = - (dg x) BS.

Combining Eqs. (7) and (10), we get the scaled, flux control-coefficient [1],

C M m ≜ E M J m ( δ J m δ E M ) = δ m M − ∑ ℓ = 1 N ∑ i = 1 N ( ν ℓ m ϕ + m − κ ℓ m ϕ − m ) B ℓ i ( κ i M − ν i M ) ( ϕ + M − ϕ − M ) ( ϕ + m − ϕ − m ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiWaaaqaaKqbakabdoeadnaaDaaabaGaemyta0eabaGaemyBa0gaaaGcbaqcfaOaeSixIaeakeaajuaGdaWcaaqaaiabdweafnaaCaaabeqaaiabd2eanbaaaeaacqWGkbGsdaahaaqabeaacqWGTbqBaaaaamaabmaabaWaaSaaaeaacqaH0oazcqWGkbGsdaahaaqabeaacqWGTbqBaaaabaGaeqiTdqMaemyrau0aaWbaaeqabaGaemyta0eaaaaaaiaawIcacaGLPaaaaOqaaaqaaiabg2da9aqaaiabes7aKnaaBaaaleaacqWGTbqBcqWGnbqtaeqaaOGaeyOeI0scfa4aaSaaaeaadaaeWaqaamaaqadabaGaeiikaGIaeqyVd42aa0baaeaacqWItecBaeaacqWGTbqBaaGaeqy1dy2aa0baaeaacqGHRaWkaeaacqWGTbqBaaGaeyOeI0IaeqOUdS2aa0baaeaacqWItecBaeaacqWGTbqBaaGaeqy1dy2aa0baaeaacqGHsislaeaacqWGTbqBaaGaeiykaKIaemOqai0aaSbaaeaacqWItecBcqWGPbqAaeqaaiabcIcaOiabeQ7aRnaaDaaabaGaemyAaKgabaGaemyta0eaaiabgkHiTiabe27aUnaaDaaabaGaemyAaKgabaGaemyta0eaaiabcMcaPiabcIcaOiabew9aMnaaDaaabaGaey4kaScabaGaemyta0eaaiabgkHiTiabew9aMnaaDaaabaGaeyOeI0cabaGaemyta0eaaiabcMcaPaqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6eaobGaeyyeIuoaaeaacqWItecBcqGH9aqpcqaIXaqmaeaacqWGobGtaiabggHiLdaabaGaeiikaGIaeqy1dy2aa0baaeaacqGHRaWkaeaacqWGTbqBaaGaeyOeI0Iaeqy1dy2aa0baaeaacqGHsislaeaacqWGTbqBaaGaeiykaKcaaiabcYcaSaaaaaa@8EC6@
(12)

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 RTBS (dg J). The unscaled, flux, control-coefficient matrix is C '= (dg J) C (dg J)-1 = I - RTBS.

Summation and connectivity theorems

Among ϵ, ε, C, and C ^ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafm4qamKbaKaaaaa@2CF4@ , we have the following relationships. First, note the relationships

C ^ = − ( S ( dg  J ) ϵ ) − 1 S ( dg  J ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafC4qamKbaKaacqGH9aqpcqGHsislcqGGOaakcqWHtbWucqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWHkbGscqGGPaqktqvzynutnfgDOLeDHXwAJbqegmwBTLwmWaacemGae8x9diRaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqqGGaaicqWHtbWucqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWHkbGscqGGPaqkcqGGSaalaaa@4D78@
(13)
C = I + ϵ C ^ , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeC4qamKaeyypa0JaeCysaKKaey4kaSYeuvgwd1utHrhAjrxySL2yaeHbJ1wBPfdmaGabdiab=v=aYkqbhoeadzaajaGaeiilaWcaaa@3C30@
(14)

between the scaled coefficients, and similarly,

C ^ ′ = − ( S ϵ ′ ) − 1 S , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafC4qamKbaKGbauaacqGH9aqpcqGHsislcqGGOaakcqWHtbWutqvzynutnfgDOLeDHXwAJbqegmwBTLwmWaacemGaf8x9diRbauaacqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabhofatjabcYcaSaaa@404C@
(15)
C ′ = I + ϵ C ^ ′ , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafC4qamKbauaacqGH9aqpcqWHjbqscqGHRaWktqvzynutnfgDOLeDHXwAJbqegmwBTLwmWaacemGae8x9diRafC4qamKbaKGbauaacqGGSaalaaa@3C47@
(16)

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 = -Je= 0, we have

C ^ 1 M × 1 = 0 N × 1 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafC4qamKbaKaacqWHXaqmdaWgaaWcbaGaemyta0Kaey41aqRaeGymaedabeaakiabg2da9iabhcdaWmaaBaaaleaacqWGobGtcqGHxdaTcqaIXaqmaeqaaOGaeiilaWcaaa@39CA@
(17)

C1M × 1= 1M × 1,

where 1i×jand 0i×jdenote the i × j matricies of ones and zeros, respectively. Eqs. (17) and (18) are known as the summation theorems for the scaled control-coefficients. For the unscaled control-coefficients, we have the generalized summation theorems, C ' K = K and C ^ ′ K = 0 N × ( M − N ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafC4qamKbaKGbauaacqWHlbWscqGH9aqpcqWHWaamdaWgaaWcbaGaemOta4Kaey41aqRaeiikaGIaemyta0KaeyOeI0IaemOta4KaeiykaKcabeaaaaa@3868@ , 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 = 0N×(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 steady-state 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 control-coefficients. On the other hand, if there are no clamped concentrations and a NESS is sustained by injection fluxes, the steady-state fluxes are homogeneous functions of enzyme concentrations with order 0. A similar argument exists for steady-state concentrations as a function of enzyme concentrations, leading to the summation theorem for the concentration control-coefficients.

Third, BSRT= I ⇒ (RTBS)(RTB) = RTB, which yields the connectivity theorems for the scaled coefficients:

C ^ ϵ = − I , C ^ ε = − B ( dg  B ) − 1 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaaiqbhoeadzaajaWeuvgwd1utHrhAjrxySL2yaeHbJ1wBPfdmaGabdiab=v=aYkabg2da9iabgkHiTiabhMeajjabcYcaSaqaaiqbhoeadzaajaaccmGae4xTduMaeyypa0JaeyOeI0IaeCOqaiKaeiikaGIaeeizaqMaee4zaCMaeeiiaaIaeCOqaiKaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqGGSaalaaaaaa@4A2C@
(19)

C ϵ = 0M×N,   C ε = 0M×N,

and the generalized connectivity theorems for the unscaled coefficients [23]:

C ^ ′ ϵ ′ = − I , C ^ ′ ε ′ = − ( dg  x ) B ( dg  B ) − 1 ( dg  x ) − 1 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaaiqbhoeadzaajyaafaWeuvgwd1utHrhAjrxySL2yaeHbJ1wBPfdmaGabdiqb=v=aYAaafaGaeyypa0JaeyOeI0IaeCysaKKaeiilaWcabaGafC4qamKbaKGbauaaiiWacuGF1oqzgaqbaiabg2da9iabgkHiTiabcIcaOiabbsgaKjabbEgaNjabbccaGiabhIha4jabcMcaPiabbccaGiabhkeacHqabiab9bcaGiabcIcaOiabbsgaKjabbEgaNjabbccaGiabhkeacjabcMcaPmaaCaaaleqabaGaeyOeI0IaeGymaedaaOGaeiikaGIaeeizaqMaee4zaCMaeeiiaaIaeCiEaGNaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqGGSaalaaaaaa@5B36@
(21)

C ' ϵ ' = 0M×N,   C ' ε ' = 0M×N.

Finally, by combining the generalized summation and connectivity theorems, we get

( C ′ C ^ ′ ) ( K ϵ ′ ) = ( K 0 M × N 0 N × ( M − N ) − I N × N ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaeWaaeaafaqabeGabaaabaGafC4qamKbauaaaeaacuWHdbWqgaqcgaqbaaaaaiaawIcacaGLPaaadaqadaqaauaabeqabiaaaeaacqWHlbWsaeaatqvzynutnfgDOLeDHXwAJbqegmwBTLwmWaacemGaf8x9diRbauaaaaaacaGLOaGaayzkaaGaeyypa0ZaaeWaaeaafaqaaeGacaaabaGaeC4saSeabaGaeCimaaZaaSbaaSqaaiabd2eanjabgEna0kabd6eaobqabaaakeaacqWHWaamdaWgaaWcbaGaemOta4Kaey41aqRaeiikaGIaemyta0KaeyOeI0IaemOta4KaeiykaKcabeaaaOqaaiabgkHiTiabhMeajnaaBaaaleaacqWGobGtcqGHxdaTcqWGobGtaeqaaaaaaOGaayjkaiaawMcaaiabcYcaSaaa@56C7@
(23)

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 control-coefficients can be calculated.

Using the steady-state elasticity-coefficients, we have

( C ′ C ^ ′ ) ( K ε ′ ) = ( K 0 N × ( M − N ) 0 M × N − ( dg  x ) B ( dg  B ) − 1 ( dg  x ) − 1 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaeWaaeaafaqabeGabaaabaGafC4qamKbauaaaeaacuWHdbWqgaqcgaqbaaaaaiaawIcacaGLPaaadaqadaqaauaabeqabiaaaeaacqWHlbWsaeaaiiWacuWF1oqzgaqbaaaaaiaawIcacaGLPaaacqGH9aqpdaqadaqaauaabaqaceaaaeaacqWHlbWsaeaacqWHWaamdaWgaaWcbaGaemOta4Kaey41aqRaeiikaGIaemyta0KaeyOeI0IaemOta4KaeiykaKcabeaaaaGccqqGGaaifaqabeGabaaabaGaeCimaaZaaSbaaSqaaiabd2eanjabgEna0kabd6eaobqabaaakeaacqGHsislcqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWH4baEcqGGPaqkcqqGGaaicqWHcbGqieqacqGFGaaicqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWHcbGqcqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabcIcaOiabbsgaKjabbEgaNjabbccaGiabhIha4jabcMcaPmaaCaaaleqabaGaeyOeI0IaeGymaedaaaaaaOGaayjkaiaawMcaaiabc6caUaaa@63D5@
(24)

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]:

J = ϕ + ( 1 − e Δ μ / k B T ) = ϕ − ( e − Δ μ / k B T − 1 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaOKaeyypa0Jaeqy1dy2aaSbaaSqaaiabgUcaRaqabaGcdaqadaqaaiabigdaXiabgkHiTiabdwgaLnaaCaaaleqabaGaeuiLdqKaeqiVd0Maei4la8Iaem4AaS2aaSbaaWqaaiabdkeacbqabaWccqWGubavaaaakiaawIcacaGLPaaacqGH9aqpcqaHvpGzdaWgaaWcbaGaeyOeI0cabeaakmaabmaabaGaemyzau2aaWbaaSqabeaacqGHsislcqqHuoarcqaH8oqBcqGGVaWlcqWGRbWAdaWgaaadbaGaemOqaieabeaaliabdsfaubaakiabgkHiTiabigdaXaGaayjkaiaawMcaaiabc6caUaaa@5076@
(25)

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 elasticity-coefficients for chemical potentials [25]

( ϵ Δ μ ) n m ≜ x n ∗ Δ μ m ( δ Δ μ m δ x n ) = k B T ( κ n m − ν n m ) Δ μ m MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaeWaaeaatqvzynutnfgDOLeDHXwAJbqegmwBTLwmWaaceiGae8x9di7aaWbaaSqabeaacqqHuoarcqaH8oqBaaaakiaawIcacaGLPaaadaqhaaWcbaGaemOBa4gabaGaemyBa0gaaOGaeSixIaucfa4aaSaaaeaacqWG4baEdaqhaaqaaiabd6gaUbqaaiabgEHiQaaaaeaacqqHuoarcqaH8oqBdaahaaqabeaacqWGTbqBaaaaamaabmaabaWaaSaaaeaacqaH0oazcqqHuoarcqaH8oqBdaahaaqabeaacqWGTbqBaaaabaGaeqiTdqMaemiEaG3aaSbaaeaacqWGUbGBaeqaaaaaaiaawIcacaGLPaaacqGH9aqpcqWGRbWAdaWgaaqaaiabdkeacbqabaGaemivaq1aaSaaaeaadaqadaqaaiabeQ7aRnaaDaaabaGaemOBa4gabaGaemyBa0gaaiabgkHiTiabe27aUnaaDaaabaGaemOBa4gabaGaemyBa0gaaaGaayjkaiaawMcaaaqaaiabfs5aejabeY7aTnaaCaaabeqaaiabd2gaTbaaaaaaaa@6A16@
(26)

as elements of the matrix ϵΔμ= k B T (dg Δ μ)-1 ST, the scaled, steady-state elasticity-coefficients for chemical potentials

( ε Δ μ ) n m ≜ x n ∗ Δ μ m ( δ Δ μ m δ x n ) = k B T ∑ ℓ = 1 N ( κ ℓ m − ν ℓ m ) B ℓ n Δ μ m B n n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaeWaaeaacqaH1oqzdaahaaWcbeqaaiabfs5aejabeY7aTbaaaOGaayjkaiaawMcaamaaDaaaleaacqWGUbGBaeaacqWGTbqBaaGccqWICjcqjuaGdaWcaaqaaiabdIha4naaDaaabaGaemOBa4gabaGaey4fIOcaaaqaaiabfs5aejabeY7aTnaaCaaabeqaaiabd2gaTbaaaaWaaeWaaeaadaWcaaqaaiabes7aKjabfs5aejabeY7aTnaaCaaabeqaaiabd2gaTbaaaeaacqaH0oazcqWG4baEdaWgaaqaaiabd6gaUbqabaaaaaGaayjkaiaawMcaaiabg2da9iabdUgaRnaaBaaabaGaemOqaieabeaacqWGubavdaWcaaqaamaaqadabaWaaeWaaeaacqaH6oWAdaqhaaqaaiabloriSbqaaiabd2gaTbaacqGHsislcqaH9oGBdaqhaaqaaiabloriSbqaaiabd2gaTbaaaiaawIcacaGLPaaaaeaacqWItecBcqGH9aqpcqaIXaqmaeaacqWGobGtaiabggHiLdGaemOqai0aaSbaaeaacqWItecBcqWGUbGBaeqaaaqaaiabfs5aejabeY7aTnaaCaaabeqaaiabd2gaTbaacqWGcbGqdaWgaaqaaiabd6gaUjabd6gaUbqabaaaaaaa@6F58@
(27)

as elements of the matrix εΔμ= k B T (dg Δ μ)-1 STB (dg B)-1, and the scaled, potential control-coefficients [25]

( C Δ μ ) n m ≜ E n Δ μ m ( δ Δ μ m δ E n ) = − k B T ∑ ℓ = 1 N ∑ i = 1 N ( κ ℓ m − ν ℓ m ) B ℓ i ( κ i n − ν i n ) ( ϕ + n − ϕ − n ) Δ μ m MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiWaaaqaamaabmaabaGaem4qam0aaWbaaSqabeaacqqHuoarcqaH8oqBaaaakiaawIcacaGLPaaadaqhaaWcbaGaemOBa4gabaGaemyBa0gaaaGcbaGaeSixIaeajuaGbaWaaSaaaeaacqWGfbqrdaahaaqabeaacqWGUbGBaaaabaGaeuiLdqKaeqiVd02aaWbaaeqabaGaemyBa0gaaaaadaqadaqaamaalaaabaGaeqiTdqMaeuiLdqKaeqiVd02aaWbaaeqabaGaemyBa0gaaaqaaiabes7aKjabdweafnaaCaaabeqaaiabd6gaUbaaaaaacaGLOaGaayzkaaaakeaaaeaacqGH9aqpaeaacqGHsislcqWGRbWAdaWgaaWcbaGaemOqaieabeaakiabdsfauLqbaoaalaaabaWaaabmaeaadaaeWaqaaiabcIcaOiabeQ7aRnaaDaaabaGaeS4eHWgabaGaemyBa0gaaiabgkHiTiabe27aUnaaDaaabaGaeS4eHWgabaGaemyBa0gaaiabcMcaPiabdkeacnaaBaaabaGaeS4eHWMaemyAaKgabeaacqGGOaakcqaH6oWAdaqhaaqaaiabdMgaPbqaaiabd6gaUbaacqGHsislcqaH9oGBdaqhaaqaaiabdMgaPbqaaiabd6gaUbaacqGGPaqkcqGGOaakcqaHvpGzdaqhaaqaaiabgUcaRaqaaiabd6gaUbaacqGHsislcqaHvpGzdaqhaaqaaiabgkHiTaqaaiabd6gaUbaacqGGPaqkaeaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGobGtaiabggHiLdaabaGaeS4eHWMaeyypa0JaeGymaedabaGaemOta4eacqGHris5aaqaaiabfs5aejabeY7aTnaaCaaabeqaaiabd2gaTbaaaaaaaaaa@8815@
(28)

as elements of the matrix CΔμ= -k B T (dg Δ μ)-1 STBS (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

( ϵ Δ μ ) ′ = ( dg  Δ μ ) ϵ Δ μ ( dg  x ) − 1 = k B T S T ( dg  x ) − 1 , ( ε Δ μ ) ′ = ( dg  Δ μ ) ε Δ μ ( dg  x ) − 1 = k B T S T B ( dg  B ) − 1 ( dg  x ) − 1 , ( C Δ μ ) ′ = ( dg  Δ μ ) C Δ μ ( dg  J ) − 1 = − k B T S T B S . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabmqaaaqaaiabcIcaOmbvLH1qn1uy0Hws0fgBPngaryWyT1wAXadaiqWacqWF1pGSdaahaaWcbeqaaiabfs5aejabeY7aTbaakiqbcMcaPyaafaGaeyypa0JaeiikaGIaeeizaqMaee4zaCMaeeiiaaccceGae4hLdqKaeqiVd0MaeiykaKIae8x9di7aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccqqGGaaicqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWH4baEcqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabg2da9iabdUgaRnaaBaaaleaacqWGcbGqaeqaaOGaemivaqLaeC4uam1aaWbaaSqabeaacqWGubavaaGccqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWH4baEcqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabcYcaSaqaaiabcIcaOGGadiab9v7aLnaaCaaaleqabaGaeuiLdqKaeqiVd0gaaOGafiykaKIbauaacqGH9aqpcqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqGFuoarcqaH8oqBcqGGPaqkcqqF1oqzdaahaaWcbeqaaiabfs5aejabeY7aTbaakiabbccaGiabcIcaOiabbsgaKjabbEgaNjabbccaGiabhIha4jabcMcaPmaaCaaaleqabaGaeyOeI0IaeGymaedaaOGaeyypa0Jaem4AaS2aaSbaaSqaaiabdkeacbqabaGccqWGubavcqWHtbWudaahaaWcbeqaaiabdsfaubaakiabhkeacjabbccaGiabcIcaOiabbsgaKjabbEgaNjabbccaGiabhkeacjabcMcaPmaaCaaaleqabaGaeyOeI0IaeGymaedaaOGaeiikaGIaeeizaqMaee4zaCMaeeiiaaIaeCiEaGNaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqGGSaalaeaacqGGOaakcqWHdbWqdaahaaWcbeqaaiabfs5aejabeY7aTbaakiqbcMcaPyaafaGaeyypa0JaeiikaGIaeeizaqMaee4zaCMaeeiiaaIae4hLdqKaeqiVd0MaeiykaKIaeC4qam0aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccqqGGaaicqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWHkbGscqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabg2da9iabgkHiTiabdUgaRnaaBaaaleaacqWGcbGqaeqaaOGaemivaqLaeC4uam1aaWbaaSqabeaacqWGubavaaGccqWHcbGqcqWHtbWucqGGUaGlaaaaaa@C64C@

Note the appearance of the matrices ST, STB, and STBS in the local and steady-state elasticity-coefficients and the control-coefficient for the potentials, as compared to the matricies RT, RTB, and RTBS in the coefficients for the fluxes. Furthermore, note the relationships

C Δ μ = ϵ Δ μ C ^ , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeC4qam0aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccqGH9aqptqvzynutnfgDOLeDHXwAJbqegmwBTLwmWaacemGae8x9di7aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccuWHdbWqgaqcaiabcYcaSaaa@40D5@
(29)
( C Δ μ ) ′ = ( ϵ Δ μ ) ′ C ^ ′ . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbaGae8hkaGIaeC4qam0aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccuGGPaqkgaqbaiabg2da9iabcIcaOmbvLH1qn1uy0Hws0fgBPngaryWyT1wAXadaiqWacqGF1pGSdaahaaWcbeqaaiabfs5aejabeY7aTbaakiqb=LcaPyaafaGafC4qamKbaKGbauaacqWFUaGlaaa@445F@
(30)

Eqs. (19) and (29) in combination yield the summation [25] and connectivity theorems for the scaled coefficients, under consideration of the steady-state condition,

CΔμ1M×1= 0M×1,

CΔμϵ = -ϵΔμ,   CΔμε= -εΔμ.

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

(CΔμ)' K = 0M×(M-N),

(CΔμ)' ϵ '= - (ϵΔμ)',   (CΔμ)' ε '= - (εΔμ)'.

Eq. (33) is a manifestation of Kirchhoff's loop law [15, 16].

By combining the generalized summation and connectivity theorems, we get

( ( C Δ μ ) ′ C ^ ′ ) ( K ϵ ′ ) = ( 0 M × ( M − N ) − ( ϵ Δ μ ) ′ 0 N × ( M − N ) − I ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaeWaaeaafaqabeGabaaabaGaeiikaGIaeC4qam0aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccuGGPaqkgaqbaaqaaiqbhoeadzaajyaafaaaaaGaayjkaiaawMcaamaabmaabaqbaeqabeGaaaqaaiabhUealbqaambvLH1qn1uy0Hws0fgBPngaryWyT1wAXadaiqWacuWF1pGSgaqbaaaaaiaawIcacaGLPaaacqGH9aqpdaqadaqaauaabaqaciaaaeaacqWHWaamdaWgaaWcbaGaemyta0Kaey41aqRaeiikaGIaemyta0KaeyOeI0IaemOta4KaeiykaKcabeaaaOqaaiabgkHiTiabcIcaOiab=v=aYoaaCaaaleqabaGaeuiLdqKaeqiVd0gaaOGafiykaKIbauaaaeaacqWHWaamdaWgaaWcbaGaemOta4Kaey41aqRaeiikaGIaemyta0KaeyOeI0IaemOta4KaeiykaKcabeaaaOqaaiabgkHiTiabhMeajbaaaiaawIcacaGLPaaacqGGUaGlaaa@6205@
(35)

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 steady-state elasticity-coefficients, we have

( ( C Δ μ ) ′ C ^ ′ ) ( K ε ′ ) = ( 0 M × ( M − N ) − ( ε Δ μ ) ′ 0 N × ( M − N ) − ( dg  x ) B ( dg  B ) − 1 ( dg  x ) − 1 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaeWaaeaafaqabeGabaaabaGaeiikaGIaeC4qam0aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccuGGPaqkgaqbaaqaaiqbhoeadzaajyaafaaaaaGaayjkaiaawMcaamaabmaabaqbaeqabeGaaaqaaiabhUealbqaaGGadiqb=v7aLzaafaaaaaGaayjkaiaawMcaaiabg2da9maabmaabaqbaeaabiGaaaqaaiabhcdaWmaaBaaaleaacqWGnbqtcqGHxdaTcqGGOaakcqWGnbqtcqGHsislcqWGobGtcqGGPaqkaeqaaaGcbaGaeyOeI0IaeiikaGIae8xTdu2aaWbaaSqabeaacqqHuoarcqaH8oqBaaGccuGGPaqkgaqbaaqaaiabhcdaWmaaBaaaleaacqWGobGtcqGHxdaTcqGGOaakcqWGnbqtcqGHsislcqWGobGtcqGGPaqkaeqaaaGcbaGaeyOeI0IaeiikaGIaeeizaqMaee4zaCMaeeiiaaIaeCiEaGNaeiykaKIaeeiiaaIaeCOqaiecbeGae4hiaaIaeiikaGIaeeizaqMaee4zaCMaeeiiaaIaeCOqaiKaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWH4baEcqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaaaaaakiaawIcacaGLPaaacqGGUaGlaaa@7248@
(36)

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 { x i ∗ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdMgaPbqaaiabgEHiQaaaaaa@2FC5@ } fluctuates in response to a random stochastic perturbation to the concentration level of an external species x N ∈ I MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbiGae8xsaKeaaa@2CF7@ . 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),

δ x i / x i ∗ = B i N δ x N / B N N x N ∗ . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiTdqMaemiEaG3aaSbaaSqaaiabdMgaPbqabaGccqGGVaWlcqWG4baEdaqhaaWcbaGaemyAaKgabaGaey4fIOcaaOGaeyypa0JaemOqai0aaSbaaSqaaiabdMgaPjabd6eaobqabaGccqaH0oazcqWG4baEdaWgaaWcbaGaemOta4eabeaakiabc+caViabdkeacnaaBaaaleaacqWGobGtcqWGobGtaeqaaOGaemiEaG3aa0baaSqaaiabd6eaobqaaiabgEHiQaaakiabc6caUaaa@480F@

From the time series of {δ x i } with an uncorrelated perturbation on δ x N , we have

〈 δ x i δ x j 〉 〈 ( δ x i ) 2 〉 〈 ( δ x j ) 2 〉 = 1 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqGHPms4cqaH0oazcqWG4baEdaWgaaqaaiabdMgaPbqabaGaeqiTdqMaemiEaG3aaSbaaeaacqWGQbGAaeqaaiabgQYiXdqaamaakaaabaGaeyykJeUaeiikaGIaeqiTdqMaemiEaG3aaSbaaeaacqWGPbqAaeqaaiabcMcaPmaaCaaabeqaaiabikdaYaaacqGHQms8cqGHPms4cqGGOaakcqaH0oazcqWG4baEdaWgaaqaaiabdQgaQbqabaGaeiykaKYaaWbaaeqabaGaeGOmaidaaiabgQYiXdqabaaaaOGaeyypa0JaeGymaedaaa@5152@
(37)

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

d d t ( δ x N ) = − k d δ x N + a ξ ( t ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGKbazaeaacqWGKbazcqWG0baDaaGaeiikaGIaeqiTdqMaemiEaG3aaSbaaeaacqWGobGtaeqaaiabcMcaPiabg2da9iabgkHiTiabdUgaRnaaBaaabaGaemizaqgabeaacqaH0oazcqWG4baEdaWgaaqaaiabd6eaobqabaGaey4kaSIaemyyaeMaeqOVdGNaeiikaGIaemiDaqNaeiykaKIaeiilaWcaaa@4805@
(38)

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 Ornstein-Uhlenbeck (OU) process with auto-covariance function

g ( τ ) ≜ 〈 δ x N ( t ) δ x N ( t + τ ) 〉 = a 2 2 k d e − k d | τ | . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4zaCMaeiikaGIaeqiXdqNaeiykaKIaeSixIaKaeyykJeUaeqiTdqMaemiEaG3aaSbaaSqaaiabd6eaobqabaGccqGGOaakcqWG0baDcqGGPaqkcqaH0oazcqWG4baEdaWgaaWcbaGaemOta4eabeaakiabcIcaOiabdsha0jabgUcaRiabes8a0jabcMcaPiabgQYiXlabg2da9KqbaoaalaaabaGaemyyae2aaWbaaeqabaGaeGOmaidaaaqaaiabikdaYiabdUgaRnaaBaaabaGaemizaqgabeaaaaGccqWGLbqzdaahaaWcbeqaaiabgkHiTiabdUgaRnaaBaaameaacqWGKbazaeqaaSGaeiiFaWNaeqiXdqNaeiiFaWhaaOGaeiOla4caaa@5A7A@
(39)

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)

d d t ( δ x i ) = ∑ j = 1 N − 1 ( A i j x j ∗ ) δ x j + A i N x N ∗ δ x N ( t ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGKbazaeaacqWGKbazcqWG0baDaaGccqGGOaakcqaH0oazcqWG4baEdaWgaaWcbaGaemyAaKgabeaakiabcMcaPiabg2da9maaqahabaWaaeWaaKqbagaadaWcaaqaaiabdgeabnaaBaaabaGaemyAaKMaemOAaOgabeaaaeaacqWG4baEdaqhaaqaaiabdQgaQbqaaiabgEHiQaaaaaaakiaawIcacaGLPaaaaSqaaiabdQgaQjabg2da9iabigdaXaqaaiabd6eaojabgkHiTiabigdaXaqdcqGHris5aOGaeqiTdqMaemiEaG3aaSbaaSqaaiabdQgaQbqabaGccqGHRaWkjuaGdaWcaaqaaiabdgeabnaaBaaabaGaemyAaKMaemOta4eabeaaaeaacqWG4baEdaqhaaqaaiabd6eaobqaaiabgEHiQaaaaaGccqaH0oazcqWG4baEdaWgaaWcbaGaemOta4eabeaakiabcIcaOiabdsha0jabcMcaPaaa@6014@
(40)

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, ...,

δ x i ( n + 1 ) = ∑ j = 1 N ( e A ˜ ˜ T ) i j δ x j ( n ) − B i N a ξ n , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiTdqMaemiEaG3aaSbaaSqaaiabdMgaPbqabaGccqGGOaakcqWGUbGBcqGHRaWkcqaIXaqmcqGGPaqkcqGH9aqpdaaeWbqaamaabmaabaGaemyzau2aaWbaaSqabeaacuWHbbqqgaacgaacaiabdsfaubaaaOGaayjkaiaawMcaamaaBaaaleaacqWGPbqAcqWGQbGAaeqaaaqaaiabdQgaQjabg2da9iabigdaXaqaaiabd6eaobqdcqGHris5aOGaeqiTdqMaemiEaG3aaSbaaSqaaiabdQgaQbqabaGccqGGOaakcqWGUbGBcqGGPaqkcqGHsislcqWGcbGqdaWgaaWcbaGaemyAaKMaemOta4eabeaakiabdggaHjabe67a4naaBaaaleaacqWGUbGBaeqaaOGaeiilaWcaaa@57CF@

in which A ˜ ˜ i j = A i j / x j ∗ − δ i N k d MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmyqaeKbaGGbaGaadaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabg2da9iabdgeabnaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaei4la8IaemiEaG3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaakiabgkHiTiabes7aKnaaBaaaleaacqWGPbqAcqWGobGtaeqaaOGaem4AaS2aaSbaaSqaaiabdsgaKbqabaaaaa@41F0@ . 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 A ˜ = { A ˜ i j = A i j / x j ∗ , i , j = 1 , 2 , ... , N − 1 } MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafCyqaeKbaGaacqGH9aqpcqGG7bWEcuWGbbqqgaacamaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeyypa0Jaemyqae0aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGGVaWlcqWG4baEdaqhaaWcbaGaemOAaOgabaGaey4fIOcaaOGaeiilaWIaemyAaKMaeiilaWIaemOAaOMaeyypa0JaeGymaeJaeiilaWIaeGOmaiJaeiilaWIaeiOla4IaeiOla4IaeiOla4IaeiilaWIaemOta4KaeyOeI0IaeGymaeJaeiyFa0haaa@4E90@ . Then the stationary solution to Eqs. (38) and (40) is

δ x i ( t ) = ∑ ℓ = 1 N − 1 ∫ 0 t ( e A ˜ ( t − η ) ) i ℓ A ℓ N x N ∗ δ x N ( η ) d η , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiTdqMaemiEaG3aaSbaaSqaaiabdMgaPbqabaGccqGGOaakcqWG0baDcqGGPaqkcqGH9aqpdaaeWbqaamaapedabaWaaeWaaeaacqWGLbqzdaahaaWcbeqaaiqbhgeabzaaiaGaeiikaGIaemiDaqNaeyOeI0Iaeq4TdGMaeiykaKcaaaGccaGLOaGaayzkaaWaaSbaaSqaaiabdMgaPjabloriSbqabaqcfa4aaSaaaeaacqWGbbqqdaWgaaqaaiabloriSjabd6eaobqabaaabaGaemiEaG3aa0baaeaacqWGobGtaeaacqGHxiIkaaaaaOGaeqiTdqMaemiEaG3aaSbaaSqaaiabd6eaobqabaGccqGGOaakcqaH3oaAcqGGPaqkcqWGKbazcqaH3oaAaSqaaiabicdaWaqaaiabdsha0bqdcqGHRiI8aaWcbaGaeS4eHWMaeyypa0JaeGymaedabaGaemOta4KaeyOeI0IaeGymaedaniabggHiLdGccqGGSaalaaa@622A@
(41)

which leads to

〈 δ x i ( t ) δ x j ( t + τ ) 〉 = ∑ ℓ , m = 1 N − 1 ∫ 0 ∞ ∫ 0 ∞ ( e A ˜ η ) i ℓ A ℓ N A m N x N ∗ x N ∗ ( e A ˜ ζ ) j m g ( τ + η − ζ ) d η d ζ . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeyykJeUaeqiTdqMaemiEaG3aaSbaaSqaaiabdMgaPbqabaGccqGGOaakcqWG0baDcqGGPaqkcqaH0oazcqWG4baEdaWgaaWcbaGaemOAaOgabeaakiabcIcaOiabdsha0jabgUcaRiabes8a0jabcMcaPiabgQYiXlabg2da9maaqahabaWaa8qmaeaadaWdXaqaamaabmaabaGaemyzau2aaWbaaSqabeaacuWHbbqqgaacaiabeE7aObaaaOGaayjkaiaawMcaamaaBaaaleaacqWGPbqAcqWItecBaeqaaKqbaoaalaaabaGaemyqae0aaSbaaeaacqWItecBcqWGobGtaeqaaiabdgeabnaaBaaabaGaemyBa0MaemOta4eabeaaaeaacqWG4baEdaqhaaqaaiabd6eaobqaaiabgEHiQaaacqWG4baEdaqhaaqaaiabd6eaobqaaiabgEHiQaaaaaGcdaqadaqaaiabdwgaLnaaCaaaleqabaGafCyqaeKbaGaacqaH2oGEaaaakiaawIcacaGLPaaadaWgaaWcbaGaemOAaOMaemyBa0gabeaakiabdEgaNjabcIcaOiabes8a0jabgUcaRiabeE7aOjabgkHiTiabeA7a6jabcMcaPiabdsgaKjabeE7aOjabdsgaKjabeA7a6bWcbaGaeGimaadabaGaeyOhIukaniabgUIiYdaaleaacqaIWaamaeaacqGHEisPa0Gaey4kIipakiabc6caUaWcbaGaeS4eHWMaeiilaWIaemyBa0Maeyypa0JaeGymaedabaGaemOta4KaeyOeI0IaeGymaedaniabggHiLdaaaa@886B@
(42)

Two limits of Eq. (42) are particularly interesting. If the correlation time of g(τ) in Eq. (39) is ≫ all the relaxation times of A ˜ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafCyqaeKbaGaaaaa@2CF3@ , then we have Eq. (37). On the other hand, if the correlation time for the random perturbation is ≪ all the relaxation times of A ˜ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafCyqaeKbaGaaaaa@2CF3@ , we have

〈 δ x i ( t ) δ x j ( t + τ ) 〉 = ∑ ℓ = 1 N − 1 〈 δ x i δ x ℓ 〉 ( e A ˜ τ ) j ℓ . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeyykJeUaeqiTdqMaemiEaG3aaSbaaSqaaiabdMgaPbqabaGccqGGOaakcqWG0baDcqGGPaqkcqaH0oazcqWG4baEdaWgaaWcbaGaemOAaOgabeaakiabcIcaOiabdsha0jabgUcaRiabes8a0jabcMcaPiabgQYiXlabg2da9maaqahabaGaeyykJeUaeqiTdqMaemiEaG3aaSbaaSqaaiabdMgaPbqabaGccqaH0oazcqWG4baEdaWgaaWcbaGaeS4eHWgabeaakiabgQYiXdWcbaGaeS4eHWMaeyypa0JaeGymaedabaGaemOta4KaeyOeI0IaeGymaedaniabggHiLdGcdaqadaqaaiabdwgaLnaaCaaaleqabaGafCyqaeKbaGaacqaHepaDaaaakiaawIcacaGLPaaadaWgaaWcbaGaemOAaOMaeS4eHWgabeaakiabc6caUaaa@61C3@
(43)

The matrix R ( t ) = e A ˜ t MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCOuaiLaeiikaGIaemiDaqNaeiykaKIaeyypa0Jaemyzau2aaWbaaSqabeaacuWHbbqqgaacaiabdsha0baaaaa@353E@ is the fundamental solution to the linear kinetic equation (4). R ij (t) is the response curve, as a function of t, for the i th species to an impulse of the j th 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 R ′ i j ( 0 ) ≠ 0 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmOuaiLbauaadaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabcIcaOiabicdaWiabcMcaPiabgcMi5kabicdaWaaa@3551@ , whereas if i and j are not directly connected, R ′ i j ( 0 ) = 0 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmOuaiLbauaadaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabcIcaOiabicdaWiabcMcaPiabg2da9iabicdaWaaa@3490@ . Analytic solutions for completely linear systems have been derived and analyzed [29].

Conclusion

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

X 0 ⇌ k − 1 k 1 X 1 ⇌ k − 2 k 2 X 2 ⇌ k − 3 k 3 ⋯ ⇌ X n − 2 ⇌ k − ( n − 1 ) k n − 1 X n − 1 ⇌ k − n k n X n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiwaG1aaSbaaSqaaiabicdaWaqabaGcdaWfWaqaaiablYCidcWcbaGaem4AaS2aaSbaaWqaaiabgkHiTiabigdaXaqabaaaleaacqWGRbWAdaWgaaadbaGaeGymaedabeaaaaGccqWGybawdaWgaaWcbaGaeGymaedabeaakmaaxadabaGaeSiZHmialeaacqWGRbWAdaWgaaadbaGaeyOeI0IaeGOmaidabeaaaSqaaiabdUgaRnaaBaaameaacqaIYaGmaeqaaaaakiabdIfaynaaBaaaleaacqaIYaGmaeqaaOWaaCbmaeaacqWImhYGaSqaaiabdUgaRnaaBaaameaacqGHsislcqaIZaWmaeqaaaWcbaGaem4AaS2aaSbaaWqaaiabiodaZaqabaaaaOGaeS47IWKaeSiZHmOaemiwaG1aaSbaaSqaaiabd6gaUjabgkHiTiabikdaYaqabaGcdaWfWaqaaiablYCidcWcbaGaem4AaS2aaSbaaWqaaiabgkHiTiabcIcaOiabd6gaUjabgkHiTiabigdaXiabcMcaPaqabaaaleaacqWGRbWAdaWgaaadbaGaemOBa4MaeyOeI0IaeGymaedabeaaaaGccqWGybawdaWgaaWcbaGaemOBa4MaeyOeI0IaeGymaedabeaakmaaxadabaGaeSiZHmialeaacqWGRbWAdaWgaaadbaGaeyOeI0IaemOBa4gabeaaaSqaaiabdUgaRnaaBaaameaacqWGUbGBaeqaaaaakiabdIfaynaaBaaaleaacqWGUbGBaeqaaaaa@7448@
(44)

under NESS with clamped concentrations x0 and x n for X0 and X n , respectively. Analogous to a continuous energy landscape, we have the discrete energy function

E j = − k B T ln ( k 1 k 2 ⋯ k j k − 1 k − 2 ⋯ k − j ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyrau0aaSbaaSqaaiabdQgaQbqabaGccqGH9aqpcqGHsislcqWGRbWAdaWgaaWcbaGaemOqaieabeaakiabdsfaujGbcYgaSjabc6gaUnaabmaajuaGbaWaaSaaaeaacqWGRbWAdaWgaaqaaiabigdaXaqabaGaem4AaS2aaSbaaeaacqaIYaGmaeqaaiabl+UimjabdUgaRnaaBaaabaGaemOAaOgabeaaaeaacqWGRbWAdaWgaaqaaiabgkHiTiabigdaXaqabaGaem4AaS2aaSbaaeaacqGHsislcqaIYaGmaeqaaiabl+UimjabdUgaRnaaBaaabaGaeyOeI0IaemOAaOgabeaaaaaakiaawIcacaGLPaaacqGGSaalaaa@5087@

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]

J = k 1 k 2 ⋯ k n k − 1 k − 2 ⋯ k − n x 0 − x n 1 k − n + k n k − n k − ( n − 1 ) + k n k n − 1 k − n k − ( n − 1 ) k − ( n − 2 ) + … + k n k n − 1 ⋯ k 2 k − n k − ( n − 1 ) ⋯ k − 1 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaOKaeyypa0tcfa4aaSaaaeaadaWcaaqaaiabdUgaRnaaBaaabaGaeGymaedabeaacqWGRbWAdaWgaaqaaiabikdaYaqabaGaeS47IWKaem4AaS2aaSbaaeaacqWGUbGBaeqaaaqaaiabdUgaRnaaBaaabaGaeyOeI0IaeGymaedabeaacqWGRbWAdaWgaaqaaiabgkHiTiabikdaYaqabaGaeS47IWKaem4AaS2aaSbaaeaacqGHsislcqWGUbGBaeqaaaaacqWG4baEdaWgaaqaaiabicdaWaqabaGaeyOeI0IaemiEaG3aaSbaaeaacqWGUbGBaeqaaaqaamaalaaabaGaeGymaedabaGaem4AaS2aaSbaaeaacqGHsislcqWGUbGBaeqaaaaacqGHRaWkdaWcaaqaaiabdUgaRnaaBaaabaGaemOBa4gabeaaaeaacqWGRbWAdaWgaaqaaiabgkHiTiabd6gaUbqabaGaem4AaS2aaSbaaeaacqGHsislcqGGOaakcqWGUbGBcqGHsislcqaIXaqmcqGGPaqkaeqaaaaacqGHRaWkdaWcaaqaaiabdUgaRnaaBaaabaGaemOBa4gabeaacqWGRbWAdaWgaaqaaiabd6gaUjabgkHiTiabigdaXaqabaaabaGaem4AaS2aaSbaaeaacqGHsislcqWGUbGBaeqaaiabdUgaRnaaBaaabaGaeyOeI0IaeiikaGIaemOBa4MaeyOeI0IaeGymaeJaeiykaKcabeaacqWGRbWAdaWgaaqaaiabgkHiTiabcIcaOiabd6gaUjabgkHiTiabikdaYiabcMcaPaqabaaaaiabgUcaRiablAciljabgUcaRmaalaaabaGaem4AaS2aaSbaaeaacqWGUbGBaeqaaiabdUgaRnaaBaaabaGaemOBa4MaeyOeI0IaeGymaedabeaacqWIVlctcqWGRbWAdaWgaaqaaiabikdaYaqabaaabaGaem4AaS2aaSbaaeaacqGHsislcqWGUbGBaeqaaiabdUgaRnaaBaaabaGaeyOeI0IaeiikaGIaemOBa4MaeyOeI0IaeGymaeJaeiykaKcabeaacqWIVlctcqWGRbWAdaWgaaqaaiabgkHiTiabigdaXaqabaaaaaaacqGGUaGlaaa@9B03@
(45)

Therefore, the flux control-coefficient is

C ℓ J = E ℓ J ( δ J δ E ℓ ) = k n k n − 1 ⋯ k ℓ + 1 k − n k − n ( n − 1 ) ⋯ k − ( ℓ + 1 ) k − ℓ 1 k − n + k n k − n k − ( n − 1 ) + k n k n − 1 k − n k − ( n − 1 ) k − ( n − 2 ) + … + k n k n − 1 ⋯ k 2 k − n k − ( n − 1 ) ⋯ k − 1 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiWaaaqaaiabdoeadnaaDaaaleaacqWItecBaeaacqWGkbGsaaaakeaacqGH9aqpaeaajuaGdaWcaaqaaiabdweafnaaCaaabeqaaiabloriSbaaaeaacqWGkbGsaaGcdaqadaqcfayaamaalaaabaGaeqiTdqMaemOsaOeabaGaeqiTdqMaemyrau0aaWbaaeqabaGaeS4eHWgaaaaaaOGaayjkaiaawMcaaaqaaaqaaiabg2da9aqaaKqbaoaalaaabaWaaSaaaeaacqWGRbWAdaWgaaqaaiabd6gaUbqabaGaem4AaS2aaSbaaeaacqWGUbGBcqGHsislcqaIXaqmaeqaaiabl+UimjabdUgaRnaaBaaabaGaeS4eHWMaey4kaSIaeGymaedabeaaaeaacqWGRbWAdaWgaaqaaiabgkHiTiabd6gaUbqabaGaem4AaS2aaSbaaeaacqGHsislcqWGUbGBcqGGOaakcqWGUbGBcqGHsislcqaIXaqmcqGGPaqkaeqaaiabl+UimjabdUgaRnaaBaaabaGaeyOeI0IaeiikaGIaeS4eHWMaey4kaSIaeGymaeJaeiykaKcabeaacqWGRbWAdaWgaaqaaiabgkHiTiabloriSbqabaaaaaqaamaalaaabaGaeGymaedabaGaem4AaS2aaSbaaeaacqGHsislcqWGUbGBaeqaaaaacqGHRaWkdaWcaaqaaiabdUgaRnaaBaaabaGaemOBa4gabeaaaeaacqWGRbWAdaWgaaqaaiabgkHiTiabd6gaUbqabaGaem4AaS2aaSbaaeaacqGHsislcqGGOaakcqWGUbGBcqGHsislcqaIXaqmcqGGPaqkaeqaaaaacqGHRaWkdaWcaaqaaiabdUgaRnaaBaaabaGaemOBa4gabeaacqWGRbWAdaWgaaqaaiabd6gaUjabgkHiTiabigdaXaqabaaabaGaem4AaS2aaSbaaeaacqGHsislcqWGUbGBaeqaaiabdUgaRnaaBaaabaGaeyOeI0IaeiikaGIaemOBa4MaeyOeI0IaeGymaeJaeiykaKcabeaacqWGRbWAdaWgaaqaaiabgkHiTiabcIcaOiabd6gaUjabgkHiTiabikdaYiabcMcaPaqabaaaaiabgUcaRiablAciljabgUcaRmaalaaabaGaem4AaS2aaSbaaeaacqWGUbGBaeqaaiabdUgaRnaaBaaabaGaemOBa4MaeyOeI0IaeGymaedabeaacqWIVlctcqWGRbWAdaWgaaqaaiabikdaYaqabaaabaGaem4AaS2aaSbaaeaacqGHsislcqWGUbGBaeqaaiabdUgaRnaaBaaabaGaeyOeI0IaeiikaGIaemOBa4MaeyOeI0IaeGymaeJaeiykaKcabeaacqWIVlctcqWGRbWAdaWgaaqaaiabgkHiTiabigdaXaqabaaaaaaacqGGSaalaaaaaa@B673@
(46)

which reaches its maximum at maximal E j . The biochemical conductance, on the other hand, is

| J Δ μ ℓ | = J k B T ( ln ϕ + ℓ ϕ + ℓ − J ) − 1 = J k B T ( ln J + ϕ − ℓ ϕ − ℓ ) − 1 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaqWaaKqbagaadaWcaaqaaiabdQeakbqaaiabfs5aejabeY7aTnaaBaaabaGaeS4eHWgabeaaaaaakiaawEa7caGLiWoacqGH9aqpjuaGdaWcaaqaaiabdQeakbqaaiabdUgaRnaaBaaabaGaemOqaieabeaacqWGubavaaWaaeWaaeaacyGGSbaBcqGGUbGBdaWcaaqaaiabew9aMnaaDaaabaGaey4kaScabaGaeS4eHWgaaaqaaiabew9aMnaaDaaabaGaey4kaScabaGaeS4eHWgaaiabgkHiTiabdQeakbaaaiaawIcacaGLPaaadaahaaqabeaacqGHsislcqaIXaqmaaGaeyypa0ZaaSaaaeaacqWGkbGsaeaacqWGRbWAdaWgaaqaaiabdkeacbqabaGaemivaqfaamaabmaabaGagiiBaWMaeiOBa42aaSaaaeaacqWGkbGscqGHRaWkcqaHvpGzdaqhaaqaaiabgkHiTaqaaiabloriSbaaaeaacqaHvpGzdaqhaaqaaiabgkHiTaqaaiabloriSbaaaaaacaGLOaGaayzkaaWaaWbaaeqabaGaeyOeI0IaeGymaedaaiabcYcaSaaa@6397@

which is at its minimum when ϕ + ℓ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqy1dy2aa0baaSqaaiabgUcaRaqaaiabloriSbaaaaa@2FDD@ (or ϕ − ℓ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqy1dy2aa0baaSqaaiabgkHiTaqaaiabloriSbaaaaa@2FE8@ ) reaches its minimum. For small flux, ϕ + ℓ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqy1dy2aa0baaSqaaiabgUcaRaqaaiabloriSbaaaaa@2FDD@ (and ϕ − ℓ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqy1dy2aa0baaSqaaiabgkHiTaqaaiabloriSbaaaaa@2FE8@ ) 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 constant-current (zero internal conductance) and constant-voltage (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)

h d r = − J Δ μ = k B T J ln ( 1 + J ϕ − ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiAaGMaemizaqMaemOCaiNaeyypa0JaeyOeI0IaemOsaOKaeuiLdqKaeqiVd0Maeyypa0Jaem4AaS2aaSbaaSqaaiabdkeacbqabaGccqWGubavcqWGkbGscyGGSbaBcqGGUbGBdaqadaqaaiabigdaXiabgUcaRKqbaoaalaaabaGaemOsaOeabaGaeqy1dy2aaSbaaeaacqGHsislaeqaaaaaaOGaayjkaiaawMcaaiabc6caUaaa@47FC@
(47)

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

h d r = Δ μ ϕ − ( 1 − e − Δ μ / k B T ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiAaGMaemizaqMaemOCaiNaeyypa0JaeuiLdqKaeqiVd0Maeqy1dy2aaSbaaSqaaiabgkHiTaqabaGcdaqadaqaaiabigdaXiabgkHiTiabdwgaLnaaCaaaleqabaGaeyOeI0IaeuiLdqKaeqiVd0Maei4la8Iaem4AaS2aaSbaaWqaaiabdkeacbqabaWccqWGubavaaaakiaawIcacaGLPaaacqGGUaGlaaa@45DF@
(48)

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 I2 R and U2/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

J = α 0 ( x + a 1 ) ( x + a 2 ) ⋯ ( x + a K ) β 0 ( x + b 1 ) ( x + b 2 ) ⋯ ( x + b L ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaOKaeyypa0tcfa4aaSaaaeaacqaHXoqydaWgaaqaaiabicdaWaqabaGaeiikaGIaemiEaGNaey4kaSIaemyyae2aaSbaaeaacqaIXaqmaeqaaiabcMcaPiabcIcaOiabdIha4jabgUcaRiabdggaHnaaBaaabaGaeGOmaidabeaacqGGPaqkcqWIVlctcqGGOaakcqWG4baEcqGHRaWkcqWGHbqydaWgaaqaaiabdUealbqabaGaeiykaKcabaGaeqOSdi2aaSbaaeaacqaIWaamaeqaaiabcIcaOiabdIha4jabgUcaRiabdkgaInaaBaaabaGaeGymaedabeaacqGGPaqkcqGGOaakcqWG4baEcqGHRaWkcqWGIbGydaWgaaqaaiabikdaYaqabaGaeiykaKIaeS47IWKaeiikaGIaemiEaGNaey4kaSIaemOyai2aaSbaaeaacqWGmbataeqaaiabcMcaPaaacqGGSaalaaa@5FDE@
(49)

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 elasticity-coefficient for the flux J can be written as

ε = x x + a 1 + x x + a 2 + … + x x + a K − x x + b 1 − x x + b 2 − … − x x + b L , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyTduMaeyypa0tcfa4aaSaaaeaacqWG4baEaeaacqWG4baEcqGHRaWkcqWGHbqydaWgaaqaaiabigdaXaqabaaaaiabgUcaRmaalaaabaGaemiEaGhabaGaemiEaGNaey4kaSIaemyyae2aaSbaaeaacqaIYaGmaeqaaaaacqGHRaWkcqWIMaYscqGHRaWkdaWcaaqaaiabdIha4bqaaiabdIha4jabgUcaRiabdggaHnaaBaaabaGaem4saSeabeaaaaGaeyOeI0YaaSaaaeaacqWG4baEaeaacqWG4baEcqGHRaWkcqWGIbGydaWgaaqaaiabigdaXaqabaaaaiabgkHiTmaalaaabaGaemiEaGhabaGaemiEaGNaey4kaSIaemOyai2aaSbaaeaacqaIYaGmaeqaaaaacqGHsislcqWIMaYscqGHsisldaWcaaqaaiabdIha4bqaaiabdIha4jabgUcaRiabdkgaInaaBaaabaGaemitaWeabeaaaaGaeiilaWcaaa@5EC6@
(50)

which also has the standard form of a Pade approximation [35]. Furthermore, BSA proposed [22] to approximate (49) by a power-law expression

J ≈ C x 1 ε 1 x 2 ε 2 ⋯ x N ε N , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaOKaeyisISRaem4qamKaemiEaG3aa0baaSqaaiabigdaXaqaaiabew7aLnaaBaaameaacqaIXaqmaeqaaaaakiabdIha4naaDaaaleaacqaIYaGmaeaacqaH1oqzdaWgaaadbaGaeGOmaidabeaaaaGccqWIVlctcqWG4baEdaqhaaWcbaGaemOta4eabaGaeqyTdu2aaSbaaWqaaiabd6eaobqabaaaaOGaeiilaWcaaa@4368@
(51)

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.

Methods

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 KTK is invertible because K contains linearly independent columns, and because all columns of K are in the right nullspace of S, the matrix

( S K T ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaeWaaeaafaqaaeGabaaabaGaeC4uamfabaGaeC4saS0aaWbaaSqabeaaieGacqWFubavaaaaaaGccaGLOaGaayzkaaaaaa@317D@

is also invertible [23]. Therefore, we have

( K ε ′ ) ( ( K T K ) − 1 K T ( I − R T B S ) ( dg  x ) ( dg  B ) S ) = ( S K T ) − 1 ( S K T ) ( K ε ′ ) ( ( K T K ) − 1 K T ( I − R T B S ) ( dg  x ) ( dg  B ) S ) = ( S K T ) − 1 ( 0 N × ( M − N ) S ε ′ K T K K T ε ′ ) ( ( K T K ) − 1 K T ( I − R T B S ) ( dg  x ) ( dg  B ) S ) = ( S K T ) − 1 ( S K T ) = I . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabqqaaaaabaWaaeWaaeaafaqabeqacaaabaGaeC4saSeabaaccmGaf8xTduMbauaaaaaacaGLOaGaayzkaaWaaeWaaeaafaqabeGabaaabaGaeiikaGIaeC4saS0aaWbaaSqabeaacqWGubavaaGccqWHlbWscqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabhUealnaaCaaaleqabaGaemivaqfaaOGaeiikaGIaeCysaKKaeyOeI0IaeCOuai1aaWbaaSqabeaacqWGubavaaGccqWHcbGqcqWHtbWucqGGPaqkaeaacqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWH4baEcqGGPaqkcqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWHcbGqcqGGPaqkcqWHtbWuaaaacaGLOaGaayzkaaaabaGaeyypa0ZaaeWaaeaafaqaaeGabaaabaGaeC4uamfabaGaeC4saS0aaWbaaSqabeaaieGacqGFubavaaaaaaGccaGLOaGaayzkaaWaaWbaaSqabeaacqGHsislcqaIXaqmaaGcdaqadaqaauaabaqaceaaaeaacqWHtbWuaeaacqWHlbWsdaahaaWcbeqaaiab+rfaubaaaaaakiaawIcacaGLPaaadaqadaqaauaabeqabiaaaeaacqWHlbWsaeaacuaH1oqzgaqbaaaaaiaawIcacaGLPaaadaqadaqaauaabeqaceaaaeaacqGGOaakcqWHlbWsdaahaaWcbeqaaiabdsfaubaakiabhUealjabcMcaPmaaCaaaleqabaGaeyOeI0IaeGymaedaaOGaeC4saS0aaWbaaSqabeaacqWGubavaaGccqGGOaakcqWHjbqscqGHsislcqWHsbGudaahaaWcbeqaaiabdsfaubaakiabhkeacjabhofatjabcMcaPaqaaiabcIcaOiabbsgaKjabbEgaNjabbccaGiabhIha4jabcMcaPiabcIcaOiabbsgaKjabbEgaNjabbccaGiabhkeacjabcMcaPiabhofatbaaaiaawIcacaGLPaaaaeaacqGH9aqpdaqadaqaauaabaqaceaaaeaacqWHtbWuaeaacqWHlbWsdaahaaWcbeqaaiab+rfaubaaaaaakiaawIcacaGLPaaadaahaaWcbeqaaiabgkHiTiabigdaXaaakmaabmaabaqbaeaabiGaaaqaaiabhcdaWmaaBaaaleaacqWGobGtcqGHxdaTcqGGOaakcqWGnbqtcqGHsislcqWGobGtcqGGPaqkaeqaaaGcbaGaeC4uamLaf8xTduMbauaaaeaacqWHlbWsdaahaaWcbeqaaiabdsfaubaakiabhUealbqaaiabhUealnaaCaaaleqabaGaemivaqfaaOGaf8xTduMbauaaaaaacaGLOaGaayzkaaWaaeWaaeaafaqabeGabaaabaGaeiikaGIaeC4saS0aaWbaaSqabeaacqWGubavaaGccqWHlbWscqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabhUealnaaCaaaleqabaGaemivaqfaaOGaeiikaGIaeCysaKKaeyOeI0IaeCOuai1aaWbaaSqabeaacqWGubavaaGccqWHcbGqcqWHtbWucqGGPaqkaeaacqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWH4baEcqGGPaqkcqGGOaakcqqGKbazcqqGNbWzcqqGGaaicqWHcbGqcqGGPaqkcqWHtbWuaaaacaGLOaGaayzkaaaabaGaeyypa0ZaaeWaaeaafaqaaeGabaaabaGaeC4uamfabaGaeC4saS0aaWbaaSqabeaacqGFubavaaaaaaGccaGLOaGaayzkaaWaaWbaaSqabeaacqGHsislcqaIXaqmaaGcdaqadaqaauaabaqaceaaaeaacqWHtbWuaeaacqWHlbWsdaahaaWcbeqaaiab+rfaubaaaaaakiaawIcacaGLPaaacqGH9aqpcqWHjbqscqGGUaGlaaaaaa@D822@

Abbreviations

BCT:

Biochemical Circuit Theory

BSA:

Biochemical Systems Analysis

CMC:

Correlation Metric Construction

FBA:

Flux Balance Analysis

MCA:

Metabolic Control Analysis

NESS:

Nonequilibrium Steady-State

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 = SRT

S :

N ×M stoichiometric matrix

x :

N-dimensional vector of species' concentrations

J :

M-dimensional vector of reaction fluxes

Δ :

μ M-dimensional vector of reaction potentials

ϵ M ×N local:

elasticity-coefficient matrix

ω M ×N steady-state:

elasticity-coefficient matrix

C ^ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbeGaf83qamKbaKaaaaa@2CFA@ :

N × M concentration, control-coefficient matrix

C M ×M flux:

control-coefficient matrix

CΔμM ×M biochemical potential:

control-coefficient 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 steady-state

::

' Distinguishes unscaled coefficients from scaled coefficients

References

  1. Fell DA, Sauro HM: Metabolic control and its analysis. Eur J Biochem. 1985, 148: 555-561. 10.1111/j.1432-1033.1985.tb08876.x

    Article  CAS  PubMed  Google Scholar 

  2. Kholodenko BN, Westerhoff HV: The macroworld versus the microworld of biochemical regulation and control. TIBS. 1995, 20: 52-54.

    CAS  PubMed  Google Scholar 

  3. Edwards JS, Palsson BØ: The E. Coli MG1655 in silico metabolic geotype: its definition, characteristics, and capabilities. Proc Natl Acad Sci USA. 2000, 97: 5528-5533. 10.1073/pnas.97.10.5528

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  4. Schilling CH, Palsson BØ: The underlying pathway structure of biochemical reaction networks. Proc Natl Acad Sci USA. 1998, 95: 4193-4198. 10.1073/pnas.95.8.4193

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  5. Arkin A, Ross J: Statistical construction of chemical reaction mechanisms from measured time-series. J Phys Chem. 1995, 99: 970-979. 10.1021/j100003a020.

    Article  CAS  Google Scholar 

  6. Samoilov M, Arkin A, Ross J: On the deduction of chemical reaction pathways from measurements of time series of concentrations. Chaos. 2001, 11: 108-114. 10.1063/1.1336499

    Article  CAS  PubMed  Google Scholar 

  7. Qian H: Phosphorylation energy hypothesis: open chemical systems and their biological functions. Ann Rev Phys Chem. 2007, 58: 113-142. 10.1146/annurev.physchem.58.032806.104550.

    Article  CAS  Google Scholar 

  8. Qian H: Open-system nonequilibrium steady-state: statistical thermodynamics, fluctuations and chemical oscillations. J Phys Chem B. 2006, 110: 15063-15074. 10.1021/jp061858z

    Article  CAS  PubMed  Google Scholar 

  9. Wyman J: The turning wheel: a study in steady state. Proc Natl Acad Sci USA. 1975, 72: 3983-3987. 10.1073/pnas.72.10.3983

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  10. Hill TL: Free Energy Transduction in Biology: The Steady-State Kinetic and Thermodynamic Formalism. 1977, New York: Academic Press

    Google Scholar 

  11. Hill TL: Free Energy Transduction and Biochemical Cycle Kinetics. 1989, New York: Springer-Verlag

    Book  Google Scholar 

  12. Qian H: The Mathematical theory of molecular motor movement and chemomechanical energy transduction. J Math Chem. 2000, 27: 219-234. 10.1023/A:1026428320489.

    Article  CAS  Google Scholar 

  13. Qian H: Cycle kinetics, steady-state thermodynamics and motors – a paradigm for living matter physics. J Phys: Condens Matter. 2005, 17: S3783-S3794. 10.1088/0953-8984/17/47/010.

    CAS  Google Scholar 

  14. Murray JD: Mathematical Biology. 1991, New York: Springer, 2

    Google Scholar 

  15. Beard DA, Liang SD, Qian H: Energy balance for analysis of complex metabolic networks. Biophys J. 2002, 83: 79-86.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  16. Qian H, Beard DA, Liang SD: Stoichiometric network theory for nonequilibrium biochemical systems. Eur J Biochem. 2003, 270: 415-421. 10.1046/j.1432-1033.2003.03357.x

    Article  CAS  PubMed  Google Scholar 

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

    Chapter  Google Scholar 

  18. Segel IH: Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady-state Enzyme Systems. 1975, New York: John Wiley & Sons

    Google Scholar 

  19. Alberty RA: Degrees of freedom in biochemical reaction systems at specific pH and pMg. J Phys Chem. 1992, 96: 9614-9621. 10.1021/j100203a012.

    Article  CAS  Google Scholar 

  20. Poland D: On the stability of mass action reactions in open systems. J Chem Phys. 1991, 94: 4427-4439. 10.1063/1.460633.

    Article  CAS  Google Scholar 

  21. Savageau MA: Biochemical system analysis I. J Theor Biol. 1969, 25: 365-369. 10.1016/S0022-5193(69)80026-3

    Article  CAS  PubMed  Google Scholar 

  22. Savageau MA: Biochemical system analysis II. J Theor Biol. 1969, 25: 370-379. 10.1016/S0022-5193(69)80027-5

    Article  CAS  PubMed  Google Scholar 

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

    Book  Google Scholar 

  24. Westerhoff HV, Chen YD: How do enzyme activities control metabolite concentrations?. Eur J Biochem. 1984, 142: 425-430. 10.1111/j.1432-1033.1984.tb08304.x

    Article  CAS  PubMed  Google Scholar 

  25. Westerhoff HV, van Dam K: Thermodynamics and control of biological free-energy transduction. 1987, Amsterdam: Elsevier

    Google Scholar 

  26. Giersch C: Control analysis and metabolic networks. Eur J Biochem. 1988, 174: 509-519. 10.1111/j.1432-1033.1988.tb14128.x

    Article  CAS  PubMed  Google Scholar 

  27. Reder C: Metabolic control theory: a structural approach. J Theor Biol. 1988, 135: 175-201. 10.1016/S0022-5193(88)80073-0

    Article  CAS  PubMed  Google Scholar 

  28. Beard DA, Qian H: Relationship between thermodynamic driving force and one-way fluxes in reversible processes. PLoS One. 2007, 2: e144- 10.1371/journal.pone.0000144

    Article  PubMed Central  PubMed  Google Scholar 

  29. Heuett WJ, Qian H: Grand canonical Markov model: a stochastic theory for open nonequilibrium biochemical networks. J Chem Phys. 2006, 124: 044110- 10.1063/1.2165193

    Article  PubMed  Google Scholar 

  30. Fell DA: Increasing the flux in metabolic pathways: a metabolic control analysis perspective. Biotechnol Bioeng. 1998, 58: 121-124. 10.1002/(SICI)1097-0290(19980420)58:2/3<121::AID-BIT2>3.0.CO;2-N

    Article  CAS  PubMed  Google Scholar 

  31. Derrida B: Velocity and diffusion constant of a periodic one-dimensional hopping model. J Stat Phys. 1983, 31: 433-450. 10.1007/BF01019492.

    Article  Google Scholar 

  32. Fisher ME, Kolomeisky AB: The force exerted by a molecular motor. Proc Natl Acad Sci USA. 1999, 96: 6597-6602. 10.1073/pnas.96.12.6597

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  33. Oster G, Wang HY: Reverse engineering a protein: the mechanochemistry of ATP synthase. Biochim Biophys Acta. 2000, 1458: 482-510. 10.1016/S0005-2728(00)00096-7

    Article  CAS  PubMed  Google Scholar 

  34. Beard DA, Qian H, Bassingthwaighte JB: Stoichiometric foundation of large-scale biochemical system analysis. Modelling in Molecular Biology. Edited by: Ciobanu G, Rozenberg G. 2004, 1-19. [Natural Computing Series]., New York: Springer

    Chapter  Google Scholar 

  35. Bender CM, Orszag SA: Advanced Mathematical Methods for Scientists and Engineers. 1999, New York: Springer-Verlag

    Book  Google Scholar 

Download references

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.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to William J Heuett.

Additional information

Authors' contributions

HQ initiated this work. All authors contributed to and wrote the manuscript together. All authors read and approved the final manuscript.

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Heuett, W.J., Beard, D.A. & Qian, H. Linear analysis near a steady-state of biochemical networks: control analysis, correlation metrics and circuit theory. BMC Syst Biol 2, 44 (2008). https://doi.org/10.1186/1752-0509-2-44

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1752-0509-2-44

Keywords