Abstract
Background
Gene networks in nanoscale are of nonlinear stochastic process. Time delays are common and substantial in these biochemical processes due to gene transcription, translation, posttranslation protein modification and diffusion. Molecular noises in gene networks come from intrinsic fluctuations, transmitted noise from upstream genes, and the global noise affecting all genes. Knowledge of molecular noise filtering and biochemical process delay compensation in gene networks is crucial to understand the signal processing in gene networks and the design of noisetolerant and delayrobust gene circuits for synthetic biology.
Results
A nonlinear stochastic dynamic model with multiple time delays is proposed for describing a gene network under process delays, intrinsic molecular fluctuations, and extrinsic molecular noises. Then, the stochastic biochemical processing scheme of gene regulatory networks for attenuating these molecular noises and compensating process delays is investigated from the nonlinear signal processing perspective. In order to improve the robust stability for delay toleration and noise filtering, a robust gene circuit for nonlinear stochastic timedelay gene networks is engineered based on the nonlinear robust H_{∞ }stochastic filtering scheme. Further, in order to avoid solving these complicated noisetolerant and delayrobust design problems, based on TakagiSugeno (TS) fuzzy timedelay model and linear matrix inequalities (LMIs) technique, a systematic gene circuit design method is proposed to simplify the design procedure.
Conclusion
The proposed gene circuit design method has much potential for application to systems biology, synthetic biology and drug design when a gene regulatory network has to be designed for improving its robust stability and filtering ability of diseaseperturbed gene network or when a synthetic gene network needs to perform robustly under process delays and molecular noises.
Background
Gene expression involves a series of molecular events, which include binding of regulators, transcription, splicing, translation, posttranslation modification and diffusion. As each of these molecular events is subject to significant biochemical process delays, intrinsic fluctuations, and extrinsic disturbances, gene expression is best viewed as a nonlinear stochastic process with multiple time delays [14]. Even in cases where population measurements are regular and reproducible, singlecell measurements often display significant heterogeneity [1]. In general, these observations suggest that the molecular events underlying cellular physiology at the nanomolar scale are easily subject to fluctuations and biochemical process delays so that we have to propose a nonlinear stochastic model with multiple time delays for gene expression and biochemistry [24]. Therefore, molecular noise processing and delay compensation in nonlinear stochastic gene network with process delays is an important topic to understanding how cells function and process information when the underlying molecular events are random with biochemical process delays. As pointed out in [1,3,4], this topic is one of the most challenging and fascinating problems for systems biologists, since it opens questions in physiology, development and evolutionary biology.
When the underlying molecular events are basically random and timedelayed, how does the physiology of the cell remain highly orchestrated and robust? Most cellular events are orderly and precisely regulated in spite of the stochastic function and process delays of gene regulatory circuits within cells [1,5,6]. Without consideration of biochemical process delay, a stochastic differential equation or the Langevin equation has recently been employed to describe the molecular fluctuation in gene networks [3,7]. Many algorithms have been developed for simulating the Langevin equation to calculate the probability density function [5,6,8,9]. The FokkerPlank equation is employed to describe the evolution of the probability function [1,10]. Most researchers analyze these stochastic models, using the MonteCarlo method, such as Gillespie algorithm, StockSim algorithm and so on to describe the evolution of biochemical networks via the discrete stochastic model. Even these modeling tools allow us to address questions concerning intracellular noise. However, if biochemical process delays have not been considered in the dynamic model of biochemical gene networks, an engineered biochemical network based on this model may lead to fluctuation, oscillation or even blowing up [11]. Therefore, biochemical process delays and molecular noises must be considered in the dynamic model to mimic the realistic cellular behaviors of biochemical gene network in cell. Recently, it is also found that most newly synthetic networks are nonfunctioning and need tuning owing to process delays, parameter fluctuations, due to thermal fluctuation, gene expression noise, mutation and extrinsic noises, due to changing extra cellular environments, interactions with cellular context [12]. Therefore, a systematic design method for noisetolerant and delayrobust gene networks is an important topic in systems biology and synthetic biology for an engineered gene network to work properly in host cell [13].
In recent years, many researchers have found it useful to invoke analogies from signal processing when investigating molecular noise. In this situation, a biochemical pathway is viewed as an analogy filter and classified in terms of its frequency response from the signal processing perspective. In terms of signal processing, the biochemical pathways function as lowpass filters, as they transduce lowfrequency signals and attenuate highfrequency signals. In addition to intrinsic chemical damping, negative feedback [14,15], integral feedback [16], and many other simple mechanisms such as redundancy mechanisms are found to attenuate molecular noise in biochemical systems. It is also found that the effect of molecular noise is also amplified in some sense by autocatalytic mechanisms (positive feedback) to give rise to population heterogeneity (and so diversity). Even some of the elementary mechanisms for molecular noise attenuation and amplification enumerated above seem simple and identifiable. However, elementary mechanisms typically do not function in isolation; rather they interact in complex gene networks involving multiple feedback loops. Although it is straightforward to understand how a single feedback loop shapes molecular noise, it is far more difficult to understand composite behavior of multiple mechanisms interconnected in complex architecture of gene networks [1,8,9]. These molecular noiserelated problems are called molecular noise filtering problems of gene network by molecular biologists [13,1719]. The analysis of molecular noise attenuation and amplification of gene regulatory network without process delays has been recently discussed from the stochastic system point of view [17]. At present, though theoretical and computational tools exist for analyzing the filtering properties of a given network, as pointed out by Rao et al. in [1], no good theory exists for identifying all possible systematic mechanisms that generate robust networks to compensate biochemical process delays and to attenuate molecular noises simultaneously.
It is clear that complex gene networks are able to function reliably despite inherent noise attributable to molecular fluctuations and unavoidable time delays due to biochemical process. Robustness has been hypothesized as an intrinsic property of intracellular networks. Although robustness is often studied independent of noise and time delay, the two problems are related. When studying robustness, the typical question is how sensitive the behavior of a gene network is to kinetic parameter perturbations in the model [13]. As these parameter perturbations are subject to molecular fluctuations, a molecular noiseresistant network is likely to be robust. Nevertheless, a gene network that is insensitive to kinetic parameter perturbations may still be sensitive to extrinsic molecular noises. A robust filtering of biochemical networks under parameter perturbations has been discussed based on the global linearization technique in [20]. However, this method is only with Ssystems and without consideration of process delays. In [21], based on Ssystem model [13], a robust circuit design is proposed for nonlinear deterministic biochemical networks from steady state approach and without considering process delay. The robust stabilization method for gene network against molecular noises has been discussed in [22] without consideration of timedelays in biochemical processes.
Time delays are common and substantial in biochemical process. They can protect biochemical systems against transient loss of input signal. Delay for proofreading is widely employed to increase the fidelity of recognition to provide security against misrecognition. Further delay can filter the nonbeneficial pulses [11]. However, timedelays may play a negative role in stability of gene networks. In [23,24], the concepts of sensitive edge and robust edge are used to analyze the robustness and fragility of synchronization. Then the authors want to increase the robustness of synchronization by its associated feedback loops to protect against attacks to the complex networks. Further, the robust synchronization designs of complex timedelayed networks are also discussed in [25,26] to increase the robustness of synchronization to tolerate timedelays and protect against attacks to the networks. In this study, the random parameter fluctuations are modeled as statedependent intrinsic molecular fluctuations of the stochastic gene network model in which multipletime delays of biochemical process is also considered to mimic the process delays in gene regulation processes. Both the upstream molecular noise and the global molecular noise affecting all genes are modeled as extrinsic molecular noises. The robust molecular noise attenuation problem will be systematically discussed for gene networks with multiple process delays from the nonlinear H_{∞ }stochastic signalprocessing perspective.
The H_{∞ }filtering theory has been developed for state estimation of nonlinear stochastic signal processing to investigate the noise attenuation problem via minimizing the worstcase effect of noise on the filtering error [27,28] and the H_{∞ }control theory has been developed to robustly stabilize the nonlinear stochastic system under parameter perturbations [29,30]. However, the time delay of the stochastic process is not considered. Recently, various timedelay control designs for nonlinear systems and fuzzy systems have been discussed in [25,26,31,32]. In this study, the robust nonlinear stochastic H_{∞ }filtering theory is applied to discuss the robust stability and noise filtering problems of nonlinear stochastic gene networks under both multipletime delays due to biochemical process delays and stochastic intrinsic molecular noise due to parameter fluctuations. Furthermore, the filtering ability of extrinsic molecular noises will be also investigated at each gene of gene network from the H_{∞ }filtering point of view. For the biotechnological purpose or drug design purpose, if robust stability is required to tolerate a prescribed range of intrinsic parameter variations and the filtering ability of gene network to filter extrinsic molecular noises from the environment has to be improved, a robust gene circuit design scheme needs to be developed for gene networks from the nonlinear noise filtering perspective to achieve the design purpose via transfection and transformation biotechnologies. For recent metabolic engineering and future synthetic biology [7], if a synthetic gene regulatory network is required to work near a desired equilibrium point in an environment with process delays, intrinsic molecular fluctuations, and extrinsic molecular noises, a robust control circuit design is necessary to improve the robust stability and molecular noise filtering ability of the synthetic gene network to achieve its required performance. This robust filtering design method will be potential for robust gene circuit design in future, from which gene therapy and drug design could be developed. However, the drawback of these nonlinear stochastic filtering approaches needs to solve a nonlinear differential HamiltonJacobi inequality (HJI) for the robust stabilization and filtering design of gene networks. In general, we still have no analytic solution or numerical solution for HJI, except for some simple cases. In this study, to avoid solving HJI, a fuzzy interpolation method is employed to interpolate several local linear gene networks at different operation points to approach the nonlinear gene network to simplify the design produce of robust stabilization and filtering of gene networks so that the robust gene circuit design could be easily designed and implemented.
In recent years, TakagiSugeno (TS) fuzzy systems have efficiently interpolated several local linear systems via fuzzy bases to approximate a nonlinear system [33,34]. Therefore, a TS fuzzy stochastic gene network is employed to approximate a nonlinear stochastic gene network with time delays by interpolating several local linear stochastic gene networks with time delays at different operation points of the nonlinear stochastic gene network. In this study, based on fuzzy interpolation approximation, the robust stabilization problem for tolerating process delays and intrinsic molecular fluctuations and H_{∞ }filtering design problem for attenuating extrinsic molecular noises of gene networks can be efficiently solved via a set of linear matrix inequalities (LMIs) developed from a set of fuzzy local linear stochastic systems. These LMIs could be easily solved via LMI toolbox in Matlab [35] to simplify the design procedure. Furthermore, based on fuzzy approximation method, the circuit design procedure for robust stability and noise filtering of a nonlinear stochastic gene network with process delay could be simplified from the linear gene network design point of view. In this situation, more biological robust stabilization and filtering insight could be investigated from the linear system point of view. We found that if the eigenvalues of these local linear gene networks are in the far lefthand side of the complex domain (i.e., with more negative real parts or more stable), then the gene network will tolerate more process delays and intrinsic molecular fluctuations, and could also filter more extrinsic molecular noises to achieve the design purpose of gene networks. The robust stabilization and H_{∞ }molecular noise filtering circuit design could be achieved by shifting the eigenvalues of fuzzy linear gene networks to the far lefthand side of the complex domain through engineering appropriate negative feedback loops via the proposed robust gene circuit design method for nonlinear stochastic timedelayed gene networks.
In general, the proposed robust filter design in nonlinear stochastic gene network is different from the conventional robust fuzzy H_{∞ }filter design for engineering control or signal processing systems [2731]. In the conventional fuzzy H_{∞ }filter design [2729], a dynamic fuzzy estimator is employed to estimate system state variables from the noisy measurement data, which needs a very complicated computation for state estimation at any time instant. In this study, only some gene circuits are implemented (or embedded) directly in gene network to achieve robust stabilization to tolerate process delay and intrinsic molecular fluctuations and to filter the extrinsic disturbances. Therefore, the proposed robust circuit design is different from the conventional control and filter designs in control engineering and signal processing in [2731]. Gene circuits are implemented by inserting the transcription factor (TF) binding site of regulatory gene product, i.e. the binding site of protein produced by regulatory gene, to the promoter region of regulated gene by the techniques of transfection and transformation [7,14]. The kinetic parameters of gene circuits of molecular noise filter are proportional to the length of binding sites of TF inserted to the promoter region of regulated genes and will be specified by the designer. In order to simplify the design procedure, only a few gene circuits are implemented. In our design procedure, the TS fuzzy approximation method is only employed to simplify the specification of kinetic parameters of gene circuit of molecular noise filter, i.e., the specification procedure of lengths of inserted binding sites of TF produced by regulatory genes in gene circuit design procedure. Recent experimental advances in sequencing, genetic engineering and bionanotechnology will make this engineered circuit implementation feasible in near future [7,3641].
Because disease may perturb the normal gene regulatory network through genetic perturbation and/or by pathological environmental molecular noise such as infection agents or chemical carcinogens [42], in future synthetic biology [7,12,19] and systems biology [13,36,42], we could improve the robust stability and molecular noise filtering ability of gene network for drug design and gene therapy by the proposed gene circuit design method via transfection and transformation biotechnologies, so that a designed gene network could be protected from the influence of process delays, intrinsic molecular fluctuations and environmental molecular noises, and thus work more reliably. Finally, two design examples in silico are given to illustrate the design procedure and to validate the performance of the proposed robust H_{∞ }gene circuit design method for nonlinear stochastic gene networks under process delays, intrinsic molecular fluctuations and extrinsic molecular noises.
Results
Firstly, a nonlinear stochastic timedelayed model is developed to mimic the realistic dynamic behavior of a gene network under process delays, random intrinsic and extrinsic molecular noises. Then the robust stability and the filtering ability of nonlinear stochastic timedelayed gene network are discussed from the H_{∞ }signal processing perspective. Further, in order to improve the robust stability and filtering ability to tolerate more parameter fluctuations, process delays and to attenuate much environmental molecular noises, based on fuzzy approximation and LMI technique, a systematic design method is proposed for nonlinear stochastic timedelayed gene networks. Finally, a simple design procedure is developed for the proposed robust gene circuit design method. These results are described in the following subsections in detail.
Robust Stability and Filtering Ability Analysis for Nonlinear Stochastic TimeDelayed Gene Networks
For the convenience of illustration, the following linear gene network with process delay is introduced first,
x(t) = x_{0}(t)∀t ∈ [τ, 0], where x(t) = [x_{1}(t),...,x_{N}(t)]^{T }denotes the vector of mRNA concentrations of N genes; the maximum time delay τ = max{τ_{k}, k ∈ [1, m]}; A_{0 }and A_{k }denote the realtime and delaytime kinetic interaction matrices among these genes, respectively. For example, x(t) is the mRNA expression vector of N genes; the diagonal component a_{0,ii }= λ_{i }of A_{0 }denotes the degradation of mRNA of ith gene with decay rate λ_{i}; the i, j component a_{k,ij }of A_{k }denotes the regulatory interaction from gene j to gene i to activate or repress gene i with time delay τ_{k}. The delay time τ_{k }may be due to the process time of gene transcription, translation, posttranslation protein modification and diffusion, which are needed for regulatory genes to product proteins (e.g. transcription factors TFs) and then convey them to their target genes (see Fig. 1). In general, in real biochemical systems, these biochemical process delays may be very long and can not be neglected in biochemical dynamic models, especially for practical gene network design applications.
Figure 1. A linear gene regulatory network of N genes under process delays, intrinsic molecular fluctuation Δa_{k,ij}n_{k}(t) and extrinsic molecular noise v(t).
Suppose the linear gene network in (1) suffers from intrinsic molecular fluctuations so that a_{k,ij }→ a_{k,ij }+ Δa_{k,ij}n_{k}(t), k = 0,1,...,m, where Δa_{k,ij }denotes the amplitude of the stochastic kinetic fluctuation and n_{k}(t) is a white noise with zero mean and unit variance, i.e., Δa_{k,ij }denotes the deterministic part of fluctuations and n_{k}(t) absorbs the stochastic property of intrinsic molecular fluctuations (see Fig. 1). Intrinsic molecular fluctuations are mainly attributed to thermal fluctuation and random molecular events in the transcription, splicing, and translation processes of gene expression [13]. The covariances of stochastic intrinsic molecular fluctuations Δa_{k,ij}n_{k}(t), k = 0, 1,...,m, are as follows
for k = 0, 1,...,m, where denotes the delta function as follows
According to the above analysis, the gene network under stochastic molecular fluctuations can be represented as
where the fluctuation matrices B_{k }are given by
If a_{k,ij }is free of fluctuation, then Δa_{k,ij }= 0 in B_{k}. By the Itô stochastic differential equation, the stochastic gene network equation in (3) is equivalent to the following stochastic process [27,29,31]
where W_{k}(t) is a standard Wiener process or Brownian motion with dW_{k}(t) = n_{k}(t)dt.
Actually, in real gene networks, the dynamic gene regulatory equations are always nonlinear. In this situation, the linear dynamic gene network regulatory equation in (4) should be modified as the following Langevin nonlinear stochastic equation with multiple time delays [1,3,10]
where the first two terms on the righthand side denote the nominal nonlinear interactions of gene network and the last two terms denote the effect of intrinsic molecular fluctuations on the gene network, which are statedependent and will influence the stability of nominal gene network.
Remark 1 If the Langevin equation in (5) is linearized at an operation point, it should be a linear stochastic gene network as (4).
For simplicity, the linear stochastic gene network in (4) and nonlinear stochastic gene network in (5) will be reformulated, respectively, as follows
and
where τ_{0 }= 0 and τ_{k }> 0, k = 1,...,m, denote the corresponding time delays.
We say that the stochastic gene network is asymptotically stable in probability at the equilibrium point x = 0 if there exists a Lyapunov (powerlike) function V(x(t)) > 0, with V(0) = 0, in the neighborhood of the equilibrium point such that the expectation of the derivative of V(x(t)), , i.e., the power of the gene network is decreasing. In general, there exist many equilibrium points for nonlinear gene networks in (7). Gene networks perform their biochemical function within some local region of an equilibrium point, which is called the phenotype of the gene network. For the convenience of design, only the robust stability of the equilibrium point at x = 0 is discussed for the nonlinear stochastic gene network in (7). If the interested equilibrium point (steady state) of the nonlinear timedelayed gene network is not at x = 0, the origin should be shifted to the equilibrium point or the steady state. In this situation, the robust stabilization and filtering ability at the origin is equivalent to the stabilization and filtering ability at the interested equilibrium point of the gene network. According to the stochastic gene network model in (7) and the nonlinear stochastic stability theory in [27,30,31], the robust stability and molecular noise filtering ability of stochastic timedelayed gene networks are analyzed in the following propositions.
Proposition 1 Assume there exists a positive Lyapunov function V(x(t)) > 0 and V (0) = 0 satisfying the following secondorder HamiltonJacobi inequality (HJI)
for all nonzero x(t) ∈ R^{N}, then the equilibrium point x = 0 of the nonlinear stochastic gene network in (7) is asymptotically stable in probability, i.e., the effects of the time delays and intrinsic molecular fluctuations could be tolerated by the gene network.
Proof. See Appendix A. ■
Remark 2 The last diffusion term in (8) is from the intrinsic molecular noises in (7). If the gene network is free of intrinsic molecular fluctuations, then the asymptotic stability of the equilibrium point is only to check the existence of V(x(t)) > 0 in the following inequality without diffusion terms
Proposition 2 For the linear stochastic gene network with process delays in (6), we could choose the Lyapunov function as [35]
for P = P^{T }> 0 and > 0, i = 1,...,m. Then the robust stability problem in Proposition 1 is reduced to check whether or not the existence of symmetric matrices P > 0 and Q_{k }> 0, k = 1,...,m, in the following LMI
Proof. See Appendix B. ■
If the linear gene network in (6) is free of intrinsic molecular fluctuations, then the asymptotic stability is only to check whether or not P > 0 and Q_{k }> 0, k = 1,..., m, exist in the following LMI
Obviously, the condition for the existence of P > 0 and Q_{k }> 0, k = 1,...,m, in inequality (10) is stricter than that in inequality (11) because A_{0 }should be more negative than that in inequality (11), i.e., the eigenvalues of A_{0 }should be in the farther lefthand side of the complex domain (the real parts of all eigenvalues of A_{0 }should be more negative) or more robustly stable in order to tolerate time delays and intrinsic molecular fluctuations.
In general, a gene network in vivo also suffers from extrinsic molecular noises such as the transmitted noise from upstream genes and the global noise affecting all genes. Therefore, the nonlinear stochastic gene regulatory equation with process delays in (7) should be modified in the following to mimic the realistic dynamic behavior in its host cell.
where signal vector denotes the extrinsic molecular noises and G denotes the coupling matrix between the extrinsic molecular noises and the gene network. In general, the intrinsic molecular fluctuations are statedependent as shown in (6)–(7). Therefore, they will influence the stability of gene network and provide an important topic for discussing the robust stability of gene network at an equilibrium point of interest, i.e. to discuss how much the gene network could tolerate the intrinsic molecular fluctuation without blowing away from the equilibrium point. Another important topic is how to quantify the ability to filter extrinsic molecular noises v(t) in the gene network, i.e. the attenuation analysis of v(t) on x(t) by the gene network. Such information is crucial for understanding the signalprocessing scheme in gene networks, from which the design of molecular noisetolerant and time delaycompensation gene circuits can be developed to achieve a robust H_{∞ }noise filtering design for biotechnological purposes.
Therefore, we want to investigate the influence of extrinsic molecular noises v(t) on the gene i, i.e. x_{i}(t), which is crucial to the design of noisetolerant gene circuit in synthetic biology [7]. Let us denote the choice of a gene of interest as
where the row vector C = [0, 0,...0, 1, 0, 0, 0, 0], i.e., all elements of C are zero except the ith element being 1, if the ith gene x_{i}(t) is the gene of interest in molecular noise filtering. Let us denote the effect of extrinsic molecular noises v(t) on a gene of interest as follows
or
for all bounded energy noises v(t). This is the socalled H_{∞ }noise filtering problem [27], which addresses the effect of extrinsic molecular noises on the gene of interest from the system gain (L_{2}gain) point of view. If the extrinsic molecular noises v(t) are deterministic signals, then the expectation E on v(t) in (14) could be neglected.
Remark 3
(i) The physical meaning of (14) is that the effect of extrinsic molecular noises v(t) on z_{o}(t) should be less than or equal to γ. If γ < 1, then extrinsic molecular noises v(t) are filtered at gene i by the gene network. If γ > 1, then it means that extrinsic molecular noises v(t) are amplified at gene i, i.e., gene i is much influenced by external molecular noises v(t) and may be a weak site of the gene network.
(ii) In the conventional control system design and signal processing [2732], the H_{∞ }noise filtering criterion in (14) has been employed to design a filter to robustly estimate the state variables from the noisy measurement output, i.e., to specify the filter design parameters to efficiently attenuate the effect of uncertain external disturbance on the state estimation error of the filter [27,28]. In this study, we only discuss the effect of external disturbance on the genes of nonlinear timedelayed gene network and then engineer some gene circuits to improve the noise filtering ability of the nonlinear timedelayed gene network to attenuate the effect of external disturbance. Therefore, there are some differences between the proposed noise filtering circuit design in gene network and the conventional H_{∞ }filter design to robustly estimate state variables in a signal processing system under external disturbance.
(iii) If we want to discuss the effect of extrinsic molecular noises v(t) on the whole gene network, then C in (13) should be an identity matrix, i.e. C = I and z_{o}(t) = x(t).
(iv) If the initial condition x(0) ≠ 0, then the filtering ability in (14) should be modified as follows [27,28]
for some Lyapunov function V (x(0)), i.e., the energy due to initial condition x(0) should be considered in the H_{∞ }filtering performance.
Then we get the following molecular noise filtering result for the nonlinear stochastic gene network with process delays, intrinsic molecular fluctuations, and extrinsic molecular noises in (12).
Proposition 3 If there exists a positive function V (x(t)) > 0 with V (0) = 0 solving the following HJI
then the influence of molecular noises v(t) on the gene network is less than γ or equivalently, the molecular noise filtering ability γ in (14) or (15) is achieved in the nonlinear stochastic gene network under process delays, extrinsic molecular noises, and intrinsic molecular fluctuations.
Proof. See Appendix C. ■
Therefore, the optimal extrinsic molecular noise filtering ability of gene network at the gene of interest is obtained by solving the following constrained optimization
i.e., the molecular noise filtering ability γ_{0 }of v(t) on gene i could be evaluated by the optimal H_{∞ }filtering ability of gene network via solving the constrained optimization in (17).
Remark 4
(i) From Propositions 1 and 3, we gain more insight into robust stability and molecular noise filtering of nonlinear stochastic timedelayed gene network from the nonlinear stochastic H_{∞ }filtering perspective. The solution for HJI in (16) is stricter than that for HJI in (8), because of two extra terms for the robust filtering of extrinsic molecular noises in (16). i.e., Proposition 3 not only guarantees robust stability against process delay and intrinsic molecular fluctuations, but also achieves a prescribed noise filtering ability γ against extrinsic molecular noises.
(ii) In general, a nonlinear stochastic gene network has many equilibrium points (i.e. phenotypes). We are only concerned about the robust stability of one equilibrium point of interest, i.e. the equilibrium point at the steady state of gene networks. Propositions 1 and 3 are true only at the equilibrium point (one phenotype) of the origin, i.e. x_{e }= 0. If we are interested in an equilibrium point x_{e }≠ 0 (i.e. another phenotype), for the convenience of investigation, the origin should be shifted to x_{e}, i.e., the stochastic system in (12) should be modified as
where x'(t) = x(t)  x_{e}. Therefore, when we are interested in an equilibrium point x_{e }≠ 0 of gene network, the variable x(t) in HJIs (8) and (16) should be replaced by x'(t) + x_{e }so that the origin of the coordinate should be shifted to x_{e}. This shift is a onetoone correspondence [43]. Therefore, instead of studying the behavior of gene network in the neighborhood of x_{e}, one can equivalently study the behavior in the neighborhood of the origin x'(t) = 0 in (18).
In general, it is still very difficult to solve the secondorder HJI in (8) or (16) with V (x(t)) > 0 and V (0) = 0 to guarantee the robust stability of a nonlinear stochastic gene network under process delay and intrinsic molecular fluctuations or to solve the constrained minimization in (17) to get the noise filtering ability γ_{0 }on extrinsic molecular noises at any gene of the gene network, especially for complex gene networks. In this situation, the TS fuzzy model is employed to interpolate several linear gene regulatory networks at different operation points to efficiently and globally approximate the stochastic gene networks in (7) and (12). According to the TS fuzzy model, both the robust stability of a gene network to tolerate process delays and intrinsic molecular fluctuations and the filtering ability of attenuating extrinsic molecular noises of the gene network could be efficiently investigated and calculated via fuzzy interpolation of several local linear gene networks. Finally, if the gene regulatory network cannot tolerate process delays and intrinsic parameter variations and achieve a prescribed noise filtering level γ <γ_{0 }for some biotechnological or drug design purpose, a robust gene circuit control design is developed from the fuzzy robust stabilization and H_{∞ }filtering theory to achieve the robust stability and to improve the filtering ability of the gene network. The proposed robust H_{∞ }filtering design method for gene networks is also useful for synthetic biology in the near future, if a synthetic gene network wants to perform reliably under process delays, intrinsic molecular fluctuations, and extrinsic molecular noises.
Robust H_{∞ }Gene Circuit Design of Nonlinear Gene Network Under Process Delays and Molecular Noises via Fuzzy Methodology
If the linear stochastic gene network in (6) is not robust enough and cannot tolerate the molecular fluctuations , which are statedependent and will affect the stability of gene network. i.e., the perturbative gene network in equation (6) becomes unstable, then a gene circuit design is necessary for the gene network. We want to engineer some feedback gene circuits to robustly stabilize the perturbative gene network as follows
where , are the circuit kinetic parameter matrices of gene circuits to be designed. The element k_{0,ij }of K_{0 }denotes the circuit parameter to be specified for the engineered gene circuit between gene j and gene i via transfection and transformation biotechnologies (see Fig. 1). The gene circuit from gene j to gene i can be implemented by inserting the binding site of gene product j (i.e. the protein of gene j) into the promoter region of gene i so that the protein of gene j could bind this inserted binding site to act as a transcription factor (TF) to regulate the gene expression of gene i. By inserting strong (long) or weak (short) binding site, we can get a large or small circuit kinetic parameter k_{0,ij}. The inserting of a TF binding site into promoter region or the deleting of a TF binding site from the promoter region can be easily done by using a highly efficient phagebased homologous recombination system, called recombineering [37,38]. Furthermore, the change of decay rates in the diagonal terms of A_{0 }by the diagonal terms of K_{0 }can be implemented by the mechanisms and controls of mRNA degradation through elongating or shortening the 3' polyadenylate tail of mRNA [39,40,4446]. These powerful biotechnologies have been employed to engineer large segments of genomic DNA to generate transgenic and knockout constructs. However, for simplicity and feasibility, only a few feasible gene circuits in K_{0}x(t) are considered for gene circuit design in gene network, i.e., the circuit kinetic parameter matrix K_{0 }has with only a few elements k_{0,ij}. The design principle of K_{k}, k = 1,...,m, is similar. A more detailed design procedure will be given in the design example in the sequel. Therefore, in the gene circuit design procedure, we engineer some feasible gene circuits k_{k,ij}x_{j}(t  τ_{k}) for a gene network with adequate k_{k,ij }to achieve robust stabilization and filtering design. In the situation, the specification of k_{k,ij }for the embedded gene circuits is different from the convenient control or filter design method, which is used to calculate control input or state estimation from state variables or output signals. However, the closedform solutions for K_{k}, k = 0, 1,...,m, are diffcult to be obtained. Therefore, some searching algorithms for components of K_{k}, k = 0, 1,...,m, are needed.
Following the robust stability inequality in (10), we could get the following result.
Proposition 4 The designed gene network in (20) is robustly stable if we could specify K_{k}, k = 0, 1,...,m, for the designed gene circuits such that there exist the positive definite symmetric matrices P > 0 and Q_{k }> 0, k = 1,...,m, solving the following inequality
where and = A_{k }+ K_{k}, i.e., the robust gene circuit design becomes how to specify K_{k}, k = 0, 1,...,m, such that there exist the positive symmetric matrices P and Q_{k}, k = 1,...,m, in LMI (21).
Proof. The proof is trivial. ■
Remark 5 In order to make sure the positive concentrations of gene regulation network in (19) or (20), the choice of designed kinetic parameter K_{k }to solve LMI in (21) should guarantee the positive concentrations of gene networks. In general, according to the positive orthant stabilization theory [35], the guaranty for the positive concentration of x(t) is that the offdiagonal entries of = A_{k }+ K_{k }are all nonnegative. Therefore, the specification of K_{k }to solve LMI should be subjected to this constraint. If these concentration constraint can not be achieved, changes of gene circuits through other loop should be tried.
If the nonlinear stochastic timedelayed gene network in equation (7) cannot tolerate molecular fluctuations , then we want to design the gene circuits in the following to robustly stabilize the perturbative gene network
where g_{k}(x(t  τ_{k})), k = 0, 1,...,m, denote the nonlinear feedback circuits and the matrices K_{k}, k = 0, 1,...,m, denote their corresponding circuit kinetic parameters to be specified to robustly stabilize the perturbative gene network. In general, for the nonlinear stochastic gene network with time delays in (22), it is not easy to obtain a systematic design method for K_{k}, k = 0,1,...,m, to achieve the robust gene circuit design.
The fuzzy dynamic model with time delay has been widely employed to interpolate several local linear timedelayed dynamic models to efficiently approximate a nonlinear timedelayed dynamic system. The TS fuzzy dynamic model is described by fuzzy IfThen rules and employed here to deal efficiently with the robust H_{∞ }filtering problem for gene networks under process delay, intrinsic molecular fluctuations and extrinsic molecular noises in (22). The ith rule of the fuzzy model for nonlinear stochastic timedelayed gene network in (22) is proposed as the following form [33,34]:
for i = 1, 2,...,L. F_{ij }is the fuzzy set; A_{k,i}, B_{k,i}, and G_{k,i }(k = 0, 1,...,m) are known constant matrices; L is the number of IfThen rules; z(t) = [z_{1}(t),...,z_{g}(t)]^{T }are the premise variables; g is the number of premise variables. If all state variables x(t) are used as premise variables, then z(t) = x(t) and g = N. The physical meaning of the fuzzy rule i is that if premise variables z_{1}(t), z_{2}(t),...,z_{g}(t) are with the fuzzy sets F_{i1}, F_{i2},...,F_{ig}, then the nonlinear stochastic gene network in (22) can be represented by interpolating the linearized systems in (23). The fuzzy stochastic system in (23) is inferred as follows [33,34]
where
F_{ij}(z_{j}(t)) is the grade of membership of z_{j}(t) in F_{ij }or the possibility function of z_{j}(t) in F_{ij }and μ_{i}(z), for i = 1, 2,...,L, are called the fuzzy bases. The denominator is only for normalization so that the total sum of fuzzy bases . The physical meaning of (24) is that the fuzzy stochastic system interpolates L local linear stochastic systems through nonlinear bases μ_{i}(z) to approximate the nonlinear stochastic gene network in (22). In this situation, the Langevin stochastic equation in (22) can be represented by the fuzzy interpolatory gene network as follows
where ; Δf_{k}(x), Δg_{k}(x), and Δh_{k}(x) denote fuzzy approximation errors as follows
There are many methods for finding A_{k,i }and B_{k,i}, k = 0, 1,...,m, i = 1, 2,...,L, for fuzzy model identification [33]. We could use fuzzy toolbox in Matlab to find A_{k,i}, B_{k,i}, and G_{k,i }easily. After finding A_{k,i}, B_{k,i}, and G_{k,i}, k = 0, 1,...,m, i = 1, 2,...,L, we could easily find the bounds of fuzzy approximation errors Δf_{k}(x), Δh_{k}(x) and Δg_{k}(x) as follows
for some positive constants a_{k}, b_{k}, c_{k}, k = 0, 1,...,m, where , e.g. if k = 0 then with τ_{0 }= 0 in (7).
According to the above fuzzy approximation method, we get the following result.
Proposition 5 The nonlinear stochastic gene network with time delays in (22) is robustly stabilizable by gene circuits K_{k}g_{k}(x(t  τ_{k})), k = 0, 1,...,m, if there exist the positive definite symmetric matrices P > 0 and Q_{k }> 0, k = 1,...,m, solving the following LMIs
for all i = 1,⋯,L, and P <βI, where
Therefore, the robust stabilization problem for the nonlinear stochastic timedelayed gene network in (22) becomes how to specify K_{k}, k = 0, 1,...,m, in LMIs (27) such that there exist the positive definite symmetric common matrices P > 0 and Q_{k }> 0, k = 1,...,m, i.e., the gene circuits K_{k}g_{k}(x(t  τ_{k})), k = 0, 1,...,m, could robustly stabilize the nonlinear stochastic timedelayed gene network, and the equilibrium point x = 0 of gene network (22) is asymptotically stable in probability under process delays and intrinsic molecular fluctuations.
Proof. See Appendix D. ■
Remark 6 In the conventional fuzzy design [33,34], we need to design a complicated fuzzy controller or fuzzy filter for engineering systems. In this paper, the fuzzy approximation in (25) is only to simplify the gene circuit design procedure by solving K_{k}, k = 0, 1,...,m, via LMIs in (27) instead of solving K_{k}, k = 0, 1,...,m, via HJI directly. So there are some differences between the proposed fuzzy gene circuit design method and the conventional fuzzy control and filter design methods.
Remark 7 For the nonlinear stochastic gene network with time delays, it is diffcult to construct a Lyapunov function V (x(t)) to satisfy the HJI in (8) or (16). However, based on fuzzy approximation method, the nonlinear stochastic gene network with time delays is approximated by interpolating several local linear systems with time delays at different operation points. The advantage of the proposed fuzzy approach is that the circuit design procedure could be simplified from the linear timedelayed system design point of view. Thus the following standard Lyapunov function for linear system with multiple timedelays can be employed for the robust stabilization of nonlinear stochastic gene network [32,35]
Therefore, we can avoid the diffculty of constructing a complex Lyapunov function for (16) and can employ the above Lyapunov function to obtain a systematic design procedure in the molecular circuit design.
Remark 8 If the gene circuit terms cannot be separated as (22) in gene circuit design and are merged into the terms f_{k}(x(t  τ_{k})), k = 0, 1,...,m, as follows
then (27) should be changed as follows
for all i = 1,⋯,L, and P <βI, where
and other matrices are same as in (27). A_{k,i}(K_{k}), k = 0, 1,...,m, denote the system matrices A_{k,i}, k = 0, 1,...,m, containing the designed circuit parameters K_{k}, k = 0, 1,...,m, as elements, respectively.
After investigating the robust stabilization design of nonlinear stochastic gene network under process delays and intrinsic molecular fluctuations by the fuzzy approximation method, in order to avoid solving the nonlinear constrained optimization for molecular noise filtering problem in (17), the extrinsic molecular noise filtering ability improvement problem of gene network could be also treated by the fuzzy interpolation method. Therefore, the robust filtering circuit design problem becomes how to engineer the gene circuits K_{k}g_{k}(x(t  τ_{k})), k = 0, 1,..., m, not only to tolerate intrinsic molecular fluctuations but also to attenuate the effect of extrinsic molecular noises v(t) on z_{o}(t) to a prescribed noise filtering level γ as (14) or (15). Especially, when the prescribed noise filtering level γ is below the optimal filtering ability γ_{0 }in (17) of the gene network, the robust filtering circuit design is needed to improve the molecular noise filtering ability. Since it is not easy to specify K_{k}, k = 0, 1,...,m, to achieve robust filtering design directly for the nonlinear stochastic gene network in (22), the TS fuzzy approximation in the following equation is employed to treat the robust filtering design problem as follows
Then we get the following result.
Proposition 6 For the stochastic gene network with time delay in (29), if we could specify K_{k}, i = 1,...,m for the designed gene circuits such that there exist the positive definite symmetric matrices P and Q_{k}, k = 1,...,m, solving the following LMIs
for i = 1, 2,⋯,L, and P <βI where
and the other matrices are defined in (27), then the effect of extrinsic molecular noises v(t) on the gene of interest is less than γ.
Proof. See Appendix E. ■
Therefore, according to fuzzy approximation method, the optimal H_{∞ }noise filtering design problem for nonlinear stochastic timedelayed gene network in (29) could be solved by the following constrained optimization problem
After solving K_{k}, k = 0, 1,...,m, from LMIs in (31) for a prescribed noise filtering level γ or solving K_{k}, k = 0, 1,...,m, from (32) for the optimal H_{∞ }noise filtering ability γ_{0}, we can design gene circuits K_{k}g_{k}(x(t  τ_{k})), k = 0, 1,...,m, in the nonlinear stochastic timedelayed gene network (22) to achieve a desired molecular noise filtering ability or an optimal H_{∞ }noise filtering ability of gene network, respectively. Since the gene circuit design could improve the noise filtering ability, the γ_{0} in (32) should be less than the γ_{O }in (17) without gene circuit design.
Remark 9 In general, the optimization problem in (32) is called the eigenvalues problem, which can be efficiently solved by the Matlab LMI toolbox [35].
Remark 10 If the gene control circuit terms cannot be separated as (22) and are merged into f_{k}(x(t  τ_{k})), k = 0, 1,...,m, as follows
then LMIs in (31) should be modified as follows
for i = 1, 2,⋯,L, and P <βI, where
and other matrices are same as in (27) and (31). A_{k,i}(K_{k}), k = 0, 1,...,m, denote the matrices A_{k,i }whose entries contain circuit kinetic parameters K_{k}. In this situation, for a prescribed noise filtering ability γ, the robust circuit design is to specify K_{k}, k = 0, 1,...,m, such that there exist P > 0 and Q_{k }> 0 for LMIs in (35).
Remark 11 If the prescribed attenuation level γ is small in H_{∞ }noise filtering design, more extrinsic molecular noises are eliminated by gene network. However, it is more diffcult to solve LMIs in (31) or (35) and more control or filtering effort is needed. This is a trade off for a designer between a good filtering ability (small γ) and a design diffculty in the robust gene circuit design.
Remark 12 Unlike the conventional fuzzy control designs [30,33,34,47] and fuzzy filter designs [20,27,28], the proposed fuzzy interpolation method is only to simplify the circuit design procedure so that we could solve K_{k}, k = 0, 1,...,m, easier via the LMI scheme.
According to the above analysis, a robust circuit design procedure using fuzzy interpolation method for stochastic gene regulatory networks are summarized as follows
Design Procedure
Step 1: Model the nonlinear timedelayed gene network as (12) and shift the interested equilibrium point to origin as (18).
Step 2: Design some feasible feedback control circuits K_{0}g_{0}(x(t)) and K_{k}g_{k}(x(t  τ_{k})), k = 1,...,m.
Step 3: Approximate the nonlinear stochastic gene network by fuzzy system as (29).
Step 4: Estimate the fuzzy approximation error bounds a_{k}, b_{k}, and c_{k}, k = 0,1,...,m, in (26), and give a prescribed noise filtering ability γ.
Step 5: Specify the feasible circuit kinetic parameters K_{k}, k = 0,1,...,m, by solving LMIs in (31) for the prescribed noise filtering ability γ or solving the constrained optimization in (32) for the optimal noise filtering ability γ_{0}.
Computational Design Examples in Silico
Two computational design examples using the developed robust gene circuit design method are presented in silico to illustrate the design procedure and to validate the performance of the proposed molecular circuit design methods.
Example 1:
Consider a benchmark genetic regulatory network as shown in Fig. 2, which is a typical gene regulatory network describing the gene, mRNA and protein interactions [36]. These genes are regulated by other genes and then expressed through transcription and translation to obtain their products, i.e. proteins. Then these proteins could be as the TFs of other genes to regulate the expressions of other genes after some process delays. We consider only the mRNA abundances x_{1}, x_{2}, x_{3}, and x_{4}. Suppose the gene network suffers from stochastic intrinsic molecular fluctuations and extrinsic molecular noises v(t) on gene x_{2 }and gene x_{4}. The stochastic gene network under process delays, intrinsic molecular fluctuations and extrinsic molecular noises can be represented as follows [36]
Here λ_{1}, λ_{2}, λ_{3 }and λ_{4 }are the firstorder rate constants of the degradation of x_{1}, x_{2}, x_{3 }and x_{4}, respectively. τ_{k }denotes the process delay of gene regulation due to transcription, transfection, posttranslation modification and transportation of TF. The Hill term describes the sigmoid formation of x_{2 }activated by x_{4 }with time delay τ_{1}, maximal rate V_{2}, dissociation constant Λ_{2 }and Hill coefficient n_{4}. The inhibition by x_{3 }is expressed by the term (Λ_{I3 }+ (t  τ_{3})). The formation of x_{3 }is modeled with Hill expression that points to a threshold of the formation of x_{3 }depending on the concentrations of x_{1 }and x_{2}. V_{3 }and Λ_{3 }are maximal rate and dissociation constant, respectively, and n_{12 }is the Hill coefficient. The production of x_{4 }depends on the maximal rate V_{4 }and on the inhibition by x_{3}. The parameters are chosen as follows [36]
The gene expression profiles of stochastic gene network with process delays, intrinsic and extrinsic molecular noises in (36)–(39) are given in Fig. 3. It is seen that there are large fluctuations in these gene expression profiles due to process delays and molecular noises. Suppose the proposed gene circuit design method is employed to attenuate these molecular fluctuations to achieve a desired steady state quickly. In general, engineering a gene circuit to a gene network may change its steady state because of nonlinear inherence. Since changing the steady state of a gene network may destroy the normal function of the gene network, it is important to make sure that a designed gene circuit does not change the steady state (phenotype) of a gene network while it can compensate process delays and attenuate molecular noises. Therefore, in order to keep the steady state of a designed gene network the same as the desired steady state, the gene network is designed as follows
Figure 3. The gene expression profiles of nonlinear stochastic timedelayed gene network without engineered gene circuit.
where k_{1}, k_{2 }are the circuit kinetic parameters at the genes x_{2 }and x_{3}. From the simulation result of the nominal gene regulatory network, we can know that the equilibrium point of the nominal system is at [x_{e1}, x_{e2}, x_{e3}, x_{e4}] = [1.0000, 0.5899, 1.0557, 0.5739] which is the desired steady state (phenotype) and should be shifted to the origin as shown in (18) in the design procedure. Suppose we want to specify the engineered circuit parameters k_{1 }and k_{2 }such that the prescribed noise filtering ability γ = 0.9 can be achieved for the molecular noise filtering of perturbative gene regulatory network in (40)–(43). We must find the operative points of the gene network to construct the fuzzy model at first. The operative points of the state x_{1 }are located at (0.1, 1). For the other states x_{2}, x_{3}, and x_{4}, the operative points are located at (0, 1), (0.5, 1.5), and (0, 1), respectively. Then we can create two membership functions for every state at the operative points, and the number of fuzzy rules is L = 16. The bounds of the fuzzy approximation errors are estimated as a_{0 }= 0, a_{1 }= 0.001, a_{2 }= 0.0524, b_{0 }= 0, b_{1 }= 0, b_{2 }= 0, c_{0 }= 0, c_{1 }= 0, c_{2 }= 0 at (26). By the robust molecular noise filtering circuit design in Proposition 4, in order to guarantee the positive definite symmetric matrices of P, Q_{1 }and Q_{2}, the engineered circuit kinetic parameters k_{1 }and k_{2 }should be in the range k_{1 }∈ (0, 1] and k_{2 }∈ [1, 20]. In the case k_{1 }= 0.1, k_{2 }= 10, we get
The molecular noise filtering results without engineered circuit and with two engineered circuit k_{1 }= 0.1 and k_{2 }= 10 are simulated in Fig. 3 and Fig. 4, respectively. Obviously, the large molecular fluctuations in Fig. 3 are significantly attenuated by the proposed gene circuits as shown in Fig. 4 and the expression profiles of four genes asymptotically achieve the desired states in probability. From the simulation, the gene network without engineered gene circuit design has much molecular fluctuations due to process delays and molecular noises. However, the designed gene network by the proposed gene circuit design can robustly achieve the desired equilibrium point x_{e }= [1.0000, 0.5899, 1.0557, 0.5739] and the molecular noise filtering is also improved by the proposed gene circuit design. For confirmation, the noise filtering ability of the designed gene network is estimated as follows
Figure 4. The gene expression profiles of nonlinear stochastic timedelayed gene network with engineered gene circuit. Four step functions in figures denote the desired steady states of gene expressions x_{1}x_{4}.
and the filtering ability of stochastic gene network in (36)–(39) without gene circuit design is given by
where z_{o}(t) = Cx(t) in which C is an identity matrix. Therefore, the prescribed filtering ability γ is achieved by the proposed robust gene circuit design. From the design example in silico, it is obvious that the proposed gene circuit design could improve the molecular noise filtering ability of the nonlinear stochastic gene network.
Recently, the gene circuit design can be implemented by using a highly efficient phagebased homologous recombination system, called recombineering [37,38]. This powerful technology has been used to engineer large segments of genomic DNA to generate transgenic and knockout construct to insert or delete the TF binding sites in the promoter region of regulated genes to increase or decrease the expression of regulated genes. Therefore, k_{1 }= 0.1 of the first term in the right hand side of (41) could be achieved by deleting 90% of TF binding sites of gene product X_{4 }of gene x_{4 }from the promoter region of gene x_{2 }through knockout construct technique of recombineering. Similarly, k_{2 }= 10 of first term in the right hand side of (42) could be implemented by inserting 10 times of TF binding sites of complex protein X_{1}X_{2 }of gene x_{1 }and x_{2 }to the promoter region of gene x_{3 }to increase its gene expression.
As for the implementation of two mRNA decay terms k_{1}λ_{2}x_{2 }and k_{2}λ_{3}x_{3 }in (41) and (42), respectively, it has been shown that mRNA decay in eukaryotic cells can be achieved by shortening the 3' polyadenylate tail found on eukaryotic mRNAs (referred to as deadenylation), which primarily triggers decapping, leading to 5' to 3' exonucleolysis. Alternatively, removal of 3' polyadenylate tail can expose the mRNA to 3' to 5' degradation [39,40,4446]. Therefore, by elongating the 3' polyadenylate tail of mRNA of gene x_{2}, we can get a small kinetic decay parameter k_{1 }of k_{1}λ_{2}x_{2 }in (41) and by shortening the 3' polyadenylate tail of gene x_{3}, we can get a large kinetic decay parameter k_{2 }of k_{2}λ_{3}x_{3 }in (42).
Example 2:
Consider a dynamic system for the regulation of induction in the lac operon [48]. This stochastic system consists of five nonlinear differential equations with multiple time delays due to transcription and translation processes (see Fig. 5). Therefore, this gene regulatory system in vivo could be described by the following nonlinear stochastic timedelayed system [48].
Figure 5. The schematic diagram of the lactose operon regulatory system.
where x_{1}(t) is the dynamic of mRNA production; x_{2}(t) is the dynamic of βgalactosidase; the dynamic of allolactose is described by x_{3}(t); the lactose dynamic is described by x_{4}(t); the permease dynamic is described by x_{5}(t). When the the glucose available for cellular metabolism is absent, the external lactose L_{e }is transported to the cell by the permease x_{5}(t). Then, by the enzyme βgalactosidase x_{2}(t), the intracellular lactose x_{4}(t) is broken down into glucose, galactose, and allolactose x_{3}(t). Finally, the allolactose feeds back to bind with the lactose repressor and enables the transcription process to produce the mRNA production x_{1}(t). The mRNA production x_{1}(t) from DNA via transcription needs a delayed time τ_{M }for RNA polymerase to transverse the three structural genes, and the βgalactosidase production x_{2}(t) through mRNA translation requires a delayed time τ_{B}. The delayed time τ_{P }is the translation time between mRNA and permease. The delay τ_{P }is the sum of the βgalactosidase and premise translation times based on the assumption that permease production can not start until βgalactosidase production is complete. Δ_{i}, i = 1,...,5 denote the corresponding parameter variations. The detailed process for the lactose operon regulatory system is described in Fig. 5 and refers to the literature [48]. The parameters for the model are given in Table 1. The initial values are chosen as x_{1}(0) = 6.26 × 10^{4}, x_{2}(0) = 0, x_{3}(0) = 3.80 × 10^{1}, x_{4}(0) = 3.72 × 10^{1}, and x_{5}(0) = 1.49 × 10^{2}. The desired steady state (the interested equilibrium point) is at x_{e}(t) = [1.2359 × 10^{3}, 8.3645 × 10^{4}, 5.7666 × 10^{1}, 4.1583 × 10^{1}, 1.1664 × 10^{2}]^{T}.
Table 1. The parameters for the lactose operon regulatory system.
If the gene regulatory system is free of intrinsic molecular fluctuations and extrinsic disturbances, i.e., W(t) ≡ 0 and v(t) ≡ 0, the timeprofiles of regulatory system are shown in Fig. 6. We can know the system can work right. However, in the realistic cellular environment, this gene regulatory system may suffer from the intrinsic molecular fluctuations and extrinsic disturbances. In this example, an extrinsic disturbance v(t) affects the lactose dynamic x_{4}(t), and the allolactose dynamic x_{3}(t) and the lactose dynamic x_{4}(t) also suffer from some molecular fluctuations. By computational simulation, the time profiles of the stochastic regulatory system are shown in Fig. 7. Obviously, this system can not work properly. Therefore, we need to redesign the some parameters of some gene circuits to ensure the system can resist the influences of intrinsic molecular fluctuations and extrinsic disturbances. From the lactose dynamic x_{4}(t) which is critically destroyed, we choose parameters L_{e}, α_{L}, , and γ_{L }to redesign the gene regulatory system.
Figure 6. The timeprofiles of the lactose operon regulatory system in example 2 which is free of the intrinsic molecular fluctuations and extrinsic disturbances.
Figure 7. The timeprofiles of the stochastic lactose operon regulatory system in example 2 without circuit design.
By the same procedure as example 1, the we could choose the three membership functions for every state x_{i}(t) to approximate nonlinear functions of the gene network. Therefore, we have L = 3^{5 }= 243 fuzzy rules to approximate the gene regulatory system. The approximation errors in (26) could be estimated as a_{0 }= 4.2355 × 10^{3}, a_{1 }= 1.2136 × 10^{9}, a_{2 }= 4.8344 × 10^{4}, a_{3 }= 7.7390 × 10^{17}, b_{0 }= 2.7578 × 10^{7}, b_{1 }= 0, b_{2 }= 0, and b_{3 }= 0. Note that the parameters of the fuzzy system are omitted, because the amount of the parameters is huge. Using LMI toolbox in Matlab to solve the optimization problem in (32), we can obtain the optimal H_{∞ }filtering ability γ_{0 }= 0.3535 and the positive definite matrices P, Q_{k}, k = 1...3 in the following.
The design parameters of this regulation system are obtained as L_{e }= 8.0 × 10^{2}, α_{L }= 2880.3, = 8460, and γ_{L }= 0.002.
From the simulation results in Fig. 7 and Fig. 8, the molecular noise filtering ability is also improved significantly by the proposed gene circuit design. For confirmation, the filtering ability of stochastic regulatory systems of lac operon without gene circuit design is given by
Figure 8. The timeprofiles of the stochastic lactose operon regulatory system in example 2 with circuit design.
and the noise filtering ability of the designed regulatory system is estimated as follows
The conservative result of noise filtering ability is mainly due to the comservative procedure in solving LMIs [35].
Conclusion
In the study, a gene network with process delays, intrinsic molecular fluctuations and extrinsic molecular noises is modeled as a nonlinear stochastic timedelayed system. Then we propose a stochastic gene circuit design method for the improvement of robust stability and molecular noise filtering ability of nonlinear gene network to tolerate timedelays and to attenuate molecular noises via LMI technique and fuzzy interpolation scheme. The TS fuzzy system can approach the nonlinear stochastic gene network with time delays via the interpolation of several local linear timedelayed systems to avoid solving HJI in nonlinear robust stabilization and filtering design problems. Therefore, the robust circuit design procedure for the nonlinear stochastic timedelay gene network can be simplified by specifying circuit kinetic parameters to satisfy a set of LMIs, which could be efficiently solved by the LMI toolbox in Matlab. Unlike the conventional trialanderror, the proposed design method provides a systematic method for robust gene circuit design of nonlinear stochastic timedelayed gene networks. Because the microarray data become popular, the construction of a dynamic model from microarray data for a gene regulatory network becomes possible. Furthermore, the experimental advances in transfection and transformation biotechnologies make gene circuit implementation easier. Therefore, the proposed gene circuit design methods have much potential for application to systems biology, synthetic biology and drug design when a gene regulatory network has to be designed for improving its robust stability and filtering ability of diseaseperturbed gene networks or when a synthetic gene network needs to perform reliably around a desired equilibrium point despite of process delays, intrinsic molecular fluctuations and extrinsic molecular noises in host cell. Finally, a benchmark design example is also given in silico to illustrate the design procedure and to validate the proposed robust circuit design method in nonlinear stochastic timedelayed gene regulatory networks.
In the future, we will focus on the development of some more general design methods for robust synthetic biologic networks under parameter fluctuations, timedelays and environmental molecular noise in host cell. After some design specifications are given beforehand, for example, the magnitudes of parameter variations to be tolerated, the possible delays to be compensated, the filtering ability to attenuate the molecular noises, the feasible ranges of kinetic parameters to be designed and the desired steady states to be achieved, we want to develop some systematic design methods for a synthetic biologic network to meet these design specifications and achieve the design objective.
Appendix
Before the proof of propositions, the following fact is necessary.
Fact 1 ([35]):
for vectors X, Y, a constant α > 0 and a positive definite matrix P > 0 with appropriate dimensions.
Appendix A Proof of Proposition 1
By the stochastic Lyapunov stability theorem, we choose a Lyapunov function V(x(t)) > 0 such that . Using Itô formula [49], we have
Then
If the inequality in (8) holds, then , i.e., the nonlinear stochastic gene network in (7) is asymptotically stable in probability at x = 0.
Appendix B Proof of Proposition 2
For the linear stochastic gene network in (6), if we choose
where P = P^{T }> 0 and , then
where = [x(t), x(t  τ_{1}),..., x(t  τ_{m})]^{T }and . If the LMI in (10) holds, then , i.e., the linear stochastic gene network in (6) is asymptotically stable in probability.
Appendix C Proof of Proposition 3
Consider the following equivalent equation
where V(x(t)) > 0.
By the Itô formula, we get [49]
Substituting the above equation into (C.1), by the fact that V(x(∞)) ≥ 0 and W_{k}(t) is a zero mean Wiener process and independent of x(t), we get
By Fact 1, we have
Therefore, we can obtain
By the inequality in (16), we get
Obviously, the filtering ability in (15) is achieved. Because the noise filtering ability is defined as the minimum effect of v(t) on the gene, so the noise filtering ability could be achieved by minimizing γ under the constraint of (16).
Appendix D Proof of Proposition 5
Let us choose a Lyapunov function for gene network (25) as
where P = P^{T }> 0 and . By Itô formula [49] and gene network (25), we get
By Fact 1, we have
and
for k = 0, 1,...,m, where α_{k, 1}, α_{k, 2}, and α_{k, 3}, are positive constants.
Suppose the following inequality holds
where β is a positive constant. By (D.2)–(D.6) and (26), we get
where = [x(t), x(tτ_{1}),⋯,x(tτ_{m})]^{T }and other matrices are defined in (27). If the following inequalities held
for i = 1,...,L, then . So the gene network (25) is asymptotically stable in probability.
Appendix E Proof of Proposition 6
Consider the following equivalent equation
where for some P = P^{T }> 0 and .
By the Itô differential equation [49] and (29), we get
Substituting the above equation into (E.1), by the fact that and W_{k}(t) is a zero mean Wiener process and independent of x(t), we get
By (D.2)–(D.6), (26) and (by Fact 1), we get
By the the following inequalities,
we get
Obviously, the noise filtering ability in (15) is achieved. Furthermore, by Shur complement [35], the inequalities in (E.5) are equivalent to the LMIs in (31).
Appendix F The design parameters of example 1
The parameters of the fuzzy model for the genetic regulatory network in (36)–(39) are listed in the following. Because the genetic regulatory network in (36)–(39) have two different time delays τ_{1 }= 1 and τ_{2 }= τ_{3 }= 2, we can separate the system to f_{0}(x(t)), f_{1}(x(t  τ_{1})), and f_{2}(x(t  τ_{2})). Therefore, we can obtain the the parameters of fuzzy model in (24) as follows
Because the function f(x(t)) is linear, the matrices A_{0,i }are the same.
Authors' contributions
BSC gives the topic and derives some results, and YTC gives some proofs of the results and performs the example simulations.
Acknowledgements
The research was supported by the National Science Council of Taiwan under Contract NSC 952221E007196MY3.
References

