Multicriteria global optimization for biocircuit design

Abstract

Background

One of the challenges in Synthetic Biology is to design circuits with increasing levels of complexity. While circuits in Biology are complex and subject to natural tradeoffs, most synthetic circuits are simple in terms of the number of regulatory regions, and have been designed to meet a single design criterion.

Results

In this contribution we introduce a multiobjective formulation for the design of biocircuits. We set up the basis for an advanced optimization tool for the modular and systematic design of biocircuits capable of handling high levels of complexity and multiple design criteria. Our methodology combines the efficiency of global Mixed Integer Nonlinear Programming solvers with multiobjective optimization techniques. Through a number of examples we show the capability of the method to generate non intuitive designs with a desired functionality setting up a priori the desired level of complexity.

Conclusions

The methodology presented here can be used for biocircuit design and also to explore and identify different design principles for synthetic gene circuits. The presence of more than one competing objective provides a realistic design setting where every solution represents an optimal trade-off between different criteria.

Background

A hallmark of Synthetic Biology is, quoting Arkin, the ambition to formalize the process of designing cellular systems in the way that traditional engineering disciplines have formalized design and manufacture, so that complex behaviours can be achieved for practical ends[1]. In formalizing the design process, as it is the case in more traditional engineering disciplines, mathematical modeling and optimization play a central role.

Over the past ten years, many advances have been achieved in the field, from the first bacterial toggle switches [2] and biological oscillators [3], to the recent mammalian cell to cell communication devices [4]. In a so called first wave of Synthetic Biology basic elements and small biological modules were successfully implemented and characterized. One of the challenges of the second wave in progress is the integration of modules to create circuits of increasing complexity [5]. However, as reported by Purnick and Weiss [5], the level of complexity achieved in synthetic circuits, measured by the number of regulatory regions, is relatively low. While circuits in Biology are complex, subject to natural tradeoffs and playing multiple roles [6], most synthetic designs are simple and perform a single task. Reported reasons for the current limited complexity in synthetic circuits include too simplistic engineering design principles [5], difficulty to independently control multiple cellular processes in parallel [7] and increasing problems to construct and test circuits as they get larger [8]. Efforts are necessary to overcome these difficulties and, quoting Lu et al.[9], advancing Synthetic Biology to the realm of higher-order networks with programmable functionality and real world applicability. In parallel, new computational tools need to be developed to support these efforts [10].

In this contribution, our goal is to set up the basis of an advanced optimization tool for the modular and systematic design of biocircuits capable of handling high levels of complexity and multiple design criteria.

Modular design requires the previous definition of standardized functional objects and interfaces [11]. From the foundations of Synthetic Biology, efforts have been held in order to characterize standard biological parts, i.e. DNA sequences encoding a function that can be assembled with other standard parts. The abstraction hierarchy proposed by Endy [12] classifies standard parts in three different layers: parts, which are defined as sequences with basic biological functions (like for example DNA-binding proteins), devices which are combinations of parts with a particular function and systems which are combinations of devices. An emerging catalogue of standard parts is available at the registry supported by the BioBricks Foundation [13].

Systematic design relies on mathematical models describing the circuit dynamics. In this regard, modular modeling tools are advancing to facilitate the mathematical representation of biological parts and their combinations [14], providing the description of the reactions taking place inside the different parts and the interfaces to connect them. Inspired by the BioBrick registry of standard parts, Marchisio and Stelling [15] developed a formal modeling framework based on the ordinary differential equations (ODE) formalism which permits modular model composition and has been recently extended for the modeling of more complex eukaryotic systems [14]. Some remarkable advances have been also achieved regarding synthetic biology computer aided design tools [16]. The systematic design of circuits combining components or parts from a list or library can be formulated as an optimization problem [16]-[18] where the circuit model structure is manipulable through decision variables, and the desired behaviour of the circuit is encoded in the objective function to optimize. This results in Mixed Integer Nonlinear Problems (MINLP) whose solution is challenging due to the simultaneous presence of binary variables and constraints in form of ODEs. Dasika and Maranas [17] developed an optimization framework for the design of biocircuits, based on the circuit modeling formulation by Hasty [19] and a multistart local outer approximation method for the optimization. A number of design problems were successfully solved within this framework including a circuit with inducer specific response, a genetic decoder and a concentration band detector.

