Abstract
Keywords
Introduction
Although the molecular bases of cellular processes are already known, the behaviours of gene regulatory networks remain poorly understand due to the complexity of their components as well as their multiple specific interactions. The main goal of the nascent field of synthetic biology is to design and construct biological systems that present a desired behaviour.1,2 Synthetic biology is foreseen with important applications in biotechnology and medicine, and with significant contributions to a better understanding of the function of complex biological systems.
3
Recently, several novel approaches to study the oscillatory interactions between genes and proteins of synthetic gene network have arisen.4–6 Two important designs of such networks advanced the study. One of them is the synthetic genetic toggle switch
6
and the other is called the repressilator. In the present study, we focus on the design of a repressilator. The genetic network study of a repressilator was proposed by Elowitz and Leibler,
4
consisting of three genes repressing each other in a closed chain (see Fig. 1). The repressilator can induce oscillation of the intracellular levels of three proteins encoded by the sequence of plasmid which is hosted by the unicellular bacterial microorganism

Construction of the repressilator network in the host cell,
Although the repressilator can show oscillations in its protein concentrations, the design of a repressilator with desired amplitude, frequency and phase is a challenging problem in gene oscillator design. Moreover, the effects of thermal dynamic noise,7,8 mutation, cell division, undefined interactions with the cellular context, 2 and changing extracellular environments can perturb oscillatory properties such as amplitude, frequency and phase. These stochastic behaviours are the consequence not only of current technological limitations but also of internal molecular fluctuations and external noise in the environment of the host cell. Therefore, designing a reliable gene oscillator with a desired amplitude, frequency and phase under these intrinsic parameter fluctuations and extrinsic disturbances on the host cell is another important topic for the design of robust synthetic gene oscillators.4,5,9,10 Robustness to intrinsic and extrinsic noises limits the range of possible circuits to only a few designs that can function properly in the synthetic gene network. 11 Ko et al have proposed a constrained optimization method for system identification of biochemical network. 12 Batt et al 13 have proposed a piecewise-multi-affine method for a desired steady input/output of synthetic gene network and Chen et al14–16 and Lin et al 17 have proposed a robust design to achieve a desired steady state of synthetic gene network. However, these desired steady state design methods can not be applied to transient oscillatory behaviours for a robust gene oscillator design problem. Indeed, these oscillatory behaviours could be represented as limit cycles in phase planes, which are more complex than steady state behaviours with a fixed equilibrium point. Therefore, how to design a robust oscillator around a desired limit cycle is more difficult than a robust synthetic gene network design with a desired constant steady state (i.e. to a desired equilibrium point) by the conventional system stabilization methods.14–16 The reason is that engineering a synthetic gene network with desired oscillatory behaviours is a tracking design problem, and more effort is needed to resolve it. To achieve the design objective of robust gene oscillator, we propose a genetic algorithm approach to imitate the adaptive design mechanisms via natural selection in the real evolutionary process.
Unlike the conventional trial and error methods in gene oscillator designs, a systematic approach is developed for the robust design of oscillatory gene networks with desired amplitude, period, and phase that is robust against internal parameter fluctuations and external disturbances. At first, the desired oscillation is modeled by periodic reference signals with desired amplitudes, frequencies, and phases which are prescribed by the designer beforehand. Then, the robust gene oscillator design problem is transformed to a stochastic optimal tracking design problem, i.e. to specify the kinetic parameters of a gene oscillator to optimally track a reference periodic signal with desired amplitude, frequency and phase despite intrinsic fluctuation and external noise. Finally, a genetic algorithm (GA) is employed to mimic natural selection in an evolutionary process, which is faster than the adaptive design rules by natural selection in the real evolutionary process. 18 This GA-based algorithm enables us to select adequate kinetic parameters to achieve optimal tracking of a repressilator on the host cell (see Fig. 2). GAs have been extensively applied in solving global optimization searching problems.19–26 Since GAs are parallel global search techniques that simulate natural genetic operators,20–22 they are useful when the closed-form optimization technique can not be applied. Further, because the GA-based design algorithm can simultaneously evaluate many points in the kinetic parameter space, it is more likely to converge toward the globally optimal tracking solution. The proposed GA-based design algorithm is inspired by the mechanics of natural selection to a population of binary strings encoding the parameter space. At each generation, it explores different areas of the parameter space, and then directs the search to regions where there is a high probability of finding improved tracking performance (fitness). By working with a population of tracking solutions, the design algorithm can search in effect for many local maximums of the fitness function, and thereby increase the likelihood of finding the global maximum of fitness function (or the optimal tracking performance). Global optimal tracking can be achieved via a number of genetic operators, e.g. reproduction, mutation, and crossover.

