Abstract
A grand challenge in synthetic biology is to push the design of biomolecular circuits from purely genetic constructs towards systems that interface different levels of the cellular machinery, including signalling networks and metabolic pathways. In this paper, we focus on a genetic circuit for feedback regulation of unbranched metabolic pathways. The objective of this feedback system is to dampen the effect of flux perturbations caused by changes in cellular demands or by engineered pathways consuming metabolic intermediates. We consider a mathematical model for a control circuit with an operon architecture, whereby the expression of all pathway enzymes is transcriptionally repressed by the metabolic product. We address the existence and stability of the steady state, the dynamic response of the network under perturbations, and their dependence on common tuneable knobs such as the promoter characteristic and ribosome binding site (RBS) strengths. Our analysis reveals tradeoffs between the steady state of the enzymes and the intermediates, together with a separation principle between promoter and RBS design. We show that enzymatic saturation imposes limits on the parameter design space, which must be satisfied to prevent metabolite accumulation and guarantee the stability of the network. The use of promoters with a broad dynamic range and a small leaky expression enlarges the design space. Simulation results with realistic parameter values also suggest that the control circuit can effectively upregulate enzyme production to compensate flux perturbations.
1. Introduction
Synthetic biology aims at engineering cellular systems to perform customized and programmable biological functions. The seminal works published in 2000 [1,2] kickstarted the development of a wide range of gene circuits with prescribed functions, including bacterial logic gates [3], mechanisms for programmed celltocell communication [4] and lightresponsive modules [5]. This progress has recently been followed by the socalled ‘second wave’ of synthetic biology [6], which aims at scaling up the designs from individual genetic modules to whole cellular systems that operate across different layers of cellular regulation, including signalling networks and metabolic pathways [7,8].
One of the most prominent applications of synthetic biology is the manipulation of bacterial metabolism for chemical production in sectors such as energy, biomedicine and food technology [6]. Effective control of metabolism hinges on the ability to upregulate or downregulate pathways in response to changes in the intracellular conditions, cell requirements or environmental perturbations [9]. These requirements call for dynamic control strategies that can modulate enzyme expression in a metabolitedependent fashion [10,11]. One of the key bottlenecks in this respect is our limited understanding of how genetic design knobs modulate the metabolic responses.
The goal of this paper is to reveal new insights into the design limitations and tradeoffs arising from the interplay between gene circuits and metabolic pathways. To that end, we analyse a dynamic model for a feedback system comprising nonlinear kinetic equations for the metabolic species, together with productdependent enzyme expression controlled by a synthetic gene circuit. We focus on the existence and stability of the steady state, the dynamic response of the network under perturbations and the dependence of these on the design knobs of the synthetic gene circuit.
Two landmark implementations of engineered genetic–metabolic circuits are the genetic control of lycopene production [12] and the metabolic oscillator described in Fung et al. [13]. These works were followed up by the recent study by Zhang et al. [14], whereby the authors reported the first successful implementation of a genetic control circuit to increase biofuel production. In a way akin to manmade technological systems, the use of feedback control plays a pivotal role in ‘robustifying’ pathway dynamics under changing environmental conditions, celltocell variability and biochemical noise. Despite the ubiquity of control engineering methods [15], only a few works have rigorously addressed the problem of genetic feedback design on the basis of mathematical models. Notably, Anesiadis et al. [16] demonstrated the use of a genetic toggle switch [2] as an ON–OFF controller for metabolism, whereas Dunlop et al. [17] explored different genetic control architectures for biofuel production.
From a control engineering standpoint, catalytic enzymes act as inputs to a metabolic pathway in order to drive the metabolite dynamics (i.e. the outputs). The pathway outputs are then sensed by metaboliteresponsive molecules that can modulate enzyme expression levels (e.g. transcription factors (TFs) or riboswitches [18]). In the control engineering jargon, this feedback system can be seen as a ‘plant’ (i.e. the pathway to be controlled), and a ‘controller’ (i.e. the gene regulatory circuit controlling the expression of the catalytic enzymes); see figure 1a. The design of the genetic controller must then account for two complementary control objectives: firstly, it must dynamically adjust pathway activity to match the cellular demand for product and sustain the homeostatic balance of native cellular processes. Secondly, a common strategy in metabolic engineering is to modify host microbes by expressing heterologous enzymes that convert metabolic intermediates into a chemical of interest [19]. The consumption of intermediates diverts part of the flux allocated to the host native processes (figure 1b), and, therefore, the controller must also alleviate the impact of these engineered pathways on the native flux.
In this paper, we study an unbranched metabolic pathway under transcriptional repression from the product. The synthetic circuit consists of an operon encoding all the catalytic enzymes that is repressed by a productresponsive TF (§2). The operon feedback architecture mimics natural circuits enabling cellular adaptations to environmental perturbations (e.g. in bacterial amino acid metabolism [20] and nutrient uptake [21]). To maintain a general analysis, we do not specify the kinetics of the metabolic model, but rather work with a generic class of enzyme turnover rates satisfying mild assumptions. These are satisfied by a wide range of saturable enzyme kinetics, including Michaelis–Menten kinetics and cooperative behaviour described by sigmoidal kinetics [22]. We parameterize the genetic model in terms of the promoter characteristic and the ribosome binding site (RBS) strengths, which are typical design elements used as tuneable knobs in synthetic biology applications. As with the enzyme kinetics, we do not fix the shape of the promoter characteristic, but rather consider a generic class of repressive functions that account, in particular, for the standard Hill equation model for transcriptional repression [23].
Model analysis revealed that enzymatic saturation and promoter leaky expression limit the RBS strength design space (§3.1). These constraints must be satisfied to guarantee the existence of an equilibrium point, to prevent the accumulation of metabolites and to ensure the stability of the network under small perturbations. The feasible set for the RBS strengths depends critically on the promoter leakiness and substrate availability. Within the feasible set, RBS strengths may be used to finetune the balance between the intermediate metabolite levels and the gene expression burden imposed on the host cell. We also obtained analytical formulae for the modes of the feedback system; these showed that the operon architecture leads to slow fixed modes, and suggests a separation principle between the effect of RBS strengths and the promoter characteristic (§3.2).
We also show that engineered pathways consuming an intermediate add further constraints to the RBS strengths design space, which can be relaxed by using promoters with a high dynamic range and small leakiness (§4.1). We performed numerical simulations of the model with physiologically realistic parameters in Escherichia coli (§§3.2 and 4.2). The simulations show that the control circuit can effectively upregulate enzyme production to compensate an increase in the cell's native demand for product and the impact of engineered pathways. These also suggest that, in terms of both flux and product homeostasis, the synthetic circuit always outperforms an uncontrolled pathway (i.e. with constant enzyme levels), thus highlighting the advantages of using a dynamic feedback control strategy.
2. Unbranched pathway under transcriptional feedback regulation
We consider an unbranched metabolic pathway as in figure 2b, where s_{0} denotes the concentration of substrate, s_{1} and s_{2} are intermediate metabolites and s_{3} is the metabolic product. The metabolic reactions occur at a rate v_{i} (each one catalysed by an enzyme with concentration e_{i}) and d denotes the rate of product consumption by the cell. The metabolic genes are encoded in a single operon controlled by a productresponsive TF that represses enzyme expression. This kind of transcriptional feedback is common, for example, in bacterial nutrient uptake systems (e.g. the lactose operon [21]) and amino acid metabolism (e.g. the tryptophan operon [20]).
2.1. Metabolic pathway
The network in figure 2 exchanges mass with the environment and/or other networks in the cell. The model accounts for this interaction via the input substrate s_{0} and the product consumption rate d. We are interested in biologically meaningful phenotypes, and, therefore, we assume that s_{0} is constant to ensure that, the network can reach a nonzero steady state [24]. Note that, if the substrate decays in time, the network eventually reaches a zero equilibrium, whereby the substrate, intermediate metabolites and product are fully depleted. The constant substrate assumption is also suitable for scenarios where s_{0} is an extracellular substrate pool shared by a lowdensity cell population (so that the effects of celltocell competition are negligible).
In a pathway with n reactions and n metabolites, the rate of change of metabolite concentrations can be described by 2.1for i = 1, 2, … , n and v_{n+1} = d(s_{n}). This model arises from the mass balance between the reactions that produce and consume s_{i}, and the enzyme kinetics are included in the reaction rates v_{i}(s_{i−1}, e_{i}). To keep a general analysis, we will not presuppose a specific form for the enzyme kinetics. Instead, we will generically assume that the metabolic reaction rates are linear in the enzyme concentrations [22] 2.2where is the enzyme turnover rate (i.e. the reaction rate per unit of enzyme concentration) satisfying g_{i}(0) = 0. We will also assume that the enzyme turnover rates are increasing and saturable functions of the metabolite concentrations, so that 2.3and 2.4Assumptions (2.2)–(2.4) account for a broad class of saturable enzyme kinetics, including both irreversible Michaelis–Menten and Hill equation kinetics [22,25].
The rate of product consumption d(s_{n}) is typically modelled as a saturable function of Michaelis–Menten type [20], but, for the sake of generality, we will consider a generic saturable function d(s_{n}) satisfying d(0) = 0, and (cf. assumptions (2.3) and (2.4) for the turnover rates). The cellular demand for product depends on the concentration of a productcatalysing enzyme (which is not explicitly modelled in (2.1)). For typical consumption kinetics such as the Michaelis–Menten or Hill equation, the maximal consumption rate d_{max} is proportional to the enzyme concentration, and therefore in our model we can describe changes in cellular demand as changes in the parameter d_{max}.
2.2. Synthetic gene circuit
In an operon architecture all the enzymes are under the control of a single promoter (for the multipromoter case see [26]), and therefore we model the expression of catalytic enzymes as 2.5for i = 1, 2, …, n. This model comes from the balance between protein synthesis and degradation. We consider a firstorder degradation process with kinetic constants γ_{i}, which accounts for the aggregate effect of degradation and dilution by cell growth. A common strategy in synthetic biology is to control protein degradation by adding a degradation tag to the gene sequence [27], and thus we assume that all enzymes are tagged and degraded at the same rate, i.e. γ_{i} = γ. In the model (2.5), we have parameterized enzyme expression in terms of the promoter characteristic and RBS strengths, both of which are common design elements in synthetic gene circuits (figure 3a):