In this work, we advance the optimization-based design of biocircuits with two contributions: increasing the computation efficiency in order to handle higher levels of complexity and introducing multiple criteria in the design. To this purpose, we first introduce a set of global MINLP solvers that reduce drastically the computation time for the monoobjective design problem in comparison with other published methods. Then we formulate a general multiobjective optimization framework that combines the efficiency of the global MINLP solvers with the ability to tackle multiple design criteria. The inducer specific response circuit design by Dasika and Maranas [17] is used to illustrate the efficiency of the MINLP methods presented and further reformulated with additional design criteria to discuss the advantages of a multiobjective formulation in the design of genetic circuits.

Methods

Global stochastic MINLP solvers for biocircuit design

Optimization based design of biocircuits requires the integration of tools for modular modeling, simulation and optimization. As reported in the Background section, modular tools for modeling in Synthetic Biology are advancing fast as well as repositories of biological parts. Searching for a generic optimization framework, the methods presented next do not bound to a specific modeling tool, but accommodate to any ODE based modeling framework such that the circuit's model structure can be obtained from the starting list of parts by giving values to a set of integer variables.

The design problem consists of finding the best solution or solutions among the set of all possible alternatives according to a number of criteria. In this first part, we focus on problems with one unique design objective. Under these assumptions, the design of biocircuits can be formulated as a Mixed Integer Nonlinear Programming Problem [17],[18], where the model structure can be encoded by integer variables and the constraints are the dynamics of the system in form of ODEs. Tunable kinetic parameters are real decision variables in the optimization model. For a complete formulation we refer to [17], where the single objective MINLP problem is formalized for a particular modeling framework [19]. Next, our focus is on the computational challenges of the resultant MINLP, since some features inherent to biological circuit models make it particularly difficult to solve.

In first instance, the dynamics of biocircuits are highly nonlinear, and the resultant optimization problem is non convex and multi-modal. In this type of problems, local methods lead to suboptimal solutions (unless we start close to the global optimum). A number of approaches have been proposed in previous works to find the global optimum in monoobjective biocircuit design. Dasika and Maranas [17] implemented a multistart local outer approximation algorithm where a convergence sequence of upper and lower bounds to the original problem is generated and a local optimum solution is identified at each iteration. In this way, a local deterministic search is performed from several points. Rodrigo et al.[18] use a stochastic metaheuristic based on simulated annealing [20],[21] and Huynh et al.[22] apply a global deterministic optimization method to a linear approximation of the model around a steady state.

On the other hand, the design of gene circuits involves in general large search spaces that combine a high number of integer variables with the presence of real variables. Our first goal is to provide global optimization methods that efficiently solve monoobjective design problems of medium/high complexity. Global deterministic methods ensure convergence to the global optimum within a desired tolerance, but the computational burden is in general very high for non convex systems with large search spaces. Therefore, we have decided to employ global stochastic methods, which offer no guarantee of convergence to the global minimum in a finite number of iterations but showed excellent results solving complex process optimization problems in reasonable computation time [23].

In this work, we use three different global stochastic methods: mixed-integer tabu search (MITSMITS) [24], mixed-integer ant colony optimization (ACOmiACOmi) [25] and the enhanced scatter search eSSeSS described in [23]. The three methods are actually hybrid, since the stochastic global search is combined with the local mixed-integer sequential quadratic programming (MISQPMISQP) developed by [26]. These methods have been shown to be efficient metaheuristics in solving complex-process optimization problems from different fields, providing a good compromise between diversification (exploration by global search) and intensification (local search).

MITSMITS uses a combinatorial component, based on Tabu Search [27], to guide the search into promising areas, where the local solver is activated to precisely approximate local minima. Exler et al.[24] made use of MITSMITS to solve complex integrated design problems where other state of the art solvers failed, including a wastewater plant for nitrogen removal and the well known Tennessee Eastman Process. ACOmiACOmi extends ant colony optimization meta-heuristic [28] to handle mixed integer search domains. Schlueter et al.[25] showed the efficiency of this method for a number of engineering benchmark problems with high levels of non-convexity. Finally, eSSeSS is an enhanced version of the scatter search for mixed integer search domain. Egea et al.[23] proved the efficiency of the method for solving complex-process models through a set of engineering benchmarks, where eSSeSS performed well even in cases in which standard local search methods failed to locate the global solution.

