- Methodology article
- Open Access

# (Im)Perfect robustness and adaptation of metabolic networks subject to metabolic and gene-expression regulation: marrying control engineering with metabolic control analysis

- Fei He
^{1, 3}, - Vincent Fromion
^{2}and - Hans V Westerhoff
^{1, 4, 5}Email author

**7**:131

https://doi.org/10.1186/1752-0509-7-131

© He et al.; licensee BioMed Central Ltd. 2013

**Received:**8 July 2013**Accepted:**12 November 2013**Published:**21 November 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 gene-expression 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 gene-expression 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 gene-expression 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 gene-expression 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.

## Keywords

- Metabolic control analysis
- Control engineering
- Transcriptional regulation
- Synthetic biology
- Robustness

## 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 gene-expression alterations, and (iii) covalent modification via signal transduction. The first is termed metabolic (or enzymatic) regulation. The second is known as gene-expression (mediated) regulation and the third as signal-transduction (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 above-mentioned similarities, the third type is addressed implicitly. Until now, significant progress has been made on the modelling of genome-scale metabolic networks in microorganisms integrating metabolic and gene-expression regulation[7, 8]. The steady-state properties of a number of representative metabolic regulatory mechanisms, such as end-product 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 Cornish-Bowden[11–13]. In order to take gene-expression 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 gene-expression 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 feed-forward 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 gene-expression regulation. In this study, the integration of metabolic and gene-expression 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 so-called ‘manipulated variable’ is non-robust (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 so-called ‘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. so-called ‘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 true-to-life 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 end-product 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, gene-expression 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 steady-state 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 kinetics-based 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 gene-expression 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

*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 gene-expression 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.

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 gene-expression. *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. co-factors 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 gene-expression regulation in general, i.e. includes regulation at the level of transcription, translation and post-translational 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.

*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 steady-state regime exists, and it is the unique solution of the following equation:

*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 gene-expression regulation, respectively. Likewise,

*f*

_{ n }(

*x*

_{ n }) is usually a monotonically increasing function of

*x*

_{ n }. As shown by Equation (2), the steady-state 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 re-defined 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

*supply*and

*demand*blocks as shown in Figure 3.

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

*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 steady-state 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

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

*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

*supply*module, such that${c}_{1}^{{J}_{1}}=1$ and

*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 end-product 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 state-steady properties of the aforementioned full regulatory system, but only for the case of metabolic regulation. When including the gene-expression regulation, a more complete interpretation can be achieved by using the hierarchical control analysis and a new ‘hierarchical supply-demand’ theory as investigated in the next subsection.

#### Gene-expression and metabolic regulation: hierarchical supply–demand theory

Westerhoff and coworkers have developed hierarchical control analysis (HCA), an extension of MCA that can take gene-expression 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.

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

*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 gene-expression regulation:

**The list of symbols and definitions**

Symbol | Definition & comments |
---|---|

${C}_{i}^{f}$ or ${c}_{i}^{f}$ | Metabolic control coefficients as defined in (3). |

${H}_{i}^{f}$ | Hierarchical control coefficient. Its mathematical definition is in the same form as metabolic control coefficients in (3), but the system under study can be more general. MCA only studies the control in a metabolic pathway or a signal transduction cascade, but not their combination. HCA investigates the control in a hierarchical regulatory network with interactions at different levels, i.e. metabolic, signal transduction, and gene-expression. |

${\Re}_{i}^{f}$ | MCA (or HCA) based robustness coefficient. |

${\Re}_{i}^{f}\equiv \frac{1}{\left(\frac{\partial lnf}{\partial ln{e}_{i}}\right)}\equiv \frac{1}{{C}_{i}^{f}}\phantom{\rule{0.5em}{0ex}}\mathit{or}\phantom{\rule{0.5em}{0ex}}\frac{1}{{H}_{i}^{f}}$. | |

${F}_{i}^{f}$ | Fragility coefficient. It is the inverse of the robustness coefficient and identical to the control coefficient. ${F}_{i}^{f}\equiv \frac{\partial lnf}{\partial ln{e}_{i}}\equiv \frac{1}{{\Re}_{i}^{f}}\equiv {C}_{i}^{f}\phantom{\rule{0.5em}{0ex}}\mathit{or}\phantom{\rule{0.5em}{0ex}}{H}_{i}^{f}.$ |

${\u03f5}_{{x}_{i}}^{{v}_{j}}$ | Elasticity coefficient defined in MCA. It denotes the immediate influences of metabolite |

