## Abstract

Recent progress in synthetic biology allows the construction of dynamic control circuits for metabolic engineering. This technology promises to overcome many challenges encountered in traditional pathway engineering, thanks to their ability to self-regulate gene expression in response to bioreactor perturbations. The central components in these control circuits are metabolite biosensors that read out pathway signals and actuate enzyme expression. However, the construction of metabolite biosensors is laborious and currently a major bottleneck for strain design. Here we present a general method for biosensor design based on multiobjective optimization. Our approach produces libraries of biosensors that optimally trade-off production flux against the genetic burden on the host. We explore properties of control architectures built in the literature, and identify their advantages and caveats in terms of performance and robustness to growth conditions or leaky promoters. We demonstrate the optimality of a control circuit for glucaric acid production in *Escherichia coli*, which has been shown to increase titer by 2.5-fold as compared to static designs. Our results lay the groundwork for the automated design of control circuits for pathway engineering, with applications in the food, energy and pharmaceutical sectors.

## I. INTRODUCTION

Metabolic engineering has led to a wealth of chemicals produced with microbial hosts. Classic pathway engineering employs combinations of heterologous expression and knock-downs of native enzymes of the host, so as to redirect metabolic flux toward production pathways^{1,2}. Progress in synthetic biology and genetic engineering now allow the construction of control circuits that sense pathway activity and actuate enzyme expression in real time. These feedback systems are an attractive strategy to build pathways that self-tune to fermentation conditions and thus overcome many of the challenges encountered in pathway engineering. For example, it is generally challenging to find the right enzyme expression levels that maximise production and, at the same time, avoid excessive genetic burden on the host, flux bottlenecks, or the accumulation of toxic intermediates. Dynamic control circuits can resolve these challenges by self-adjusting enzyme expression in response to metabolic signals and continuously controlling pathway activity in response to intracellular or bioreactor perturbations.

Almost a decade ago, Zhang and colleagues reported the first successful implementation of such a control circuit, which was designed to improve production of fatty acid ethyl esters in *Escherichia coli*^{3}. Since then there has been a sharp increase in the number of pathways and organisms where this technology has been deployed, and numerous reviews have discussed the opportunities and challenges offered by dynamic pathway engineering^{2,4–8}. Control circuits have been built into diverse production pathways, including fatty acids^{9}, carotenoids^{10}, and other high-value chemicals^{11,12}, using industrially-relevant hosts such as *Saccharomyces cerevisiae*^{13,14} and *Bacillus subtilis*^{15,16}.

Current implementations of dynamic control differ substantially on their components and architectures, but generally speaking, control circuits require biosensors to read out metabolite concentrations and up-/down-regulate enzyme expression accordingly. Metabolite biosensors are thus central to dynamic pathway engineering, and as a result substantial efforts have been made toward their characterisation^{17}. For example, a number of works have focused on the tuneability of metabolite-responsive transcription factors (TFs) using computational methods^{18} or large screens of biosensor dose-response curves^{19,20}. Metabolite-responsive TFs have become the most widespread mechanisms for sensing because they can be tuned with many techniques, including promoter engineering^{18,21} and protein engineering^{22,23}, and can be repurposed across species and pathways^{24,25}.

Despite enormous progress in biosensor design, their construction is costly, labour intensive, and significantly slows down strain design. In particular, the dose-response curve of metabolite-responsive TFs depends on a complex interplay among various processes, such as metabolite-TF binding, TF binding to its target promoter site, and the expression dynamics of the TF itself. Such processes are poorly characterised in all but a handful of well studied TFs (e.g. the LacI repressor). Therefore, prior to implementation it is unclear how the shape of the dose-response curve affects pathway production, which ultimately leads to substantial efforts for screening large libraries of biosensors.

The design of dynamic control systems requires a quantitative understanding of how dose-response curves relate to pathway activity. Such analyses can be done on the basis of kinetic models based on ordinary differential equations (ODEs). Works in this space have primarily focused on the impact of parameters on metabolic flux using dynamical systems^{26,27}, model simulation^{28–;30}, and control theory^{31–33}. Other works have looked at how various control architectures, i.e. combinations of negative and positive feedback loops, affect the steady state response of engineered pathways, for example using exhaustive sampling of the design space^{34}, or through bespoke mathematical analyses^{35}. Most methods developed so far have been tailored *ad hoc* to specific pathways, and to date there are no generally applicable methods for dynamic pathway engineering.