In this contribution, we evaluate the efficiency of these methods in the context of Synthetic Biology and in particular for the systematic design of genetic circuits. For illustrative purposes we chose a representative design example from Ref. [17], with one single design objective.

Starting from a list of components, the goal is to build a circuit with a specific response upon stimulation by two different inducers. There are eight different promoter elements (denoted by P ₁…P ₈): P l a c ₁, P l a c ₂, P l a c ₃, P l a c ₄, P _λ , P t e t ₁, P t e t ₂ Para and four transcripts (denoted by R ₁…R ₄): tetR, lacI, cI, and araC. The inducers of interest are IPTG and aTc. The dynamic model of the overall reaction network is constituted by a set of ordinary differential equations of the form:

(1)

where E _j is the expression term for the transcripts, K _{j decay} z _j is the degradation rate and V _j is the production/consumption rate of z _j due to other reactions. The expression rates for the transcripts are known and they read:

(2)

where v _ji is the rate of production of R _j from P _i , and Y _ij is a binary variable such that:

The structure of the model is given by a 8×4 superstructure matrix Y containing the 32 binary variables of the model. We define the vector of binary variables y as the vector obtained by converting the matrix Y to a vector by columns. The tunable parameters are contained in a vector of real variables denoted by x.

As mentioned, the goal is to achieve a specific response upon induction. Namely, the steady state level of LacI must be high upon aTc and low upon IPTG induction whereas the steady state level of tetR must be low upon aTc and high upon IPTG induction. This design goal is encoded in the following objective function Z to be maximized:

(3)

where the maximum value Z=1 is achieved for .

The design problem is formulated as a MINLP where the decision variables are contained in the vectors y and x and the objective function to maximize is Z in (3), subject to the system's dynamics (1). The following constraint on the maximum number of active pairs (M _max ) is also imposed:

(4)

thus limiting the complexity of the circuit.

First we use the original formulation of the problem by Dasika and Maranas [17] with a maximum of two promoter-transcript pairs, and compare the performance of the methods with the published results. Afterwards, we gradually increase the network complexity to evaluate how the methods proposed scale with the increasing problem size. The results obtained are included in the Results and discussion section.

Multiobjective framework for automatic biocircuit design

In traditional engineering disciplines design problems are often multicriteria, where a number of design objectives are conflicting (typically production and cost) since we cannot increase one without decreasing the other. Problems with multiple and conflicting design criteria do not have a unique optimal solution, but a trade-off front between the competing objectives, also known as Pareto optimal front of solutions.

In biological systems, trade-offs between robustness, fragility, performance, and resource demands have been conjectured [6],[29]-[32]. We know that living organisms allocate limited resources to various competing traits, and arising tradeoffs are central to evolutionary biology. Furthermore molecular pathways have been shown in many cases to play diverse and complex roles. However, de novo engineered circuits have been designed to perform a single task, and optimization based designs in Synthetic Biology have been formulated as problems with a single objective.

In this contribution we propose a multiobjective optimization framework for the design of biocircuits. In first instance, the design is formulated as a multicriteria optimization problem with a number of conflicting objectives and then a multiobjective optimization strategy is implemented to find the Pareto optimal set of solutions.

In order to mathematically define the multiobjective design problem, let us first introduce the following vectors: is the vector of state variables coding for the levels of all the species involved in the circuit (we will denote its time derivative by ); is the vector of continuous variables containing a set of tunable parameters; is the vector of integer variables determining the circuit model structure; is the vector of fixed parameters and for i=1,…,s is the set of conflicting objectives, where one subset of objectives encodes aspects related to the performance of the circuit and a second subset encodes aspects related to robustness and/or cost.

The design of a biocircuit can be formulated as finding a vector of continuous variables and a vector of integer variables which minimize the vector J of s objective functions:

(5a)

subject to:

the circuit dynamics in the form of ODEs or differential algebraic equations (DAEs) with the state variables z and additional parameters k:

(5b)

additional requirements in the form of equality and inequality constraints:

(5c)

(5d)

upper and lower bounds for the real and integer decision variables:

(5e)

(5f)

In order to evaluate the solutions of the multiobjective optimization problem, we need to introduce the notion of Pareto optimality [33]. Given two pairs (x ^*,y ^*), (x ^**,y ^**), we say that the vector J(x ^*,y ^*) dominates J(x ^**,y ^**) if J(x ^*,y ^*)…J(x ^**,y ^**) for all i=1,…,s with at least one strict inequality. A feasible circuit defined by (x ^*,y ^*) is a Pareto optimal solution of the multiobjective optimization problem if it is not dominated by other feasible circuits. The set of all Pareto optimal solutions is known as Pareto front.

Computing the Pareto optimal set is a very challenging task in the context of complex biocircuit design. On the one hand, as indicated previously, high complexity imply large search spaces, and on the other hand the expected Pareto front is discrete and possibly non-convex, due to the high nonlinearity of the biocircuits dynamics and the existence of discrete decision variables.

There are a number of approaches to solve multiobjective optimization problems (MOPs) [34]. Evolutionary approaches [35] allow to compute an approximation of the entire Pareto front in one single run, but require large population sizes and consequently a high computational effort for the systems with the complexity we want to tackle. Scalar approaches consist in transforming the MOP into one or more single objective problems, and include among others the well known weighted sum approach, Normal Boundary Intersection (NBI) and ε-constraint methods [33].

In the weighted sum approach, weights must be changed in order to generate different solutions in the Pareto front and the performance depends on the choice of the weighting coefficients, which is in general not straightforward. This method cannot find solutions in concave parts of the Pareto front.

NBI first builds a plane in the objective space which contains all convex combinations of the individual minima, denoted as convex hull of individual minima (CHIM) and then constructs normal lines to this plane. The MOP is reformulated as to maximize the distance from a point on the CHIM along the normal through this point. When dealing with integer variables, there may not exist a feasible solution on the selected normal to the CHIM, and therefore NBI at least in its original formulation has limited applicability for discrete Pareto fronts.

In the ε-constraint strategy [33], the MOP is reduced to a number of MINLP, where each MINLP is obtained by minimising one of the objectives and converting the rest of criteria to inequality constraints. Different solutions can be obtained by changing the upper bounds on the objectives not minimised. This methodology has two important advantages for the design of complex biocircuits: the methodology works well for discrete and non-convex Pareto fronts and, in addition, it allows exploiting the MINLP solvers introduced in the previous section, which solve efficiently the resultant MINLPs at a reasonable cost.

Next, we describe the ε-constraint strategy implemented in this work. The proposed optimization process is composed of the following steps (for simplicity and without loss of generality we have considered two objective functions J ₁ and J ₂):

Search for the optima of each of the individual objectives, i.e.:

Compute the individual objective bounds as:

where low and upper bars denote lower and upper bounds respectively.

Select the objective function to be minimized, denoted in what follows as the primary objective (without loss of generality let us take J ₁ as the primary objective).

For the non-minimized objective J ₂, generate a vector

such that , and ε ₁<ε ₂<…<ε _m .

Solve the MINLP:

subject to:

for k=1,…,m-1 by means of one of the MINLP solvers introduced in the previous section.

Evaluate the solutions obtained and construct the Pareto front with the non dominated optimal ones.

Continuing with the example introduced in the previous section where the goal was to find a circuit with a specific response upon induction, we introduce now an additional design criterion. As mentioned, in the original formulation the design objective was unique and given by Eq. (3). Here we consider the protein production cost as an additional objective to minimize, competing with the circuit performance. This criterion has been suggested as a design principle by several authors [6],[36]. The cost of protein production C is encoded in an objective function that, considering the mass balance equations (1) takes the form:

(6)

where T is the final time.

