 Methodology article
 Open Access
 Published:
(Im)Perfect robustness and adaptation of metabolic networks subject to metabolic and geneexpression regulation: marrying control engineering with metabolic control analysis
BMC Systems Biologyvolume 7, Article number: 131 (2013)
Abstract
Background
Metabolic control analysis (MCA) and supply–demand theory have led to appreciable understanding of the systems properties of metabolic networks that are subject exclusively to metabolic regulation. Supply–demand theory has not yet considered geneexpression regulation explicitly whilst a variant of MCA, i.e. Hierarchical Control Analysis (HCA), has done so. Existing analyses based on control engineering approaches have not been very explicit about whether metabolic or geneexpression regulation would be involved, but designed different ways in which regulation could be organized, with the potential of causing adaptation to be perfect.
Results
This study integrates control engineering and classical MCA augmented with supply–demand theory and HCA. Because geneexpression regulation involves time integration, it is identified as a natural instantiation of the ‘integral control’ (or near integral control) known in control engineering. This study then focuses on robustness against and adaptation to perturbations of process activities in the network, which could result from environmental perturbations, mutations or slow noise. It is shown however that this type of ‘integral control’ should rarely be expected to lead to the ‘perfect adaptation’: although the geneexpression regulation increases the robustness of important metabolite concentrations, it rarely makes them infinitely robust. For perfect adaptation to occur, the protein degradation reactions should be zero order in the concentration of the protein, which may be rare biologically for cells growing steadily.
Conclusions
A proposed new framework integrating the methodologies of control engineering and metabolic and hierarchical control analysis, improves the understanding of biological systems that are regulated both metabolically and by gene expression. In particular, the new approach enables one to address the issue whether the intracellular biochemical networks that have been and are being identified by genomics and systems biology, correspond to the ‘perfect’ regulatory structures designed by control engineering visàvis optimal functions such as robustness. To the extent that they are not, the analyses suggest how they may become so and this in turn should facilitate synthetic biology and metabolic engineering.
Background
With the development of quantitative functional genomics approaches, it has become possible to analyse the cellular adaptation of cell physiology to altered environmental conditions experimentally, by monitoring changes in fluxes, metabolites, proteins or mRNAs. Such adaptations tend to occur at multiple regulatory levels if not simultaneously, then subsequently, depending on the time scales of observation[1–3]. In principle, an adaptive change in the rate of an enzyme (or flux) can be mediated by changes in (i) the concentration of metabolites (e.g. substrates, products and effectors) with direct, cooperative and allosteric effects on the activity of the enzyme[4], (ii) changes in the concentration of the enzyme through geneexpression alterations, and (iii) covalent modification via signal transduction. The first is termed metabolic (or enzymatic) regulation. The second is known as geneexpression (mediated) regulation and the third as signaltransduction (mediated) regulation. Because of similar properties, the latter two types of regulation have been considered together under the term ‘hierarchical regulation’[2, 5, 6]. Although in this paper only the former two types of adaptive changes will be discussed explicitly, because of the abovementioned similarities, the third type is addressed implicitly. Until now, significant progress has been made on the modelling of genomescale metabolic networks in microorganisms integrating metabolic and geneexpression regulation[7, 8]. The steadystate properties of a number of representative metabolic regulatory mechanisms, such as endproduct inhibition, have been investigated substantially both in terms of metabolic control analysis (MCA)[9, 10] and by the supply–demand theory championed by Hofmeyr and CornishBowden[11–13]. In order to take geneexpression regulation into account, hierarchical control analysis (HCA)[14, 15] has been developed as an extension to MCA, but it has not yet been linked up with the supply–demand theory. Developing such a link would seem useful as in quantitative experimental studies geneexpression regulation turned out to be as important as metabolic regulation[1, 2, 5, 16].
The adaptive changes of reaction rates through metabolic and genetic regulation are usually due to feedback and/or feedforward mechanisms. In biology, there is a perception that evolutionary optimization has made these mechanisms perfect. If this were so, this would suggest that such mechanisms might be identical to ‘perfect’ regulatory mechanisms designed by control engineering[17]. Indeed, Csete and Doyle[18] have suggested that such a convergent evolution of engineering and biology may have occurred. In particular, they came with an integral control structure containing both an actuator unit (corresponding to an integrator) and a controller/sensor unit. They showed that this regulatory structure would lead to a phenomenon called perfect adaptation and then proposed that such structures should be common to biology. In systems biology contexts, several biochemical processes have been discussed in terms of their control system structures. For example, robust perfect adaptation in bacterial chemotaxis signalling system, in mammalian iron and calcium homeostasis, and in yeast osmoregulation, have been interpreted as integral feedback control systems[19–22], without however proving that they corresponded precisely to the very same regulatory topology or even performance. A recent study identified the three different types of control structures used in control engineering, i.e. proportional, integral, and derivative control, in the regulation of energy metabolism[23]. With the exception of[22], the above work focused only on metabolic regulation, whereas[22] did not compare metabolic regulation with geneexpression regulation. In this study, the integration of metabolic and geneexpression regulation plus the integration between Metabolic Control Analysis and Control Engineering will be investigated.
Control engineering has examined which network structures may make adaptation of a network upon a sustained perturbation of a network component, ‘perfect’. Perfection was defined as the phenomenon that some important system variables (known as ‘controlled variables’) should be completely robust to the perturbations, i.e. with steady states values unaffected by the perturbations. Such perfect robustness can be achieved when a time integrator is applied to any variation of the controlled variable (or system error). This control feature is known as ‘integral control’. Through this time integral, the network would continue to change until the controlled variable is restored completely to its initial value. Because there must be some compensation for the perturbation, a different system variable then has to move away from its initial state. This socalled ‘manipulated variable’ is nonrobust (fragile) to the perturbation, but enables the controlled variable to be robust.
If the control action is proportional to the variation of the perturbed variable itself, or a function thereof that is zero when that variation equals zero, the ultimate deviation of the controlled variable from its value before the perturbation, will be nonzero. This is the socalled ‘proportional control’ of control engineering. Perturbations may also result in a sustained oscillation of the controlled variable and to prevent this from happening, the third type of control focused on by control engineering can be useful, i.e. socalled ‘derivative control’ , which will not be discussed in this paper, but has been exemplified in reference[23].
As mentioned above, the mechanism of integral control is often referred to as ‘perfect adaptation’. Other authors have referred to similar network behaviour that was not based on the same integrative mechanism by the same phrase of perfect adaptation. One such case is that all steps in a metabolic pathway are regulated identically, i.e. their activities being modulated by the same factor. Tyson et al.[24] referred to this as perfect adaptation, but the mechanism hinges on precise regulation of various steps, we would suggest to refer to this as ‘perfect regulation’ since the adaptation part is not crucial. Kacser and Acerenza[25] called this the universal method for metabolic engineering. Fell and Thomas[26] proposed that this may be a common motif in biological regulation and Adamczyk et al.[27] elaborated it into the stealthy engineering principle. This paper will not discuss this perfect regulation mechanism, but focus on the robust perfect adaptation mechanism operating through integral control loops.
In this work, we shall try to bridge two rather unconnected approaches in analysing regulation of network properties. The one is that of control engineering which has devised networks structures that lead to perfect or imperfect adaptation. The other is that of biochemistry with MCA and supply–demand theory, as well as truetolife examples of intracellular biochemical networks involving both metabolism and gene expression. We shall focus on pathways synthesizing precursors for macromolecule synthesis (proteins, nucleic acids) in which that precursor often inhibits an enzyme early in the pathway, both directly and through gene expression. Such endproduct regulatory structures allow for some simplifications[28]. This makes them suitable for illustrating our relatively simple conclusions that are however valid more generally. We shall hypothesize that because a time integration of protein synthesis is involved, geneexpression regulation should be a prime example of integral control, whilst metabolic regulation is our candidate for the role of proportional control. We shall then interpret both these steadystate robustness properties and the control properties in terms of a new hierarchical supply–demand framework.
Methods
Kinetic description and classical control analysis
In this section, we demonstrate that the unique steady state of a metabolic network under regulations can be analysed by both the kineticsbased analysis and by metabolic (or hierarchical) control analysis. A hierarchical supply–demand theory linking hierarchical control analysis with classical supply–demand analysis, is developed for when geneexpression regulation is active. As a result the steady state properties of a metabolic network subject to various regulatory mechanisms can be analysed within a unified theoretical framework.
Basic regulatory architecture and kinetic analysis
The overall regulatory behaviour of a pathway can be decomposed into a number of elementary structures. Feedback inhibition by endproduct has been reported for quite a few metabolic pathways, particularly in anabolism[29]. Another regulatory motif is the feedforward activation of downstream enzymes (see Appendix A and[30]). In addition, in many metabolic pathways one or a few reactions are product insensitive. The reactions catalyzed by hexokinase, phosphofructokinase, and pyruvate kinase in mammalian glycolysis[31] and several steps in the central carbon metabolism in B. subtilis[28], constitute examples. The activities and concentrations of some of these enzymes are regulated by allosteric effectors, covalent modification or transcription. In this study we take a linear pathway with metabolic and geneexpression regulation of the first reaction through the end metabolite as the example of choice[28] (Figure 1). The first reaction (catalyzed by enzyme E_{1}) is assumed to be insensitive to its immediate product. With this example we will be able to illustrate the essence of the principles we are after.
The following differential equations describe this endproduct regulation pathway:
Here, we assume that the concentration of the substrate x_{1} is not influenced by the pathway and that only E_{1} is regulated through geneexpression. x_{ i } is the concentration of i^{th} metabolite, and f_{ i } describes the kinetics of the i^{th} reaction. Parameter p corresponds to other factors that could affect the activity of the first enzyme, e.g. cofactors or external metabolic modulators. Such effects on other enzymes are not addressed here. The gene expression function g is here assumed to depend on the ultimate metabolite only, the latter acting on the synthesis of the first enzyme. More realistic situations involving the dynamics of mRNA will be discussed in later sections. This paper is relevant for geneexpression regulation in general, i.e. includes regulation at the level of transcription, translation and posttranslational modification, but our examples will mostly deal with only one type of these at a time and mainly consider transcription regulation. k_{ED} is the degradation rate constant of the first enzyme. In fast growing organisms and for stable proteins k_{ED} may merely represent the dilution effect due to cell growth and division (i.e. k_{ED} = μ, μ denoting the specific growth rate)[5, 28], but in other cases it will depend on proteolysis, which will be discussed later.
By definition of the steadystate in living cells[32], x_{1}(t) and p(t) are constants at steady state, i.e.${x}_{1}\left(t\right)={\overline{x}}_{1}$ and$p\left(t\right)=\overline{p}$. Similarly, the concentrations of the enzymes are constants and equal to${\overline{E}}_{2},\dots ,{\overline{E}}_{n}$ respectively. Accordingly, if such a constant steadystate regime exists, and it is the unique solution of the following equation:
which only depends on the first and the end enzyme features and is only a function of the end product concentration x_{ n }. Usually, both${f}_{1}\left({\overline{x}}_{1},{x}_{n},\overline{p}\right)$ and g(x_{ n }) are monotonically decreasing functions of x_{ n }, which describe the negative, metabolic and geneexpression regulation, respectively. Likewise, f_{ n }(x_{ n }) is usually a monotonically increasing function of x_{ n }. As shown by Equation (2), the steadystate regimen corresponds to the intersection of the two functions, and is unique due to the monotonic characteristics of f_{1}, f_{ n } and g, as illustrated in Figure 2.
Alternatively, if the i^{th} reaction (i > 1) is product insensitive but the first reaction is not, then f_{1} is redefined and it also becomes a function of x_{2}, while the kinetic function of the i^{th} step only depends on x_{ i }. It can be proven that at steady state, x_{2} then only depends on functions f_{2}, …, f_{ i } and is independent of f_{i + 1}, …, f_{n  1}. This conclusion is both theoretically attractive and practically useful, because the steady state properties of an otherwise complex metabolic pathway may only depend on a limited number of enzyme features[33]. This is a case where the complexity of a pathway is limited; its flux and the concentrations of the upstream metabolites are only controlled by the properties of the upstream enzymes and the corresponding genes.
We will now examine how metabolic control analysis and supply–demand theory deal with these types of metabolic control structures.
Metabolic regulation: MCA and the supply–demand theory
Metabolic Control Analysis (MCA) has mostly dealt with the steady state properties of the metabolic part of schemas such as the ones mentioned in the previous subsection. Modular MCA[34] has divided metabolic networks of this type into modules with relatively autonomous activities, connected through wellidentified metabolites. In the supply–demand theory, Hofmeyr and colleagues[12, 13] have used a similar simplification to demonstrate much of the essence of the regulation of cell function. According to the latter, the metabolic part of the endproduct pathway in Figure 1 can be partitioned into a twostep linear pathway with supply and demand blocks as shown in Figure 3.
The concentration and flux control coefficients of metabolic control analysis measure the steady state change in the concentration of a metabolite X and flux J, respectively, in response to a change in the activity of a process i. That change in activity may be effected by changing a parameter p_{ i } (e.g. concentration of an inhibitor) that is specific for step i:
The v_{ i } is an activity rather than the rate of reaction i at steady state. When the parameter p_{ i } is changed, v_{ i } is not equal to the resultant change of the steadystate rate of reaction i; it is the change in that rate only if all the other variables that affect that rate been kept constants. The concentration control coefficients of the aforementioned supply–demand system can then be defined as${C}_{s}^{{x}_{n}}=\partial ln\left({x}_{n}\right)/\partial ln{v}_{\mathrm{supply}}$ and${C}_{d}^{{x}_{n}}=\partial ln\left({x}_{n}\right)/\partial ln{v}_{\mathrm{demand}}$, and the flux control coefficients,${C}_{s}^{J}$ and${C}_{d}^{J}$ similarly. These control coefficients describe the control exercised by a specific reaction or enzyme on the overall system variables or fluxes, while the ‘local’ regulatory properties of individual enzymes are quantified by the elasticity coefficients, such as
${\u03f5}_{{x}_{n}}^{s}$ and${\u03f5}_{{x}_{n}}^{d}$ measures how the reaction rates of the supply and demand blocks are influenced by the end product concentration x_{ n }. Because the rates are those of blocks rather than of single enzymes, these elasticities are named ‘overall elasticities’[32]. They can be further expressed into more, though not quite elemental elasticity coefficients:
where${\u03f5}_{{x}_{n}}^{{v}_{1}}$ denotes the total elasticity of the first reaction with respect to x_{ n } through the metabolic regulation. The lowercase flux control coefficients c now refer to the local control within the supply module only (if seen as if in isolation with x_{ n } fixed). They can be obtained from summation and connectivity theorems as control coefficients of the module (see Appendix B). When assuming the first reaction is product insensitive, the supply modular elasticity in (5) can be further simplified as
In such a case, the first enzyme has the full control of flux (i.e. ratelimiting step) within the supply module, such that${c}_{1}^{{J}_{1}}=1$ and
Using the summation and connectivity theorems in MCA, all the four concentration and flux control coefficients can be expressed in terms of the elasticity coefficients (see Appendix B). For the concentration control of the supply step and the flux control of the demand step this becomes, using (5) and (7):
These equations show that if the first reaction has complete flux control within the supply module, the control of steady state endproduct concentration and flux are only functionally depend on two elasticity coefficients, i.e. the ones that correspond to the first and last reactions. When the feedback is very strong, i.e.$\left{\u03f5}_{{x}_{n}}^{{v}_{1}}\right\gg \left{\u03f5}_{{x}_{n}}^{{v}_{n}}\right$, the control of the demand flux${C}_{d}^{J}$ is close to 1. The MCA and supply–demand simplification thereof discussed in this section partially explain the statesteady properties of the aforementioned full regulatory system, but only for the case of metabolic regulation. When including the geneexpression regulation, a more complete interpretation can be achieved by using the hierarchical control analysis and a new ‘hierarchical supplydemand’ theory as investigated in the next subsection.
Geneexpression and metabolic regulation: hierarchical supply–demand theory
Westerhoff and coworkers have developed hierarchical control analysis (HCA), an extension of MCA that can take geneexpression regulation and signal transduction into account[14, 26]. We shall here implement this by extending the meaning of the elasticity coefficients in the previous subsection to include regulation through gene expressions.
When considering the metabolic part of the network alone, i.e. if gene expression were always the same, the control on the concentration of intermediate X in a supply–demand system follows (8). If the roles of the synthesis and degradation of metabolic enzymes are considered explicitly, as illustrated in Figure 4, HCA has to be introduced, and the corresponding hierarchical control coefficient becomes:
Capital H is here used for the hierarchical control coefficients as defined in Table 1.$*{\u03f5}_{X}^{s}$ is an “overall” elasticity coefficient, including a classical ‘direct elasticity’ only related with metabolic responses (i.e.${\u03f5}_{X}^{s}$ similar to the${\u03f5}_{{x}_{n}}^{s}$ defined in the MCA in (8)) and an ‘indirect elasticity’ due to geneexpression regulation:
The lower case c is used for ‘metabolic control coefficients’ , i.e. control coefficients that only take the local network (metabolic, or gene expression but not their combination) into account.${\u03f5}_{{E}_{s}}^{s}$ is often equal to 1, i.e. when the rate of the reaction in isolation is proportional to the concentration of the enzyme catalyzing it.
Using metabolic control analysis for the gene expression part of the network, the control coefficient of the protein synthesis reaction with respect to the concentration of the protein synthesized is:
Combining the above expressions, one can express the hierarchical coefficient quantifying the control exerted by the supply enzyme on the concentration of the metabolic intermediate X in terms of all the elasticity coefficients in the network:
The terms in parentheses are usually positive. The equation shows that the control by supply (i) decreases with the absolute magnitudes of the elasticities with respect to X of the supply, of the demand, and of the protein synthesis, but (ii) increases for increasing elasticities of the protein synthesis and degradation reactions with respect to the concentration of the enzyme. The equation also shows that for finite nonzero magnitudes of the elasticities, the hierarchical control coefficients for control by supply may be decreased by elasticities in the geneexpression network, but is usually not brought down all the way to zero. The same applies to the control by demand, which is equal to minus the control by supply.
Now let us recall the endproduct module with both metabolic and geneexpression regulatory feedbacks of Figure 1. As the end product (x_{ n }) regulates the first reaction through both metabolic regulation and geneexpression regulation, the corresponding ‘overall’ elasticity (see Table 1) of the first reaction can be expressed as:
where${\u03f5}_{{x}_{n}}^{{v}_{1}}$ denotes the elasticity through metabolic regulation,${c}_{{v}_{\mathrm{Trans}}}^{{E}_{1}}$ denotes the control of gene expression (i.e. transcription and translation) on the concentration of the first enzyme[15, 35]. By replacing the${\u03f5}_{{x}_{n}}^{{v}_{1}}$ in (7) with the more complete expression (13), i.e. using${\u03f5}_{{x}_{n}}^{s}=*{\u03f5}_{xn}^{v1}$, the hierarchical control coefficient${H}_{s}^{{x}_{n}}$ can be expressed into elasticity coefficients, in a similar manner as (12). Since both metabolic and geneexpression regulation constitute negative feedbacks,${\u03f5}_{{x}_{n}}^{s}$ is negative and becomes more negative with increasing x_{ n } due to increasing product inhibition.${\u03f5}_{{x}_{n}}^{d}$ is positive and decreases asymptotically to zero with increasing x_{ n }. Therefore, the relationship between the reaction rate and end product concentration can be described in terms of the elasticities as depicted in Figure 5. Figure 5 explains the unique steadystate results obtained from classical steady state analysis (see Figure 2).
An illustration of the supply–demand relationship similar to Figure 5 has been presented in[13]. However, here we extend the interpretation of the system and corresponding elasticity coefficients to the more general case that includes geneexpression regulation. As an extension to the classical supply–demand theory, the analysis given in this section can be named the ‘hierarchical supplydemand’ theory.
Control engineering
The discipline of Control Engineering first identifies a socalled controlled variable, which it sees as the output of the system. In metabolic biochemistry, output often relates to a flux, but can also be the concentration of an important metabolite in the pathway. Control engineering next examines the various categories of mechanism that may contribute to the capability of the network of maintaining the controlled variable close to its original steady state value when the system is subject to a sustained perturbation. The ‘error (function)’ is the deviation (δX) of the value of the controlled variable (X) from its value before the perturbation, or the difference to a reference signal r. The network ‘adapts’ to the perturbation of the controlled variable, i.e. to the error function, in a socalled ‘control action’. RNA polymerase plus the ribosomes that together translate changes in the concentration of metabolites to changes in gene expression, or direct metabolic regulation of the activity of an enzyme correspond to such control actions. The output of the control action (or of the ‘controller’) is often named the manipulated variable. In metabolism, reaction rates v(X, E) are variables manipulated either by the concentration of metabolites (f(X)) or by the concentration of the enzyme that catalyses the reaction concentration (E). A mechanical control system often includes an actuator that converts the control signal into some kind of mechanical motion. For a biochemical system discussed in this study, this may correspond to the enzyme catalysing the reaction synthesizing or degrading X. Usually there exists a sensor measures the controlled variable and translates its error function into the input signal of the controller. In a geneexpression regulation, the transcription factor can be regarded as the sensor.
The three most widely used categories of control are the proportional, integral, and derivative (PID) control mechanisms[17]. They differ depending on whether the control system’s response is a function of the ‘error function’ itself, the time integral thereof or the time derivative thereof, respectively. In systems biology literature, the proportional control mechanism has already been referred to in terms of metabolic regulation[19, 20, 23] (e.g. feedback inhibition). However, when Control Engineering discusses proportional control mechanisms, response is proportional to the error function. In actual biochemistry, enzyme activity is rarely a linear function of the concentrations of metabolites X, which includes the enzyme’s substrate, its product and allosteric modifiers. MCA accommodates this nonlinearity by allowing the elasticity coefficient to differ from 1. Metabolic regulation by the ‘error function’, is part of the nonlinear dynamics of the process or system, i.e. f(X), both conceptually and in the mathematical modelling. It would seem therefore that the proportional control of Control Engineering can be nonlinear in biochemical networks.
Integral control action through the accumulation of molecules in the metabolic process has also been reported[19, 20, 23]. Here the systems response should be a function not of the error function itself but of the integral of that error function. In the present study we examine geneexpression regulation from this point of view, since protein synthesis requires time integration and depend on the error function, and because changes in protein concentration directly affect the rate of the reaction the protein may catalyzes. We may expect that this integral control somehow corresponds to the ‘indirect elasticities’ of HCA. Whether indeed geneexpression regulation corresponds to an exact (or ideal) integral control mechanism will be further discussed in the Results and Discussion section. Whether there exists a derivative control action and whether it relates to a specific type of regulation in a metabolic pathway will not be investigated in this paper. Reference[23] already identified an example.
Considering the dynamics of a metabolic system, we can write the time dependence of the concentrations of its metabolites as:
N is the stoichiometry matrix. E is a diagonal matrix with the concentrations of the enzymes that catalyse the various reactions along its diagonal. f(X) is a vector function of the concentrations of the metabolites and kinetic parameter values. The regulated metabolic pathway can be described in terms of a closedloop feedback control system, as indicated in Figure 6, in which k_{ p }, k_{ i } and k_{ d } are the PID control parameters. We here note that Figure 6 and subsequent figures refrain from biochemical detail. This is because the analysis in the present paper aims at obtaining a set of conclusions with general significance. Being specific in the schemes we use as illustrations would detract from this aim.
When considering a sustained perturbation γ∙δp (e.g. change in a parameter p) and denoting by δ the (small) deviation from the steady state prior to this perturbation, the time dependent variation in the metabolite concentrations may be observed:
The subscript ss refers to the steady state values. By substituting the time integration of gene expression, i.e.${k}_{i}\cdot {\displaystyle {\int}_{0}^{t}q\left(\tau \right)}\cdot \mathit{d\tau}$, for δE, and assuming the proportional and derivative actions a part of the metabolic process (14),
This describes the overall dynamics of a closedloop metabolic system under perturbation as a sum of three terms. The first term of these is a nonlinear function of (or in first order proportional to) the perturbation of the controlled variable (i.e. the error function) δX. This term describes all the direct elasticities, including nonregulatory system kinetics such as substrate and product effects, and (other) metabolic regulation such as allosteric activation. The second term corresponds to a time integral of a function of the perturbation of the controlled variable δX, and can also depend on other system variables as discussed in the Results and Discussion section.
We shall now examine whether these two terms in the equation (16), correspond to the proportional and integral control loops of Control Engineering.
Results and discussion
A simple example of combined metabolic and geneexpression regulation of an important intracellular process: ATP (energy) metabolism
Let us consider the simple example given in Figure 7, i.e. a twostep pathway with ATP and ADP as combined intermediate and with the expression of the gene encoding the first enzyme E increasing in proportion to the concentration of ADP. The ‘moiety conservation sum’ C is the sum of the concentrations of ATP and ADP and a constant here (i.e. C = [ATP] + [ADP]) because the reactions only convert the one into the other. The metabolic regulation addresses the interplay between the supply and demand processes (s and d).
The dynamics of ADP and enzyme E are here assumed to follow mass action kinetics:
The degradation of the enzyme is here written as the sum of two terms, which will serve to emphasize the implication of this degradation to be independent (for k_{ b } = 0) or dependent (for k_{0} = 0; see below) on the enzyme concentration. The first order degradation rate reflects the assumption that there is a rather unspecific protease activity for which the particular enzyme E we are considering here is a minority substrate. The zero order degradation would reflect a case where there is a specific protease system for enzyme E (e.g. an ubiquitination followed by a generic protease) that is saturated by the already high concentration of the enzyme relative to the K_{ M } of the ubiquitin transferase. In (17) all rate constants are considered nonnegative and k_{0} = 0 whenever E is nonpositive. The closedloop control system structure of the pathway can be represented as in Figure 8.
The control system diagram suggests that the ADP concentration is the controlled variable and the enzyme concentration E a manipulated variable in the geneexpression control loop. The zero order degradation rate k_{0} can be treated as a reference signal to the system. The metabolic regulation is included as a part of the ADP kinetic process. By considering a perturbation of k_{ d } from its steady state value (i.e. δk_{ d }), and reformulating the kinetics of ADP and E (see Appendix C), we have
Comparing this to a general closedloop control system (16), with ADP for X, we recognize on the righthand side first a proportional response term, then an integral response term, and then the perturbation term. The proportional response corresponds to the direct ‘elasticity’ of the supply and demand reactions with respect to the error function δ[ADP], which is a metabolic and instantaneous regulation. The integral response is related to the protein synthesis and degradation and thus to geneexpression regulation.
When the degradation of enzyme is zero order in terms of E (i.e. when k_{ b } = 0), the geneexpression regulation becomes an ideal integral control loop, and the metabolic network can exhibit robust perfect adaptation to the external or parametric perturbations. This can be understood by requiring (18) to be valid in a steady state, i.e. with time independent values for [ADP] and E. Because the time integral reaches to infinity this requires that the argument of the integral, i.e. k_{ a } ⋅ δ[ADP]  k_{ b } ⋅ δE, must equal zero. The mechanism for this perfect adaptation is that after the perturbation the concentration of the enzyme will vary until it makes the time dependences equal zero by itself, forgoing the more usual process in metabolic regulation that changes in the controlled variable arrange for the steady state to be reattained in the presence of the sustained perturbation.
For the more general case where the enzyme degradation may depend on the enzyme’s concentration, the (hierarchical) control of the enzyme level by the demand reaction can be expressed in terms of the kinetic parameters and the steadystate ADP concentration (see Appendix C):
The flux control exercised by the demand reaction is quantified by:
Both the control of enzyme level and the control of demand flux by the perturbation equal 1 minus a hyperbolic function of k_{ b }. For the ideal integral control scenario of k_{ b } = 0, the enzyme concentration E tracks the activity of the pathway degrading ATP perfectly, i.e.${H}_{{k}_{d}}^{E}=1$. More importantly, the pathway flux perfectly tracks the perturbation in the demand flux and${H}_{{k}_{d}}^{J}=1$. This is the case of robust perfect adaptation. For other cases when k_{ b } ≠ 0, the adaptation of the pathway to the perturbation will not be perfect and both control coefficients are smaller than 1.
Also the robustness coefficient[36] of the ADP concentrations visàvis perturbations in the demand reaction can be expressed in terms of kinetic constants and the concentration of ADP (see Table 1 for definition):
Only if k_{ b } = 0, k_{ d } = 0, or [ATP]_{ ss } =0, the ADP and ATP are perfectly robust (${\Re}_{{k}_{d}}^{\left[\mathit{ADP}\right]}=\infty $) versus perturbations. The fragility[36] of the ADP concentration visàvis perturbation in the demand reaction, which is the inverse of the robustness, can be quantified by the concentration control coefficient for the concentration of ADP with respect to the degradation process. It reads as:
This fragility is a hyperbolic function of the first order degradation rate constant of the enzyme and hence zero when that degradation is zeroorder (see Figure 9 for illustration). The fragility has the ATP/(ADP + ATP) ratio as its maximum value (hence the minimum robustness equals the (ADP + ATP)/ATP ratio). Half maximum fragility is attained for:
This means that the fragility may be low for a substantial magnitude of the first order rate constant of protein degradation if the rate constant for protein synthesis is also high. The control coefficients for the enzyme level${H}_{{k}_{d}}^{E}$ and the demand flux${H}_{{k}_{d}}^{J}$ attain a maximum of 1 and a value of ½ when the fragility of ADP is half maximal.
These conclusions can be generalized somewhat by directly implementing HCA and hierarchical supply–demand result given in (12). For the above model the elasticity coefficients assume the following magnitudes:
For the robustness of the ADP concentration visàvis increased demand one then finds:
This shows that robustness is infinite if the enzyme degradation reaction is zero order (i.e.${\u03f5}_{E}^{b}=0$), and that the robustness becomes smaller with increasing order (elasticity) of this reaction.
The most important conclusion here is that there is no discontinuity in the ability of integral control loops to lead to good adaptation. The closer the degradation of the enzyme that enables the adaptation is to a zeroorder reaction, the stronger its tracking of the perturbation in the demand flux, and the higher the robustness of the variable that is to be kept homeostatic. Important perhaps is the phenomenon that the robustness is not determined by the magnitude of the degradation reaction but by its kinetic order (i.e. elasticity of effective Hill coefficient). The corresponding conclusions pertain to the tracking of the demand by the enzyme level E and the control of the pathway flux by the demand reaction.
A further issue in control engineering is the robustness of systems versus perturbations at various frequencies. In engineering, an airplane wing has to be robust to variations of air pressures at high frequencies, as well at low frequencies. In order to achieve this combined robustness, different control loops may have to be put in place simultaneously, although a tradeoff limits what one can do[18]. In systems biology, this can be illustrated for the endproduct feedback regulation in Figure 1. If the flux demand of the pathway increases rapidly, the concentration of the end product decreases rapidly and as a result of the direct allosteric product inhibition effect, the activity of the first enzyme will increase quickly too. This metabolic control of enzyme activity is a fast actuator of the system. However, if there is a further increase in the flux demand, the first enzyme may ‘lose’ its control capacity since its activity may be approaching its maximum capacity (k_{cat}). At this stage, the system may then undergo a second ‘adaptation’ through gene expression which should be expected to be slower because the cell has to produce enzyme, but still ultimately lead to an increase in the concentration of the first enzyme. This increase should then decrease the direct metabolic stimulation of the catalytic activity of the enzyme discussed above. In this sense, the regulation of the first enzyme of the pathway is bifunctional in dynamic terms[18]: The metabolic regulation rapidly rejects high frequency perturbations but possibly with small amplitude or capability, while the geneexpression regulation is slow to adapt, but may be able to reject very large constant perturbations.
Conditions for integral and pseudointegral control: endproduct pathway
In this section, we recall the endproduct module example (of which the simple ATP metabolism example is a special case) to analyse geneexpression regulation in a more general metabolic pathway. In particular, the conditions will be identified for which the geneexpression regulation constitutes ideal integral control or pseudointegral control will be discussed. A simple but representative endproduct example is given in Figure 10, i.e. a linear pathway with three metabolites, where for all the three enzymes the gene expression is regulated by the last metabolite. The kinetic model of this example with all the parameters is provided in[5].
Gene transcription and translation are modelled here explicitly by further including the dynamic function of mRNA as
Here g_{Trnl}(R) = k_{Trnl} ⋅ R is a function of mRNA concentration R. k_{ED} essentially consists of three parts. One term is due to dilution, which is proportional to the specific cell growth rate μ. The other terms correspond to proteolysis, as below:
The final term denotes the zero order proteolysis as would be caused by proteases that are saturated with the protein of interest. We here assume that the specific growth rate μ is independent of the activity of the pathway under study, an assumption that is sometimes but not always realistic. Since after multiplying with E the last term is independent of the protein concentration, it is convenient to move this term into the protein synthesis function g_{Trnl}(R). Hence, the new protein degradation rate can be defined as,
and the protein dynamics can be rewritten as,
In the exponential growth phase or if proteolysis is first order (i.e. k_{ED} ≠ 0), the above pathway example corresponds to a pseudo or nonintegral control scenario. The control structure of the regulatory system is then given by Figure 11. The dynamics of the ‘sensor’ is decomposed here by addressing both transcription and the translation through mRNA. At steady state, often g_{Trnl}(R)_{ ss } = k_{ED} ⋅ E_{ ss } ≠ 0. Therefore, after perturbation of a system parameter (i.e. kinetic constants), the new steady state values of x_{3}, R, and g_{Trnl}(R) will no longer be the same as the old steady state values, which indicates that then the regulatory system does not achieve perfect adaptation. However, when k_{TnD} is very small, nearperfect adaptation behaviour should be observed.
An ideal integral control scenario will happen only when cell enters a stationary phase and there only exists zero order proteolysis, i.e.${k}_{\mathrm{ED}}=\mu +{k}_{\text{proteolysis}}^{1}=0$. In such case, at steady state, g_{Trnl}(R)_{ ss } ≡ 0 due to integral control, and R and x_{3} will also always keep the same constant values at steady state, no matter how the system parameters are perturbed. The adaptation of the system is then perfect in the sense of making the system properties R and x_{3} robust against the perturbations. The functions g_{Trnl}(·) and g_{Trsc}(·) can also be multivariate, i.e. contain other system variables. Only if those variables do not depend functionally on the protein concentration (E), the adaptation visàvis perturbations of the system kinetic parameters will remain perfect.
The ideal integral control scenario that we considered above led to the perfect robustness of certain system variables with respect to certain (external or parametric) perturbations. A second aspect of the perfect adaptation scenario is the perfect tracking by a second system variable of perturbations. Perfect tracking means that the relative change in the variable is identical to the relative change in the perturbing parameter. If the parameter is the activity of a process, then this means that the corresponding control coefficient is equal to 1. The perfect tracking of references and the perfect robustness of controlled variables to perturbations are two aspects of integral control systems, i.e. the two features can be observed simultaneously for the same system. A specific pathway then shows perfect robustness of a system variable visàvis multiple perturbations (or parameters) but perfect tracking with respect to only a limited set of parameters (e.g. a reference signal r).
For the aforementioned example, the zero order proteolysis can be regarded as a ‘reference’ signal$r={k}_{\text{proteolysis}}^{0}$. Effectively perturbation of this rate constant corresponds to an external perturbation of protein synthesis. The effect is that even at an ideal integral control case, i.e. when k_{ED} = 0, the steady state concentration of mRNA will not become 0. Rather, it will always ‘track’ the reference r (i.e. R = r/k_{Trnl} ≠ 0) whenever${\dot{E}}_{\mathit{ss}}=0$. The reference tracking control system structure is included in Figure 11 with reference r referred to by a dashed arrow, and with g_{Trnl}(R) then representing k_{Trnl} ⋅ R.
Practical concerns and assumptions on degradation rate constant k_{ED}
In general,${k}_{\mathrm{ED}}=\mu +{k}_{\text{proteolysis}}^{1}>0$. It may seem that the condition of integral control can be approached, when either i) k_{ED} < < α with α a very small positive value, or ii) k_{ED} < < k_{Trnl} (or k_{ED} · E < < g_{Trnl}(R)). In the latter case, for the level of protein to remain bounded, there should be a background degradation rate of the protein independent of the concentration of that protein. This would be so if:
During the exponential growth phase of bacteria such as E. coli less than 1% of a protein may be degraded during a cell division cycle[37]. Consequently the major term in the protein degradation is the dilution term μ and this term is generally very small. Indeed, by definition, e^{μ ⋅ T} = 2, where T is the doubling time of the bacteria and then μ = ln2/T. Practically, the smallest value for T is close to 20 minutes[38], which corresponding to the fastest growth rate of E. coli:
Such a growth rate in microbes such as E. coli and yeast, which is fast for organisms but slow at the time scale of RNA and protein synthesis, produces a small effective degradation rate constant for the proteins. Below we shall see whether such a small degradation rate constant suffices to produce nearperfect adaptation behaviour in practice, which would be interpreted as a quasiperfect integral control system.
Simulation study
In this section, both the integral and the nonintegral control scenarios are simulated based on the example given in Figure 10 with all the kinetic parameters given in[5]. Three different systems are considered, with k_{ED} = 0, 0.2, and 0.4, respectively. All three systems are simulated from the same initial condition. The concentration changes of mRNA, E and x_{3} are shown in Figure 12. After a period of time (e.g. 30 seconds) the three systems reach different steady states. Then at 50 seconds we perturb one system parameter (i.e. the k_{cat} of the third reaction) by 20% for all three cases. After a while the three systems reach new steady states. After perturbation only the system with zero order protein degradation (i.e. k_{ED} = 0), returns to the same steady state values for mRNA and x_{3} as those before perturbation, which indicates perfect adaptation: this is an ideal integral control system. In this case the manipulated variable, the enzyme level E, varies strongly. Its adaptation enables x_{3} and mRNA to be completely robust visàvis the perturbations of the k_{cat} of the third reaction. In the other two cases (i.e. k_{ED} = 0.2 or 0.4), the change in enzyme level is smaller, but the new steady state values of mRNA and x_{3} deviate from their previous steadystate values and adaptation is not perfect. However, the deviations are not large; only 3% for x_{3} and 4% for mRNA when k_{ED} = 0.2, and 4% for x_{3} and 8% for mRNA when k_{ED} = 0.4.
Now let us consider the reference tracking scenario. The responses of enzyme and mRNA to a 20% perturbation in the protein stability (r) at t = 50 seconds, are given in Figure 13 for three different rate constants of enzyme degradation k_{ED}. Only when k_{ED} = 0 the mRNA concentration tracked the reference value with zero steady state ‘deviation’. This indicates the existence of a perfect integral action of the feedback regulatory system. When k_{ED} was not equal to zero (i.e. k_{ED} = 0.2 or k_{ED} = 0.4), the mRNA response did not track the reference signal, indicating that in these two cases the controller of the system was not an ideal integral controller, although it changed less than did the reference signal.
HCA and an hierarchical supply–demand interpretation
This simulation example can be represented by a hierarchical supply–demand structure such as in Figure 4. To obtain ideal integral control, both protein synthesis and protein degradation should be independent of the concentration of the protein that is being degraded (i.e.$\dot{E}={g}_{\mathrm{Trnl}}\left(R\right)$). Since this implies that protein degradation is zero order in protein concentration:
So that the hierarchical control coefficient of the metabolite concentration becomes
Here${\u03f5}_{{E}_{s}}^{a}={\u03f5}_{E}^{{v}_{\mathrm{Trnl}}}$,${\u03f5}_{{E}_{s}}^{b}={\u03f5}_{E}^{{v}_{\mathrm{ED}}}$,${H}_{s}^{X}={H}_{3}^{{x}_{3}}$ for the biosynthetic pathway example. Because the hierarchical control by supply and demand must add to zero (due to the concentration control summation law[15]), also the control by demand on the metabolite concentration becomes precisely equal to zero in the zero order protein degradation case.
Feedforward activation: a case study of a leucine biosynthetic pathway
In previous sections, either a metabolic intermediate (ADP) or a penultimate product (x_{ n }) inhibited or repressed upstream enzymes. In this section a different regulatory structure is investigated, one in which a metabolite activates downstream enzymes through geneexpression. A simplified mathematical model describing the leucine biosynthetic pathway in Saccharomyces cerevisiae[39] is used to demonstrate that the analysis integrating hierarchical supply–demand theory and control engineering continues to apply. The pathway converts pyruvate to leucine by the sequential reactions described in Figure 14. There are two major regulatory mechanisms in the pathway. One is a metabolic feedback inhibition of Leu4 and Leu9 by leucine, which is an endproduct module similar to the ones discussed above. The other is the transcriptional (geneexpression) activation of downstream enzymes Leu1 (E_{1}) and Leu2 (E_{2}) by αIPM (I_{1}) (through transcription factor Leu3), which we shall call initialproduct modules (see Appendix A and Figure 17). The model predictions fit the experimental data and all the parameter values have been estimated and provided in[39].
The hierarchical control coefficient quantifying the control of supply enzymes E_{s} with respect to the concentration of I_{1} (αIPM) can be expressed into the various elasticity coefficients (see also Figure 15),
and its value can be computed from the simulation of the pathway model by perturbing E_{ u } and recording the changes in the steady state value of I_{1}. This led to${H}_{s}^{{I}_{1}}=0.12$. The hierarchical control coefficient is small but not zero, which indicates that the regulatory network can react and adapt to this external perturbation although the adaptation is not ‘perfect’.
The closedloop control system structure of the downstream geneexpression activation is the same as that of the upstream inhibition case described in Figure 11 (with endproduct metabolic regulation dominating process dynamics). It is still a negative feedback control rather than a feedforward control in the control context: For a negative feedback control system, a change (increase/decrease) in some controlled variable will result an opposite change (decrease/increase) in the operation of the process itself in such a way as to reduce changes. For the endproduct module given in (1), when E_{1} increases the concentration of the end product x_{ n } also increases, whilst as x_{n} increases E_{1} would decreases as a result of negative feedback inhibition (through metabolic or geneexpression regulation). Hence this is a negative feedback control system, since when E_{1} increases the system attempts to reduce such an increase. For case of the positive activation of downstream enzyme by a metabolite upstream of that enzyme (i.e. initialproduct module given in (34)), when E_{1} increases the concentration of x_{1} decreases since it is a downstream enzyme, while E_{1} would also decreases as x_{1} decreasing because of the positive feedforward activation. This gives the same result as the endproduct case. Herewith, both two systems are negative feedback control system in the control context and both can be represented by the same feedback control structure as given in Figure 11.
In this model, the protein degradation rates purely depend on the dilution effects, and the estimated growth rate of μ = 0.0058 min^{1} = 9.7 × 10^{5} s^{1} is very small as compared to metabolic turnover times and may imply the existence of a quasiintegral control scenario. For testing, a simulation study is provided where a 50% increase on the k_{cat} of the first reaction is applied at 800 minutes as an environmental perturbation, the calculated concentration changes of the intermediate metabolites are shown in Figure 16. This perturbation could be due to an allosteric activation by a substance added to the system. The concentrations of both αIPM and βIPM reach new steady states. The differences between the new and the old steadystate values are around 7% and 20% respectively. Again, the adaptation is imperfect; the two concentrations are robust but not infinitely so.
Conclusions
In this paper, two existing approaches to the analysis of the robustness and adaptation of networks have been integrated: control engineering and MCA. The former designs control structures for networks that lead to optimal behaviour, for instance in terms of ‘perfect adaptation’ leading to infinite robustness. The latter quantifies the extents to which processes in a network determine fluxes and concentrations and identifies the molecular interactions that determine the corresponding distributions of control. Two extensions of MCA have also been integrated into our analysis, i.e. HCA adding the possibility to analyse regulation through gene expression, and supply–demand theory, greatly simplifying the analysis towards understanding of the essence. We also integrated the two latter approaches into a novel hierarchical supply–demand theory. The steady state properties of exemplary metabolic pathways served as test cases; they were analysed in terms of the robustness of their steady state properties. They included a pathway of free energy transduction, a pathway with feedback inhibition and repression, and a pathway with feedforward regulation. Most substrates for the synthesis of macromolecules such as proteins and nucleic acids are the end product of such pathways, making this analysis important for the understanding for the control of cell growth.
We then used the resulting framework to address the question whether metabolic pathways regulated by both metabolic interactions and through gene expression, come close to the ideal control structures designed in control engineering. We focused on the control structures leading to socalled ‘perfect adaptation’ , as defined by (i) complete robustness of the concentration of the pathway’s end product towards perturbations in supply and demand, (ii) perfect tracking of perturbations by variables involved in the adaptation. Control engineering distinguishes between proportional, derivative and integral control and showed that of these only integral control loops produce perfect adaptation. In a hierarchical control analysis, in a hierarchical supply–demand analysis, as well as in a number of computer simulations, we showed that such perfect adaptation (and perfect tracking) should not be expected for the usual geneexpression regulation pathways, even though they seem to engage in integral control. Although that integral control increased the robustness of the concentration of the pathway product visàvis perturbation in the activities of metabolic enzymes, that robustness was not perfect, except for the singular case where the metabolic enzymes would be infinitely stable. For complete steady states, such infinite stability should not be considered realistic.
We expected that for cases where proteins are highly, adaptation of the anabolic networks studied should be high even if not quite perfect. This was however not much observed. HCA can show that this is not to be expected either. The perfectness of the adaptation of the networks in complete steady state should not be a function of the magnitude of any firstorder protein degradation rate constant, but rather a function of the order of the protein degradation reaction, or to be precise, of its elasticity coefficient: the protein degradation should be zero order in protein concentration for the adaptation to become perfect.
Because gene expression regulation involves a time integral that is a function of changes at the level of the metabolic pathway, we suggest that at least when applied to biological systems, control engineering is extended so as to explicitly include such ‘integral control’ loops even if they do not lead to perfect adaptation. By adding the subtlety of HCA one can then analyse how perfect the adaptation furnished by the integral control loop actually is. We propose a simplification of the controlengineering concept of ‘proportional control’ to the more general case where the adaptation is any function of the displacement of the controlled variable but neither a function of the time integral nor the time derivative of that displacement. That function could be more complex than the proportionality used in control engineering and thereby be realistic for biochemical networks, thereby making control engineering much more useful for the life sciences. Perhaps the terms ‘proportional control’ should then be replaced by ‘direct or metabolic control’ , referring to the direct nature of the interactions.
Metabolic networks typically have more than one metabolic intermediate and when the network is perturbed by affecting a process activity, many metabolite concentrations tend to change. If it is of particular interest to maintain one of these as constant as possible whereas variations in other metabolite concentrations are less detrimental to biological function, then the former may be designated as the ‘controlled variable’ and the latter as ‘manipulated variables’ in the control engineering analysis. Because many enzymes do not serve a function other than through the reaction they catalyse, they may be the more obvious ‘manipulated variables’. On the other hand it may well be that some metabolite concentrations serve as ‘manipulated variables’ with the sole function of providing near perfect control loops. cAMP might be an example.
Integral control loops do not require geneexpression regulation. Also in exclusively metabolic networks, integral control may arise. An example would be a linear pathway with the penultimate metabolite affecting the first reaction of the pathway, whilst its degradation rate is independent of its concentration. If the first step of the pathway is then perturbed, the concentration of the penultimate metabolite will be a function of the time integral of the perturbation of the first metabolite concentration in the pathway, and it may effect perfect adaptation. On the other hand, if the degradation rate of that penultimate metabolite were first order, then that concentration would be a direct function of the concentration of the first metabolite of the pathway, and the regulation would turn into proportional rather than integral control. In hindsight, reference[33] was an early example of this.
As shown at length in the present paper, geneexpression regulation does not always lead to perfect adaptation. Indeed, if protein degradation is first order, then the deviation in the concentration of the enzyme may well be proportional to the perturbation in the controlled variable and one effectively obtains ‘proportional control’ through geneexpression regulation. We conclude that whether control loops correspond to the integral or proportional category of control engineering depends on the elasticity coefficients (orders) of the reactions involved with respect to the controlled variable, rather than on time integration being involved. The extension of control engineering with hierarchical control analysis that was initiated here, may well provide the subtlety that helps analyse the complex networks that mankind is confronted with today, both in the life sciences and in economics and environmental sciences. It may also help design new and better networks, if only for synthetic biology and biotechnology.
We intend the present work to serve as a beginning of a development where multiple principles of control engineering may be compared with achievements by biological evolution. One example where the present work may be extended to is the case where rather than that one piece of DNA encodes three metabolic enzymes (as in Figure 8), one metabolic enzyme catalyses reactions in three metabolic modules in the metabolic network: a multifunctional enzyme. The gene expression of that enzyme may be regulated by metabolites in all three modules. This is a type of network structure that control engineering may come up with as serving a function of coordination. And it will be of interest whether this type of network may serve a similar function in Biology.
Appendices