Here we present a method for biosensor design that can be employed in a wide range of metabolic engineering problems. Using multiobjective optimization^{36–38}, we show that metabolite-responsive TFs can be designed to optimally tradeoff competing objectives. We focus on designs that weigh production flux against the genetic burden imposed by expression of heterologous enzymes. In principle, high flux can be achieved with strong enzyme expression. But in many cases, production performance is counteracted by the negative impact of heterologous expression on the host. Heterologous expression draws molecular resources away from native genes, such as RNA polymerases, σ-factors, tRNAs and ribosomes. This phenomenon is generally known as “genetic burden”^{39} and can lead to homeostatic imbalances that impair growth and ultimately limit production^{40,41}. This causes production flux and genetic burden to conflict with each other, in the sense that higher production flux entails stronger enzyme expression and thus higher burden.

To resolve the conflict between production and burden, we propose a method to find libraries of biosensors that trade-off both objectives. As shown in the diagram in Fig. 1A, each biosensor dose-response curve (left) corresponds to a point in the performance space (right), defined by the values of two objective functions to be jointly minimised. The boundary of such point cloud is the *Pareto front*, i.e. the curve containing all optimal trade-off solutions whereby no objective function can be further reduced without increasing the value of the other^{36}. Therefore, points along the Pareto front correspond to specific dose-response curves from which designers can choose biosensors to match a desired trade-off between production and burden. In Fig. 1B we illustrate the concept behind our proposed pipeline. Starting from a kinetic model for a pathway of interest plus a list of candidate control architectures and toxicity constraints, our method produces libraries of biosensor dose-response curves that optimally navigate the flux vs. burden space.

We show the effectiveness of our approach in a simple model for a branched pathway that contains many features found in metabolic engineering problems, such as nonlinear kinetics, toxic intermediates, and a limited carbon supply. Using this exemplar system, we explore the advantages and caveats of various control architectures implemented in the literature so far. We identified a number of trade-offs between optimality and circuit robustness to fluctuations in carbon sources and leaky promoters. We further tested our method in a more complex pathway for production of glucaric acid in *E. coli*, a high-value precursor for a number of applications. Based on published data, we built a detailed kinetic model for the pathway and computed optimal biosensors for different control architectures based on negative feedback. Our results show that a recently implemented control architecture^{11} out-performs alternative implementations of negative feedback, thus demonstrating the untapped potential of multiobjective optimization for dynamic pathway engineering.

## II. RESULTS

### A. Multiobjetive optimization of pathway dynamics

We focus on heterologous pathways described by kinetic models of the form:
where *x* and *e* are vectors of metabolite and enzyme concentrations, respectively. The first-order term represents the dilution effect on molecular concentrations caused by cell growth with rate constant λ. The function *f(x)* describes the kinetic model in terms of mass balance equations, whereas the term *u(x)* models the rate of expression of pathway enzymes.

In the case of static pathways, heterologous enzymes are expressed at constitutive rates so that *u(x)* in Eq. (1) is constant and independent of pathway intermediates. Conversely, for pathways under dynamic control, the expression rates of pathway enzymes are made dependent on metabolite concentrations. In such case the term *u(x)* describes the doseresponse curve of a biosensor that employs metabolite readouts to control enzyme expression. Different control architectures can thus be implemented with combinations of activation and repression feedback loops, encoded in the specific shape of the dose-response curve.

To illustrate the utility of our approach, we first studied a prototypical pathway that contains many of the features encountered in metabolic engineering problems. As shown in Fig. 2A, we consider the production of a metabolite *x*_{1} from a native intermediate *x*_{0} in central carbon metabolism. The two-step production pathway is catalyzed by heterologous enzymes E_{1} and E_{2}. We model each reaction using Michaelis-Menten kinetics and assume that a carbon source is taken up at a constant rate *V _{in}*; we have included the detailed model in the Methods section and all model parameters can be found in Supplementary Data.

To gain an initial understanding of the design trade-offs of this model system, we first optimized its enzyme expression levels in open-loop, i.e. in the absence of dynamic control. This corresponds to the traditional approach in static pathway engineering, whereby heterologous enzymes are expressed at constant levels. In this case, the concentrations of both heterologous enzymes follow the ODE model:
where *a _{i}* is the rate of constitutive expression. To model the trade-offs between production and burden, we define two objective functions to be minimised:

where *u _{i}(t)* is the temporal expression rate of the i

^{th}pathway enzyme, and

*T*is an optimization horizon after activation of heterologous expression at

*t*= 0. In the open-loop case, the enzyme expression rates are

*u*=

_{i}(t)*a*as in Eq. (2). The chosen objective functions, illustrated in Fig. 2B, quantify the production performance (

_{i}*J*

_{1}) and the genetic burden (

*J*

_{2}). The objective

*J*

_{1}describes the loss in pathway production, so that low values of

*J*

_{1}imply a high production flux

*v*

_{2}(

*t*) close to the carbon influx

*V*

_{in}. The objective

*J*

_{2}, conversely, corresponds to the total concentration of heterologous protein expressed during the temporal horizon. Note that for convenience, we have defined performance as a

*production loss*, so that both objectives need to be jointly minimised.

We computed the enzyme expression rates ai as solutions to the multiobjective optimization problem of the form min (*J*_{1}, *J*_{2}); the full statement of the optimization problem is shown in Table IA. We imposed an upper bound on the promoter strengths *a _{i}* < 1 × 10

^{−3}

*μ*Ms

^{−1}to model the limited biosynthetic capacity of the host; this constraint corresponds to maximal enzyme concentrations of ~ 5

*μ*M in fast growing

*E. coli*with doubling time of 26 min. We additionally introduced a constraint on the intermediate along the whole optimization horizon (

*x*180 ≤

_{1}(t)*μ*M), to model cytotoxic effects commonly encountered in pathway engineering

^{34}.

We solved the optimization problem using a genetic algorithm detailed in the Methods section; all solutions can be found in Supplementary Data. The resulting Pareto front, shown in Fig. 2C, has three distinct regimes: (i) low production and light burden, (ii) high production and heavy burden, and (iii) an optimal trade-off regime between the two. These regimes are in agreement with the intuition that stronger enzyme expression levels (shown in insets of Fig. 2C) lead to higher production flux at the expense of an increased cost in protein expression. The convexity of the Pareto front indicates that the optimization problem is well-posed, in the sense that both objective functions oppose each other across the whole space of optimal expression levels.

### B. Optimal biosensors for negative feedback control

We next turned our attention to the design of dynamic control circuits for the general branched pathway of previous section. As shown in Fig. 3A, we consider three control architectures in which the biosensor detects the concentration of intermediate (*x _{1}*) and controls enzyme expression accordingly. For generality, we assume that the biosensor is implemented via a metabolite-responsive transcription factor (TF), which is currently the most widespread strategy for implementation of metabolite biosensors

^{2}.

In this case, as in Eq. (1) the expression of heterologous enzymes follows the ODE model:
where *u _{i}(x_{1})* is the dose-response curves of the biosensor with respect to the the concentration of the intermediate. In this case, there are a total 3

^{2}− 1 = 8 possible control architectures, because each of the two enzymes can be subject to three modes of control (positive, negative, and no control) and we exclude the open-loop case studied in the previous section. Binding between the TF and the intermediate

*x*

_{1}causes a reduction in the pool of regulators that can bind to the operator region of enzymatic genes. Therefore, the metabolite exerts negative regulation on the TF by molecular sequestration, a mechanism commonly found in many biosensors employed in the literature so far, including BetI

^{42}, FadR

^{3}, IpsA

^{11}, FapR

^{9}, LacI

^{23}, and others.

As shown in Fig. 3A, we focused on three different implementations of negative feedback: a transcriptional repressor that downregulates upstream enzymes (upstream repression), a transcriptional activator that upregulates downstream enzymes (downstream activation), or a TF with dual regulatory function that combines both control strategies (dual control). Negative feedback is known to offer various benefits, including robustness^{26}, accelerated response times^{43}, and noise control^{35}. The three architectures have been implemented in several metabolic control systems in *E. coli*. For example, an upstream repression circuit was employed for speeding up the production of fatty acids using the regulator FadR as biosensor^{43}. The downstream activation circuit has been recently built for production of glucaric acid with IpsA as a metabolite-responsive biosensor^{11}. The dual control architecture was employed for the production of malonyl-CoA using FapR as a dual transcriptional regulator^{9}.

By optimizing the biosensor dose-response curves *u _{i}*(

*x*) for each control architecture, we computed libraries of biosensors that optimally navigate the production-burden space. To this end, we modelled the dose-response curves as sigmoids: where

