Bayesian Optimization for Design of Multiscale Biological Circuits

Recent advances in synthetic biology have enabled the construction of molecular circuits that operate across multiple scales of cellular organization, such as gene regulation, signaling pathways, and cellular metabolism. Computational optimization can effectively aid the design process, but current methods are generally unsuited for systems with multiple temporal or concentration scales, as these are slow to simulate due to their numerical stiffness. Here, we present a machine learning method for the efficient optimization of biological circuits across scales. The method relies on Bayesian optimization, a technique commonly used to fine-tune deep neural networks, to learn the shape of a performance landscape and iteratively navigate the design space toward an optimal circuit. This strategy allows the joint optimization of both circuit architecture and parameters, and provides a feasible approach to solve a highly nonconvex optimization problem in a mixed-integer input space. We illustrate the applicability of the method on several gene circuits for controlling biosynthetic pathways with strong nonlinearities, multiple interacting scales, and using various performance objectives. The method efficiently handles large multiscale problems and enables parametric sweeps to assess circuit robustness to perturbations, serving as an efficient in silico screening method prior to experimental implementation.


S2 Toy Pathway
The first pathway under consideration is a toy pathway with two pathway intermediates and enzymes under transcriptional control (1 ). The substrate represses a metabolite-responsive transcription factor (MRTF) which can act as an activator or repressor for the transcription of enzymes 1 and 2. (see Figure 1A). The influx from native precursors and draw from pathway metabolism are modeled as constant rates. While there is no direct biological correlate for this pathway, it exemplifies the most basic engineered pathway possible and thus is useful primarily as a proof-of-concept for the BayesOpt method.
The mass balance equations for the toy branched pathway are to parameterize using limited kinetic data. Additionally, dynamic pathways often express the TF constitutively, which reduces the effect of TF expression on performance (2 ). We consider three architecture forms: activation, repression, and open loop control: As a result, the toy pathway has four continuous parameters under study: two associated with each enzyme. The four optimizable parameters are constrained on the following ranges: The possible architectures are limited to three architectures which include only negative feedback loops, in addition to an open-loop control with no dynamic feedback loops (see Figure 1A). This architecture limitation was made to remove the possibility of undesirable multistability or oscillatory behavior. All parameter values are given in Table S1. Supplementary The ↵ 1 and ↵ 2 values given in Table S1 were chosen to scale the two objectives to be similar in magnitude and the overall loss to be between 0 and 1. The loss equation for the objective function is We ran each simulation of the model for a total time of 5 · 10 4 seconds to ensure the final values would be at steady state. All metabolites and enzymes had an initial concentration of 0µM except for x 0 , which had an initial concentration of 2290µM (1 ). A summary of the model details is given in Table S2. Supplementary