A.
The initialproduct module
The initialproduct feedforward regulation module given in Figure 17 can be mathematically described by the following differential equations:
The first reaction here is still assumed product insensitive and other factors p can act on the enzyme. Different from the endproduct module, the function g is assumed to be an increasing function of its argument. It has been demonstrated in[28] that if a constant steadystate regimen exists, the following simple relationship should be satisfied, with${\overline{E}}_{1}=g\left({\overline{x}}_{1}\right)/{k}_{\mathrm{ED}}$ and$\left(g\left({\overline{x}}_{1}\right)/{k}_{\mathrm{ED}}\right)\cdot {f}_{1}\left({\overline{x}}_{1}\right)={v}_{1}$. Here, it is assumed the intermediate reactions of the metabolic pathways do not ‘saturate’.
B. Calculation of global and local control coefficients in a supply–demand system
According to the summation and connectivity laws, for a supply–demand system, we have
By solving the above four equations, the ‘global’ concentration and flux control coefficients with respect to the supply and demand steps can be derived as
The expressions of the ‘local’ flux control coefficients given in (5) can be obtained by solving the following summation and connectivity laws with respect to the local linear pathway within the supply module as given in Figure 3.
C. Steady state analysis of the ATP metabolism example
The steady state values of ADP and E before the perturbation are determined when d[ADP]/dt = 0 and dE/dt = 0:
and
It is assumed that cell function requires ADP concentration to be at a certain level (i.e. [ADP]_{ ss }) and that in the absence of the perturbation; the cell has adjusted E_{ ss } and the rate constant k_{0} (and perhaps k_{ d } and k_{ s }) to meet this requirement. We assume that the cell will do this also at different values of k_{ b }. The implication is that if different values of k_{ b } are considered, k_{ a } is also different as defined by (36).
Considering a perturbation of k_{ d } from its steady state value, one finds for the time dependence of the variation in the enzyme level:
and for the time dependence of the variation in the level of ADP:
By integrating the time dependence of the change in enzyme level δE in (37) into (38) one finds (18). By substituting the steady state condition of E_{ ss } in (35),
For that change in ADP level to be time independent, the integrand in the integral control term should equal zero around steady state, so that:
By using this expression to eliminate the change in enzyme level from (39) for the time dependence of the change in ADP level and set the latter to zero, the robustness coefficient in (21) can be obtained. Similarly, by using (40) to eliminate the change in ADP level from (39), the control of enzyme level (19) can be obtained.
References
 1.