_{1}*a*is the promoter strength,

_{i}*θ*is a regulatory threshold, and

_{i}*h*is a Hill coefficient to model cooperative effects. Note that the sigmoidal dose-response curves in Eq. (5) are effective models that lump together various molecular processes, including metabolite-TF and TF-DNA binding. As a result, the parameters of the dose-response curve depend on a combination of experimentally tunable parameters such as binding affinities, operator sequence and others

_{i}^{18,20}. For simplicity, we fixed

*h*because it is particularly challenging to tune experimentally; we have chosen

_{i}*h*= 2 to capture the effect of cooperativity on the basis that many biosensors require the TF to dimerise before binding to the operator site of the promoter.

_{i}To optimize the biosensor parameters *a _{i}* and

*θ*, for each architecture we solved a multiobjective optimization problem of the form min (

_{i}*J*) subject to constraints on maximal enzyme expression (

_{1}, J_{2}*a*), biosensor thresholds (

_{i}*θ*), and a toxicity constraint on the intermediate

_{i}*x*as in the previous case. The full statement of the optimization problem is shown in Table IB, and all solutions are given in Supplementary Data.

_{1}As shown in Fig. 3B, to quantitatively compare the Pareto fronts of each circuit we computed their minimum distance to the origin (indicative of the optimal trade-off achieved by each architecture) and the minimal production loss (representing the ability to maximise the production flux). Overall, we observe that the three circuits achieve similar Pareto fronts with a nearly identical optimal trade-off. Yet we found that the minimal production loss is ~30% higher in the downstream activation circuit, which means that the best production achievable by downstream activation is substantially worse than the other two circuits. Taken together, the analysis suggests that upstream repression and dual control are similarly strong in terms of performance. However, this also suggests that the extra cost of implementing the more complex dual control system does not offer advantages in terms of the trade-offs considered in this paper.

Detailed inspection of the biosensor libraries reveals that architectures require different designs to optimize performance (Fig. 3C). For example, both upstream repression and downstream activation require increasing promoter strengths (*a _{i}*) to increase production. Yet the optimal promoter regulatory thresholds (

*θ*) exhibit opposite trends. In the upstream repression circuit, increasing the regulatory threshold helps to increase production, but the downstream activation circuit requires a lowered threshold to increase production. In the case of the dual control architecture (Fig. 3C), we observe a strong asymmetry between promoters. The upstream promoter needs to be at least 5-fold stronger than the downstream one, and their regulatory thresholds follow opposite trends, similarly as in the first two architectures. Altogether these results suggest that strategies for fine-tuning the trade-off between production and burden are highly dependent on the chosen control architecture.

_{i}### C. Robustness and sensitivity of the control circuits

The optimization results in Fig. 3A suggest that the Pareto fronts of the three negative feedback circuits are similar to the open-loop case in Fig. 2C. While this seems counter-intuitive, from control theory it is known that open-loop control can indeed outperform feedback systems, at the expense of their robustness to perturbations^{44}. Open-loop solutions are fragile in the sense that their performance can significantly degrade in face of perturbations or inaccuracies in the model of the system to be controlled. We hence sought to identify salient properties of the circuits in terms of their robustness to perturbations commonly encountered in applications: a) fluctuations in growth conditions, and b) leaky expression of promoters.

#### a. Robustness to growth conditions

To compare architectures in terms of their robustness to carbon supply, we challenged the optimal circuits with perturbations to the influx *V*_{in}. We first simulated each optimal circuit to steady state under a constant *V*_{in}, and then introduced a step perturbation of a given size. As shown in Fig. 4A, we obtained a relation between changes to Vin and the resulting change in production flux *V*_{2} for each circuit architecture and each combination of optimized parameters. Ideally, the production flux should be insensitive to growth conditions; this would translate into a flat curve in Δ *V*_{2} vs Δ *V*_{in} in Fig.4A. At the other end, fragile designs would produce a 45 ° line in the Δ*V*_{2} > vs Δ *V*_{in} space, meaning that the circuit is completely unable to compensate for perturbations in carbon influx. The latter case is what happens in static pathways, whereby carbon flux perturbations are relayed directly onto the flux of production pathways. In contrast, by sensing changes in the concentration of the intermediate (*x _{1}*), the feedback circuits adjust enzyme expression so as to keep the production flux close to its level prior to the perturbation.