Robust synthetic gene oscillator design process. Robust synthetic gene oscillator design process based on stochastic optimal reference tracking via GA searching. Based on the tracking error, the GA-based algorithm can select the design parameter vector
After specifying kinetic parameters in feasible ranges by the proposed GA-based design algorithm for robust optimal tracking design, the robust synthetic gene oscillator could be realized with the recently advanced biological techniques. Combinatorial promoters with multiple TF binding sites, or operators, can assist in the programming gene expression to carry out the designed optimal transcription rates.5,27–29 Another characteristic method of protease modification is to control the degradation rates of protein by fusing ssrA- tagged proteins with corresponding protein, which can reduce the protein degradation of corresponding gene.5,30–32 Recently, a simple method to select adequate biological parts or devices from biological device datasheets (or libraries) to construct a gene network with desired kinetic parameters and decay rates has been proposed.33–35 In this way, synthetic biologists can increase the efficiency of gene circuit design through registries of biological parts from standard datasheets so that these biological parts or devices with desired parameters can be efficiently assembled into a desired synthetic gene oscillator.33–35
Gene oscillators have many useful applications. For instance, gene oscillators can be applied to control the dosage of drugs, e.g. melatonin can be released at night to aid sleeping. 36 Moreover, oscillators are also essential for many biochemical networks which require synchronization among circuit elements.4,37 Further, since all these techniques can be effectively applied to any synthetic gene network design, in the future there should be many potential applications of robust synthetic gene networks.
Finally, a design example is given to describe the design procedure for a desired synthetic gene oscillator using the proposed GA-based design algorithm and to confirm its robust performance under intrinsic parameter fluctuation and extrinsic disturbance.
Methods
Stochastic model for repressilator under intrinsic and extrinsic molecular noises
In the previous repressilator design, the repressilator consists of two plasmids (see Fig. 1), 4 one of which is a plasmid containing three in-chain repressor genes. The other plasmid consists of the reporter gene which encodes the green fluorescence protein (GFP). Because the GFP sequence is coupled to a promoter corresponding to one of proteins in the repressilator, if this protein is produced, it will repress the production of the GFP. That is to say, the oscillation expressions of the repressilator encoded proteins would be presented by the expression level of reporter gene. Hence, the oscillations of the system can be detected by measuring the fluorescence emitted by the cells.
In the repressilator shown in Figure 1, the first repressor protein,
where
In this model, the network behaviour depends on the transcription rate of repressor concentration, the translation rates and decay rates of protein and mRNA. Depending on the values of these parameters, the network may be stable, chaotic or leading to sustained limit-cycle oscillations. Oscillations are favored by gene regulatory networks with strong promoters containing an efficient ribosome-binding site, tight transcriptional repression (low ‘leakiness’), cooperative repression characteristics, and comparable protein and mRNA decay rates.
4
A further obstacle to the design of oscillatory biochemical networks is internal uncertainty, e.g. the thermal fluctuation and the stochastic effects due to the small number of particles involved, characterized as the fluctuations of parameters, and external disturbance on the host cell from the environment. These intrinsic parameter fluctuations and extrinsic molecular noises also may lead sustained oscillations to stable steady states or chaos. Although synthetic oscillators are much simpler than the real biological oscillations, at present these synthetic oscillators still can not work reliably for a long time and need further tuning before application.1,2 It is still difficult to systematically design a synthetic gene oscillator with desired amplitude, frequency and phase specified beforehand by the user. In practical applications, a robust synthetic gene oscillator with the desired amplitude, frequency and phase under intrinsic and extrinsic molecular noises is more useful. More efforts are still needed to achieve this kind of robust synthetic gene oscillator
Therefore, a robust synthetic oscillator network with desired amplitude, frequency and phase is more appealing for synthetic biologist. Before further discussion on the robust design of synthetic biological oscillators, a stochastic model for synthetic biological oscillator with intrinsic fluctuations and extrinsic disturbances
where ∆
Suppose the parametric fluctuations are stochastic as follows
where
i.e.
Substituting (9) into (3) to (8), we get the following stochastic synthetic oscillator
A more general form of stochastic system for synthetic biological oscillator (10) under intrinsic parameter fluctuations and external disturbances in the context of the host cell can be represented by the following nonlinear equation
where
Recent developed technique in synthetic biology may allow tuning promoter, ribosome binding and the protein degradation with relative ease and precision.5,27–32 Hence, we may select biological parts or devices with desired parameters to engineer gene circuits from biological device datasheets (or libraries) in future.33–35 Based on these biotechnologies or biological device datasheets, we want to design the transcription rates
Based on the analysis above, design specifications of robust synthetic gene oscillators are given as follows
Give a desired oscillation with the following desired amplitudes, frequencies, and phases as follows a
The sinusoidal signal in (12) is given only for the convenience of the desired periodic signal. Actually, it can be any periodic signal to be designed for the synthetic genetic oscillator. For example, it could be any signal
i.e.
Give the standard deviations
Give the ranges of feasible design parameters
Then our design objective is to search for
If the above mean-square tracking errors could be minimized by some design parameters
Robust synthetic genetic oscillator design via GA-based design algorithm
It is generally not easy to get a closed-form solution to solve the optimal tracking design problem in (14) for a nonlinear stochastic oscillation system (11) to meet the robust synthetic oscillator design specifications (i)–(iii). In this study, the genetic algorithm (GA) is employed to mimic natural selection in the evolutionary process of a gene oscillator but with a faster evolutionary computation. Genetic algorithms are stochastic optimization algorithms that are originally motivated by the mechanisms of natural selection and evolutionary genetics.20–22 The underlying principles of GAs and mathematical frameworks were presented in Holland's pioneering work, Adaptation in Natural and Artificial Systems. 22 GAs have been proven to be efficient in many areas 21 and more details about GAs can be found in. 20
GAs are powerful searching algorithms based on the mechanics of natural genetic and are inherently parallel because they simultaneously evaluate many points in the parameter space (search space). Therefore, they are very suitable for the robust optimal design of synthetic gene oscillators. In the optimal tracking design problem (14), let us denote the cost function
Our objective is to search for a set of design parameters
where the constants

