Design of synthetic genetic oscillators using evolutionary optimization.

Efforts have been made to establish computer models of genetic oscillation. We have developed a real structured genetic algorithm (RSGA) which combines advantages of the traditional real genetic algorithm (RGA) with those of the structured genetic algorithm (SGA) and applies it as an optimization strategy for genetic oscillator design. For the generalized design, our proposed approach fulfils all types of genes by minimizing the order of oscillator while searching for the optimal network parameters. The design approach is shown to be capable of yielding genetic oscillators with a simpler structure while possessing satisfactory oscillating behavior. In silico experiments show effectiveness of the proposed algorithm to genetic oscillator design. In particular, it is shown that the proposed approach performs better than the traditional GAs in the sense that a cheaper structure of genetic oscillators can be obtained.

Introduction

Building a biological circuit by using artificial biologic components has been regarded as a significant tendency in synthetic biology, which is an interdisciplinary application in molecular biology and engineering. Information of the bio-computing process is based on the concentrations of biochemical molecules, like voltage signal in an electrical circuit. The ultimate goal is to construct systems which integrate very large scale bio-circuits and bio-chips, similar to those in the field of electrical engineering. Inspired by electronic logic elements, several recent studies have stated the possibility of realizing a genetic circuit, which is foreseen to have significant contributions to cancer medicine as well as to a variety of bio-energy sources.1,2

For the design of bio-systems, several optimization algorithms were taken into consideration. Among such algorithms, genetic algorithms (GAs) (see, eg, Holland, 1975)3 have been widely applied to engineering optimization problems. The methods were firstly introduced by Holland3 in 1970s. Conventional genetic algorithms emphasize binary coding in chromosomes. However, binary genetic algorithms (BGAs) require excessive computing time when dealing with high-dimensional problems, and the premature convergence of solutions often occurs. To compensate for the disadvantages of BGAs, real genetic algorithms (RGAs) have made changes in floating point coding of chromosomes, and are proven to have advancement on computing speed and precision.4 Furthermore, structured genetic algorithms (SGAs)5 and hierarchical genetic algorithm (HGAs)6 have been proposed to address the premature convergence issue in solving optimization problems. Combining the advantages of RGA and SGA, a novel real structured genetic algorithm (RSGA)7,8 was proposed to solve complex multi-objective optimization problems. This approach exhibits advantages to simultaneously deal with the parameter and structure optimization problems based on a specific structural mapping operator used to determine appropriate numeric values of effective parameter genes in individuals.

Among a variety of problems in the bio-systems design, a problem concerning design of a genetic network with desired dynamic behaviors has attracted much attention from many researchers. In the pioneering work of Elowitz and Leibler,9 a genetic ring oscillator, known as a repressilator, consisting of three genes: tetR, lacI, λcI was proposed. The three genes repress each other like a ring oscillation circuit comprising three BioLogic “NOT” gates, and induce oscillating light in the green fluorescent protein. The transcriptional regulation that generates the oscillatory behavior can be modeled based on several factors, including the transcription rates on repressor concentration, the decay rates of the protein, and the translation rates of mRNA.9 Chen and Chen10 proposed a GA-based approach to achieve a robust design of oscillatory behaviors in genetic networks.

Design of a genetic oscillator with desired amplitude, phase, and frequency could be thought as a tracking design problem in automatic control engineering. However, the stochastic intrinsic dynamic noise and extrinsic disturbances can perturb the oscillatory behaviors, and preservation of robustness to these noise perturbations has become an important issue.2 The existing GA-based design searches for an optimal pair of parameters including the decay rates of protein concentrations and the transcription rates of mRNA. However, this approach can only solve the parameter optimization problem. In the presently described work, we combine the advantages of RGA and SGA to simultaneously determine the optimal structure and the optimal parameters for synthetic biological oscillators. In particular, designing the optimal transcription rates of mRNA and the decay rates of protein concentrations to track the desired sinusoidal signals is considered. Moreover, this paper proposes two methods to improve the tolerance of the intrinsic and extrinsic noises. The first method is based on the increase of searching undefined parameters to extend the dimension of the design parameters. The second method is to adjust hill coefficients to bring changes in protein concentration output. These approaches not only extend the design freedom but also simplify the structure of the oscillator module.