S3 Glucaric Acid Production
The second pathway describes the synthesis of glucaric acid, a key precursor for a number of applications (3 ). This pathway has been implemented using dynamic control architectures which led to a 2.5-fold increase in product titer over static metabolic engineering methods (4 ). Glucose is converted to glucose-6-phosphate, which can either be drawn into native metabolism and converted, reversibly, into fructose-6-phosphate, or into the engineered pathway (see Figure 2A). The pathway requires three foreign enzymes: inositoal-3-phosphate synthetase, or Ino1, from Saccharomyces cerevisiae, myoinositol oxidase (MIOX), from Mus musculus, and uronate dehydrogenase (Udh) from Pseudomonas syringae. In addition to being exported from the cell and acting allosterically on MIOX, myoinositol can sequester the transcription factor IpsA, a dual transcriptional regulator from Corynebacterium glutamicum. IpsA can then act on the transcription of Ino1 or MIOX. The glucaric acid pathway provides a more complex, real-world application that builds on the toy model and includes reversible and allosteric reactions.
The glucaric acid model was chosen for its increased complexity. While maintaining the same number of possible architectures as the toy model, this model incorporates allosteric control and reversible reactions, two more complex biological interactions. We specify the model equations and describe a sample optimization. We then experimented with perturbing the growth conditions and kinetic parameters to investigate model robustness.
The glucaric acid pathway is based on the one designed in (1 ). The mass balance equations for the glucaric acid production pathway (see Figure 2A) are where Ino1 and MIOX are the enzymes in the pathway and g6p, f6p, and MI are the substrates. The enzyme parameters k cat , k eq , k m, x and k m, y are all fixed kinetic parameters specific to each enzyme. The effective substrate activation constant k cat, eff has two additional activation kinetic parameters k a and a which must be specified: The function u(x, k, ✓) describes the genetic control topology at the enzyme's promoter.
There are three options for this functional form: activation, repression, and no control as in the toy model (see Equation S2).
Similarly to the toy pathway, there are four continuous parameters: The four decision variables are constrained on the following ranges: The possible architectures are limited to those only containing negative feedback loops. Similarly to the toy model, we term these three architectures upstream repression, downstream activation, and dual control based on their mode of action (see Figure 2A). All kinetic parameters are given in Table S3. Supplementary The ↵ 1 and ↵ 2 values in the objective function were chosen to scale the two objectives to be similar in magnitude and the overall loss to be between 0 and 1 (see Table S3). The loss equation for the objective function is (S11) A summary of the model details is given in Table S4. We ran each simulation of the model where FFA is the concentration of free fatty acids and tesA is the concentration of thioesterase enzyme (tesA). The cellular growth rate is set for all architectures at 3.85 · 10 4 mM/s, which is calculated from an E. coli doubling time of 26 minutes. The parameter k tesA is the kinetic parameter of the thioesterase reaction. There is one free optimisable parameter in this model, the tesA promoter constant r lac .
The negative gene loop architecture expresses a transcriptional repressor on the same promoter as tesA, creating a negative feedback loop on enzyme production. The mass balance equations for this architecture are (S13) The variable tetR is the concentration of the repressor expressed on the same gene as tesA.
While the other models allow multiple modes of transcriptional control (activation, repression, no control), this model only allows for activation due to the type of transcription factor used (5 ). The relative expression strengths r tl and r tl, tetR are free parameters.
The negative metabolic loop architecture includes a transcription factor which acts on the tesA promoter. The mass balance equations for this architecture are where the free parameters of this model are k i and r fl 0 , which control the shape of the transcription factor dose-response curve. Unlike the other free model parameters, k i varies on the range 0 to 0.12.
Finally, the layered negative feedback loop uses tetR as an intermediate repressor. The transcription factor is repressed by the product and in turn represses the repressor TetR, which can repress the expression of the enzyme tesA. This negative feedback loop is described by the following equations: (S15) The repressor TetR is produced on the AR2 promoter, which has its own strength parameter r ar2 . The free parameters on this model are r ar2 and r tl .
Unless otherwise stated, all continuous model parameters are restricted to the range from 10 · 10 11 to 10 · 10 8 . We ran each simulation of the model to 5 · 10 4 s to ensure steady state values had been reached. All pathway components started with a concentration of 0mM. Two objective functions were implemented for this pathway: a production-burden function and a speed-accuracy function. The loss equation for the production-burden objective function is the sum of the production loss and pathway cost: where J prod is the production loss. The ↵ 1 and ↵ 2 values in the objective function were chosen to scale the two objectives to be similar in magnitude and the overall loss to be between 0 and 1 (see Table S3). The production loss is defined as The pathway cost J cost varies for each architecture: (S18) The second speed-accuracy objective function has two terms, percent overshoot and normalized rise time: The percent overshoot J os is the percent difference between the maximum FFA concentration