${\u03f5}_{{x}_{i}}^{s}$ or ${\u03f5}_{{x}_{i}}^{d}$ | Elasticity coefficient of metabolite |

$*{\u03f5}_{{x}_{i}}^{{v}_{j}}$, $*{\u03f5}_{{x}_{i}}^{s}$ (or$*{\u03f5}_{{x}_{i}}^{d}$) | The overall elasticity (see[32]) for a reaction step under both metabolic and gene expression regulation, or the overall elasticity in a hierarchical supply or demand module. |

| Reaction rate of transcription or translation |

| Reaction rate of mRNA or protein degradation |

| mRNA or protein degradation rate constant |

${k}_{\text{proteolysis}}^{i}$ | Rate constant of |

| The cells’ specific growth (division) rate |

| Concentration of enzyme |

| Concentration of intermediate metabolite |

| Concentration of mRNA |

| Reference signal. |

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.

*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 non-zero magnitudes of the elasticities, the hierarchical control coefficients for control by supply may be decreased by elasticities in the gene-expression 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.

*x*

_{ n }) regulates the first reaction through both metabolic regulation and gene-expression regulation, the corresponding ‘overall’ elasticity (see Table 1) of the first reaction can be expressed as:

*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 steady-state 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 gene-expression regulation. As an extension to the classical supply–demand theory, the analysis given in this section can be named the ‘hierarchical supply-demand’ theory.

### Control engineering

The discipline of Control Engineering first identifies a so-called *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 so-called ‘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 gene-expression 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 gene-expression 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 gene-expression 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.

*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 closed-loop 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.

*γ∙δ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:

*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 closed-loop 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 non-regulatory 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 gene-expression regulation of an important intracellular process: ATP (energy) metabolism

*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*).

*E*are here assumed to follow mass action kinetics:

*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 non-negative and

*k*

_{0}= 0 whenever

*E*is non-positive. The closed-loop control system structure of the pathway can be represented as in Figure 8.

*E*a manipulated variable in the gene-expression 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 closed-loop control system (16), with ADP for *X*, we recognize on the right-hand 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 gene-expression regulation.

When the degradation of enzyme is zero order in terms of *E* (i.e. when *k*_{
b
} = 0), the gene-expression 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 re-attained in the presence of the sustained perturbation.

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.

*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):

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

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 zero-order 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 trade-off limits what one can do[18]. In systems biology, this can be illustrated for the end-product 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 bi-functional in dynamic terms[18]: The metabolic regulation rapidly rejects high frequency perturbations but possibly with small amplitude or capability, while the gene-expression regulation is slow to adapt, but may be able to reject very large constant perturbations.

### Conditions for integral and pseudo-integral control: end-product pathway

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

*μ*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,

*k*

_{ED}≠ 0), the above pathway example corresponds to a pseudo- or non-integral 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, near-perfect 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}

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

*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 near-perfect adaptation behaviour in practice, which would be interpreted as a quasi-perfect integral control system.

#### Simulation study

*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 steady-state 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.

*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

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.

### Feed-forward activation: a case study of a leucine biosynthetic pathway

*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 gene-expression. 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 end-product module similar to the ones discussed above. The other is the transcriptional (gene-expression) activation of downstream enzymes Leu1 (

*E*

_{1}) and Leu2 (

*E*

_{2}) by αIPM (

*I*

_{1}) (through transcription factor Leu3), which we shall call initial-product 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].

*E*

_{s}with respect to the concentration of

*I*

_{1}(αIPM) can be expressed into the various elasticity coefficients (see also Figure 15),

*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 closed-loop control system structure of the downstream gene-expression activation is the same as that of the upstream inhibition case described in Figure 11 (with end-product metabolic regulation dominating process dynamics). It is still a negative feedback control rather than a feed-forward 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 end-product 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 gene-expression 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. initial-product 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 feed-forward activation. This gives the same result as the end-product 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.

*μ*= 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 quasi-integral 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 steady-state 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 feed-forward 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 so-called ‘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 gene-expression 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 first-order 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 control-engineering 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 gene-expression 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, gene-expression 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 gene-expression 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 initial-product module

The first reaction here is still assumed product insensitive and other factors *p* can act on the enzyme. Different from the end-product module, the function *g* is assumed to be an increasing function of its argument. It has been demonstrated in[28] that if a constant steady-state 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

*supply–demand*system, we have

*supply*and

*demand*steps can be derived as

*supply*module as given in Figure 3.

## C. Steady state analysis of the ATP metabolism example

