Synthetic gene circuits for metabolic control: design trade-offs and constraints

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 trade-offs 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.


Introduction
Synthetic biology aims at engineering cellular systems to perform customized and programmable biological functions. The seminal works published in 2000 [1,2] kick-started the development of a wide range of gene circuits with prescribed functions, including bacterial logic gates [3], mechanisms for programmed cell-to-cell communication [4] and light-responsive modules [5]. This progress has recently been followed by the so-called '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 metabolite-dependent 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 trade-offs 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 geneticmetabolic 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 man-made technological systems, the use of feedback control plays a pivotal role in 'robustifying' pathway dynamics under changing environmental conditions, cell-to-cell 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 metabolite-responsive 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 product-responsive 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 fine-tune 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.

Unbranched pathway under transcriptional feedback regulation
We  rsif.royalsocietypublishing.org J R Soc Interface 10: 20120671 single operon controlled by a product-responsive 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]).

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 non-zero 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 low-density cell population (so that the effects of cell-to-cell competition are negligible).
In a pathway with n reactions and n metabolites, the rate of change of metabolite concentrations can be described by ð2:1Þ for 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 i21 , 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] v i ðs iÀ1 ; e i Þ ¼ g i ðs iÀ1 Þe i ; ð2:2Þ where g i ðÁÞ 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  [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, d 0 ¼ @dðs n Þ=@s n . 0 and lim sn!1 dðs n Þ ¼ d max (cf. assumptions (2.3) and (2.4) for the turnover rates). The cellular demand for product depends on the concentration of a product-catalysing 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 .

Synthetic gene circuit
In an operon architecture all the enzymes are under the control of a single promoter (for the multi-promoter case see [26]), and therefore we model the expression of catalytic enzymes as Àg i e i ; ð2:5Þ for i ¼ 1, 2, . . ., n. This model comes from the balance between protein synthesis and degradation. We consider a first-order degradation process with kinetic constants g 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. g i ¼ g. 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 sðÁÞ 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 sðÁÞ 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 s(0) ¼ 1, s 0 ¼ @s=@s n , 0 and lim sn!1 sðs n Þ ¼ 0.
Promoters are typically described in terms of their tightness (k 0 ) and strength (k 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 m,  rsif.royalsocietypublishing.org J R Soc Interface 10: 20120671 [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 . . . ; n: ð2:7Þ In 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.

Trade-offs 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 steady-state metabolite concentrations, enzyme concentrations and reaction rates as s i , e i and v i , respectively. We first note that the steady-state enzyme concentrations can be obtained by setting _ e i ¼ 0 in (2.7), leading to At steady state, the product consumption rate dð s n Þ determines the metabolic flux of the network by the relation v 1 ¼ dð s n Þ, and therefore the steady-state product concentration must satisfy dð s n Þ ¼ g 1 ðs 0 Þ e 1 . Combining this expression with (3.1) for the first enzyme, we obtain an implicit equation for the steady-state concentration of the product The equilibrium concentrations of the intermediates can be obtained by setting g i ð s iÀ1 Þ e i ¼ g 1 ðs 0 Þ e 1 for i ! 2, ; for all i ! 2: ð3:3Þ The solution of the implicit equation where Fð s n Þ ¼ ðd 0 ð s n Þ=g 1 ðs 0 Þ À ðb 1 k 1 =gÞs 0 ð s n ÞÞ À1 . We note that since d 0 ð s n Þ . 0 and s 0 ð s n Þ , 0, the function FðÁÞ is positive and therefore the sensitivities in (3.4) -(3.6) are also positive. We thus conclude that an increase in the promoter parameters (k 0 and k 1 ) or RBS strength b 1 will lead to a higher steady-state product concentration, which in turn translates into a higher flux. Moreover, since the genes cannot be transcribed at a rate beyond (k 0 þ k 1 ), from (3.2), we observe that the flux is constrained by the promoter parameters according to 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 ) which defines a constraint on the RBS strength of the first enzyme. Since both sides of (3.2) are monotonic in s n , the solution is unique. By equation (3.1), the existence of s n also guarantees the existence of the steady-state enzyme concentrations.
Secondly, since the enzyme turnover rates saturate atĝ i , equation (3.3) has a finite solution provided that b i b 1 . 2) can be seen as the intersection of two curves, The intersection does not exist when condition (3.9) fails.
low RBS ratio  rsif.royalsocietypublishing.org J R Soc Interface 10: 20120671 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 k 0 , together with the substrate availability s 0 and the enzyme saturationĝ i ), 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 leak-less promoter (i.e. k 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).

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 _ x ¼ FðxÞ and compute its Jacobian matrix ( J ¼ @F=@x), then trajectories starting in a small vicinity of the steady state x can be approximated as . . , q, are the q distinct eigenvalues of J evaluated at x ¼ x, q i is the algebraic multiplicity of l i and the coefficients a ij depend on the initial conditions and the eigenvectors of the Jacobian. The terms t jÀ1 e Àlit 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 l fixed , l RBSi and l prom . The system has the following: -(n 2 1) stable eigenvalues at l fixed ¼ 2 g , 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 2 1) modes of the form e 2t/g , t e 2t/g , . . . , t n22 e 2t/g . The enzyme degradation rates g are inversely proportional to their half-lives, which are in turn much longer than metabolic time scales (enzymatic half-lives 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 2 1) stable eigenvalues at l RBSi ¼ À e i g 0 i ð s iÀ1 Þ , 0 for i ¼ 2; 3; . . . ; n; ð3:12Þ with g i ð s iÀ1 Þ ¼ @g i =@s iÀ1 j s iÀ1 . 0. Since the steady-state concentration of the enzyme e i ; i ! 2, and the intermediate s i ; i , n, 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 fine-tune the feedback mode associated with l RBSi . -Two stable eigenvalues at with d 0 ¼ d 0 ð s n Þ . 0 and s 0 ¼ s 0 ð s n Þ , 0. Unlike l fixed and l RBSi , these two eigenvalues depend on the steady-state product concentration, and therefore they can be finetuned through the promoter characteristic (see equation (3.2)). To study the dependence of l prom on the promoter design parameters, we computed them for a pathway with realistic parameter values. In figure 6a, we show the steady-state values of the product, flux and first enzyme level for a wide span of promoter dynamic range m. 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 dð s n Þ ¼ g 1 ðs 0 Þ e 1 , the flux corresponds to a scaled version of the concentration of the first enzyme.
In figure 6b, we plot the location of the promoter-dependent eigenvalues l prom in the complex plane. These indicate that, in the case of weak promoters, the eigenvalues l prom lie on the real axis, becoming complex only for a sufficiently broad dynamic range m. For strong promoters, the real part <fl prom g 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 S-shaped 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 steady-state equation in (3.2), a higher d max inevitably leads to a higher flux and a lower steady-state product concentration (see also figure 4a). Before the perturbation, the network is in steady state with a pre-stimulus flux d pre ¼ 19.5 mM min 21 . 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 half-maximal 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 post-stimulus 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 pre-stimulus value. The real and imaginary parts in (b) are normalized to the degradation rate g. We considered a pathway with two metabolites and two enzymes with Michaelis -Menten kinetics (g i (s i21 ) ¼ k cat i s i21 /(K M i þ s i21 )) with k cat i ¼ 32s 21 and K M i ¼ 4.7 mM. These are representative values for PRA isomerase (extracted from the BRENDA database [29], EC number 5.3.1.24), a transcriptionally regulated enzyme in the tryptophan pathway of E. coli. We took the enzyme degradation rate as 2 Â 10 24 s 21 (half-life 1 h), and used a substrate concentration s 0 ¼ K M1 (so that g 1 (s 0 ) is at half-saturation). We modelled the product consumption as a Michaelis -Menten function d(s n ) ¼ d max s n /(K d þ s n ) with d max ¼ 24.96 mM min 21 and K d ¼ 0.2 mM, both taken from an experimentally validated model for the tryptophan pathway [20]. Promoter tightness was fixed to k 0 ¼ 0.03 nM min 21  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 dð s unc n Þ ¼ v unc 1 ; ð3:15Þ and therefore any increase in d max translates into a lower product concentration s unc n , which in turn depends only on the kinetic parameters of the consumption rate d(s n ). In contrast, the feedback-controlled 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 pre-stimulus 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 (m ¼ 10) leads to a drop in product of approximately 47 per cent, whereas a strong promoter (m ¼ 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.

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 optimization-based 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.

Trade-offs and constraints in the design of the RBS strengths
To account for an engineered pathway consuming the intermediate s ' at a constant rate d ext , we include d ext as a consumption rate in the ODE for s ' Þe n À dðs n Þ; . . . ; n: The system is shown in figure 8a, and in this case the steady-state equation for the product and the first enzyme is where the left-hand side is a shifted version of the one in (3.2). For the intermediates before the branch point, the steady-state concentration is given by the same equation as in (3.3) whereas for the intermediates after the branch point, we have a modified equation stronger and tighter promoter For equation (4.4) to have a solution, in principle, we need ðg 1 ðs 0 Þb 1 =b i À d ext = e i Þ ,ĝ i , 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.
is 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) which reveals that promoters with a broad dynamic range and small leakage enlarge the RBS design space (see figure 8c).

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.
As in the case without the branch, the expression in (4.9) indicates that all enzymes are upregulated by an identical fold-factor 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 e t ðtÞ ¼ e pre i chosen to match the flux d pre of the controlled case). In terms of both flux and product concentration, we observe that the feedback-controlled 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 steady-state 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 post-stimulus flux is given by and hence the post-stimulus flux scales linearly with the rate of the branch. In figure 9b, we show this linear dependence together with the feedback-regulated 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.

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 product-responsive 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][35][36][37][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. Model-based 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 half-lives. 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][40][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 steady-state equations reveal a trade-off between the steady-state 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 high-enzyme/low-intermediate or low-enzyme/high-intermediate 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 closed-form 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 n-step pathway, we found that only two depend on the promoter characteristic, whereas further (n 2 1) modes depend exclusively on the RBS ratio. The remaining (n 2 1) modes correspond to the enzyme halflives and are independent of the promoter characteristic and RBS strengths. Since enzyme half-lives are considerably slower than metabolic time constants (even with the use of rsif.royalsocietypublishing.org J R Soc Interface 10: 20120671 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 feedback-controlled 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 genome-wide 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], metabolite-responsive riboswitches [18,46,47] and natural TFs (for a comprehensive catalogue of natural metabolite-responsive 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 metabolite-responsive 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 multi-promoter 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 steady-state equations and the fact that most of them are decoupled from each other. However, in the case of an n-step pathway with reversible reactions, the steady-state 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 enzyme-only 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 fine-tuning, can be costly and time-consuming. Our work provides a first step towards understanding the fundamental limitations and trade-offs that must be addressed at the design stage, potentially facilitating the implementation using a model-guided rationale.