Rao CV, D MW, Arkin AP: Control, exploitation and tolerance of intracellular noise.
Nature 2002, 420(14):231237. PubMed Abstract  Publisher Full Text

Berg OG: A model for the statistical fluctuations of protein numbers in a microbial population.
J Theor Biol 1978, 71(4):587603. PubMed Abstract

McAdams HH, Arkin AP: Stochastic mechanisms in gene expression.
Proc Natl Acad Sci USA 1997, 94:814819. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hasty J, McMillen D, Isaacs F, Collins JJ, et al.: Computational studies of gene regulatory networks: in numero molecular biology.
Nat Rev Genet 2001, 2(4):268279. PubMed Abstract  Publisher Full Text

Thattai M, Oudenaarden AV: Intrinsic noise in gene regulatory networks.
Proc Natl Acad Sci USA 2001, 98:86148619. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kepler TS, Elston TC: Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations.
Biophys J 2001, 81:31163136. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hasty J, McMillen D, Collins JJ: Engineered gene circuits.
Nature 2002, 420(14):224230. PubMed Abstract  Publisher Full Text

Ozbudak EM, Thattai M, Kurtser I, Grossman AD, Oudenaarden AV: Regulation of noise in the expression of a single gene.
Nature Genet 2002, 31:6973. PubMed Abstract  Publisher Full Text

Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell.
Science 2002, 297:11831186. PubMed Abstract  Publisher Full Text

Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems.