We observed substantial variations in the response of the three control circuits to perturbations in the carbon influx (see Supplementary Figure 1). For example, the downstream activation circuit fails to compensate for drops in carbon influx, while the other two architectures display a compensatory response across both negative and positive perturbations. To quantitatively explore such differential responses to perturbations, we defined the robustness score:
Where and are relative changesin the both steady state fluxes, and is the optimal steady state production flux achieved with a constant . The second term in Eq. (6) is the slope of the Δ*V*_{2} vs Δ *V*_{in} curves when Δ*V*_{in} = 0 as marked in Fig. 4A; all the computed curves for each control circuits can be found in Supplementary Figure 1.

Under the *ø* metric, robust circuits would score *ø* ≈ 1 while fragile designs would score as *ø* ≈ 0 We thus computed the robustness metric for every optimal design along the three Pareto fronts in Fig. 3B. The results, shown in Fig. 4B, suggest an inverse relation between production and robustness to carbor influx; the three circuits display poor robustness (*ø* ≈ 0) in the high production regime of the Pareto front, and high robustness *ø* ≈ 1 in the low production regime. Moreover, we found that upstream repression generally outperforms the other architecture in terms of robustness to growth conditions, and this happens even in the optimal trade-off area of the Pareto front (highlighted in Fig. 4B). The downstream acti vation circuit appears to be particularly fragile, with near zero low production regime. Dual control, on the other hand, displays an intermediate level of robustness across many production levels, but as in the previous section, given its extra complexity and cost, it does not provide obvious benefits over upstream repression.

#### b. Sensitivity to promoter leakiness

Regulated promoters typically express a residual amount of protein in the absence of a regulator. This is commonly known as leakiness and is generally detrimental for circuit design. Since the mechanisms behind leaky expression depend on the specific regulatory mechanism of the promoter, leakiness is poorly understood and makes it difficult to predict system performance.

To quantify the dependencies between leakiness and optimal performance, we re-simulated each circuit in Fig. 3C with a modified model for enzyme expression:
where represents the promoter leakiness as a fraction of the optimized promoter strength and u* _{i}*(

*x*) is the dose-response curve of the biosensor assuming a non leaky promoter. To quantify the robustness of each architecture to the leakiness α, we computed the change in the production loss: whereis the optimal production loss computed for nonleaky promoters, and is the value of the production loss obtained with a leaky promoter. Note that for convenience we have defined Δ

_{1}*J*as the negative change in loss. Intuitively, since leakiness increases the total amount of enzyme expressed, we expect a higher genetic burden accompanied by an increase in production, i.e. a decrease in the production loss

_{1}*J*and Δ

_{1}*J*> 0. On the contrary, designs where leakiness causes a drop in production would lead to Δ

_{1}*J*< 0, an undesirable scenario where the extra genetic burden is accompanied by a lower production flux. This concept is illustrated in Fig. 4C, where perturbations to the Pareto front can shift designs to the left (more production, Δ

_{1}*J*> 0) or to the right (less production, Δ

_{1}*J*< 0).

_{1}We found stark differences in the Δ*J _{1}* scores for each architecture (Fig. 4D). Both upstream repression and dual control display regimes with lowered production flux (Δ

*J*< 0), despite the additional availability of enzymes caused by leaky expression. Moreover, this phenomenon appears in the high production regime in both cases and across all tested levels of leakiness. This suggests that leaky expression generally leads to sub-optimal production in architectures that include an upstream repression component. Upstream repression can prevent an increase in the production flux due to the interplay between an increased enzyme abundance and the stronger repression caused by elevated concentration of the intermediate metabolite. In contrast, we found the downstream architecture displays negligible Δ

_{1}*J*scores across all optimal designs and all tested levels of leakiness; the highest value for down-stream activation in Fig. 4D is ΔJ2 ≈ 2%. This result suggests that downstream activation produces a pathway flux that is extremely robust to promoter leakiness.

_{1}### D. Control circuits for production of glucaric acid

To showcase the power of our approach, we optimized control circuits for the production of glucaric acid in *Escherichia coli*. Glucaric acid is a high-value precursor to several key products such as polymers^{45} and therapeutics^{46}. The seminal work of Moon *et al*^{47} built a pathway that synthesises glucaric acid from glucose-6-phosphate (g6p) in upper glycolysis. As shown in Fig. 5A, the four-step pathway requires the native *E. coli* enzyme inositol monophosphatase (SuhB), plus three heterologous enzymes: inositol-3-phosphate synthase (Ino1) from *Saccharomyces cerevisiae*, *myo*-inositol oxygenase (MIOX) from *Mus musculus*, and uronate dehydrogenase (Udh) from *Pseudomonas syringae*.