Rossell S, Van Der Weijden CC, Lindenbergh A, van Tuijl A, Francke C, Bakker BM, Westerhoff HV: Unraveling the complexity of flux regulation: a new method demonstrated for nutrient starvation in Saccharomyces cerevisiae. Proc Natl Acad Sci. 2006, 103: 21662171. 10.1073/pnas.0509831103.
 2.
Ter Kuile BH, Westerhoff HV: Transcriptome meets metabolome: hierarchical and metabolic regulation of the glycolytic pathway. FEBS Letters. 2001, 500 (3): 169171. 10.1016/S00145793(01)026138.
 3.
Gerosa L, Sauer U: Regulation and control of metabolic fluxes in microbes. Curr Opin Microbiol. 2011, 22 (4): 566575.
 4.
Holzer H, Duntze W: Metabolic regulation by chemical modification of enzymes. Annu Rew Biochem. 1971, 40: 345374.
 5.
Bruggeman FJ, de Haan J, Hardin H, Bouwman J, Rossell S, van Eunen K, Bakker BM, Westerhoff HV: Timedependent hierarchical regulation analysis: deciphering cellular adaptation. IEE Syst Biol. 2006, 153 (5): 318322. 10.1049/ipsyb:20060027.
 6.
Bruggeman FJ, Snoep JL, Westerhoff HV: Control, responses and modularity of cellular regulatory networks: a control analysis perspective. IEE Syst Biol. 2008, 2 (6): 397410.
 7.