— Promoter characteristic. It describes the regulatory effect of the TF on gene transcription. The function depends on the specific molecular mechanisms underlying the product–TF and TF–promoter interactions. In order to keep a generic description of the regulatory effect, and to parameterize the model in terms of experimentally accessible design parameters, we opt for a phenomenological description of the promoter characteristic. We therefore consider a function that depends directly on the product concentration and represents the net effect of the product on the transcription rates. Gene transcription under the action of a repressible promoter is typically modelled using Hill functions, but to keep the analysis general we consider generic transcription repression functions satisfying σ(0) = 1, and .
Promoters are typically described in terms of their tightness (κ^{0}) and strength (κ^{1}); see figure 3b. The tightness refers to the level of baseline transcription (i.e. under full repression by the product), whereas the strength is the gap between the ON and OFF transcription levels. The promoter strength is quantified in terms of the dynamic range μ, 2.6Note that, since the promoter strength is always positive, the dynamic range satisfies μ ≥ 1, and we can model the uncontrolled case (i.e. pure constitutive expression without regulation) by taking μ = 1.

— Ribosome binding site strengths. RBSs are mRNA sequences that are bound by the ribosomes to initiate translation [10]. The translation rate of the enzymes can then be modified by choosing RBS sequences with different affinities to ribosome binding [28]. We model the effect of the RBS strengths on the enzyme expression rate via the parameters b_{i}.
With the above assumptions and definitions, we can write the complete model for the feedback system as 2.7In the remainder of the paper, we focus on the existence and stability of the steady state of the model (2.7), its response to perturbations, and its behaviour as a function of the promoter characteristic and the RBS strengths.
3. Circuit design for cellular demands
3.1. Tradeoffs and constraints in the design of ribosome binding site strengths
The operon circuit must be able to sustain a metabolic flux that feeds the product into the downstream native processes of the host. In this section, we show how this essential requirement translates into constraints on the RBS strength design space.
We will denote the steadystate metabolite concentrations, enzyme concentrations and reaction rates as , and , respectively. We first note that the steadystate enzyme concentrations can be obtained by setting in (2.7), leading to 3.1At steady state, the product consumption rate determines the metabolic flux of the network by the relation , and therefore the steadystate product concentration must satisfy . Combining this expression with (3.1) for the first enzyme, we obtain an implicit equation for the steadystate concentration of the product 3.2
The equilibrium concentrations of the intermediates can be obtained by setting for i ≥ 2, 3.3The solution of the implicit equation (3.3) gives the steadystate concentrations of the intermediates as a function of the RBS ratio b_{i}/b_{1}. Note that because the enzyme turnover rates g_{i} are monotonically increasing functions, increasing the b_{i}/b_{1} ratio leads to a lower steadystate concentration of the intermediate s_{i−1}. From equations (3.2) and (3.3), we can infer how the different tuning knobs affect the steady state of the network:

— Effect of the promoter characteristic. From the steadystate equation in (3.2), we can calculate (see appendix A.1 for a detailed derivation) the sensitivity of the product concentration to changes in promoter tightness κ^{0}, promoter strength κ^{1} or RBS strength b_{1}, 3.4 3.5 3.6where . We note that since and , the function is positive and therefore the sensitivities in (3.4)–(3.6) are also positive. We thus conclude that an increase in the promoter parameters (κ^{0} and κ^{1}) or RBS strength b_{1} will lead to a higher steadystate product concentration, which in turn translates into a higher flux. Moreover, since the genes cannot be transcribed at a rate beyond (κ^{0}+κ^{1}), from (3.2), we observe that the flux is constrained by the promoter parameters according to 3.7

— Effect of the RBS strengths. Using (3.1), we can write the steady state of the downstream enzymes as 3.8which indicates that a higher b_{i}/b_{1} ratio leads to a higher concentration for the ith enzyme. Taken together, equations (3.3) and (3.8) indicate that the concentrations of enzymes and intermediates can both be adjusted by tuning the RBS ratio b_{i}/b_{1}. Comparing the dependencies of (3.3) and (3.8) on the RBS ratio reveals a design tradeoff between enzyme expression and the intermediate metabolite concentrations (figure 5a): low b_{i}/b_{1} ratios lead to low enzyme expression levels at the expense of high concentrations for the intermediates. Conversely, high b_{i}/b_{1} ratios tend to increase enzyme expression (and therefore the gene expression burden on the host cell) in favour of low concentrations for the intermediates.
In the above discussion, we have implicitly assumed that a solution to equations (3.2) and (3.3) exists. However, because of the saturable characteristic of the product consumption rate (d) and enzyme kinetics (g_{i}), both equations may lack a solution. Firstly, the solution of (3.2) can be computed as the intersection of the two curves shown in figure 4. From these plots, we can see that an intersection exists only when d_{max}/g_{1}(s_{0} ) > b_{1}κ^{0}/γ, or equivalently 3.9which defines a constraint on the RBS strength of the first enzyme. Since both sides of (3.2) are monotonic in , the solution is unique. By equation (3.1), the existence of also guarantees the existence of the steadystate enzyme concentrations.
Secondly, since the enzyme turnover rates saturate at , equation (3.3) has a finite solution provided that 3.10which defines a constraint on the b_{i}/b_{1} ratio. Taken together, conditions (3.9) and (3.10) define a feasible region in the RBS strengths design space that prevents the accumulation of intermediates and product (figure 5b). If the condition in (3.9) is not satisfied, the substrate will be consumed at a higher rate than the maximal product consumption, and therefore the design will lead to an infinite accumulation of the product. Likewise, violation of at least one of the bounds in (3.10) will cause enzymatic saturation and lead to infinite accumulation of an intermediate.
Conditions (3.9) and (3.10) link together genetic and metabolic parameters (the RBS strengths b_{i} and promoter tightness κ^{0}, together with the substrate availability s_{0} and the enzyme saturation ), and therefore they shed light on how the design constraints appear due to the interplay between metabolic and enzyme expression dynamics. In figure 5c,d, we illustrate the effect of promoter tightness and substrate availability on the feasible region for the RBS strengths. Tighter promoters relax condition (3.9) and therefore enlarge the feasible region (figure 5c). In the limit case of a perfect leakless promoter (i.e. κ^{0} = 0), condition (3.9) does not limit the RBS strength of the first enzyme. Conversely, by conditions (3.9) and (3.10), a higher substrate tends to tighten the feasible region (figure 5d).
3.2. Adaptation to changes in cellular demand
One of the purposes of the genetic feedback circuit is to sustain pathway operation under changes in the cellular demand for product. From a control engineering standpoint, a change in cellular demand can be seen as a perturbation signal acting on the network. A useful approach to study dynamical systems under perturbations consists in examining their linear approximation around their equilibrium points. If we write the model (2.7) as and compute its Jacobian matrix (), then trajectories starting in a small vicinity of the steady state can be approximated as 3.11where λ_{i}, i = 1, 2, … , q, are the q distinct eigenvalues of J evaluated at , q_{i} is the algebraic multiplicity of λ_{i} and the coefficients a_{ij} depend on the initial conditions and the eigenvectors of the Jacobian. The terms in (3.11), known as feedback modes, provide a local approximation of the trajectories around the equilibrium point.
In the case of the feedback system in (2.7), we can exploit the structure of the Jacobian matrix to obtain analytic expressions for its eigenvalues in terms of the design knobs of the gene circuit (see appendix A.2 for details). We found that the 2n eigenvalues can be classified into three categories λ_{fixed}, λ_{RBSi} and λ_{prom}. The system has the following:

— (n−1) stable eigenvalues at λ_{fixed} = −γ < 0. These eigenvalues are independent of the circuit design parameters, and therefore they lead to fixed modes, which can be adjusted only by changing the degradation rate (e.g. with various degradation tags). They cannot be suppressed or changed by tuning the circuit design knobs, and, from (3.11), we see that they translate into (n − 1) modes of the form e^{−t/γ}, t e^{−t/γ}, … , t^{n−2} e^{−t/γ}. The enzyme degradation rates γ are inversely proportional to their halflives, which are in turn much longer than metabolic time scales (enzymatic halflives are of the order of minutes to hours, whereas metabolic time scales are typically milliseconds to seconds [22]). Therefore, depending on the initial conditions the network can potentially display very slow transients, and this appears to be aggravated in long pathways.

— (n−1) stable eigenvalues at 3.12with . Since the steadystate concentration of the enzyme , and the intermediate , depend only on the corresponding RBS ratio b_{i}/b_{1} (see equations (3.3) and (3.8)), this ratio can be used to independently finetune the feedback mode associated with λ_{RBSi}.

— Two stable eigenvalues at 3.13with and . Unlike λ_{fixed} and λ_{RBSi}, these two eigenvalues depend on the steadystate product concentration, and therefore they can be finetuned through the promoter characteristic (see equation (3.2)). To study the dependence of λ_{prom} on the promoter design parameters, we computed them for a pathway with realistic parameter values. In figure 6a, we show the steadystate values of the product, flux and first enzyme level for a wide span of promoter dynamic range μ. We observe that strong promoters tend to increase pathway flux, in agreement with the sensitivity equation previously derived in (3.5). We also see that, as shown by the steadystate relation , the flux corresponds to a scaled version of the concentration of the first enzyme. In figure 6b, we plot the location of the promoterdependent eigenvalues λ_{prom} in the complex plane. These indicate that, in the case of weak promoters, the eigenvalues λ_{prom} lie on the real axis, becoming complex only for a sufficiently broad dynamic range μ. For strong promoters, the real part becomes closer to the imaginary axis, potentially leading to slow transients. Moreover, since stronger promoters lead to a higher flux, the eigenvalues in figure 6b suggest that flux maximization may entail a reduction in the response speed.
To illustrate the dynamic response of a pathway under the control of the transcriptional control circuit, we simulated the network under a change in the cell demand for product (see figure 7a). We modelled a change in the cell demand as a slow Sshaped temporal increase in the maximal product consumption rate d_{max} (see the inset in figure 7a). This describes, for example, cases in which the demand increase is due to native processes upregulating the enzyme that metabolizes the product. Note that, from the steadystate equation in (3.2), a higher d_{max} inevitably leads to a higher flux and a lower steadystate product concentration (see also figure 4a). Before the perturbation, the network is in steady state with a prestimulus flux d^{pre} = 19.5 µM min^{−1}. We considered a Michaelis–Menten consumption rate of the form d(s_{n}) = d_{max}s_{n}/(K_{d} + s_{n}), with K_{d} being the product concentration needed for halfmaximal consumption. Upon the increase in d_{max} at t = 50 min, we observe that the promoter responds to the drop in product concentration and upregulates enzyme expression so as to drive the pathway to a new poststimulus flux d^{post} that is approximately 40 per cent higher than d^{pre}, and a product concentration that is approximately 20 per cent lower than its prestimulus value. Using equations (3.1) and (3.2), we can compute the enzyme upregulation factor as 3.14which is equivalent to the relative change in pathway flux. The dynamic upregulation of enzyme expression can be seen in the lower panel of figure 7a, where we can also verify that the upregulation factor is approximately per cent as predicted in (3.14). Note that as a consequence of the operon architecture, all the enzymes are upregulated by the same foldfactor. This factor depends on the pre and poststimulus fluxes, which by (3.2) depend only on the promoter design and first RBS strength.
As a way of comparison, in figure 7a, we have also simulated the response of a pathway without feedback regulation (i.e. with constant enzyme levels chosen to match the flux of the controlled case d_{pre}). The uncontrolled pathway is unable to increase the flux, and we observe that this leads to a considerable decrease in product concentration (approx. 70% reduction). In the uncontrolled case, the rate of substrate uptake is fixed to v_{1}^{unc} so that the equilibrium product satisfies 3.15and therefore any increase in d_{max} translates into a lower product concentration , which in turn depends only on the kinetic parameters of the consumption rate d(s_{n}). In contrast, the feedbackcontrolled pathway can partly compensate the drop in product by dynamically upregulating enzyme expression, substantially outperforming the uncontrolled case.
We study the performance of the control circuit in more detail in figure 7b, where we show the drop in product concentration relative to the prestimulus level as a function of the change in flux and dynamic range of the promoter. We observe that stronger promoters can significantly improve the compensation of the drop in product concentration (a perfect compensation would correspond to a flat curve at 0% in figure 7b). For example, under a 50 per cent increase in the pathway flux, a mild promoter (μ = 10) leads to a drop in product of approximately 47 per cent, whereas a strong promoter (μ = 100) can bring down the drop in product to approximately 20 per cent (the latter corresponds to the design simulated in figure 7a). As predicted by the upper bound in (3.7), the flux is limited by the promoter strength, and therefore weak promoters do not allow for large increases in flux (as a consequence, the domain of the curves in figure 7a decreases with decreasing promoter strength); for example, for the weakest promoter tested, the flux could not be increased beyond approximately 10 per cent.
4. Circuit design for compensation of flux perturbations
A common strategy in metabolic engineering is to modify bacteria by expressing heterologous enzymes that convert natural metabolic intermediates into a compound of interest [19]. The target compound is synthesized by ‘branching out’ a specific intermediate from a natural pathway, and therefore part of the metabolic flux needed to sustain the host native processes is redirected to the production of the foreign chemical. The choice of a good branching point (i.e. one that does not lead to lethal metabolic imbalances for the host) is a major problem typically addressed with the aid of optimizationbased computational tools [31,32]. In this section, we turn our attention to the effect of a perturbation in the native flux as a consequence of branching out from an intermediate metabolite.
4.1. Tradeoffs and constraints in the design of the RBS strengths
To account for an engineered pathway consuming the intermediate at a constant rate d_{ext}, we include d_{ext} as a consumption rate in the ODE for 4.1The system is shown in figure 8a, and in this case the steadystate equation for the product and the first enzyme is 4.2where the lefthand side is a shifted version of the one in (3.2). For the intermediates before the branch point, the steadystate concentration is given by the same equation as in (3.3) 4.3whereas for the intermediates after the branch point, we have a modified equation 4.4with . From these steadystate equations, we observe similar properties as in the case without a branch. The promoter characteristic and the first RBS strength determine the metabolic flux, whereas the RBS ratio b_{i}/b_{1} can be used to finetune the balance between enzyme expression and the concentrations of the intermediates. In addition, in this case, we see that the intermediates after the branch point also depend on the promoter characteristic.
Using similar arguments to those in figure 4, we find that a solution to (4.2) exists if 4.5and 4.6For equation (4.4) to have a solution, in principle, we need , but this condition is less stringent than the one previously derived in (3.10) for the case d_{ext} = 0. Since the design must also prevent the accumulation of intermediates in the absence of perturbations, we conclude that (3.10), i.e. 4.7is sufficient for the existence of all the intermediates. The inequalities in (4.5)–(4.7) define the feasible region for the RBS design space under a perturbation consuming one of the intermediates (see figure 8b). As in the case without the branch (figure 5b), the limits (4.5) and (4.7) prevent the accumulation of the product and intermediates, respectively. The condition in (4.6), however, adds a new type of constraint to the design space: it guarantees that the synthetic gene circuit can upregulate enzyme expression strongly enough to cope with the flux through the branch, hence preventing the depletion of the product. This new constraint also depends on the promoter dynamic range, which was absent in the case without a branch. From (4.5) and (4.6), we can compute the gap between the upper and lower bounds for the first RBS strength (see figure 8b) 4.8which reveals that promoters with a broad dynamic range and small leakage enlarge the RBS design space (see figure 8c).
4.2. Adaptation to a flux perturbation
To illustrate the effect of an engineered branch on the dynamic response of the feedback system, we simulated a network with two metabolites and two enzymes under a flux perturbation that consumes the intermediate s_{1} (see figure 9a). Before the perturbation, the network is in steady state with a native flux d^{pre} = 19.5 µM min^{−1}. We modelled the engineered branch as an Sshaped increasing rate d^{ext}(t) (see the inset in figure 9a). Upon the activation of the branch, induced at t = 50 min, the synthetic operon circuit upregulates enzyme expression by approximately 45 per cent to drive the pathway to a new native flux d^{post}. Using equation (3.1), together with the pre and poststimulus steadystate equations ((3.2) and (4.2)), we find that the enzymes are upregulated by the factor 4.9
As in the case without the branch, the expression in (4.9) indicates that all enzymes are upregulated by an identical foldfactor that depends on the promoter design and the first RBS strength.
In figure 9a, we have also simulated the response of a pathway without feedback regulation (i.e. with constant enzyme levels chosen to match the flux d^{pre} of the controlled case). In terms of both flux and product concentration, we observe that the feedbackcontrolled network displays a dramatic improvement compared with the uncontrolled case: the operon circuit reduces the loss in native flux from 50 per cent to approximately 5 per cent, whereas the decrease in steadystate product concentration is brought down from approximately 82 per cent to 20 per cent. In the uncontrolled case, the rate of substrate uptake is fixed to v_{1}^{unc} and therefore the poststimulus flux is given by 4.10and hence the poststimulus flux scales linearly with the rate of the branch. In figure 9b, we show this linear dependence together with the feedbackregulated case for a wide span of the promoter dynamic range. We observe that the feedback control circuit outperforms the uncontrolled case even with promoters with a narrow dynamic range, and that this improvement can be achieved with a relatively low enzyme upregulation factor.
5. Discussion and outlook
In this paper we have presented a detailed analysis of a synthetic gene circuit designed to dynamically control metabolic pathways. The goal of this feedback control system is twofold: to adjust pathway activity so as to match the cell demand for product, and to dampen flux perturbations that divert the native flux to the synthesis of foreign molecules. The control strategy relies on encoding the metabolic genes in a single operon repressed by a productresponsive TF. The TF can sense a drop in product concentration and upregulate enzyme expression to bring the pathway close to its homeostatic levels.
Since the seminal operon paper [33], the interaction between the genetic machinery and metabolism has been extensively studied in the context of natural systems. These studies typically focus on understanding how observed phenotypes emerge from the genetic–metabolic cross talk [34–38], and a number of detailed mechanistic models for operon regulation have been developed (e.g. [20,21]). The goal in synthetic biology, however, is to design regulatory circuits for controlling metabolism in a customized fashion. Modelbased design therefore requires mathematical descriptions that are explicitly parameterized in terms of the design knobs that can be manipulated in synthetic biology applications. Consequently, we have used a gene expression model that is deliberately not mechanistic, and instead describes the genetic feedback in terms of tuneable parameters such as the promoter's dynamic range, RBS strengths and protein halflives. This approach has proved to be adequate to explore the genetic design space and to quantify the impact of the promoter characteristic and RBS strengths on the system response.
A typical complication in engineered pathways is that enzymatic saturation may cause intermediates to accumulate in prohibitively large concentrations, thus affecting the viability of the host due to toxic effects [11]. Metabolite accumulation arises when the steady state lies beyond the saturation limit of a catalytic step, and available models for pathways under transcriptional regulation [20,39–41] have generally overlooked the impact of enzyme saturation on the existence of a metabolic steady state. In our aim to carry out a general analysis, we have used a metabolic model that accounts for a whole class of saturable enzyme kinetics under mild assumptions. By explicitly accounting for enzyme saturation, we characterized a feasible set for the design parameters which ensures that the steady state lies within the saturation limits. The feasible set also guarantees the local stability of the network, and we found that the constraints on the RBS strengths can be relaxed with the use of promoters with a high dynamic range and small leakiness. The geometry of the feasible set depends on a combination of genetic and kinetic parameters, thus highlighting the emergence of design constraints as a consequence of the interplay between the genetic and metabolic subsystems.
The steadystate equations reveal a tradeoff between the steadystate enzyme expression levels and the concentration of intermediates: the enzyme concentrations are inversely proportional to the concentration of the intermediate they catalyse. We found that a critical parameter is the RBS ratio, i.e. the relative strength of an RBS with respect to the strength of the first one in the operon, which can be used to finetune the circuit between highenzyme/lowintermediate or lowenzyme/highintermediate designs.
The two considered design knobs, promoter characteristic and RBS strengths, seem to have decoupled roles in the steady state and transient behaviour of the network. The promoter characteristic together with the first RBS determine the steady state of the product and the first enzyme. A strong promoter and a strong RBS for the first enzyme can be used to increase the pathway flux, but this may come at the expense of slow modes in the transient response. In the absence of an engineered pathway consuming an intermediate, the remaining RBS strengths can be used to independently adjust the concentrations of the intermediates and the remaining enzymes. In the case of consumption of an intermediate, however, this design rule applies only to the metabolites upstream of the consumed intermediate, i.e. the steady state of the downstream metabolites depends on a combination of the RBS strength, promoter characteristic and the size of the perturbation.
The closedform expressions for the transient modes of the feedback system show further evidence of the separation principle between promoter and RBS design. From the 2n modes of an nstep pathway, we found that only two depend on the promoter characteristic, whereas further (n − 1) modes depend exclusively on the RBS ratio. The remaining (n − 1) modes correspond to the enzyme halflives and are independent of the promoter characteristic and RBS strengths. Since enzyme halflives are considerably slower than metabolic time constants (even with the use of protein degradation tags), the system dynamics can be dominated by slow transients.
We ran numerical simulations that demonstrate the potential of the proposed control strategy. Using physiologically realistic parameter values for E. coli, the synthetic operon control circuit can dramatically compensate the loss in flux by sensing the drop in product concentration and subsequently upregulating the enzyme concentrations. The feedbackcontrolled pathway outperforms the uncontrolled one even when weak promoters are used, thus underscoring the tremendous advantage of taking a feedback approach to metabolic control.
In this work we focused on a control circuit with an operon architecture, a choice inspired by the fact that operons are one of the building blocks in genomewide bacterial networks [42]. The ubiquity of natural metabolic pathways under operon regulation [43] makes them a reasonable choice as template architectures for engineered circuits. In addition, the main difficulty in building genetic–metabolic systems is to find suitable regulatory molecules to interface a metabolite of interest with the genetic machinery. Some of the available alternatives are engineered promoters [44,45], metaboliteresponsive riboswitches [18,46,47] and natural TFs (for a comprehensive catalogue of natural metaboliteresponsive TFs see Zhang et al. [14, supp. table 5]). In this respect, an operon architecture stands as a simple yet effective topology, as it requires only one metaboliteresponsive TF. More complex architectures can certainly add more flexibility to the design, but this will probably come at the expense of more intricate relationships between the design parameters and the metabolic response. For example, the use of multipromoter circuits allows for independent tuning of the enzyme upregulation factor, but at the same time the pathway may display sustained oscillations if the characteristics of the different promoters are not carefully designed [26].
We should point out that the derived design constraints guarantee the existence and stability of the metabolic steady state, and thus they are only baselines for the correct functioning of the genetic control circuit. In most applications, the design must also account for more demanding objectives such as maximization of flux, minimization of energy expenditure, or a combination of these. Since these objectives may conflict with each other, selecting an appropriate combination of circuit parameters requires the use of multiobjective optimization methods within the feasibility sets derived here (see, for example, figures 5 and 8). Optimization routines can therefore be used to single out the parameter values that lead to an acceptable compromise between mutually colliding objectives; see Banga [48] for a review of a number optimization methods available.
As a consequence of a compromise between model complexity and the generality of the analytic results, our results have two main limitations. Firstly, we have restricted the analysis to pathways with irreversible reactions, and, secondly, our results are limited to unbranched pathways operating in isolation of the remaining metabolism of the host cell. Enzymatic reactions are inherently reversible processes and, although many biosynthetic reactions operate in a regime where the forward reaction is much more likely to occur than its backward counterpart [22], their reversibility cannot always be neglected [49]. In our case, the use of irreversible reactions is an important simplification that allowed the derivation of intuitive and easy to interpret relations between the network parameters and its steady state. Other instances where the analysis of irreversible pathways led to new insights into biological design principles include, for example, the works in [38,43]. Our derivation of the design constraints on the promoter parameters and RBS strengths relies on the structure of the steadystate equations and the fact that most of them are decoupled from each other. However, in the case of an nstep pathway with reversible reactions, the steadystate equations form a system of 2n coupled algebraic equations. These equations may admit an analytic solution for specific enzyme kinetics (see Heinrich & Klipp [50] for the solution in the case of linear and Michaelis–Menten kinetics with constant enzyme concentrations), but its extension to transcriptionally controlled enzymes and general reversible kinetics is cumbersome and lies outside the scope of our paper.
A possible workaround to deal with reversible kinetics is to exploit the natural timescale separation between enzyme expression and metabolic reactions. In this approach, the metabolite trajectories are assumed to evolve in a much faster time scale than the enzyme concentrations. This allows us to approximate the metabolite concentrations as algebraic functions of the enzymes, leading to an enzymeonly ODE model subject to the algebraic relations between metabolites and enzymes. We have previously used such an approach in the case of ON–OFF promoters [36] (i.e. promoters that are either fully active or inactive, without intermediate levels of gene expression), and future work will focus on its use with graded promoters such as those considered here. Another advantage of the time scale separation is that it may allow for the analysis of pathways with more complex stoichiometries. This is of enormous relevance in practical applications, as the cross talk between the controlled pathway and the rest of the host metabolism is likely to have a detrimental impact on the performance of the feedback control system.
We are exploring a number of extensions to this work, aiming primarily at the use of alternative feedback topologies and at quantifying the impact of biochemical noise on the pathway performance. The implementation of genetic–metabolic circuits, let alone parameter finetuning, can be costly and timeconsuming. Our work provides a first step towards understanding the fundamental limitations and tradeoffs that must be addressed at the design stage, potentially facilitating the implementation using a modelguided rationale.
Acknowledgements
We thank Dr John Heap for fruitful discussions and all the members of the Control Engineering for Synthetic Biology group at Imperial College London. We also thank the anonymous referees, whose detailed and constructive feedback considerably improved this paper.
Appendix A
A.1. Derivation of the sensitivities
To obtain the sensitivities in (3.4)–(3.6), we differentiate the steadystate equation in (3.2) with respect to the parameter of interest, and then use the chain rule. For example, to obtain (3.4) we differentiate (3.2) with respect to κ^{0} A1and then solve for to obtain A2where . The remaining sensitivities and can be calculated in an analogous manner.
A.2. Local stability analysis
Here we show that, under the conditions (3.9) and (3.10), the network (2.7) has a locally stable steady state. Moreover, its Jacobian has:
— n−1 stable eigenvalues at λ_{fixed} =−γ;
— n−1 stable eigenvalues at , ; and
— and two stable eigenvalues at
A3
We first write the vector of metabolite and enzyme concentrations as s and e, respectively, so that the model (2.7) can be written as A4
The conditions in (3.9) and (3.10) guarantee the existence of the steady state. The entries of the Jacobian matrix are A5 A6 A7 A8and the coefficients a_{i} are defined as A9and A10Note that because g_{i} and d are nondecreasing, it follows that a_{i} ≥ 0. Using (3.3), we can write the terms in terms of the RBS ratio A11We can therefore write the Jacobian as a block matrix A12where the four blocks are matrices A13 A14 A15
The characteristic polynomial of J (i.e. p(λ) = det( J − λI)) is A16where the product J_{12}J_{21} is A17and q = g_{1}(s_{0} )σ′κ^{1}b_{1}. From the structure of J_{11} and the product J_{12}J_{21}, we can carry one (λ + γ) term into the determinant in (A 16) and then into the last column of its argument. This leads to A18The factor is a polynomial given by A19We therefore conclude that the Jacobian J has (n−1) stable eigenvalues at λ=−γ < 0, (n−1) stable eigenvalues at λ=−a_{i} < 0, for i = 2, 3, …, n, and two eigenvalues at A20These last two eigenvalues are also stable because and imply that the quadratic polynomial A21has only positive coefficients and therefore its roots satisfy .
 Received August 20, 2012.
 Accepted September 20, 2012.
© 2012 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.