The inverse relation between the cost function
i.e. the fitness function
The GA-based deign algorithm is an iterative procedure to search for
Generate a population (i.e.
Calculate the fitness in (16) for each string in the population.
Create offspring by GA operators (i.e. reproduction, crossover, and mutation)
Evaluate the new strings and calculate the fitness of natural selection for each string.
If the searching goal is achieved, or an allowable generation is attained, stop and return; else go to (iii).
The convergence of genetic searching algorithm employed for our design problem has been discussed from the viewpoint of schema or similarity template scheme.
20
It can be guaranteed that the optimal fitness function
Chromosome coding and decoding
Since GAs work with a population of binary strings, not the parameters themselves, for simplicity and convenience, binary coding is used in this article. With the binary coding method to transform the phenotype space into a genotype space, the design parameters (phenotype)
where the upper and lower bounds
The resolution (
If (
Reproduction
Reproduction is based on the principle of survival of the better fitness. The fitness of the
Once the strings are reproduced or copied for possible use in the next generation, they are reproduced in a mating pool where they await the action of the other two operators, crossover and mutation (see Fig. 4).

Flow Chart for the design procedure via GA.
Crossover
By the second operator, the strings exchange information via probabilistic decision. Crossover provides a mechanism for strings to mix and match their desirable qualities through a random process. After reproduction, simple crossover proceeds in three steps. First, two newly reproduced strings are selected from the mating pool produced by reproduction. Second, a position along the two strings is selected uniformly at random. This is illustrated below where two binary coded strings, (
The third step is to exchange all characters following the crossing site. For example, the two strings (
Although crossover uses random choice, it should not be thought of as a random walk through the search space. When combined with reproduction, it is an effective means of exchanging information and combining portions of high-quality solutions.
Mutation
Reproduction and crossover give GAs most of their search power. The third operator, mutation, enhances an ability of GAs to find a near-optimal solution. Mutation is the occasional alternation of a value at a particular string position, an insurance policy against the permanent loss of any simple bit, and it is applied with a low probability such that it is chosen so that on average one string in the population is mutated. For example,
In the case of binary coding, the mutation operator simply flips the state of a bit from 0 to 1 at the 9th code or vice versa. Mutation should be used sparingly because it is a random search operator, and with high mutation rates, the algorithm could become a little more than a random search.
The convergence of a genetic search algorithm employed in our design problem can be shown by schema or a similarity template theorem 20 to the maximization of fitness, b i.e. the optimal oscillation tracking for robust parameter design in (14) can be achieved under intrinsic parameter fluctuations and extrinsic disturbances in the design specifications (i)–(iii).
In addition, in the genetic algorithm, the elitist strategy can be incorporated. This strategy directly copies the best chromosome from the old population into the next population to prevent losing the best solutions in the succeeding iteration to improve the genetic algorithm performance.
GAs are more suitable to the iterative synthetic gene oscillator design problem than other major searching methods such as gradient-based algorithms and random searching algorithms for the following reasons. First, the searching space may be very large. Second, the performance surface does not require a differentiability assumption with respect to changes in kinetic parameters of the repressilator. Hence, the gradient-based searching algorithms that depend on the existence of derivatives are inefficient. Third, the likely fit terms are less likely to be destroyed under a genetic operator, thereby often leading to faster convergence. Similarly, other synthetic gene networks with desired output responses can be obtained by the same design procedure.
Results
In order to illustrate the design procedure of the proposed robust synthetic gene oscillators, the following example with numerical simulation is given for the description of the design procedure.
Consider the synthetic gene oscillator example with the same parameters as,
4
except that

Time-responses of protein concentrations, (a) The nominal repressilator time-response with
Suppose the biochemical regulatory network is affected by the four random intrinsic parameter fluctuations from

where
and the external noises
Under the intrinsic fluctuations and extrinsic noises, the nominal and the noise-corrupted protein time-responses of a synthetic gene network with
In order to solve the robust optimal tracking design problem of synthetic gene oscillator via the proposed GA-based design algorithm, we set the GA operators as follows: first, we use the roulette wheel selection to increase the selecting efficiency of the populations which have higher fitness score; second, the crossover rate is 0.8; third, the chromosome mutates uniformly with the mutation rate 0.05; and fourth, we set the constants

Convergence of fitness value the best fitness score and best cost value evolve during the generations, whose relationship is shown in Figure 3 The vibrations of the best cost value and the best fitness score come from the stochastic intrinsic parameter fluctuations and extrinsic noises, which fluctuate in each generation and directly affect the reliability of the synthetic gene network.
Based on the design parameters via the proposed GA-based design method, the time responses of robust synthetic gene oscillator under intrinsic parameter fluctuations and extrinsic noises are shown in Figure 7. Through the robust GA-based design method, we can obtain desired oscillatory behaviour in this repressilator system. In Figure 7(b), the repressilator system with GA optimal solutions shows the robust desired characteristics of oscillation under intrinsic parameter fluctuations and extrinsic disturbances. Although there are still some discrepancies between the desired oscillation signals and the protein concentrations of repressilator, mainly due to the intrinsic parameter fluctuations and extrinsic disturbances, these results are much better than the synthetic design in, 4 as shown in Fig 5b. Clearly, the proposed robust synthetic gene oscillator design method has potential for practical applications in future.

Time-response of the synthetic gene oscillator via the proposed GA-based design method solution. (a) Time-responses of these three proteins. (b) Time-response tracking of each protein and its reference. Under the parameter fluctuations and environmental noises, the designed repressilator can maintain its characteristics of oscillation and function properly. There are still some discrepancies between the desired reference signals and the protein concentrations of the repressilator, which are mainly due to parameter perturbations and environmental noises.
In this
Comprehensive datasheets are used for quantitative descriptions of devices in many standardized engineering disciplines. A synthetic gene network designer can quickly and easily select some desired devices from biological device datasheets to meets their design requirements of a system.
34
Thus, via the help of engineering theory and experience, a conceived system could be constructed by a set of devices with standard characteristics, which are typically reported on datasheets and are common across a wide range of devices type, such as sensors, logic elements and actuators. Recently, biological datasheets have been set as standards for the characterization, manufacture and sharing of information on modular biological devices to promote a more efficient, predictable and design-driven genetic engineering science.33–34 Because datasheets of biological parts or devices embody engineering standards for synthetic biology,
34
a good device standard should show sufficient information about biological parts or devices to allow the design of synthetic gene networks with optimal parameters. Datasheets contain a formal set of context-dependent, input-output behaviours, tolerances, requirements and other details about a particular biological part or device.33–34 Since parameters
But there are many uncertainties about the behaviour of synthetic oscillators. For example, the cellular functions from devices will fluctuate and there are also many parasitic and unpredictable uncertainties among components as well as on the host cell. Since the transcriptional rate,
Discussion
By using the GA-based design method along with Matlab, we can easily solve the design parameters for this optimal reference tracking problem of a robust synthetic gene oscillator under intrinsic and extrinsic noises. However, there are still some disadvantages in the GA method. First, this method requires a great deal of time for the coding and decoding process in the natural selection if the number of design parameters increases. Fortunately, there are many advanced GA methods, like Hybrid Genetic Algorithm (HGA)23,26,38 or the combination with Simulated Annealing Algorithm (SA),24,25 that can save time and increase the probability of finding the global optimum solution. Second, the solution may be only a near-optimum due to limitations of GA method, for example, the limitation of finite bit length
Despite these disadvantages of GA methods, their primary advantage is that the highly nonlinear constrained minimization problem in (14) can be solved for a robust synthetic gene oscillator, which does not have a closed-form solution. To avoid finding a local optimal solution, the proposed GA-based design method can help approach the global optimal solution by the ‘mutation’ and ‘crossover’ processes. Even though the GA method does not always find the global optimal solution, its solution is often close to the optimum, whereas other conventional searching algorithm can only obtain a local optimal solution. By the property of mimicking natural selection in the GA method, most optimal solutions are not reproducible in the repeated biological simulations. This is not surprising because the GA searching process contains not only the different initial conditions but also the different random mutations and crossovers, as in real world evolutionary processes. For example, in some
From an engineering point of view, when we synthesize a prescribed biological oscillator as the repressilator, its function could suffer interference from the intrinsic fluctuations and environmental noises that affect the host cell. These fluctuations and noises will corrupt the synthetic gene oscillator so that it can not achieve the desired behaviour. In this study, we proposed a design procedure using the GA method which mimics the natural selection in the evolutionary process of the real world to optimize the desired reference tracking of synthetic gene oscillator and to tolerate parameter fluctuations and external disturbances on the host cell. In this respect, our design is a rapid selection scheme. This can save the evolution time for optimal selection in the revolutionary process for increasing the robust oscillation characteristics and for improving the reliability of a synthetic gene network. Thus, the time responses in Figure 7b compared with the time response in Figure 5b show that the robustly designed repressilator can efficiently eliminate the effect of uncertainties due to effects of intrinsic parameter fluctuation and the extrinsic noise on the oscillation.
Clearly, the proposed GA-based design method provides a systematic design method for a robust synthetic gene oscillator with desired amplitude, frequency and phase in a host cell with intrinsic parameter fluctuations and external disturbances. Therefore, combined with the recently advanced synthetic techniques such as promoter library, ssrA- tagged protein or the Biobrick assembly standard devices in biological device datasheets, the proposed design method has good potential for practical applications of robust synthetic genetic oscillators in future.
Recently, the synchronization problems of coupled biochemical oscillations have been widely studied.42–45 This is an important topic of synthetic gene oscillators for practical applications. Therefore, the robust synchronization design problem of a large number of coupling synthetic gene oscillators under intrinsic fluctuations and external disturbances will be our future work.
Conclusions
This study proposes a simple but efficient robust synthetic gene oscillator design method via a genetic algorithm. To mimic the natural selection in evolution in order to select adequate design parameters for obtaining a robust synthetic gene oscillator with desired amplitude, frequency and phase under intrinsic parametric fluctuations and extrinsic disturbances on the host cell, the proposed GA-based design method can search for design parameters to achieve the fitness maximization which is equivalent to the optimal tracking of desired oscillation under the effects of intrinsic and extrinsic noises on the host cell. The contributions of this study are given in the following. First, the intrinsic parametric fluctuations and environmental noises can be modelled as state-dependent noises and external disturbances of nonlinear stochastic oscillatory systems to mimic the stochastic behaviour of synthetic gene oscillators in a host cell. Second, the robust oscillator design problem can be formulated as an optimal tracking design problem and then transformed to a fitness maximization problem. Third, based on the fitness function, a GA-based design method is proposed to mimic natural selection in actual evolutionary processes to search for the design parameters of a synthetic oscillator to achieve the desired robust oscillation. The simulation results show that the robustness performance of the synthetic gene oscillator is guaranteed by the proposed design method. Therefore, the proposed GA-based design method has good potential for the practical design of robust synthetic gene oscillators. Further, it can also be extended to the robust design of other synthetic gene networks which could track their desired behaviours.
Disclosures
This manuscript has been read and approved by all authors. This paper is unique and is not under consideration by any other publication and has not been published elsewhere. The authors and peer reviewers of this paper report no conflicts of interest. The authors confirm that they have permission to reproduce any copyrighted material.