We apply the constraint strategy combined with the MINLP solvers to obtain the Pareto front for different degrees of circuit complexity. First, we set the maximum number of pairs to M _max =2, and then we increase the maximum number of pairs to evaluate how the Pareto boundary evolves, and how the methodology proposed scales with the systems size. The results obtained are included in Results and discussion section.

One interesting application of the methodology presented is to explore new topologies of medium or high order that perform a desired (complex) functionality. To illustrate this we make use of the same library of components of the previous example, but in this case searching for a circuit topology with the capability to perform adaptation, setting a priori the desired level of complexity.

Adaptation is defined as the ability of the circuit to reset itself after responding to a stimulus [37]. Here, we evaluate the levels of LacI (output) in response to a sustained stimulus of aTc (input). Ma et al.[37] assessed the ability of a circuit to adapt after a given stimulus by measuring two functional quantities encoded in two competing objectives related to the sensitivity and the precision of the system's response. On the one hand, in order to maximize adaptation after a given stimulus we need to maximize the circuit's sensitivity:

(7)

where O _peak is the level of the output (in this case LacI) at its maximum upon induction and O _t=0 is the level of the output at the steady state before induction. On the other hand, in order to maximize adaptation we need to maximize the circuit's precision, i.e. we need to minimize the following function:

(8)

where O _t=T is the level of the output at steady state reached upon induction. The search for an adaptive circuit can be formulated then as a multiobjective optimization problem where the constraints are imposed by the circuit's dynamics. In this way, it is possible to elucidate whether is it possible to construct a circuit with capacity for adaptation from the available set of components. The maximum and minimum number of allowed components can be adjusted by means of inequality constraints. The details and results of the corresponding multiobjective optimization problem are included in Results and discussionsection.

Results and discussion

Single objective global optimization design of a circuit with inducer specific response

In this subsection we present the results obtained for the monoobjective problem described in Methods section. Starting from the indicated library of transcripts and promoters (the corresponding generic circuit equations are included in the Additional file 1) we search for the circuit with best performance by maximizing Z defined in Eq. (3). We use the MINLP solvers MITSMITS, ACOmiACOmi and eSSeSS with the goal of minimizing J=-Z. We solve the optimization problem for increasing levels of complexity, i.e. for an increasing upper bound in the maximum number of pairs (M _max ). Note that, for a library with p different pairs, the number of possible circuits containing exactly M active pairs is:

(9)

According to this formula, the number of combinations nC increases with p and also with the maximum number of pairs M _max as illustrated in the Additional file 2: Figure S1 and S2. In what follows we do not modify the original library of transcripts and promoters (p=32) and evaluate the performance of the methods for different values of M _max .

For M _max =2, the three MINLP solvers, MITSMITS, ACOmiACOmi and eSSeSS, reached the same solution, the circuit with active pairs (P l a c ₁,t e t R) and (P t e t ₂,L a c I). In Figure 1, we illustrate the best circuit found together with the corresponding superstructure matrix, coefficients of the model and active pairs. The value of the objective function for the optimal circuit is J=-0.99998. This solution coincides with the one obtained by Dasika and Maranas using the outer approximation method [17].

Figure 1

Optimum of the single objective design problem from[17]with the corresponding active pairs, superestructure matrix and dynamic model equations. The full model equations can be found in the Additional file 1.

The three global MINLP methods achieve the solution in substantially less computation time than the outer approximation method and in particular MITSMITS showed the best performance for this example. Whereas the time reported to find the optimum with the outer approximation method was of 200 minutes in an Intel 3.4 GHz Xeon processor [17], MITSMITS arrived to the same solution in less than 200 seconds using an Intel 2.8 GHz Xeon, thus reducing the computational cost at least by a factor of 60. To test the algorithm, we have used as starting guess the zero vector , since the objective function value is very far from the optimum and the constraint is fulfilled. We repeat the analysis starting from different initial guesses fulfilling the constraint and the solver reaches the same solution in similar time. The corresponding convergence curves are illustrated in Additional file 2: Figure S3.

Here it is worthy of note that for the monoobjective problem there exist a number of different circuits with similar performance. In Figure 2 we include four different solutions (circuits 2 to 5) showing very good performance, with values of J below -0.95 (note that by definition the minimum value that J can reach is -1). In absence of additional design criteria, and taking into account that different sources of error limit the precision of biocircuit implementations, the selection of the best design among a set of candidates with close objective function values is rather arbitrary.