S5 P-Aminostyrene Synthesis
The final pathway we consider is the synthesis of p-aminostyrene (p-AS) in E. coli Here, the constant L is the rate of loss of p-AF through the leaky cellular membrane, and V chorismate is a constant parameter defined in Table S8. The constant ⌧ is the toxicity factor which scales concentrations to account for intermediate metabolite toxicity. The constant is the dilution rate due to cellular growth. The function f is a variation on the Michaelis- Menten equation: where: e 2 {papA, papB, papC, deaminase, LAAO, P efflux }, The promoter and mRNA transcript concentrations are modeled separately here, unlike in the previous models. The mass-balance equations for these elements are where the function u(x, k, ✓) is the sigmoidal control function defined in Equation S2. The binary parameters w i define which metabolite exerts control on each promoter: 1 papABC controlled by p-AF, 0 papABC controlled by p-ACA, 1 LAAO controlled by p-AF, 0 LAAO controlled by p-ACA, 1 P efflux controlled by p-AF, 0 P efflux controlled by p-ACS. (S23) To avoid circuits with multistable dynamics, we restrict the search to architectures without positive feedback loops. This means that we consider a total of 3 3 = 27, i.e. those where: • promoter papABC is either unregulated or repressed by one of the two intermediates (p-AF or p-ACA), • promoter LAAO is unregulated, activated by p-AF, or repressed by p-ACA, • promoter P efflux is either unregulated or activated by one of the two intermediates (p-AF or p-ACA).
The constant is the DNA duplication rate. The constant µ is the mRNA degradation rate constant, which is assumed to be constant for all mRNAs. PapA, PapB, and PapC are all expressed from the same promoter but their translation rates are variable, so their mRNAs are modeled separately. Finally, the protein folding process is modeled explicitly, with each of the 4 enzymes and the efflux pump having a folded and unfolded state. Additionally, the enzyme deaminase is expressed constitutively but its concentration rises to steady state and thus its dynamics are also modeled. The enzyme mass balance equations are Here, the constant ⇢ is the protein degradation rate and ⌘ is the protein folding rate. The function v is the translation rate equation: where T init is the transcription initiation rate, L m is the length of the mRNA, and R is the ribosome elongation rate. The objective function for the p-AS model is where pathway cost is the sum of all heterologous enzyme transcription rates across all loci and ligands, and the rate f (P efflux , pACA) is defined in Equation S21. The constant V chorismate as well as any other fixed parameter values in the model are defined in Table S8. Table S7 summarizes the model details. Supplementary

S6 Hyperparameter tuning
The tree-structured Parzen estimator method (TPE) implemented in the Hyperopt package has a hyperparameter which controls the balance of exploration and exploitation. TPE splits all samples into two distributions, a "good" and a "bad" distribution and selects the next point from the Expected Improvement computed over only the "good" samples. The parameter is the fraction of points incorporated into the "good" distribution. The default value of in the Hyperopt package is = 15, which corresponds to 15% of samples being partitioned into the "good" sample distribution and 85% in the "bad" distribution. Future samples are then drawn at the maximum of the "good" distribution. Higher values of correspond to increased exploration of the loss landscape, whereas lower values of mean increased exploitation towards local minima.
We found that regardless of the value chosen, the distribution of architectures explored was similar (see Figure S1). To give a single scalar metric, we computed the fraction of samples in each quarter taken from the majority architecture, or the architecture most commonly sampled across the optimisation. This fraction was computed as majority architecture samples optimization iterations , for each quarter of the optimization trace. In the case of the glucaric acid model, the majority architecture is dual control. The majority architecture fraction rose with each quarter as the BayesOpt routine increasingly focused on low-loss dual control parameter value combinations; however, there was no statistical difference between the different gamma values chosen. These results show that the default value is sufficient as TPE is relatively insensitive to hyperparameter tuning. This feature of TPE makes it quick to implement without significant application-specific tuning but may cause it to fail to adapt to certain specific cases. However, we did not observe TPE fail to converge in any of the systems under study.
Supplementary Figure S1: Hyperparameter Tuning The TPE hyperparameter was tuned by cloning the Hyperopt package repository and manually changing the value in the source code. Hyperparameter values of = 10, = 25, = 50, and = 90 were considered and a single optimisation of the glucaric acid pathway was run for each value. The total optimisations were split into four quarters by iteration and the percentage of sample drawn from the majority architecture (in this case, dual control).