Alon U: An Introduction to Systems Biology. Chapman & Hall/CRC, New York; 2007.

Andrianantoandro E, Basu S, Karig D, Weiss R: Synthetic biology: new engineering rules for an emerging discipline.
Mol Syst Biol 2006, 2:2006.0028. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Voit EO: Computational Analysis of Biochemical Systems: A Practical Guide for Biochemists and Molecular Biologists. Cambridge University Press, Cambridge, UK; 2000.

Smolen P, Baxter DA, Byrne JH: Modeling transcriptional control in gene networks methods, recent results, and future directions.
Bull Math Biol 2000, 62:247292. PubMed Abstract

Hasty J, Issacs F: Designer gene networks: Towards fundamental cellular control.
Chaos 2001, 11(1):207220. PubMed Abstract  Publisher Full Text

Yi TM, Huang Y, Simon MI, Doyle J: Robust perfect adaptation in bacterial chemotaxis through integral feedback control.
Proc Natl Acad Sci USA 2000, 97:46494653. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chen BS, Wang YC: On the attenuation and amplification of molecular noise in gene regulatory networks.
BMC Bioinformatics 2006, 7:114. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

McAdams HH, Arkin A: It is a noisy business! Genetic regulation at the nanomolar scale.
Trends Genet 1999, 15:6569. PubMed Abstract  Publisher Full Text