In summary, this paper attempts to develop a more efficient method than the traditional evolutionary algorithm-based approaches for solving the structure and parameter optimization problem while synthesizing the biological genetic oscillators and biological logic gates. In silico experiments show that the RSGA approach is effective in obtaining a genetic oscillator with the cheapest structure.

Real Structured Genetic Algorithm

RSGA is a combination of RGA with SGA, and is a method developed for both the optimal parameter and structure searching. The genetic operations of RSGA include reproduction, crossover, and mutation which were pioneered by Tsai, Huang, and Lin.7,8 The major difference between RSGA and SGA is that both the control genes and parameter genes in the former are real numbers (control genes in SGA are binary numbers whereas parameter genes are real numbers), which improves the mathematical mechanism due to the consistent operators of crossover and mutation, and requires less computing time.

For the current purpose in biogenetic oscillator design, a structured genetic mapping in the RSGA that varies parameters in three ways was applied. The structured mapping generates different structures and searches out the optimal solution structure according to the pre-specified fitness function. During evolution of offspring, it improves the individuals’ fitness by crossover and mutation unceasingly. Based on this feature, the RSGA may achieve the optimal structure and parameters simultaneously for application in design of synthetic biological devices if the index reflects the oscillator structure and the desired amplitude, frequency, and phase. A demonstrative structure of its representation of the individual (chromosome) is shown in Figure 2 which shows a case consisting of 6 control genes, 6 control dependent parameter genes, and some control independent parameter genes.

Figure 2

Activation function under boundary sizing.

Structured genetic mapping

The chromosomes Y = of each individual in RSGA is an order set consisting of the set of control genes c and the set of parameter gene p = with p and p representing the control dependent and control independent parameter genes respectively. The order set of control genes c and the order set of parameter genes p are both real numbers within (Rmin, Rmax). The parameter genes of p are regulated by the control genes. The RSGA structured mapping is equipped with three gene operations: activate, inactivate, and linear ratio variation. The operating function of the control genes is determined by the correlation between the boundary values of the control genes, denoted Bmax and Bmin. When the value of control gene is greater than Bmax, the corresponding parameter genes are activated (ON). If its value is less than Bmin, the corresponding parameter genes are inactivated (OFF). When the value is within (Bmin, Bmax), the parameter gene is regulated by a linear scaling factor. The functions of mapping for various statuses are explained as those shown in Figure 2.

For the current purpose, the order set of parameter genes contains the key parameters for transcription of the biogenetic oscillator which are to be determined, such as transcription rate and sensitivity of mRNA, and decay rates of protein and mRNA, etc. Variations of these parameters are alternated by the control genes during the computational evolutionary process.

The structured genetic mapping from Y to Ỹ is defined as

where the chromosome Ỹ represents an ordered set consisting of c and p̃, and

An operator ⊗ denoting the genetic switch is defined by

where j = 1, …, n, and

with ϕ denoting an empty element. Variations of Bmax and Bmin are generationally dependent and are defined by the following rules:

with

where Rmax and Rmin are the maximum and minimum boundaries, g denotes the current generation, g and g are the initial and final generations, and ΔB = kg with k being a positive constant. The current boundaries Bmax and Bmin would be shifted when g increases such that the parameter searching range extends. Through the boundary sizing technique, the RSGA is able to provide a switch function in the early stage of evolution, emphasizing structural optimization. While the boundary increases gradually with generation number, it is a linear scaling function. When gradually reaching the end of the generation, (Bmin, Bmax) replaces (Rmin, Rmax) and the working algorithm gradually focuses onto the parameter optimization.

Reproduction

In the reproduction process, the algorithm follows the usual manner utilized in evolutionary computational algorithms to create a new population for the next generation from a population in the current generation. The selection operation imitates the mechanism that describes the survival of the fittest in natural selection. The chromosomes are selected for mating, which depends on their relative fitness values, i.e., roulette wheel selection. The chromosome selection probability is given by

where f is the fitness value of the i-th member, and m is the population size, as with usual GAs. The chromosomes of high probability associate with a relative high-fitness value among the population.

Crossover

As with usual GAs, the probability of the chromosome being selected to crossover is P, in general 0.5 ≤ P ≤ 0.9. The crossover mechanism utilizes extrapolation or interpolation to generate new individuals. Initially, it operates in extrapolation. When the parameters of the offspring exceed the allowable ranges Rmax and Rmin, it then switches to interpolation. Therefore, it can avoid parameters varying over the range. The crossover operation is determined by