Feist AM, Herrgard MJ, Thiele I, Reed JL, Palsson BO: Reconstruction of biochemical networks in microbial organisms. Nat Rev Microbiol. 2009, 7 (2): 129143.
 8.
McCloskey D, Palsson BO, Feist AM: Basic and applied uses of genomescale metabolic network reconstructions of Escherichia coli. Mol Syst Biol. 2013, 9: 661
 9.
Fell DA: Understanding the control of metabolism. 1997, London and Miami: Portland Press
 10.
Kacser H, Burns JA: The control of flux. Symp Soc Exp Biol. 1973, 27: 65104.
 11.
Hofmeyr JH: Metabolic regulation: a control analytic perspective. J Bioenerg Biomembr. 1995, 27 (5): 479490. 10.1007/BF02110188.
 12.
Hofmeyr JH, CornishBowden A: Quantitative assessment of regulation in metabolic systems. Eur J Biochem. 1991, 200: 223236. 10.1111/j.14321033.1991.tb21071.x.
 13.
Hofmeyr JH, CornishBowden A: Regulating the cellular economy of supply and demand. FEBS Letters. 2000, 476: 4751. 10.1016/S00145793(00)016689.
 14.
Kahn D, Westerhoff HV: Control theory of regulatory cascades. J Theor Biol. 1991, 153 (2): 255285. 10.1016/S00225193(05)804266.
 15.