Figure 2

Alternative circuit configurations with a maximum of 2 active pairs and similar levels of performance. For every circuit we include the configuration, the superstructure matrix and the objective function value.

For M _max =3, the best solution found is the circuit 6 in Figure 2, with J=-0.999999. Again, MITSMITS showed the best performance, achieving the solution in less than 300 seconds, as it is shown in the convergence curves illustrated in Additional file 2: Figure S4.

For M _max =32, i.e. increasing the maximum level of complexity to 32 pairs (note that this is equivalent to the problem with no constraint on the number of maximum pairs), the best solution found is the circuit 7 in Figure 3 with 14 active pairs. It is important to remark that for increasing levels of complexity the number of solutions with similar values of the objective function (and consequently similar performance) also increases. As an example, we show the circuits 8 to 11 in Figure 3 with similar level of performance and rather different topologies (for space reasons we depict only the superstructure matrix for all circuits except from 7). Note also that in terms of performance, circuit 7 in Figure 3 is equivalent to circuit 6 in Figure 2. This fact leads to arbitrariness when it comes to select the best solution, and suggest the convenience of introducing additional competing criteria in order to provide more realistic design settings.

Figure 3

Circuits found with a maximum of 32 active pairs and similar levels of performance. The superstructure matrix is depicted for every circuit and the value of the objective function J is indicated. The circuit structure shown corresponds to the best circuit found (circuit 7).

Regarding solvers performance we observe that, at least for short computation times, the solution found depends on the initial guess (this dependency increases with complexity) and therefore we test every method starting from different initial guesses. Additional file 2: Figure S5 illustrates the convergence curves of MITSMITS starting from different initial guesses. Concerning the best circuits found, circuits 7, 9 and 11 in Figure 3 were obtained by MITSMITS in less than 1500 s, circuit 8 was found by ACOmiACOmi in less than 3 hours and circuit 10 was found by eSSeSS in less than 300 seconds.

Remarkably, the three methods MITSMITS, ACOmiACOmi and Ess provided solutions with objective function values below J=-0.9999 in less than 300 seconds, for all the initial guesses tested.

Multiobjective global optimization design of a circuit with inducer specific response.

Next, we introduce the protein production cost as an additional criterion to the design problem. Our primary objective is now the performance function J ₁=-Z where Z is given by Eq. (3) and the secondary objective is the cost J ₂=C, where C has been defined in Eq. (6). The problem is solved for increasing levels of complexity, applying the ε-constraint strategy.

For M _max =2 we know the solution from the previous monoobjective analysis, and the value of the cost at this optimum is . We search now the individual optimum for the secondary objective, finding the circuit with active pairs (P l a c ₁,L a c I) and (P l a c ₁,t e t R). Solutions with values of J ₁>0 are discarded. The value of the cost at the optimum is . Taking into account that the upper and lower bounds for the secondary objective function are precisely and , and with a step size of 50, we obtain six non dominated points P ₁…P ₆ corresponding to six circuits with different topologies. The Pareto front is illustrated in Figure 4. The three MINLP solvers have been used in order to compare the results, and an exhaustive search was also implemented, arriving to the same Pareto optimal front. Let us remind that the exhaustive search is possible only for low levels of complexity, since the computation time increases exponentially as the number of maximum pairs increases, as deduced from Eq. 9.

Figure 4

Pareto front for a maximum of 2 active pairs. The circuit configuration and the superstructure matrix are shown for every point of the Pareto front.

Following the same strategy, we compute the Pareto front for M _max =3. The front obtained is shown in Figure 5, and consists of four different points, labeled Q ₁…Q ₄ (note that Q ₂=P ₁).

Figure 5

Pareto front for a maximum of 3 active pairs ( Q ₁ , Q ₂ = P ₁ , Q ₃ , Q ₄ ). The Pareto front for a maximum of 2 active pairs is also shown. Blue dots represent circuits with 2 active pairs and red dots represent circuits with 3 active pairs.