Hasty J, Pradines J, Dolnik M, Collins JJ: Noisebased switches and amplifiers for gene expression.
Proc Natl Acad Sci USA 2000, 97:20752080. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chen BS, Wu WS: Robust filtering circuit design for stochastic gene networks under intrinsic and extrinsic molecular noises.
Math Biosci 2008, 211(2):342355. PubMed Abstract  Publisher Full Text

Chen BS, Wu WS, Wang YC, Li WH: On the robust circuit design schemes of biochemical networks: steady state approach.
IEEE Trans Biomedical Circuits and Systems 2007, 1(2):91104.

Chen BS, Chang YT, Wang YC: Robust H_{∞}stabilization design in gene networks under stochastic molecular noises: fuzzyinterpolation approach.
IEEE Trans Syst Man Cybern B Cybern 2008, 38(1):2542. PubMed Abstract  Publisher Full Text

Lü J, Chen G: A timevarying complex dynamical network model and its controlled synchronization criteria.

Lü J, Yu X, Chen G, Cheng D: Characterizing the synchronizability of smallworld dynamical networks.

Yu W, Cao J, Lü J: Global synchronization of linearly hybrid coupled networks with timevarying delay.

Zhang Q, Lu J, Lü J, Tse CK: Adaptive feedback synchronization of a general complex dynamical network with delayed nodes.