Snoep JL, van der Weijden CC, Andersen HW, Westerhoff HV, Jensen PR: DNA supercoiling in Escherichia coli is under tight and subtle homeostatic control, involving geneexpression and metabolic regulation of both topoisomerase I and DNA gyrase. Eur J Biochem. 2002, 269: 16621669. 10.1046/j.14321327.2002.02803.x.
 16.
DaranLapujade P, Rossell S, van Gulik WM, Luttik MA, de Groot MJ, Slijper M, Heck AJ, Daran JM, de Winde JH, Westerhoff HV, Pronk JT, Bakker BM: The flux through glycolytic enzymes in Saccharomyces cerevisiae are predominantly regulated at posttranscriptional levels. Proc Natl Acad Sci USA. 2007, 104: 1575315758. 10.1073/pnas.0707476104.
 17.
Nise NS: Control systems engineering. 2004, Hoboken, NJ, USA: John Wiley & Sons Inc, 4
 18.
Csete ME, Doyle J: Reverse engineering of biological complexity. Science. 2002, 295: 16641669. 10.1126/science.1069981.
 19.
Yi TM, Huang Y, Simon MI, Doyle J: Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc Natl Acad Sci USA. 2000, 97: 46494653. 10.1073/pnas.97.9.4649.
 20.
ElSamad H, Goff JP, Khammash M: Calcium homeostasis and parturient hypocalcemia: an integral feedback perspective. J Theor Biol. 2002, 214: 1729. 10.1006/jtbi.2001.2422.
 21.