where x and x refer to the parent individuals, x̃ and x̃ are the offspring of x and x respectively, and λ0 ∈ [0,1] is a uniform random number. Both the control and parameter genes are real numbers and share the same operation of mutation. It simplifies the mechanism of two crossover operators applied in traditional SGAs.5

Mutation

The mutation operator applies randomly chosen individuals to gain fine tuning in chromosomes. The probability of the chromosome being selected to mutate is P, in general 0.01 ≤ P ≤ 0.1. For the mutation operator, we adopt the non-uniform mutation method to change genes in a chromosome which can be realized as

and

where x ∉ {ϕ} is the j-th element of the i-th individuals in the current generation, the value h ∈ {0,1} depends on the random production, xmax (xmin) is the maximal (minimal) j-th element of the i-th individual in the current generation, and y is a scaling factor.

Dynamic probability

Inspired by the characteristic advantages that the gain of Butterworth filter in Bode plot is flat in the passband and approaches zero near the stopband, a dynamic probability is proposed to ensure that the emphasis is placed on the object’s structure first and then on its parameters:

where P is the new probability, P is the current probability, ΔP is the dynamic probability factor, g is the cut-off generation, and q is the order of the dynamic probability. Applying the new dynamic probability, the crossover and mutation probabilities become constant after reaching g. Related probability rates are defined as:

where P is new crossover probability, P is current crossover probability, P is new mutation probability, and P is current mutation probability. In general, g is selected to be the median of the total generation number. This ensures that the emphasis of the optimization will switch towards parameter optimization after reaching g. The design also ensures that the effect of crossover and mutation for parameter genes will not decay to zero.

We can sieve out good genes by evaluating the fitness values and generating elite chromosomes to the next generation. Figure 3 displays the schematic diagram for illustrating the operational process. All chromosomes in the operational process of RSGA are real numbers, thus it does not require any encoding step.

Figure 3

Schematic flowchart of real structured genetic algorithm.

Biologically Genetic Oscillator Design

This section provides a general framework for the design of genetic oscillators based on the RSGA. Several methods to improve tolerance of the intrinsic and extrinsic noises are also proposed in the current framework.

Genetic oscillator

Elowitz and Leibler9 used three transcriptional repressor systems (tetR, lacI, λcI) to build an oscillating network. The first repressor, lacI, inhibits the transcription of the second repressor, tetR, which in turn inhibits the expression of the third gene, λcI. Finally, tetR inhibits lacI’s expression. Because the green fluorescence protein is coupled to a protein in the oscillating network, the oscillations are observed by the fluorescence in the cell. That transcriptional regulation oscillation behavior can be modeled in details based on several factors including mRNA, the dependence of transcription rates on repressor concentration, the decay rates of the protein, and the translation rates of mRNA. The dynamic system can be described by the following coupled first-order differential equations:

where (i,j) ≡ (lacI, λcI), (tetR, lacI), or (λcI, tetR), m ∈ ℝ+ is the concentration of mRNA, p,p ∈ ℝ+ are concentrations of proteins in three genes, α ∈ ℝ is the transcription rate of mRNA, α0 is the leakiness of the promoter, which is usually zero for stable state, β ∈ ℝ+ denotes the ratio of the protein decay rate to the mRNA decay rate, γ ∈ ℝ+ and γ ∈ ℝ+ are decay rates of proteins and mRNA, respectively, n ∈ ℝ is the Hill coefficient, and w,j = 1,2 denotes the effect of environmental perturbations. Fundamental oscillation behavior can exhibit in the concentration of the three repressor proteins.

The three-genes oscillation model was extended to an N coupled genetic model by Strelkowa and Barahona,11 Hori, Hara and Kim,12 Zeiser, Muller and Liebscher,13 and Hori and Hara.14 That is, one can use different number of genes to synthesize an oscillator. Figure 4 shows the configuration of a generalized N-stage gene oscillator. For compactness, we represent a class of stochastic models for the N-gene oscillator with intrinsic fluctuations and extrinsic disturbances as follows:

Figure 4

N-gene oscillator.