It is of relevance that the solution Q ₄ is significantly better in terms of cost than any other and at the same time it shows a very good performance (J ₁<-0.95). The multiobjective formulation allowed in this case to find a non intuitive topology which is a very good candidate for a successful laboratory implementation. It can be deduced also from Figure 5 how a small increase in complexity from M _max =2 to M _max =3 leaded to significant improvement in the Pareto front, where Q ₃ and Q ₄ are non dominated by any of the circuits with two active pairs (P ₁…P ₆).

Finally, we compute the Pareto front for M _max =32. The circuit Q ₁ (circuit 6 in Figure 2) obtained for M _max =3 is also the best solution found for the unconstrained problem (together with the circuit 7 in Figure 3). By constraining also the minimum level of complexity by setting M _max >3 we obtain the set of non-dominated solutions depicted in Figure 6, together with the corresponding superstructure matrices. In this figure it can be seen that the multiobjective strategy employed allowed us to find points in non-convex regions of the Pareto front, as it is the case of the circuit R ₅.

Figure 6

Pareto front for a minimum of 3 active pairs and a maximum of 32 active pairs. The superstructure matrix is shown for every point of the Pareto front.

Adaptive biocircuit with predefined complexity

Now, starting from the same library of components of the previous example we search for a circuit configuration with the ability for adaptation. We assume that one of the parameters can be manipulated, in this case a kinetic parameter related to the P _tet promoter α _tet (see Additional file 1). As indicated in Methods section the adaptive capacity of the circuit is evaluated by the levels of the output protein LacI in response to a sustained stimulus of aTc, in particular by the levels at its maximum upon induction O _peak , at the steady state before induction O _t=0 and at the steady state upon induction O _t=T . Two competing objective functions are considered, the circuit's sensitivity defined by Eq. (7) and the circuit's precision measured through the formula in Eq. (8). The multiobjective MINLP problem with 32 integer and 1 real decision variables is solved with the ε-constraint strategy proposed, maximizing as a primary objective the sensitivity, i.e. minimizing -S with S defined in Eq. (7) and setting the precision as a constraint. In Figure 7A, we depict one of the solutions of the Pareto front, where P<20 with P defined in Eq. (8). As it is shown in Figure 7B, the circuit is able to adapt upon a sustained stimulus of aTc. The optimal value for the kinetic constant is also indicated.

Figure 7

Adaptive circuit found by multiobjective optimization. A) Topology and superstructure matrix. B) Circuit's response upon aTc stimulus.

Conclusions

In this work we have introduced a multiobjective formulation for the design of biocircuits. The presence of more than one competing objective provides more realistic design settings where the solution is not unique and every solution represents trade-off between different criteria.

The multiobjective optimization in the context of genetic circuit design posed a number of challenges mainly due to the inherent nonlinear nature of the gene circuit's dynamics and the large search spaces involved combining the presence of integer and real variables, which makes the expected Pareto front discrete and possibly non-convex.

In order to overcome these difficulties we made use of global optimization algorithms, showing their efficiency for the MINLP problem resultant of the monoobjective formulation of the design. Then, we provided a multiobjective optimization framework for the design of biocircuits that combines the efficiency of the global MINLP solvers with the capacity to handle multiple design criteria.

Looking for further extensions the method presented is flexible, accommodating to any ODE based modeling framework such that the circuit's model structure is obtained from the starting list of parts by giving values to a set of binary variables.

The advantages of this multiobjective formulation were shown through the design of a biocircuit with specific response upon induction. Due to the efficiency of the global solvers it was possible to obtain in reasonable times the Pareto fronts for different levels of complexity including circuits belonging to non-convex regions of the optimal set of solutions. The capacity to handle circuits with higher number of regulatory regions implies more opportunities for parameter tuning.

Through an illustrative example, we have demonstrated how using this framework we can obtain non intuitive designs to perform a desired functionality setting up a priori the desired level of complexity. This can be useful in future contributions to explore and identify different design principles for synthetic gene circuits.

Authors' contributions

IOM performed the computations. IOM and JRB contributed to the conception and design of the work. JRB supervised the project. IOM and JRB wrote the manuscript. Both authors read and approved the final manuscript.