Ni XY, Drengstig T, Ruoff P: The control of controller: molecular mechanisms for robust perfect adaptation and temperature compensation. Biophysical Journal. 2009, 97: 12441253. 10.1016/j.bpj.2009.06.030.
 22.
Muzzey D, GomezUribe CA, Mettetal JT, van Oudenaarden A: A systemslevel analysis of perfect adaptation in yeast osmoregulation. Cell. 2009, 138 (1): 160171. 10.1016/j.cell.2009.04.047.
 23.
Cloutier M, Wellstead P: The control systems structures of energy metabolism. J. R. Soc. Interface. 2010, 7: 651665. 10.1098/rsif.2009.0371. 2010
 24.
Tyson JJ, Chen KC, Novak B: Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol. 2003, 15: 221231. 10.1016/S09550674(03)000176.
 25.
Kacser H, Acerenza L: A universal method for achieving increases in metabolite production. Eur J Biochem. 1993, 216 (2): 361367. 10.1111/j.14321033.1993.tb18153.x.
 26.
Fell DA, Thomas S: Physiological control of metabolic flux: the requirement for multisite modulation. Biochem J. 1995, 311: 3539.
 27.
Adamczyk M, Westerhoff HV: Engineering of selfsustaining systems: substituting the yeast glucose transporter plus hexokinase for the Lactococcus lactis phosphotransferase system in a Lactococcus lactis network in silico. Biotech J. 2012, 7 (7): 877883. 10.1002/biot.201100314.
 28.
