Design of Synthetic Genetic Oscillators Using Evolutionary Optimization

Abstract

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 bio-logic components has been regarded as a significant tendency in synthetic biology, which is an inter-disciplinary 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) ³ have been widely applied to engineering optimization problems. The methods were firstly introduced by Holland ³ 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. ⁴ Furthermore, structured genetic algorithms (SGAs) ⁵ and hierarchical genetic algorithm (HGAs) ⁶ 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, ⁹ 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. ⁹ Chen and Chen ¹⁰ 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. ² 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.

Figure 1.

Chromosome structure of RSGA.

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 = <c, p> of each individual in RSGA is an order set consisting of the set of control genes c and the set of parameter gene p = <p_c, p> with p _c and p _u 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 (R_min, R_max). The parameter genes of p_c 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 B_max and B_max. When the value of control gene is greater than B_max, the corresponding parameter genes are activated (ON). If its value is less than B_min, the corresponding parameter genes are inactivated (OFF). When the value is within (B_min, B_max), 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 $\tilde{Y}$ is defined as $$\tilde{Y}=〈c,\tilde{p}〉$$ (1) where the chromosome $\tilde{Y}$ represents an ordered set consisting of c and $\tilde{p}$ , and $$\tilde{p}=<{\tilde{p}}_c,p_u>$$ (2)