where X = [m1p1m2p2 ... mp]T is the state vector, f(X) represents the nonlinear interactions of gene oscillation, g(X)v are intrinsic parameter fluctuations due to random intrinsic noise sources, and w = [w11w12 ... w1w2]T is the extrinsic disturbance vector. The nonlinear terms f(X) and g(X)v are defined, respectively, by

and the Hill function

where

Stable behavior is exhibited when oscillators are constructed with an odd number of repressor genes.15 In other words, protein concentration of that kind of oscillator attracts globally to a stable limit cycle. However, oscillators with even number of genes tend to be a quasi-stable periodic cycle that would diverge after a long period of time.16 When the number of genes is even, the number of repressive loops is even as well. Thus, the traditional way generally adjusts the set of nonlinear equations to ensure the number of repressive loops to be odd. This will guarantee that the bio-system insufficiently robust and attracted to the stable limit cycle.14,16 For example, a four-gene oscillator has four genes and four loops. One should change one nonlinear term of the four to be in activation form, thereby leaving the number of repressive loops at three. While applying the RSGA for oscillator design, one does not have to be concerned about the issue. Rather, the algorithm is able to determine the optimal combination of the repressors and activators.

Structured genetic mapping for genetic oscillator

A chromosome Y consists of two parts: the control gene set c and the parameter gene set p. For designing genetic oscillators, we rearranged the structured mapping operator. The first control gene controls the first and second parameter genes, which represent the decay rate of protein concentration and the transcription rate of mRNA in (16) respectively. The second control gene controls only the third parameter gene, which represents the hill coefficient. The third control gene controls the fourth and fifth parameter genes, which represent the decay rate of protein concentration and the transcription rate of mRNA. The fourth control gene controls only the sixth parameter gene, and so on.

Let the original chromosome be

where

No control independent genes are considered in (21) for briefness of demonstration. After performing the structural mapping, the new chromosome becomes

where

Consider, for example, a randomly generated chromosome given and will be transformed as follows:

Given Rmax = 5, Rmin = 0, g = 10, g = 5, k = 0.1 and the variations of Bmax and Bmin defined, respectively, by

Then, the structured genetic mapping yields

After mapping, the chromosome Y is denoted by Ỹ:

The corresponding linear ordinary differential equations with six key variables (α1,γ1,n1,α2,γ2,n2), after rearrangement of gene order, are obtained accordingly as follows

Fitness function

The objective function consisting of two parts, related to parameter and solution structure, is defined by

where ρ ∈ [0,1] is the weighting factor, J is the normalized performance index, and J is the normalized structure index.

The design objective is to search for the optimal decay rates of proteins (α), the transcription rates of mRNA (γ), and the number of genes in order for the objective function J to be minimized. Furthermore, we would like to discover a solution candidate with the simplest structure for the requirement of a compact structure.

To specify the oscillating signal we consider the reference signal denoted by

where A, ω, and ϕ are amplitude, frequency, and phase, respectively. The goal is to design α and γ that make concentrations of each protein emerging oscillating behavior. The normalized J and J are thus defined as

where T is the integration time period, J can be explained as the objective value of tracking error, J represents the objective function for gene number, and the constants N and nmax are the order of genetic oscillator and the maximum order of genetic oscillator, respectively. The smaller J, the more accurate of tracking will be. A small J implies a simpler object structure. The corresponding fitness function can be expressed as follows

Appropriately selecting α, γ and N to maximize F (α, γ, N, n), or equivalently, minimize J (α, γ, N, n) to obtain the small tracking error, is the objective of the succeeding approach. See Figure 5 for illustration.

Figure 5

RSGA based synthetic genetic oscillator design.

It is notable that oscillators with different gene orders have limited capabilities of amplitude and frequency. However, the structured genetic mapping technique makes it possible to work in different gene-number oscillator module. Therefore, the optimal order can be discovered to fit the desired oscillated amplitude, frequency, and phase. In regard to the accurate performance, we improve the selection, crossover, and mutation probability to a dynamic probability. This idea will be illustrated in the next section.

In Silico Experiments

Consider the cyclic gene regulatory network composed of 1–5 genes to track the sinusoidal wave A sin(ωt + φ) where A = 0.362, ω = 0.075π,∀i, φ = i2π/N. The parameters β = 0.3465 and γm = 0.167 were set as those by Chen and Chen.10