Goelzer A, Briki FB, MarinVerstraete I, Noirot P, Bessieres P, Aymerich S, Fromion V: Reconstruction and analysis of the genetic and metabolic regulatory networks of the central metabolism of Bacillus subtilis. BMC Systems Biology. 2008, 2: 2010.1186/17520509220.
 29.
Gerhart JC, Pardee AB: The enzymology of control by feedback inhibition. J Biol Chem. 1962, 237: 891896.
 30.
ShenOrr SS, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of Escherichia coli. Nature Genet. 2002, 31: 6468. 10.1038/ng881.
 31.
Berg JM, Tymoczko JL, Stryer L: Biochemistry. 2002, New York: W H Freeman, 5
 32.
Westerhoff HV, Van Dam K: Thermodynamics and control of biological freeenergy transduction. 1987, Amsterdam: Elsevier
 33.
Teusink B, Westerhoff HV: ‘Slave’ metabolites and enzymes. A rapid way of delineating metabolic control. Eur J Biochem. 2000, 267: 18891893. 10.1046/j.14321327.2000.01220.x.
 34.
Schuster S, Kahn D, Westerhoff HV: Modular analysis of the control of complex metabolic pathways. Biophys Chem. 1993, 48 (1): 117. 10.1016/03014622(93)80037J.
 35.