An operator ⊗ denoting the genetic switch is defined by $${\tilde{p}}_c=\left[c_i\right]\otimes {\left(p_{ij}\right)}_{n_j}\equiv \{\begin{array}{l}{p}_{ij},\text{if }{B}_{\max }\le c_i\hfill \\ p_{ij}t,\text{if }{B}_{\min }\le c_i\le B_{\max }\hfill \\ \phi ,\text{if }{c}_i\le B_{\min }\hfill \end{array}$$ (3) where j = 1, …,n_j, and $$t=\frac{c_i-B_{\min }}{B_{\max }-B_{\min }}$$ with φ denoting an empty element. Variations of B_max and B_min are generationally dependent and are defined by the following rules: $$\{\begin{array}{ll}{B}_{\max }=\hfill & B_{\min }=B_{mid},\text{if }{g}_i=g_{init}\hfill \\ B_{\max }=\hfill & R_{mid}+0.5\Delta B,\hfill \\ \hfill & B_{\min }=R_{mid}-0.5\Delta B,\text{if }0<g_i<g_{fin}\hfill \\ B_{\max }=\hfill & R_{\max },B_{\min }=R_{\min },\text{if }{g}_i=g_{fin}\hfill \end{array}$$ (4) with $$R_{mid}=\frac{1}{2}\left(R_{\max }+R_{\min }\right)$$ (5) where R_max and R_min are the maximum and minimum boundaries, g_i denotes the current generation, g_init and g_fin are the initial and final generations, and ΔB = kg_i with k being a positive constant. The current boundaries B_max and B_min would be shifted when g_i 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, (B_min, B_max) replaces (R_min,R_max) 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 $$P_r=\frac{f_i}{\sum _{j=1}^m{f}_j}$$ (6) where f_i 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_c, in general 0.5 ≤ P_c ≤ 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 R_max and R_min, it then switches to interpolation. Therefore, it can avoid parameters varying over the range. The crossover operation is determined by $$\{\begin{array}{l}{\tilde{x}}_{di}=x_{di}-\lambda \left(x_{di}-x_{mi}\right)\hfill \\ {\tilde{x}}_{mi}=x_{mi}+\lambda \left(x_{di}-x_{mi}\right)\hfill \end{array},\begin{array}{l}\text{if }{x}_{di}>R_{\max }\text{ or }{x}_{mi}<R_{\min }\hfill \end{array}$$ (7) $$\{\begin{array}{l}{\tilde{x}}_{di}=x_{di}-\lambda \left(x_{di}-x_{mi}\right)\hfill \\ {\tilde{x}}_{mi}=x_{mi}+\lambda \left(x_{di}-x_{mi}\right)\hfill \end{array},\begin{array}{l}\text{if }{R}_{\min }\le x_{mi}\le x_{di}\le R_{\max }\hfill \end{array}$$ (8) $$\lambda ={\lambda }_0\left(1-\frac{g_i}{g_{fin}}\right)$$ (9) where x_di and x_mi refer to the parent individuals, ${\tilde{x}}_{di}$ and ${\tilde{x}}_{mi}$ are the offspring of x_di and x_mi respectively, and λ₀ ∈ [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. ⁵

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_mut, in general 0.01 ≤ P_mut ≤ 0.1. For the mutation operator, we adopt the non-uniform mutation method to change genes in a chromosome which can be realized as $${\tilde{x}}_{ij}=\{\begin{array}{l}{x}_{ij}+\Delta \left(g_i,x_{ij\max }-x_{ij}\right),\text{if }h=0\hfill \\ x_{ij}-\Delta \left(g_i,x_{ij}-x_{ij\min }\right),\text{if }h=1\hfill \end{array}$$ (10) and $$\Delta \left(g_i,y\right)=y{\lambda }_0\left(1-\frac{g_i}{g_{fin}}\right)$$ (11) where x_ij ∉ {(φ) is the j-th element of the i-th individuals in the current generation, the value h ∈ {0,1} depends on the random production, x_ijmax (x_ijmin) 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: $$P_{new}=P_{cur}\frac{1}{1+{\left(\frac{g_i}{g_c}\right)}^{2q}}$$ (12) $$\Delta P=\frac{1}{1+{\left(\frac{g_i}{g_c}\right)}^{2q}}$$ (13) where P_new is the new probability, P_cur is the current probability, ΔP is the dynamic probability factor, g_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_c. Related probability rates are defined as: $$P_{crossnew}=P_{crosscur}\Delta P$$ (14) $$P_{mutnew}=P_{mutcur}\Delta P$$ (15) where P_crossnew is new crossover probability, P_crosscur is current crossover probability, P_mutnew is new mutation probability, and P_mutcur, is current mutation probability. In general, g_c 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_c. 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 Leibler ⁹ 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: $$\begin{array}{c}\frac{d{m}_i}{dt}=-{\gamma }_{m_i}{m}_i+{\alpha }_i{b}_i\left(p_j^{n_j}\right)+{\alpha }_0+w_{i1}\left(t\right),\\ \frac{d{p}_i}{dt}=\beta m_i-{\gamma }_{pi}{p}_i+w_{i2}\left(t\right)b_i\left(p_j^{n_i}\right)=\frac{1}{1+p_j^{ni}}\end{array}$$ (16) where (ij) = (lacI, λcI), (tetR, lacI), or (λcI, tetR), m_i ∈ ℝ is the concentration of mRNA, p_i, p_j ∈ ℝ are concentrations of proteins in three genes, α_i ∈ ℝ is the transcription rate of mRNA, α₀ 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, γ _pi ∈ ℝ and γ_mi ∈ ℝ are decay rates of proteins and mRNA, respectively, n_i ∈ ℝ is the Hill coefficient, and w_ij, 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, ¹¹ Hon, Hara and Kim, ¹² Zeiser, Muller and Liebscher, ¹³ and Hon and Hara. ¹⁴ 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 TV-gene oscillator with intrinsic fluctuations and extrinsic disturbances as follows: $$\dot{X}=f\left(X\right)+\sum _{i=1}^4{g}_i\left(X\right)v_i+w$$ (17) where X = [m₁ p₁ m₂ p₂ … m_N p_N]^T is the state vector, f(X) represents the nonlinear interactions of gene oscillation, g_i(X)ν_i are intrinsic parameter fluctuations due to random intrinsic noise sources, and w = [w₁₁ w₁₂ … w_N1 w_N2] ^T is the extrinsic disturbance vector. The nonlinear terms f(X) and g_i(X)ν_i are defined, respectively, by $$f\left(X\right)=\left[\begin{array}{c}-{\gamma }_{m_1}{m}_1+b_1\left(p_N\right)\\ {\beta }_1{m}_1-{\gamma }_{p_1}{p}_1\\ -{\gamma }_{m_2}{m}_2+b_2\left(p_1\right)\\ {\beta }_2{m}_2-{\gamma }_{p_2}{p}_2\\ -{\gamma }_{m_3}{m}_3+b_3\left(p_2\right)\\ {\beta }_3{m}_3-{\gamma }_{p_4}{p}_3\\ ⋮\\ -{\gamma }_{m_N}{m}_N+b_N\left(p_{N-1}\right)\\ {\beta }_N{m}_N-{\gamma }_{p_{N+1}}{p}_N\end{array}\right],g_1\left(X\right)v_1=\left[\begin{array}{c}\delta {\gamma }_{m_1}\\ 0\\ \delta {\gamma }_{m_2}\\ 0\\ \delta {\gamma }_{m_3}\\ 0\\ ⋮\\ \delta {\gamma }_{m_N}\\ 0\end{array}\right]$$ (18) $$g_2\left(X\right)v_{2=}\left[\begin{array}{c}\delta b_1\left(p_N\right)\\ 0\\ \delta b_2\left(p_1\right)\\ 0\\ \delta b_3\left(p_2\right)\\ 0\\ ⋮\\ \delta b_N\left(p_{N-1}\right)\\ 0\end{array}\right],g_3\left(X\right)v_3=\left[\begin{array}{c}0\\ \delta {\beta }_1\\ 0\\ \delta {\beta }_2\\ 0\\ \delta {\beta }_3\\ ⋮\\ 0\\ \delta {\beta }_N\end{array}\right],g_4\left(X\right)v_4=\left[\begin{array}{c}0\\ \delta {\gamma }_{p_1}\\ 0\\ \delta {\gamma }_{p_2}\\ 0\\ \delta {\gamma }_{p_3}\\ ⋮\\ 0\\ \delta {\gamma }_{p_N}\end{array}\right]$$ (19) and the Hill function $$b_i\left(p_j^{n_i}\right)=\frac{1}{1+p_j^{n_i}}$$ (20) where $$n_i\triangleq \{\begin{array}{ll}\ge 0,\hfill & \text{for repression}\hfill \\ <0,\hfill & \text{for activation}\hfill \end{array}$$

Figure 4.

N-gene oscillator.

Stable behavior is exhibited when oscillators are constructed with an odd number of repressor genes. ¹⁵ 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. ¹⁶ 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 $$Y=<c,p>=<\stackrel{C}{\overbracket{c_1\left({\alpha }_1,{\gamma }_{p1}\right)c_2\left(n_1\right)c_3\left({\alpha }_2,{\gamma }_{p2}\right)c_4\left(n_2\right)\cdots c_{2N-1}\left({\alpha }_N,{\gamma }_{pN}\right)c_{2N}\left(n_N\right),}}p>$$ (21) where $$p=\left({\alpha }_1,{\gamma }_{p1},n_1,{\alpha }_2,{\gamma }_{p2},n_2,\cdots ,{\alpha }_N,{\gamma }_{pN},n_N\right)$$ No control independent genes are considered in (21) for briefness of demonstration. After performing the structural mapping, the new chromosome becomes $$\tilde{Y}=<c,\tilde{p}>$$ where $$\tilde{p}=\left[c_i\right]\otimes \left[p_i\right]\equiv \{\begin{array}{c}{p}_i,\text{if }{B}_{\max }\le c_i\left({\alpha }_i,{\gamma }_{pi}\right)\\ p_{i}t,\text{if }{B}_{\min }\le c_i\left({\alpha }_i,{\gamma }_{pi}\right)\le B_{\max },\\ \phi ,\text{if }{c}_i\left({\alpha }_i,{\gamma }_{pi}\right)\le B_{\min }\end{array}t=\frac{c_i\left({\alpha }_i,{\gamma }_{pi}\right)-B_{\min }}{B_{\max }-B_{\min }}$$ (22)

Consider, for example, a randomly generated chromosome given and will be transformed as follows: $$\begin{array}{c}Y=<c,p>\\ =<\stackrel{C}{\overbracket{0.52\hspace{1em}1.5\hspace{1em}1.35\hspace{1em}2\hspace{1em}4.6\hspace{1em}3.8\hspace{1em}0.02\hspace{1em}3.6\hspace{1em}2.33\hspace{1em}5.1}},\stackrel{p}{\overbracket{1.58\hspace{1em}3.45\hspace{1em}4.79\hspace{1em}0.98\hspace{1em}0.11\hspace{1em}2.63\hspace{1em}3.01\hspace{1em}1.44\hspace{1em}0.7\hspace{1em}1.33\hspace{1em}3.52\hspace{1em}1.02\hspace{1em}4.4\hspace{1em}0.1\hspace{1em}-2}}>\end{array}$$

Given R_max = 5, R_min = 0, g_fin = 10, g_i = 5, k = 0.1 and the variations of B_max and B_min defined, respectively, by $$\begin{array}{c}{B}_{\max }=R_{mid}+\frac{1}{2}\Delta B=2.75\\ B_{\min }=R_{mid}-\frac{1}{2}\Delta B=2.25\end{array}$$ Then, the structured genetic mapping yields $$\tilde{p}=\left[c_i\right]\otimes \left[p_i\right]\equiv \{\begin{array}{l}{p}_i,\text{if }2.75\le c_i\hfill \\ p_{i}t,\text{if }2.25\le c_i\le 2.75,\hfill \\ \phi ,\text{if }{c}_i\le 2.25\hfill \end{array}t=2.25\frac{c_i\left({\alpha }_i,{\gamma }_{pi}\right)}{2.75-2.25}$$

After mapping, the chromosome Y is denoted by $\tilde{Y}$ : $$\begin{array}{c}\tilde{Y}=<c,c\otimes p>\\ =<\stackrel{C}{\overbracket{0.52\hspace{1em}1.5\hspace{1em}1.35\hspace{1em}2\hspace{1em}4.6\hspace{1em}3.8\hspace{1em}0.02\hspace{1em}3.6\hspace{1em}2.33\hspace{1em}5.1}},\stackrel{\tilde{p}}{\overbracket{\underset{\text{off}}{\underbrace{\phi \hspace{1em}\phi }}\hspace{1em}\underset{\text{off}}{\underbrace{\phi }}\hspace{1em}\underset{\text{off}}{\underbrace{\phi \hspace{1em}\phi }}\hspace{1em}\underset{\text{off}}{\underbrace{\phi }}\hspace{1em}\underset{\text{on}}{\underbrace{3.01\hspace{1em}1.44}}\hspace{1em}\underset{\text{on}}{\underbrace{0.7}}\hspace{1em}\underset{\text{off}}{\underbrace{\phi \hspace{1em}\phi }}\hspace{1em}\underset{\text{on}}{\underbrace{1.02}}\hspace{1em}\underset{\text{linear scaling}}{\underbrace{0.704\hspace{1em}0.016}}\hspace{1em}\underset{\text{on}}{\underbrace{-2}}}}>\end{array}$$

The corresponding linear ordinary differential equations with six key variables (α₁γ_p1n₁α₂γ_p2n₂), after rearrangement of gene order, are obtained accordingly as follows $$\{\begin{array}{c}{\dot{m}}_1=-{\gamma }_{m1}{m}_1+\frac{3.01}{1+p_2^{0.7}},\\ {\dot{p}}_1=\beta m_1-1.44p_1,\\ {\dot{m}}_2=-{\gamma }_{m1}{m}_2+\frac{0.704p_1^2}{1+p_1^2},\\ {\dot{p}}_2=\beta m_2-0.016p_2\end{array}$$

Fitness function

The objective function consisting of two parts, related to parameter and solution structure, is defined by $$J_{tot}\left({\alpha }_i,{\gamma }_{pi},N,n\right)=\underset{{\alpha }_i\in \left[{\alpha }_{ia}\hspace{1em}{\alpha }_{ib}\right],{\gamma }_{pi}\in \left[{\gamma }_{pia}\hspace{1em}{\gamma }_{pib}\right],n\in \left[n_{ia}\hspace{1em}{n}_{ib}\right],i=1,2,\cdots ,N}{\min }\left[\rho J_p+\left(1-\rho \right)J_s\right]$$ (23) where ρ ∈ [0,1] is the weighting factor, J_p is the normalized performance index, and J_s is the normalized structure index.

The design objective is to search for the optimal decay rates of proteins (α_i), the transcription rates of mRNA (γ_pi), and the number of genes in order for the objective function J_tot 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 $$r_i\left(t\right)=A_i\sin \left({\omega }_{i}t+{\phi }_i\right),i=1,2,\dots ,N$$ (24) where A_i, ω_i and φ are amplitude, frequency, and phase, respectively. The goal is to design α and γ_Pi that make concentrations of each protein emerging oscillating behavior. The normalized J_p and J_s are thus defined as $$J_p=\frac{\sum _{i=1}^N\underset{0}{\overset{T}{\int }}{\left(p_i-A_i\sin \left({\omega }_{i}t+{\phi }_i\right)\right)}^{2}dt}{T\sum _{i=1}^N{A}_i^2},J_s=\frac{N}{n_{\max }}$$ (25) where T is the integration time period, J_p can be explained as the objective value of tracking error, J_s represents the objective function for gene number, and the constants N and n_max are the order of genetic oscillator and the maximum order of genetic oscillator, respectively. The smaller J_p, the more accurate of tracking will be. A small J_s implies a simpler object structure. The corresponding fitness function can be expressed as follows $$F\left({\alpha }_i,{\gamma }_{pi},N,n\right)=1/J_{tot}\left({\alpha }_i,{\gamma }_{pi},N,n\right)$$ (26)

Appropriately selecting α_i, γ_pi and N to maximize F(α_i, γ_pi, N, n), or equivalently, minimize J(α_i, γ_pi, 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_isin(ω_it + φ_i) where A_i=0.362, ω_i = 0.075π,∀i, φ = i2π/N. The parameters β = 0.3465 and γ_mi = 0.167 were set as those by Chen and Chen. ¹⁰

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, ¹⁰ 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 γ_Pi and the transcription rate of mRNA α_i. 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 γ_mi 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 γ_mi, α_i β_i γ_Pi 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 α_i γ_pi, γ_mi, β_i n.

Figure 13.

Fitness convergence of the oscillator obtained by RSGA while extending the searching parameters to α_i, γ_pi, γ_mi, β_i 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.

Author Contributions

Conceived and designed the experiments: Y-CC, C-LL. Analyzed the data: Y-CC, C-LL. Wrote the first draft of the manuscript: Y-CC, C-LL. Contributed to the writing of the manuscript: Y-CC, C-LL, TJ. Agree with manuscript results and conclusions: Y-CC, C-LL, TJ. Jointly developed the structure and arguments for the paper: Y-CC, C-LL, TJ. Made critical revisions and approved final version: Y-CC, C-LL, TJ. All authors reviewed and approved of the final manuscript.

Competing Interests

Author(s) disclose no potential conflicts of interest.

Disclosures and Ethics

As a requirement of publication author(s) have provided to the publisher signed confirmation of compliance with legal and ethical obligations including but not limited to the following: authorship and contributorship, conflicts of interest, privacy and confidentiality and (where applicable) protection of human and animal research subjects. The authors have read and confirmed their agreement with the ICMJE authorship and conflict of interest criteria. The authors have also confirmed that this article is unique and not under consideration or published in any other publication, and that they have permission from rights holders to reproduce any copyrighted material. Any disclosures are made in this section. The external blind peer reviewers report no conflicts of interest.