In a subsequent work^{11}, a dynamic control circuit was built using the HTH-type transcription factor IpsA, a dual transcriptional regulator from *Corynebacterium glutamicum* that is sequestered by pathway intermediate *myo*-inositol (MI). The circuit was engineered to put MIOX under negative control by IpsA, thus creating a negative feedback system akin to the downstream activation circuit in Fig. 3A, and resulted in a 2.5-fold increase of glucaric acid titer compared to static designs^{11}. We employed our multiobjective optimization framework to compare the architecture implemented by Doong *et al*^{11} against alternative negative feedback architectures, exploiting IpsA as a dual transcriptional regulator to explore the three architectures as in the previous exemplar model: downstream activation of MIOX, upstream repression of Ino1, and dual control of both Ino1 and MIOX (Fig.5B).

We built a kinetic model of the glucaric acid pathway in Fig. 5A with parameters and enzyme kinetics mined from the literature and calibrated to multi-omics data^{48}. To develop a minimal model, we assume that the reactions catalyzed by SuhB and Udh occur spontaneously, so pathway dynamics can be accurately captured by a reduced model for Ino1 and MIOX alone, as shown in Fig. 5B. This assumption is based on the observation that *myo*-inositol-1-phosphate (mi1p) and glucuronic acid (GUA) were not detected in culture products^{47}, indicating that SuhB and Udh are not rate limiting. Moreover, activity of Udh is also reported to be two orders of magnitude faster than than Ino1^{47}, providing further support that Udh is not rate limiting. The model also includes the reported leakage of *myo*-inositol to the media^{47}, as well as the allosteric activation of MIOX by its own substrate MI (dashed black arrow in Fig. 5B). Further details of model equations, assumptions and parameter values are given in the Methods and the Supplementary Information.

We model the dynamics of Ino1 and MIOX in Fig. 5B through the ODEs:
where enzyme expression rates follow a lumped model
and MI is the concentration of myo-inositol sensed by the transcription factor IpsA. As in the previous example, the parameters to be optimized are the promoter strengths (*a _{i}*) and regulatory thresholds (

*θ*), and we fix the Hill coefficients to

_{i}*n*=

_{1}*n*= 1. Since there are no quantitative data on the expression levels of Ino1 and MIOX in E. coli, we define the promoter strengths as so that we can directly optimize the fold-change factors instead of the enzyme expression rates themselves.

_{2}As in the previous example, our aim is to determine control designs that maximise production whilst minimising the bur den of expressing Ino1 and MIOX. To this end, we define two objective functions to minimise:
where *T* is an optimization horizon chosen as 10 times the doubling time log2/λ We formulated a multiobjective opti mization problem to jointly minimise both objectives *J _{1}* and

*J*, as in the previous example. We additionally imposed bounds on promoter strengths as well as upper limits on the biosensor thresholds (

_{2}*θ*). The full statement of the optimization problem can be found in Table IC and all solutions are reported in the Supplementary Data.

_{i}The optimization results in Fig. 6A indicate that both objective functions mutually oppose each other and that the Pareto front is convex in all cases. However, unlike the previous example, we found that downstream activation outperforms the other two architectures: for the same level of genetic burden, the downstream activation circuit produces a higher production flux. This difference can be attributed to the pathway structure, as the glucaric acid pathway includes a number of processes absent from the general branched pathway of the previous section, including allosteric substrate activation of MIOX, export of *myo*-inositol, and a more detailed model of central carbon metabolism. Similarly to the results in Fig. 3C, here we observe that stronger promoter strengths, and so higher MIOX expression, generally improves production for the downstream activation architecture.