Westerhoff HV, Koster JG, Van Workum M, Rudd KE: On the control of gene expression. Control of metabolic processes. Edited by: CornishBowden A. 1990, New York: Plenum, 399412.
 36.
QuintonTulloch MJ, Bruggeman FJ, Snoep JL, Westerhoff HV: Tradeoff of dynamic fragility but not of robustness in metabolic pathways in silico. FEBs Journal. 2013, 280: 160173. 10.1111/febs.12057.
 37.
Roostalu J, Joers A, Luidalepp H, Kaldalu N, Tenson T: Cell division in Escherichia coli cultures monitored at single cell resolution. BMC Microbiology. 2008, 8: 6810.1186/14712180868.
 38.
VieiraSilva S, Rocha EPC: The systemic imprint of growth and its uses in ecological (meta)genomics. PLoS Genet. 2010, 6 (1): e100080810.1371/journal.pgen.1000808.
 39.
Chin CS, Chubukov V, Jolly ER, DeRisi J, Li H: Dynamics and design principles of a basic regulatory architecture controlling metabolic pathways. PLoS Biol. 2008, 6 (6): 13431356.
Acknowledgements
This study was supported by the Netherlands Organization for Scientific Research NWO through the MOSESSysMO grant, as well as by the BBSRC(BRICBBSRC)/EPSRC’s funding of the Manchester Centre for Integrative Systems Biology (DTC, BB/F003528/1, BB/C008219/1), BBSRC’s MOSESSysMO project (DTC, BB/F003528/1), various other BBSRC projects (BB/G530225/1, BB/I004696/1, BB/I017186/1, BB/I00470X/1, BB/I004688/1, BB/J500422/1, BB/J003883/1, BB/J020060/1), and the EUFP7 projects SYNPOL, ECMOAN, UNICELLSYS, ITFoM, and BioSim. FH is partially supported by EPSRC and the European Research Council (ERC) Advanced Investigator Grant. We thank an anonymous reviewer for suggestions for improvements and additions that we took on board.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
FH and HVW jointly discussed and developed the main ideas of paper, and both authors wrote the manuscript. FH developed the hierarchical supply–demand theory, proposed the control engineering analysis of the regulatory systems, and implemented the simulation studies. HVW provided the MCA and HCA analysis, developed the first ATP metabolism example, identified the three terms in control analysis, and interpreted the results biochemically. VF contributed to the kinetic analysis of the endproduct and initialproduct modules. All authors read and approved the final manuscript and HVW rewrote the final version.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Metabolic control analysis
 Control engineering
 Transcriptional regulation
 Synthetic biology
 Robustness