Zhang W, Chen BS, Tseng CS: Robust H_{∞ }filtering for nonlinear stochastic systems.

Chen BS, Tsai CL, Chen YF: Mixed H _{2}/H_{∞ }filtering design in multirate transmultiplexer systems: LMI approach.

Chen BS, Zhang W: Stochastic H _{2}/H_{∞ }control with statedependent noise.

Zhang W, Chen BS: H_{∞ }control for nonlinear stochastic systems.

Feng J, Zhang W, Chen BS: H_{∞ }control for a class of nonlinear stochastic timedelay systems.

Zhang H, Wang Y, Liu D: Delaydependent guaranteed cost control for uncertain stochastic fuzzy systems with multiple time delays.
IEEE Trans Syst Man Cybern B Cybern 2008, 38(1):126140. PubMed Abstract  Publisher Full Text

Passino KM, Yurkovich S: Fuzzy control. AddisonWesley, Menlo Park; 1998.

Chen BS, Tseng CS, Uang HJ: Robustness design of nonlinear dynamic systems via fuzzy linear control.

Boyd SP, Ghaoui LE, Feron E, Balakrishnan V: Linear matrix inequalities in system and control theory. SIAM; 1994.

Klipp E, Herwig R, Kowald A, Wierling C, Lehrach H: Systems Biology in Practice: Concepts, Implementation and Application. WileyVCH; 2005.