The optimal steady state expression levels shown in Fig. 6C suggest that production can be increased with higher Ino1 expression, reflected by the linear increase in optimal Ino1 expression and low expression of MIOX. Increases in Ino1 expression drive a higher accumulation of *myo*-inositol that results in a higher flux through the MIOX reaction. Examination of the MIOX turnover rate *υ*_{MIOX}/MIOX (i.e. the reaction per unit enzyme), shown in Fig. 6D, suggests that biosensors in the optimal trade-off region of the Pareto front tend to maximise MIOX turnover rate. Moreover, comparison of Fig. 6C–D indicates that Ino1 and MIOX levels have opposite effects on the turnover rate. Although increases in Ino1 boost the efficiency of the MIOX reaction, increasing MIOX expression itself has the converse effect and leads to lower efficiency. This is because higher amounts of MIOX reduce the steady state pool of *myo*-inositol, which in turn drives MIOX away from saturation. As a result, in the Pareto front a higher level of MIOX is needed to counteract the low turnover rate, causing the pronounced increase in burden and marginal gains in production observed at the top-left of the Pareto front in Fig. 6A. This observations are also true of the other feedback circuits (Supplementary Figure 2, Supplementary Figure 3). Altogether, these observations reiterate the complex dependencies and between system performance and design parameters. These results provide strong computational evidence that the circuit architecture implemented by Doong *et al*^{11} has substantial benefits over other implementations of negative feedback control.

## III. DISCUSSION

Dynamic control has emerged as a promising solution to some of the caveats in metabolic engineering, including the optimization of enzyme expression timing and strength, the avoidance of metabolic intermediate accumulation, and the sub-optimal performance in large fermentations. Dynamic control systems, however, are inherently more complex than static designs, as they have more degrees of freedom and their effectiveness often varies with pathway structures. Previous dynamic control systems were mostly built based on intuitive understanding of pathway features, but in general it is extremely challenging to identify suitable control architectures and parameters via purely experimental trial-and-error approaches. Computational design strategies can dramatically improve the efficiency of dynamic control systems, making them more widely and easily adoptable. However, such systematic computational tools are strongly lacking, limiting the broad use of dynamic control systems. This is in stark contrast with traditional pathway engineering, where designers can avail of powerful computational methods based on genomescale models^{49,50} and a large collection of software tools for analysis and design^{51}.

Here we have proposed the use of multiobjective optimization as a design method for dynamic pathway control. We have shown that for a given pathway model and candidate control architectures, various design objectives and toxicity constraints can be compiled into an optimization problem. The resulting problem can be solved with global optimization algorithms available in many scientific computing platforms. The optimal solutions represent different dose-response curves for a metabolite biosensor, and thus the method produces libraries of biosensors that optimally navigate the design trade-offs specified by the objective functions.

A powerful application of our method is the identification of optimal architectures among candidate circuits. In both examples considered, the general branched pathway and the glucaric acid pathway, we compared the Pareto fronts of each architecture and singled out circuits that outperform others in specific operating regimes. Moreover, our sensitivity analysis of the Pareto fronts revealed a number of previously unknown features of dynamic control circuits. For example, we observed that upstream repression is beneficial in terms of performance and robustness to growth conditions, yet fragile to promoter leakiness. Conversely, downstream activation has a poorer performance but is extremely robust to promoter leakiness. We moreover showed that downstream activation outperforms other architectures in the production pathway for glucaric acid built by Doong et al^{11}. Taken together, these results suggest that circuit design is highly specific to the pathway of interest and the metrics employed to quantify performance and robustness. Our analysis thus highlights the need for model-based approaches, as such complex dependencies between control circuits and performance are difficult to detect from intuition alone

Although the use of optimization theory in metabolic dynamics dates almost two decades, the focus has mainly confined to natural systems^{52}. For example, the seminal work by Zaslaver et al53 showed that gene regulation of amino acid biosynthesis in E. coli could be understood as the solution of a cost-benefit optimization problem. Subsequent theoretical54 and experimental work^{55,56} validated some of those early observations in other model systems, demonstrating that optimization can be a powerful tool to reverse-engineer control strategies found in nature. In the case of engineered pathways, however, the literature has largely focused on ad hoc simulation-based methods, and to date there are no general design frameworks for dynamic control. In this paper we show how multiobjective optimization can bridge this gap and provide new insights on the design of control systems for engineered pathways

Our approach offers two key advantages: the incorporation of constraints, and the resolution of design trade-offs. optimization-based design allows the inclusion of promoter or metabolite toxicity constraints directly into the design process. For pathway design, such constraints can account for e.g. upper bounds on enzyme or metabolite concentrations^{34}, biophysical bounds on the promoter parameters^{57}, or thermodynamic limits for metabolic fluxes58. Our focus on a milti objective framework moreover reveals that dynamic control can effectively be employed to find optimal trade-offs between conflicting performance objectives. We have deliberately focused on genetic burden (cost) and production flux (benefit) as objective functions, on the basis that they are general and relevant across a wide range of design scenarios. The framework can be readily extended to other objective functions relevant in applications, for example by including differential toxicity of specific enzymes, the optimization of metabolite concentrations over time, or the control of cell-to cell heterogeneity^{59,60}