*E*before the perturbation are determined when

*d*[

*ADP*]/

*dt*= 0 and

*dE*/

*dt*= 0:

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

*k*

_{ d }from its steady state value, one finds for the time dependence of the variation in the enzyme level:

*δE*in (37) into (38) one finds (18). By substituting the steady state condition of

*E*

_{ ss }in (35),

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.

## Declarations

### Acknowledgements

This study was supported by the Netherlands Organization for Scientific Research NWO through the MOSES-SysMO grant, as well as by the BBSRC(BRIC-BBSRC)/EPSRC’s funding of the Manchester Centre for Integrative Systems Biology (DTC, BB/F003528/1, BB/C008219/1), BBSRC’s MOSES-SysMO 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 EU-FP7 projects SYNPOL, EC-MOAN, 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.

## Authors’ Affiliations

## References

- 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: 2166-2171. 10.1073/pnas.0509831103.PubMedPubMed CentralView ArticleGoogle Scholar
- Ter Kuile BH, Westerhoff HV: Transcriptome meets metabolome: hierarchical and metabolic regulation of the glycolytic pathway. FEBS Letters. 2001, 500 (3): 169-171. 10.1016/S0014-5793(01)02613-8.PubMedView ArticleGoogle Scholar
- Gerosa L, Sauer U: Regulation and control of metabolic fluxes in microbes. Curr Opin Microbiol. 2011, 22 (4): 566-575.Google Scholar
- Holzer H, Duntze W: Metabolic regulation by chemical modification of enzymes. Annu Rew Biochem. 1971, 40: 345-374.View ArticleGoogle Scholar
- Bruggeman FJ, de Haan J, Hardin H, Bouwman J, Rossell S, van Eunen K, Bakker BM, Westerhoff HV: Time-dependent hierarchical regulation analysis: deciphering cellular adaptation. IEE Syst Biol. 2006, 153 (5): 318-322. 10.1049/ip-syb:20060027.View ArticleGoogle Scholar
- Bruggeman FJ, Snoep JL, Westerhoff HV: Control, responses and modularity of cellular regulatory networks: a control analysis perspective. IEE Syst Biol. 2008, 2 (6): 397-410.View ArticleGoogle Scholar
- Feist AM, Herrgard MJ, Thiele I, Reed JL, Palsson BO: Reconstruction of biochemical networks in microbial organisms. Nat Rev Microbiol. 2009, 7 (2): 129-143.PubMedPubMed CentralView ArticleGoogle Scholar
- McCloskey D, Palsson BO, Feist AM: Basic and applied uses of genome-scale metabolic network reconstructions of Escherichia coli. Mol Syst Biol. 2013, 9: 661-PubMedPubMed CentralView ArticleGoogle Scholar
- Fell DA: Understanding the control of metabolism. 1997, London and Miami: Portland PressGoogle Scholar
- Kacser H, Burns JA: The control of flux. Symp Soc Exp Biol. 1973, 27: 65-104.PubMedGoogle Scholar
- Hofmeyr JH: Metabolic regulation: a control analytic perspective. J Bioenerg Biomembr. 1995, 27 (5): 479-490. 10.1007/BF02110188.PubMedView ArticleGoogle Scholar
- Hofmeyr JH, Cornish-Bowden A: Quantitative assessment of regulation in metabolic systems. Eur J Biochem. 1991, 200: 223-236. 10.1111/j.1432-1033.1991.tb21071.x.PubMedView ArticleGoogle Scholar
- Hofmeyr JH, Cornish-Bowden A: Regulating the cellular economy of supply and demand. FEBS Letters. 2000, 476: 47-51. 10.1016/S0014-5793(00)01668-9.PubMedView ArticleGoogle Scholar
- Kahn D, Westerhoff HV: Control theory of regulatory cascades. J Theor Biol. 1991, 153 (2): 255-285. 10.1016/S0022-5193(05)80426-6.PubMedView ArticleGoogle Scholar
- 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 gene-expression and metabolic regulation of both topoisomerase I and DNA gyrase. Eur J Biochem. 2002, 269: 1662-1669. 10.1046/j.1432-1327.2002.02803.x.PubMedView ArticleGoogle Scholar
- Daran-Lapujade 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: 15753-15758. 10.1073/pnas.0707476104.PubMedPubMed CentralView ArticleGoogle Scholar
- Nise NS: Control systems engineering. 2004, Hoboken, NJ, USA: John Wiley & Sons Inc, 4Google Scholar
- Csete ME, Doyle J: Reverse engineering of biological complexity. Science. 2002, 295: 1664-1669. 10.1126/science.1069981.PubMedView ArticleGoogle Scholar
- 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: 4649-4653. 10.1073/pnas.97.9.4649.PubMedPubMed CentralView ArticleGoogle Scholar
- El-Samad H, Goff JP, Khammash M: Calcium homeostasis and parturient hypocalcemia: an integral feedback perspective. J Theor Biol. 2002, 214: 17-29. 10.1006/jtbi.2001.2422.PubMedView ArticleGoogle Scholar
- Ni XY, Drengstig T, Ruoff P: The control of controller: molecular mechanisms for robust perfect adaptation and temperature compensation. Biophysical Journal. 2009, 97: 1244-1253. 10.1016/j.bpj.2009.06.030.PubMedPubMed CentralView ArticleGoogle Scholar
- Muzzey D, Gomez-Uribe CA, Mettetal JT, van Oudenaarden A: A systems-level analysis of perfect adaptation in yeast osmoregulation. Cell. 2009, 138 (1): 160-171. 10.1016/j.cell.2009.04.047.PubMedPubMed CentralView ArticleGoogle Scholar
- Cloutier M, Wellstead P: The control systems structures of energy metabolism. J. R. Soc. Interface. 2010, 7: 651-665. 10.1098/rsif.2009.0371. 2010PubMedPubMed CentralView ArticleGoogle Scholar
- 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: 221-231. 10.1016/S0955-0674(03)00017-6.PubMedView ArticleGoogle Scholar
- Kacser H, Acerenza L: A universal method for achieving increases in metabolite production. Eur J Biochem. 1993, 216 (2): 361-367. 10.1111/j.1432-1033.1993.tb18153.x.PubMedView ArticleGoogle Scholar
- Fell DA, Thomas S: Physiological control of metabolic flux: the requirement for multisite modulation. Biochem J. 1995, 311: 35-39.PubMedPubMed CentralView ArticleGoogle Scholar
- Adamczyk M, Westerhoff HV: Engineering of self-sustaining 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): 877-883. 10.1002/biot.201100314.View ArticleGoogle Scholar
- Goelzer A, Briki FB, Marin-Verstraete 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: 20-10.1186/1752-0509-2-20.PubMedPubMed CentralView ArticleGoogle Scholar
- Gerhart JC, Pardee AB: The enzymology of control by feedback inhibition. J Biol Chem. 1962, 237: 891-896.PubMedGoogle Scholar
- Shen-Orr SS, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of Escherichia coli. Nature Genet. 2002, 31: 64-68. 10.1038/ng881.PubMedView ArticleGoogle Scholar
- Berg JM, Tymoczko JL, Stryer L: Biochemistry. 2002, New York: W H Freeman, 5Google Scholar
- Westerhoff HV, Van Dam K: Thermodynamics and control of biological free-energy transduction. 1987, Amsterdam: ElsevierGoogle Scholar
- Teusink B, Westerhoff HV: ‘Slave’ metabolites and enzymes. A rapid way of delineating metabolic control. Eur J Biochem. 2000, 267: 1889-1893. 10.1046/j.1432-1327.2000.01220.x.PubMedView ArticleGoogle Scholar
- Schuster S, Kahn D, Westerhoff HV: Modular analysis of the control of complex metabolic pathways. Biophys Chem. 1993, 48 (1): 1-17. 10.1016/0301-4622(93)80037-J.PubMedView ArticleGoogle Scholar
- Westerhoff HV, Koster JG, Van Workum M, Rudd KE: On the control of gene expression. Control of metabolic processes. Edited by: Cornish-Bowden A. 1990, New York: Plenum, 399-412.View ArticleGoogle Scholar
- Quinton-Tulloch MJ, Bruggeman FJ, Snoep JL, Westerhoff HV: Trade-off of dynamic fragility but not of robustness in metabolic pathways in silico. FEBs Journal. 2013, 280: 160-173. 10.1111/febs.12057.PubMedView ArticleGoogle Scholar
- 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: 68-10.1186/1471-2180-8-68.PubMedPubMed CentralView ArticleGoogle Scholar
- Vieira-Silva S, Rocha EPC: The systemic imprint of growth and its uses in ecological (meta)genomics. PLoS Genet. 2010, 6 (1): e1000808-10.1371/journal.pgen.1000808.PubMedPubMed CentralView ArticleGoogle Scholar
- 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): 1343-1356.View ArticleGoogle Scholar

## Copyright

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.