Coplland N, Jenkins N, Court D: Recombineering: a powerful new tool for mouse functional genomics.
Nat Rev Genet 2001, 2:76979. PubMed Abstract  Publisher Full Text

Court D, Sawitzke J, Thomason L: Genetic engineering using homologous recombination.
Ann Rev Genet 2002, 36:36188. PubMed Abstract  Publisher Full Text

Decker CJ, Parker R: A turnover pathway for stable and unstable mRNAs in yeast: evidence for a requirement for deadenylation.
Genes Dev 1993, 7:16321643. PubMed Abstract  Publisher Full Text

Decker CJ, Parker R: Mechanisms of mRNA degradation in eukaryotes.
Trends Biochem Sci 1994, 19:336340. PubMed Abstract

Voit EO: Design principles and operating principles: the yin and yang of optimal functioning.
Mathematical biosciences 2003, 182:8192. PubMed Abstract  Publisher Full Text

Hood L, Heath J, Phelps M, Lin B: Systems biology and new technologies enable predictive and preventative medicine.
Science 2004, 306:640643. PubMed Abstract  Publisher Full Text

Slotine JJE, Li W: Applied Nonlinear Control. Prentice Hall; 1991.

Beelman CA, Parker R: Degradation of mRNA in eukaryotes.
Cell 1995, 81:179183. PubMed Abstract  Publisher Full Text

Tucker M, Parker R: Mechanisms and control of mRNA decapping in Saccharomyces cerevisiae.
Annu Rev Biochem 2000, 69:571595. PubMed Abstract  Publisher Full Text

Steiger M, Parker R: Analyzing mRNA decay in Saccharomyces cerevisiae.
Methods Enzymol 2002, 351:648660. PubMed Abstract

Takagi T, Sugeno M: Fuzzy identification of systems and its applications to modeling and control.

Yildirim N, Mackey MC: Feedback regulation in the lactose operon: a mathematical modeling study and comparison with experimental data.
Biophysical Journal 2003, 84:28412851. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Has'minskii RZ: Stochastic Stability of Differential Equations. Alphen, Sweden: Sijthoff and Noordhoff; 1980.