By considering only maximization of the fitness value F with three genes, the method of conventional RGA generated an oscillator with related responses (Fig. 6) with the corresponding fitness value of F converging to 4.2 (Fig. 7).

Figure 6

Three-gene oscillator obtained by RGA.

Figure 7

Fitness convergence of the oscillator obtained by RGA.

We applied a variable number of genes into experiments by the proposed RSGA. The results obtained are summarized in Table 1. While considering maximization of F, the RSGA yielded an oscillator with only two genes and the corresponding fitness value is 4.4. The related responses are illustrated in Figure 8. Convergence of the fitness value is displayed in Figure 9. The result demonstrates superiority of the proposed approach over the conventional GA-based approach,10 in the sense that ours yields comparative tracking performance while requiring less number of genes to realizing the oscillator.

Table 1

Parameters of two-gene oscillator design using RSGA.

α1	γp1	α2	γp2	N	n
0.34	1.285	0.079	0.277	2	5

Figure 8

Two-gene oscillator obtained by RSGA.

Figure 9

Fitness convergence of the oscillator obtained by RSGA.

Robust design of genetic oscillator to tolerate intrinsic and extrinsic noises is more meaningful in vivo. To bring uncertain fluctuations into the module, we consider (17) and show responses of the repressor protein concentration in Figure 10. As is seen in the figure, the amplitude and phase are seriously affected by the intrinsic parameter fluctuations and extrinsic disturbances. This shows that existence of the intrinsic and extrinsic noises would deteriorate accuracy of the oscillation.

Figure 10

Oscillator with noise corruption.

However, the RSGA approach described above and the traditional GA-based approach only determine the decay rate of protein concentration γ and the transcription rate of mRNA α. Their capabilities are restricted because the searching spaces are not robust enough to cover extra performance requirements. The resulting design would fail to track the oscillating commands with different amplitudes and phases under noise interference.

There are two possible ways to improve the deficiency. The first is to add extra parameters for performance tuning: the decay rates of mRNA γ and the ratio of the protein decay rate to the mRNA decay rate β in (16). This extends the parameter search from two to four parameters. It is also possible to adjust the hill coefficient n to bring changes in transcription sensitivity. However, it is more difficult to realize in lab implementation. It was observed from the experimental results that the amplitude and frequency of oscillation increase as the hill coefficient increases. In contrast, a smaller hill coefficient results in smaller amplitude and frequency. As a result, if the hill coefficient is relatively small, the protein concentration will converge to a steady level. See Figure 11 for the previous design.

Figure 11

Oscillators with different hill coefficients.

Combining the two approaches, searching for four parameters γ, α, β, γ and n in each gene, as shown in Figure 12, improves the tracking performance and the corresponding fitness value f rises to 5.25. Figure 13 illustrates convergence of the fitness function. The resulting parameters are summarized in Table 2. However, the downside is that it loses the benefit of cheap network structure as the order of the oscillator is increased to three. The resulting oscillator is still of a simple structure, though.

Figure 12

Oscillator obtained by RSGA while extending the searching parameters to α, γ, γ, β, n.

Figure 13

Fitness convergence of the oscillator obtained by RSGA while extending the searching parameters to α,γ, γ, β, n.

Table 2

Improved oscillator design using RSGA.

γm1	α1	β1	γp1	γm2	α2	β2
1.328	0.395	1.621	0.541	1.621	0.541	1.295
γp2	γm3	α3	β3	γp3	N	n
0.624	1.295	0.264	0.908	0.806	3	4.05

From the numerical results presented above, it can be seen that the RSGA design approach is easy and efficient to apply while dealing with the situation of trade-off in simultaneous parameter and solution structure optimization problem.

Conclusions

RSGA approach has been applied to deal with the optimal design of synthetic genetic oscillator. The proposed approach is able to achieve the synthesized oscillator with a simplified structure and a lower number of parameters than that used by the existing evolutionary computational approach. For the noisy situation, it was observed that more tuning parameters, such as the decay rates of mRNA and the ratio of the protein decay rate to the mRNA decay rate, could be considered to improve the performance of the oscillator while clarifying the benefit of the compact network structure. Finally, a possible extension of our approach is to establish a conceptualized designing framework of steady-state combinational and sequence logic circuits, which may become a foundation while constructing Boolean computing devices.