For simplicity here we have focused on biosensors implemented via metabolite-response TFs, because they are the most widely used mechanism for metabolite sensing. A number of works have explored alternative sensing mechanisms, e.g. utilising nucleic acid aptamers^{61,62} or riboswitches^{63,64}. In general our methodology is applicable to such implementations and to other sensing mechanisms that can be described by sigmoidal dose-response curves. Other strategies, such as the use of quorum sensing to sense culture density as a proxy for growth rate^{65,66}, will require bespoke modelling frameworks once their implementation becomes more standardised7.

Dynamic pathway engineering requires the specification of control architectures as well as the parameters of the circuit components. A key step for the widespread adoption of this technology is the development of design methods that can integrate pathway dynamics with implementation constraints and performance objectives. Here we have proposed such a strategy based on multiobjective optimization that shows promise for its deployment across a range of metabolic engineering problems.

## IV. METHODS

### A. Computational models

#### a. Branched pathway

The mass balance equations for the exemplar pathway (Fig. 2A) in continuous culture are:
where *x*_{0} and *x _{1}* are the metabolite concentrations, and λ is the culture growth rate fixed at λ = 1.93 10

^{−4}s

^{−1}(corresponding to a fast doubling rate of 26min in

*E. coli*). The reaction rates are assumed to follow Michealis-Menten kinetics: where

*e*are the enzyme concentrations and (

_{i}*k*

_{cat}

*,*

_{i}*K*

_{M}

*) are kinetics parameters of each enzyme; we used representative values for the kinetic parameters*

_{i}^{67}and for simplicity assumed equal kinetics for the three enzymes. The model further assumes constitutive expression of the native enzyme

*e*

_{0}, and regulated expression of the heterologous enzymes

*e*

_{1}and

*e*

_{2}, as described in Eq. (5). All kinetic parameters can be found in the Supplementary Data.

#### b. Glucaric acid pathway

The mass balance equations for the GA pathway (Fig. 5A) in continuous culture are:

We assume a constant uptake of glucose *υ*_{PTS} and growth rate λ, with values *υ*_{PTS} = 1.34 mmol/gDCW/h = 0.1656 mM/s and A = 0.1 *h*^{-1} = 0.278 10^{−4}s^{−1} taken from the Keio multi-omics data set^{48}. Details on the model construction, assumptions, reaction kinetics and parameters can be found in the Supplementary Information and Supplementary Data.

### B. Model simulation and optimization

In all cases we first simulated pathway dynamics in absence of heterologous enzymes, and use the resulting steady state metabolite concentrations as initial conditions for the dynamic control simulations. This strategy emulates common lab protocols in which strains are first pre-cultured in absence of heterologous expression, and then inoculated into fresh media containing inducers of enzyme expression^{11}. The initial conditions for each case study can be found in the Supplementary Data.

The full statement of the optimization problems solved in Fig. 2, Fig. 3, and Fig. 6 can be found in Table I. All numerical calculations were done in Matlab 2019b. The model for the branched pathway was integrated with the ode15s stiff solver, and the glucaric acid model was integrated with the Matlab implementation of the SUNDIALS CVode solver for precision and speed^{68}. In all cases, parameter optimization was done with the multiobjective optimizer gamultiobj with parameter bounds and path constraints to account for cytotoxic effects of the intermediate. In all analyses, we computed a total of 70 Pareto points. For ease of visualisation, in Fig. 2, Fig. 3, and Fig. 6 we plotted every other solution, i.e. 35 out of 70 solutions. All optimized dose-response parameters are reported in the Supplementary Data.

## AUTHOR CONTRIBUTIONS

BKV constructed and analyzed the model of the branched pathway; AAM constructed and analyzed the model of the glucaric acid pathway. DAO designed and supervised the research; all authors contributed to analysis and discussion of results.

## ACKNOWLEDGEMENTS

BKV was funded by a Wellcome Trust ISSF3 fund awarded to DAO. AAM acknowledges funding from BBSRC Research Grant BB/M017982/1.