Abstract
Background
Noise has many important roles in cellular genetic regulatory functions at the nanomolar scale. At present, no good theory exists for identifying all possible mechanisms of genetic regulatory networks to attenuate the molecular noise to achieve regulatory ability or to amplify the molecular noise to randomize outcomes to the advantage of diversity. Therefore, the noise filtering of genetic regulatory network is an important topic for gene networks under intrinsic fluctuation and extrinsic noise.
Results
Based on stochastic dynamic regulation equation, the intrinsic fluctuation in reaction rates is modeled as a statedependent stochastic process, which will influence the stability of gene regulatory network, especially, with low concentrations of reacting species. Then the mechanisms of genetic regulatory network to attenuate or amplify extrinsic fluctuation are revealed from the nonlinear stochastic filtering point of view. Furthermore, a simple measure of attenuation level or amplification level of extrinsic noise for genetic regulatory networks is also introduced by nonlinear robust filtering method. Based on the global linearization scheme, a convenient method is introduced to measure noise attenuation or amplification for each gene of the nonlinear stochastic regulatory network by solving a set of filtering problems, which correspond to a set of linearized stochastic regulatory networks. Finally, by the proposed methods, several simulation examples of genetic regulatory networks are given to measure their robust stability under intrinsic fluctuations, and to estimate the genes' attenuation and amplification levels under extrinsic noises.
Conclusion
In this study, a stochastic nonlinear dynamic model is developed for genetic regulatory networks under intrinsic fluctuation and extrinsic noise. By the method we proposed, we could determine the robust stability under intrinsic fluctuations and identify the genes that are significantly affected by extrinsic noises, which we call the weak structure of the network. This method will be potential for robust gene circuit design in future, on which a drug design could be based.
Background
Molecular noise has been shown to play many roles in the biological functions of genetic regulatory networks, including noisedriven divergence of cell fates and population heterogeneity, noiseinduced amplification of signals, generation of errors in DNA replication leading to mutation and evolution, and maintenance of the quantitative individual of cells [1]. Other cellular processes influenced by noise include ionchannel gating [2], neural firing [3], developmental module [4,5], cytoskeleton dynamics [6] and motors [7]. Phase variation in pathogenic bacteria, where cells alternate randomly between expressing certain genes and silencing others, is thought to be a form of cultivated noise [8]. These molecularlevel noisy phenomena are deeply rooted in the statistical mechanical behavior of socalled nanoscale chemical systems, where concentrations of reacting species are extremely low and, consequently, fluctuations (noises) in the reaction rates are large [9]. Even though the molecular fluctuations leading to phase variation seem random in the individual, regulatory factors tune the variation to ensure mean levels of heterogeneity for the population, i.e., the random noises can be shown to be filtered or attenuated by the genetic regulatory networks [1].
Since molecular events in cells are subject to significant thermal fluctuations and noisy process with transcriptional control, alternative splicing, translation, diffusion and chemical modification reaction, thus gene expression is best viewed as a stochastic process. Many observations suggest that molecular events underlying cellular physiology are subject to fluctuations and have lead to the proposal of a stochastic model for gene expressions and biofunctions [913].
Noise attenuation can be considered from the signal processing perspective [1,14,15]. From this perspective, a pathway is viewed as an analog filter in terms of its frequency response. In terms of signal processing, these pathways function as lowpass filters to transduce lowfrequency signal and to attenuate highfrequency noise [14,16]. Negative feedback schemes are the most common noiseattenuation regulatory mechanism, which make a gene network robust to gene expression noise [1719]. Intrinsic chemical damping, integral feedback and redundancy have also found to be efficient noise filtering schemes in genetic regulatory systems [9,20].
Some cellular processes amplify or exploit noise in some sense, rather than just controlling or eliminating it. These processes fall into two mechanisms – one gives rise to population heterogeneity (and thus produces diversity with fitness advantage) [19], and the other uses noise to attenuate noise. Several results of isogenic heterogeneity illustrate how intrinsic molecular noise is used to generate diversity and how intrinsic molecular noise results in phenotypic variation and cellular differentiation [19,21]. How gene expression noise can lead to interesting dynamics in gene regulatory networks in found in [22]. The fim network regulates phase variation of type 1 pili uropathic E. coli . This system provides an example of how integrated regulatory modules in a network can function to both shape and filter noise, thereby creating environmentally tuned heterogeneity in a cell population [8,23]. Positive feedback can amplify the effect of noise and population heterogeneity by autocatalytic mechanisms [19,2426]. In addition to generating heterogeneous populations, cells also use noise to filter noise. Although in most systems, noise degrades a signal, when certain nonlinear effects are present, noise actually enhance a signal; for example, stochastic resonance [27].
Some elementary mechanisms for noise attenuation and amplification are tractable because they appear simple and identifiable. However, elementary mechanisms typically do not function in isolation, but rather interact in complex networks involving multiple feedback loops. Although it is straightforward to understand how a simple feedback loop shapes noise, it is far more difficult to understand the composite noise shaping behavior of multiple mechanisms interconnected in complex regulatory networks [28]. Recently, several computational models have been proposed for genetic regulatory networks that control circadian rhythms with regular oscillations in the presence of noise, using both deterministic and stochastic model [2931]. The stochastic model is able to produce regular oscillations, whereas the deterministic model can not [31], suggesting that the regulatory networks may utilize molecular fluctuations to their advantage.
Several examples have shown that noise attenuation arises from systematic properties of networks rather than a single mechanism. What specific mechanisms confer robust functionality in the presence of noise? It is clear that large, complex networks are able to function reliably despite inherent noise attributable to molecular fluctuations. Apparently, noise attenuation arises from complex network mechanisms involving multiple feedback loops [1,32].
Although theoretical and computational tools exist for analyzing the properties of a given genetic regulatory network, no good theory exists for identifying noise attenuation and amplification mechanisms of the networks. Although noise attenuation and amplification examples abound, how cells are able to manipulate biochemical noise remain unknown. By what means do regulatory networks attenuate the noise? And how and why do networks amplify noise? As pointed out by Rao et al. [1], these questions present one of the most challenging and fascinating problems for systems biologists, since they open questions in physiology, development and evolutionary biology. The answer likely resides in complex gene regulatory networks. Stochastic dynamic models are the ideal tools for such investigations, because they allow us to describe quantitatively the current states of network structure and component interactions to explore the network stochastic dynamics under intrinsic fluctuations and extrinsic noises.
Recently, a robustness measure for biochemical networks has been discussed at steady states by the singular value of the system matrix for a Ssystem model with deterministic parameter perturbations [33]. However, using this approach, it is not easy to discuss the robustness of nonlinear gene regulatory networks under stochastic intrinsic and extrinsic noises. In this study, based on a stochastic dynamic model of genetic regulatory networks with intrinsic and extrinsic noises, the robust stability to tolerate the intrinsic noise is initially revealed from the Lyapunov (energy) function point of view, and then a measure of the extrinsic noise attenuation level or amplification level of regulatory networks is developed from the H_{∞ }filtering point of view. In the last two decades, H_{∞ }control theory has been developed to efficiently attenuate the effect of extrinsic disturbance from the L_{2}gain point of view [3436]. Recently, noise attenuation of nonlinear stochastic systems has been developed for state estimation from the H_{∞ }filtering perspective [15]. These attenuation methods for extrinsic disturbances could be employed with some modifications as theoretical and computational bases for noise attenuation and amplification of gene regulatory networks. For the simplicity of illustration, a linear stochastic model with extrinsic and intrinsic noises is given at first to discuss the noise shaping of genetic regulatory network. Then the nonlinear stochastic equation is introduced to describe a general genetic regulatory network under extrinsic and intrinsic noises.
In this study, the intrinsic noises due to parametric fluctuations are modeled as state dependent noise, which will influence the stability of the gene regulatory network. The robust stability to tolerate these intrinsic parameter fluctuations by gene regulatory networks is discussed by the Lyapunov stability theory of nonlinear stochastic systems. We need to solve a HamiltonJacobi inequality (HJI) to measure the robustness of stability of nonlinear gene regulatory networks [37]. The ability to attenuate the extrinsic noises of nonlinear gene regulatory networks is measured based on the nonlinear H_{∞ }filtering theory [15]. In order to avoid solving HJI in nonlinear stochastic gene regulatory networks, based on the global linearization technique, a set of linear matrix inequalities (LMIs) [38], which can be efficiently solved by LMI Toolbox of Matlab, are employed to replace HJI. This allows us to measure the robust stability with respect to parametric fluctuations and to estimate the attenuation level or amplification level of nonlinear stochastic gene regulatory networks under intrinsic extrinsic noises.
Finally, a genetic regulatory network under intrinsic and extrinsic noises is considered for the illustration and the performance confirmation of the proposed method. According to intense simulation results, the proposed attenuation and amplification schemes could be a satisfactory method to interpret the stability robustness and noise filtering of each gene in a gene regulatory network under intrinsic fluctuations and extrinsic noises from the systems biology perspective.
Results
Stochastic system model and noise filtering
Network robust stability under intrinsic fluctuation
Initially, for the convenience of illustration, we consider only the following linear biochemical dynamics of a ngenes genetic regulatory network in Figure 1
Figure 1. The linear genetic regulatory network with intrinsic fluctuation ΔN_{ij }and extrinsic fluctuation υ_{i}(t).
where the concentration vector x(t) and stoichiometric matrix N are given by
in which x_{i}(t) denotes the concentration of the ith gene, and N_{ij }denotes the interaction between gene j and gene i.
Suppose the linear genetic regulatory network suffers intrinsic molecular fluctuations mainly due to stochastic thermal fluctuation so that stoichiometric matrix N is perturbed as N + ΔN, where the random components of the perturbation ΔN could be modeled by a stochastic Wiener process ω(t) on a probability space (Ω, p), which is a mathematical description of the socalled Brownian motion [39].
In intrinsic fluctuations, the perturbed dynamics of the nominal genetic regulatory network in equation (1) could be modeled as the following linear stochastic equation
dx(t) = Nx(t)dt + Mx(t)dω(t) (2)
This is a standard linear stochastic dynamic equation with state dependent noise. Mx(t)dω(t) denotes the term due to the intrinsic fluctuation ΔNx(t), where the stochastic property of fluctuation ΔN is extracted out to the standard Wiener process (or Brownian motion) ω(t) (which, roughly speaking, is integral of white noise) using E{ω(t)  ω(τ)^{2 }}= σ^{2 }t  τ with unit covariance σ^{2 }= 1 For example, , where M_{ij }denotes the deterministic part (amplitude) of fluctuation and n(t) is a white Gaussian noise with zero mean and unit variance to denote the stochastic part of fluctuation, i.e., the stochastic part of fluctuation is absorbed to n(t) with dω(t) = n(t)dt [15,36,37]. If some components ij of N are free of intrinsic fluctuation, then the corresponding M_{ij }should be equal to zero. And the covariance of ΔN_{ij}(t) is Cov(M_{ij}n(t), M_{ij}n(τ)) = δ (t  τ), where δ (t) is the impulse function.
In the conventional notation of engineering and system science, the stochastic dynamic equation could be represented by [39]
where n(t) ≡ dω(t)/dt denotes a normalized white Gaussian noise with unit covariance. Since ω(t) is a stochastic process, x(t) in equation (2) or (3) is also a stochastic process. Actually, the stochastic dynamic equations of genetic regulatory networks are always nonlinear. In order to meet the nonlinear stochastic regulatory networks, equation (3) should be generalized as the following Langevin equation [1]
dx(t) = N(x)dt + M(x)dω(t) (4)
where N(x) denotes the nonlinear interaction equation of nonlinear genetic regulatory network and M(x)dω(t) is due to nonlinear intrinsic fluctuation.
Remark: (i) At different operation points, the linearization of nonlinear stochastic regulatory network of equation (4) will be of the form in equation (3), i.e., at an operation point x = x_{0}, and . (ii) In this study, the white noise is normalized with unit variance, in which the covariance of stochastic noise is absorbed by M.
Since the effect of intrinsic biochemical kinetic parametric fluctuation is state dependent and will influence the stability of genetic regulatory networks, it will be discussed first. Let V(x) > 0 with V(0) = 0 denote the Lyapunov (power like) function of a stochastic genetic regulatory network. In the linear genetic regulatory network, the Lyapunov function is always chosen as V(x) = x^{T}Px for some symmetric positive definite matrix P. Then the equilibrium point x = 0 of a stochastic genetic regulatory network is stable in probability at the equilibrium point if the expectation of the derivative of V(x) is not positive [39], i.e., the total power (squares of concentrations) of the genetic regulatory network could not increase again in probability
Remark: In general, the nonlinear systems have many equilibrium points. If we are interested in the stability of the equilibrium point x_{e }≠ 0, for the convenience of discussion [35], the origin should be shifted to the equilibrium point we are interested in, i.e., x^{' }= x  x_{e}. Thus
dx'(t) = N (x' + x_{e})dt + M (x' + x_{e}) dω (t) (6)
Then the stability problem of the equilibrium point x_{e }in nonlinear stochastic network (4) is equivalent to the stability problem of x' = 0 in the nonlinear stochastic system (6).
In the linear stochastic genetic regulatory network (3), the robust stability under intrinsic fluctuation is given below
Theorem 1: The linear gene regulatory network with stochastic perturbation in equation (3) is stable in probability if the following Lyapunovtype equation
PN + N^{T}P + M^{T}PM ≤ 0 (7)
has a symmetric positive definite solution P > 0.
Proof: See Appendix 1.
Remark 1: (i) In the intrinsic noise free case, i.e., the nominal case in equation (1), the stable condition is that the matrix inequality PN + N^{T}P ≤ 0 has a symmetric positive definite solution P > 0. Obviously, the existence of a symmetric positive definite solution in (7) is more strict because the eigenvalues of system interaction matrix N should be located at the far left hand side of complex domain with large negative real values in order to overcome the extra term M^{T}PM due to intrinsic noise in equation (7). If some eigenvalues of system interaction matrix N are near the jω axis, then these modes are more easily to perturb by intrinsic molecular fluctuation across the jω axis, such that the linear genetic regulatory network becomes unstable. Therefore, the smallest distance between the locations of the eigenvalues of N to jω axis can be considered as a robustness measure of linear genetic regulatory networks. (ii) It is easy to use LMI Toolbox of Matlab to find a positive definite P to solve (7) if it exists.
For the nonlinear stochastic regulatory network in equation (4), we obtain the robust stability theory under intrinsic parametric fluctuation around the equilibrium point x_{e }= 0 as follows.
Theorem 2: The nonlinear perturbative genetic regulatory network in equation (4) is still stable in probability if the following Hamilton Jacobi inequality (HJI) has a positive solution V(x) > 0 with V(0) = 0
i.e., the nonlinear stochastic regulatory network in equation (4) is stable in probability.
Proof: See Appendix 2.
Remark 2: (i) The inequality of equation (8) is an extension of inequality in equation (7) from the linear stochastic case to the nonlinear stochastic case. (ii) In general, it is not easy to solve the nonlinear inequality in equation (8) to see whether the stability of nonlinear stochastic genetic network is guaranteed under intrinsic fluctuation. A convenient method based on the socalled global linearization [38] is proposed in the following to solve the robust stability problem of nonlinear stochastic gene regulatory networks. Suppose the global linearization of N(x) is defined as
and suppose Ω is bounded by the following polytope with n vertices
where Co{…} denotes the convex hull consisted by the vertices {…}. Therefore, all linearized systems of dx(t) = N(x)dt + M(x)dω(t) are bounded in the polytope consisting of the linearized vertices dx(t) = N_{i}x(t)dt + M_{i}x(t)dω(t), i = 1,2, … L. Then the perturbative nonlinear regulatory network is robustly stable under intrinsic noise if the following LMIs have a common positive solution P = P^{T }> 0 [38]
In general, it is very easy to find a common P > 0 to solve the above LMIs by the LMI Toolbox in Matlab if it exists. Therefore, it is more appealing to solve a P > 0 of LMIs in (10) than to solve a V(x) > 0 of HJI in equation (8) to guarantee the robust stability of nonlinear stochastic gene regulatory network under intrinsic noise, but the results of the LMI method may be more conservative.
Attenuation and amplification of extrinsic noise
After the robust stability of genetic regulatory network is guaranteed under the intrinsic biochemical parametric fluctuation, the effect of the extrinsic fluctuation on the network will be discussed. If the linear regulatory network in equation (2) also suffers the extrinsic signals υ(t) (or noises, see Figure 1) outside the network, then stochastic equation (2) is modified as follows
where H is a coupling matrix which denotes the influence of extrinsic signal on the state x(t). Z(t) denotes the concentration of some genes that we are interested in. For example, if we only want to discuss the effect of noises of ω(t) and υ(t) on gene i, i.e., x_{i}(t), then we let C = [000…1…00] i.e., every element of C is zero except 1 at the ith element. If we want to discuss the effect of noises on the whole genetic regulatory network, then we let C = I, the identity matrix.
Similarly, the nonlinear stochastic regulatory networks in equation (4) under extrinsic fluctuation should be modified as
The attenuation and amplification of extrinsic noise of stochastic regulatory networks in equations (11) and (12) will be discussed in the following theorems.
Let us denote L_{2}norm of υ(t) as
where E denotes the expectation. We denote υ(t) ∈ L_{2}, if  υ(t) _{2 }< ∞ Then the positive value ρ in the following inequality is called the effect (or gain) from the extrinsic noise υ(t) to Z(t) in the perturbative genetic regulatory network with x(0) = 0 [34,36,38]
If ρ < 1, we say the extrinsic noise υ(t) is attenuated at Z(t) by the genetic regulatory network. If ρ > 1, we say the extrinsic noise υ(t) is amplified at Z(t) by the genetic regulatory network. In this situation, ρ is called the attenuation level if ρ < 1 or amplification level if ρ > 1. If ρ = 1, we call it lossless. In the inequality (14), it is under the assumption of x(0) = 0, i.e., all signals in the gene network are driven by noises. If the initial condition is not zero, i.e., x(0) ≠ 0, then an extra term of initial condition should be added as follows [15,36]
for some positive function V(x(0)) of the initial condition x(0) .
Remark: If υ(t) is a deterministic signal from the environment, then should be changed to in (14) and (15).
Based on the analysis above, some theorems about the attenuation or amplification of extrinsic noise of stochastic genetic regulatory networks are discussed below
Theorem 3: The attenuation level ρ of linear perturbative genetic regulatory network in equation (11) is guaranteed if the following inequality has a positive definite solution P > 0
By Shur complement [38], inequality (16) is equivalent to the following LMI
Proof: See Appendix 3.
The optimal attenuation level ρ_{o }of the linear stochastic genetic regulatory network (11) can be obtained by solving the following constrained optimization
The above LMI optimization can be solved easily by decreasing ρ until no positive definite solution P > 0 for equation (17) can be found. Software packages such as LMI Optimization Toolbox in Matlab have been developed to easily solve the above LMI optimization.
Remark 3: (i) If ρ_{o }< 1 in equation (18), the extrinsic noise υ(t) is attenuated by the genetic regulatory network at Z(t) i.e., this gene is less sensitive to the environmental noise. (ii) If ρ_{o }> 1, the extrinsic noise υ(t) is amplified by the genetic regulatory network at Z(t) i.e., this gene is more sensitive to the environmental noise. (iii) If we only want to check if the genetic regulatory network has a prescribed attenuation level ρ, we just check if inequality (17) has a positive definite solution P > 0.
For the nonlinear perturbative genetic regulatory network in equation (12), the attenuation or amplification of extrinsic noise υ(t) at Z(t) is discussed in the following theorem.
Theorem 4: The attenuation level ρ of nonlinear genetic regulatory network (12) is guaranteed if the following HamiltonJacobi inequality has a positive definite solution
Proof: See Appendix 4.
The optimal attenuation level ρ_{o }of the nonlinear stochastic regulatory network can be obtained by solving the following constrained optimization
Remark 4: (i) In general, there exists no systematic method to solve the constrained optimizations in (20). It should be solved case by case by decreasing ρ until no positive solution V(x) for the HamiltonJacobi inequality in (19) can be found. (ii) An approximation solution for (20) based on the "global linearization" techniques of (9) and (10) in Remark 2 (ii) is introduced in the following. Thus, if the nonlinear perturbative genetic regulatory network (12) could be bounded in the polytope consisting of L linearized vertices [38] dx(t) = (N_{i}x(t) + H_{i}υ(t))dt + M_{i}x(t)dω(t), i = 1,2, …L, then after some rearrangements the optimal attenuation level ρ_{o }in Theorem 4 can be approximated by solving the following constrained optimization problem
This result is similar to the constrained optimization in (18) except that a set of LMI constraints, i.e., the HJI constraint in (19) is replaced by a set of LMI constraints in (21) to make the solution feasible. However, it may lead to a conservative result.
Computational simulations
To confirm the validity of the stability robustness and the noise attenuation or amplification schemes we proposed, we conducted several computational simulations of a typical genetic regulatory network. The detailed equations and parameters are shown in Appendix 5.
Consider a typical genetic regulatory network, as shown in Figure 2 [40,41]. This is a typical gene interaction system describing the gene, mRNA and protein interactions. X_{1 }is an mRNA produced from gene 1, X_{2 }is an enzyme protein that it produces, and X_{3 }is an inducer protein catalyzed by X_{2}. In addition, X_{4 }is an mRNA produced from gene 4 and X_{5 }is a regulator protein it produces. Positive feedback from the inducer protein X_{3 }and negative feedback from the regulator protein X_{5 }are assumed in the mRNA production processes of gene 1 and gene 4 [42]. The genetic regulatory network can be represented as follows
Figure 2. A typical genetic regulatory network describing the gene, mRNA and protein interactions [40].
In the nominal case, the kinetic parameters are
and the dynamic response of the genetic regulatory network is shown in Figure 3(a). X_{e }= [0.7339 0.7339 1 0.9283 0.9283]^{T }is a stable equilibrium point of the genetic regulatory network. Since we are interested in the robust stability of the equilibrium X_{e }= [0.7339 0.7339 1 0.9283 0.9283]^{T }under stochastic parameter perturbation, the origin should be shifted to the equilibrium X_{e }= [0.7339 0.7339 1 0.9283 0.9283]^{T }. Then the genetic regulatory network should be rewritten under coordinate shift (see equation (A17)) and the dynamic response of the genetic regulatory network under coordinate shift is shown in Figure 3(b). It is seen that the origin is shifted to the equilibrium point X_{e}, i.e., x' = 0 is a stable equilibrium point.
Figure 3. The dynamic responses of the genetic regulatory network in computational simulation under different stochastic perturbative cases. (a) The case without perturbation. (b) The case without perturbation but under coordinate shift. (c) The case with perturbation in (A18). (d) The case with perturbation and under coordinate shift in (A19). (e) The case with perturbation in (A20). (f) The case with perturbation and under coordinate shift in (A21). (g) The noisedriven case in (A22). (h) The signaldriven case in (A22) when υ(t) = 2+sin(10t).
Suppose the genetic regulatory parameter suffer some stochastic perturbations shown as (A 18). In order to discuss the robust stability of the stochastic regulatory network in (A18) at the equilibrium point X_{e }= [0.7339 0.7339 1 0.9283 0.9283]^{T}, we can rewrite the genetic regulatory network under such perturbations by coordinate shift as (6) (shown in (A19)). Following remark 2(ii), we can globally linearize the system and solve the LMIs in (10) to see if the perturbed stochastic system is stable or not under this kind of stochastic intrinsic noise. In the perturbative case, we can find a positive definite P of the LMIs (shown in Appendix 5), so the stochastic system (A19) is stable in probability at x' = 0 under this intrinsic noise. The dynamic responses of the perturbed stochastic systems in (A18) and (A19) are shown in Figure 3(c) and Figure 3(d), respectively, which confirms the conclusion of our proposed computational method of robust stability.
Suppose the kinetic parameters suffer another stochastic perturbations shown in (A20). In order to discuss the robust stability of the stochastic regulatory network in (A20) at the equilibrium point X_{e }= [0.7339 0.7339 1 0.9283 0.9283]^{T}, we again rewrite the genetic regulatory network under such perturbations by coordinate shift as (A21). In this stochastic noise case, we can not find a positive definite P of the LMIs, so the stability of the stochastic regulatory system is not guaranteed. This result is consistent with the dynamic responses we present in Figures 3(e) and 3(f) for the stochastic system (A20) and the shifted stochastic system in (A21), respectively.
After validating the network robust stability under intrinsic fluctuations, we want to confirm the theorem in (21) about the estimation of the attenuation or amplification level ρ_{o }of the extrinsic noise. Suppose the genetic regulatory network in (22) suffers an intrinsic stochastic kinetic parameter perturbation as (A18) and an extrinsic noise (shown in (A22)). In this case, H in (10) is an identity matrix. Let υ(t) = [υ_{1}(t) υ_{2}(t) υ_{3}(t) υ_{4}(t) υ_{5}(t)]^{T }denote the extrinsic noise vector. We are interested in the effects of extrinsic noise υ(t) on individual gene i, so we let C in (11) be [1 0 0 0 0], [0 1 0 0 0], [0 0 1 0 0], [0 0 0 1 0], [0 0 0 0 1], for gene 1, gene 2, gene 3, gene 4, gene 5, respectively. For the convenience of validation, we let υ(t) be white Gaussian noise with mean = 2 and variance = 1. By Remark 4(ii), we again globally linearize the system and solve the LMIs in (21) to calculate the attenuation (amplification) level ρ_{o}. After solving the constrained optimization problem in (21), the optimal attenuation level ρ_{o }of genes 1, 2, 3, 4, 5 are 0.5723, 0.5328, 0.5489, 0.6980, 0.6964, for X_{1}, X_{2}, X_{3}, X_{4 }and X_{5}, respectively. This means that the noise attenuation levels of these genes can not exceed these values, respectively. Obviously, extrinsic noises are all attenuated at these genes in the gene network, so these genes are robust to these noises.
The dynamic response of the noisedriven case is shown in Figure 3(g). We can compute the L_{2}norm of Z(t) in Figure 3(g) to verify the attenuation level ρ_{o }that we obtain by the computation method above. Thus, by the system simulation results, we get attenuatoion level ≈ 0.4954 < 0.5723 for gene 1, ≈ 0.5133 < 0.5328 for gene 2, ≈ 0.4916 < 0.5489 for gene 3, ≈ 0.6281 < 0.6980 for gene 4, and ≈ 0.6424 < 0.6954 for gene 5, respectively. Obviously, the inequality in (14) holds, i.e., the proposed optimal noise attenuation levels are the upper bounds of those computed by the simulation results. In this genetic regulatory network, we found that the noises at X_{1}, X_{2}, X_{3}, X_{4}, and X_{5 }are all attenuated by the network, which are validated by both the proposed optimal attenuation measure method and system simulation.
In the stochastic genetic regulatory system (A22), υ(t) can not only be of extrinsic noise but also be of extrinsic signal such as sinusoid signal. In this situation, we let υ(t) in (A22) be 2+sin(10t) to verify the attenuation levels of extrinsic signal from the environment. The dynamic response of the signaldriven case is shown in Figure 3(h). By the simulation results, we can also obtain the attenuation level ≈ 0.4895 < 0.5723 for gene 1, ≈ 0.5113 < 0.5328 for gene 2, ≈ 0.4831 < 0.5489 for gene 3, ≈ 0.6080 < 0.6980 for gene 4, and ≈ 0.6227 < 0.6964 for gene 5, respectively, which are also consistent with the optimal attenuation levels calculated by the method we proposed.
Discussion
Genetic and biochemical networks are influenced by unavoidable fluctuations which cause random perturbations of biochemical parameters. These perturbations are reflected in dynamic models as numerical changes in reaction rate constants of the system. Genetic regulatory networks also suffer from environmental disturbances. We classify these perturbations and disturbances as intrinsic and extrinsic noise which are modeled as state dependent noise and state independent noise, respectively. In this study, we propose a method to examine the network robust stability under intrinsic noise and to measure the attenuation or amplification of extrinsic noise. We found that the dynamic stability of a nonlinear perturbative genetic regulatory network under intrinsic noise is guaranteed if we can find a Lyapunov function which satisfies the Hamilton Jacobi inequality (HJI) in (8). We can also estimate the attenuation (amplification) level of the nonlinear perturbative genetic regulatory network if we can find a Lyapunov function which satisfies the HJI in (19). The robustness measurement of Ssystems in biochemical network [33] is based on the deterministic parameter perturbations at steady state. The attenuation or amplification of extrinsic noises has not been considered. This study considers the robust stability problem of a more general gene regulatory network with intrinsic and extrinsic noises.
The effect of different types of noise on biochemical networks or genetic regulatory networks is studied based on the robust stability and noise attenuation (amplification) from the stochastic system perspective. The proposed method can work on both linear and nonlinear biosystems. In addition, it can be applied to all kinds of model types. Once we model the genetic regulatory or biochemical network, we can analyze the robust stability and noise attenuation by the methods proposed here. In general, it is not easy to solve HJI in (8) to find the robust stability of genetic regulatory networks under intrinsic noise, or to solve HJI in (19) to get the attenuation level or amplification level to extrinsic noise. In order to avoid solving the HJI in a nonlinear stochastic gene regulatory network, which is not easy to solve directly, we replace the HJI by global linearization with a set of linear matrix inequalities (LMIs). We can easily solve the LMIs by Matlab LMI Toolbox. Therefore, we can analyze the effect of noises on each gene of a genetic regulatory network more efficiently.
This method provides a way to determine whether the robust stability of the genetic regulatory network is guaranteed under intrinsic fluctuations. Furthermore, we could also estimate to what extent the fluctuations can cease the stability of the network. Therefore, we can gain more insight into how the genetic regulatory network is robustly stabilized. It also provides a quantitative measurement to see to what extent the extrinsic noise is attenuated or amplified at each gene in the network. By choosing different values for C in (11) or (12) from [1 0 0…0 0] to [0 0 0…0…0 1], we can quantify the attenuation or amplification level of extrinsic noises at different genes in a genetic regulatory network. We can also compare the attenuation or amplification levels of noises at different networks by choosing C as I, the identity matrix.
From LMI in (7) or (10), in order to let P > 0, it is seen that there are two schemes for genetic regulatory networks to improve the robust stability against intrinsic noise. One scheme is to make the eigenvalues of the system matrix N as negative as possible (i.e., the eigenvalues (or system poles) of N are far to the left hand side of the complex domain) so that large parameter variation matrix M could be tolerated [9]. Adequate negative feedbacks in networks could place the eigenvalues of N to far left hand side of the complex domain to achieve this purpose. Another scheme is to make the parameter variation matrix M as small as possible so that the robust stability condition (7) or (10) is not violated. Pervasive redundancy in genetic regulatory network could achieve this purpose [9].
Similarly, from (16) or (17), in order to let P > 0, we examine the condition where N is more negative, and M and H are smaller, i.e., the eigenvalues of system matrix N are all far to the left hand side in the complex domain, the magnitude of intrinsic noise is smaller, and the coupling between the network and environment is weak. In this case, the genetic regulatory network will attenuate the extrinsic noise to a smaller level.
Our noise analysis has been confirmed by the typical genetic regulatory network with numerical simulations. From the simulations, we found that if we can find a common P > 0 of LMIs in (10), then the robust stability of the network is guaranteed. Therefore, the network can tolerate this kind of noises. Moreover, we also found that we can calculate the attenuation or amplification level of noises for each gene by solving LMIs in (18) or (21). Comparing the results of attenuation (or amplification) level by the proposed analysis method with those by the true system simulations, however, we find that there are still some discrepancies. This is because the HJI problem is replaced by solving LMIs, which are only an approximation method and may lead to a conservative result and the proposed optimal attenuation measurements are the upper bounds of those computed by the simulation results.
Once we could analyze the effects of the noises on the genetic regulatory network, we could know which genes are influenced much by noise, i.e., are more sensitive to noise. These genes are at the locations which we call the weak structure of the network. This information could provide potential therapeutic targets for treatment of disease. Understanding the effects of noise and the weakness of the network is helpful for improving the stability robustness or achieving a desired noise attenuation or amplification of the genetic regulatory network, on which a drug design could be based. In the future, we may be able to design drugs or construct some other pathways by transfection and transformation biotechnologies to improve the robustness of weak structures of gene networks, so that genetic regulatory networks or biochemical networks can be protected from the influence of the environmental extrinsic noise and the perturbative intrinsic noise.
Conclusion
In this study, a nonlinear stochastic system is developed to model a gene regulatory network under intrinsic noise and extrinsic noise. Based on the stochastic stability, the robust stability condition is developed as in Theorem 1 for linear stochastic gene regulatory network and in Theorem 2 for nonlinear gene regulatory network. By the global linearization method, the robust stability condition of nonlinear gene regulatory network is reduced by solving LMIs in (10) for the convenience of computation.
Based on the robust H_{∞ }filtering theory, the amplification or attenuation of extrinsic noise at an interested gene is to solve the minimum ρ_{o }in (18) for linear stochastic gene regulatory network and ρ_{o }in (2) for nonlinear stochastic gene regulatory network, which could be reduced to the constrained optimization problem in (21). If ρ_{o }< 1, the extrinsic noise is attenuated at that gene and if ρ_{o }> 1, the extrinsic noise is amplified at that gene. The genes with ρ_{o }< 1 are more robust than the genes with ρ_{o }> 1. If ρ_{o }of a gene is much larger than 1, it should be significantly affected by extrinsic noises and is at a weak structure of the gene network. The robust filtering analysis will be potential for robust gene circuit design in future, on which gene therapy and drug design could be based.
The future work will focus on gene circuit control design method to improve the robust stability of nonlinear gene network to tolerate the much large intrinsic noise and to achieve a prescribed attenuation level ρ_{o }for a desired extrinsic noise filtering. These robust circuit design methods are useful for the biotechnological purpose or the therapeutic purpose.
Methods
The stochastic stability of the Langevin equation (4) for nonlinear stochastic regulatory networks is discussed by the stochastic Lyapunov theory [35,39]. Suppose V(x) > 0 with V(0) = 0 denotes the Lyapunov function of the Langevin equation in (4). Then the equilibrium point x = 0 of the stochastic gene regulatory network is stable in probability under intrinsic noise if the total power could not be increase, i.e., the inequality (5) should hold. In the linear stochastic gene regulatory network in (3), the robust stability is reduced to the existence of P > 0 in the LMI of (7) in Theorem 1. By the global linearization, the robust stability condition for nonlinear stochastic gene regulatory network is reduced to the existence of a common P > 0 in the LMIs in (10), which could be easily solved by the LMI toolbox in Matlab.
The stochastic filtering of extrinsic noise υ(t) in the nonlinear stochastic gene network is considered from the nonlinear H_{∞ }filtering theory of nonlinear stochastic system [15], i.e., the effect of extrinsic noise on an interested gene Z(t) should be equal to or less than a level ρ_{o}. That means the inequality (14) or (15) should hold. Then the filtering problem for the nonlinear stochastic gene regulatory network under intrinsic noise and extrinsic noise is to solve the constrained optimization problem in (20) to obtain ρ_{o}. If ρ_{o }< 1, then the extrinsic noise is attenuated and if ρ_{o }> 1, then the extrinsic noise is amplified. In the linear stochastic gene regulatory network in (11), the filtering problem is reduced to solving the constrained optimization problem in (18). By the global linearization method, the nonlinear filtering problem in (20) is reduced to the constrained optimization problem in (21), which could be easily solved by LMI toolbox in Matlab.
Appendix
Appendix 1: proof of Theorem 1
Let us denote the Lyapunov (energy) function of linear stochastic regulatory network in equation (2) as V(x) = x(t)^{T}Px(t), for some symmetric positive definite matrix P^{T }= P > 0. Then the change of Lyapunov function is obtained by following Itô derivative [36,39]
From the inequality (7), we get
Since the change of Lyapunov (energy) function is nonpositive, the linear stochastic genetic regulatory network is stable in probability.
Appendix 2: proof of Theorem 2
Let us denote the Lyapunov function V(x) of the perturbative nonlinear genetic regulatory network. By following Itô derivative, we get [15,37]
By the inequality (8), we get
Then the nonlinear perturbative stochastic regulatory network in equation (4) is stable in probability.
Appendix 3: proof of Theorem 3
Consider the following equality,
By Itô formula for stochastic equation (11), we get
Let us choose the Lyapunov function as V(x) = x(t)^{T }Px(t). Then
Taking expectation to both sides of (A5) and substituting equation (A7) into equation (A5), we get
By the fact that
X^{T}Y + Y^{T}X ≤ ρ^{2 }XX^{T }+ ρ^{2}Y^{T}Y, for any ρ^{2 }> 0 (A9)
and E{V(x(∞))} ≥ 0, we get
By the inequality (16), we get
or
Appendix 4: proof of Theorem 4
Consider the following equality,
By Itô formula for stochastic equation (12), we get
Taking expectation to both sides of (A13) and substituting equation (A14) into equation (A13), we get
By the fact E{V(x(∞))} ≥ 0 and the in equality (A9), we get
By the inequality (19), we get
or
Appendix 5: details of computational simulations
The dynamic response of the typical genetic regulatory network in (22) is shown in Figure 3(a), where we can find that X_{e }= [0.7339 0.7339 1 0.9283 0.9283]^{T }is a stable equilibrium point of the genetic regulatory network. Since we are interested in the robust stability of the equilibrium under stochastic parameter perturbation, the origin should be shifted to the equilibrium. Then the genetic regulatory network in (22) under coordinate shift can be rewritten as follows
and the dynamic response of the genetic regulatory network under coordinate shift is shown in Figure 3(b).
Suppose the kinetic parameters k_{i}, i = 1,…,10, suffer some stochastic perturbations k_{i }+ Δk_{i}, i = 1, …,10, respectively, where Δk_{1 }= 1n(t), Δk_{2 }= 2n(t), Δk_{3 }= 2n(t), Δk_{4 }= 0.3n(t), Δk_{5 }= 0.8 n(t), Δk_{6 }= 1.4n(t), Δk_{7 }= 1.2(t), Δk_{8 }= 2n(t), Δk_{9 }= 1.8n(t) and n(t) denotes the standard zero mean white Gaussian noise with unit variance. Then the genetic regulatory network under such perturbations becomes the following nonlinear stochastic system
We can rewrite the genetic regulatory network under such perturbations by coordinate shift as (6) as follows
In the perturbative case, we can use global linearization method and solve the LMIs in (21), Then the solution P of the LMIs in (21) is found as
where the eigenvalues of P are 0.074151, 0.15722, 0.20491, 0.22885, 0.50216, which are all positive, so the stochastic system (A19) is stable in probability at x' = 0 under this intrinsic noise. The dynamic responses of the perturbed stochastic systems in (A18) and (A19) are shown in Figure 3(c) and Figure 3(d), respectively, which confirms the conclusion of our proposed computational method of robust stability.
If the kinetic parameters suffer the following stochastic perturbations Δk_{1 }= 4n(t), Δk_{2 }= 2n(t), Δk_{3 }= 6n(t), Δk_{4 }= 7n(t), Δk_{5 }= 5n(t), Δk_{6 }= 2n(t), Δk_{7 }= 6n(t), Δk_{8 }= 2.5n(t) Δk_{9 }= 8n(t), Δk_{10 }=  6n(t) the network of the stochastic parameter pertubative case becomes the following nonlinear stochastic system
We again rewrite the genetic regulatory network under such perturbations by coordinate shift as follows
In this stochastic noise case, we can not find a positive definite P of the LMIs in (21), so the stability of the stochastic regulatory system is not guaranteed. This result is consistent with the dynamic responses we present in Figures 3(e) and 3(f) for the stochastic system (A20) and the shifted stochastic system in (A21), respectively.
Suppose the genetic regulatory network in (22) suffers an intrinsic stochastic kinetic parameter perturbation as (A18) and an extrinsic noise υ(t), which is white Gaussian noise with mean = 2 and variance = 1, i.e.,
We can calculate the attenuation (amplification) level ρ_{o }of the system by solving the LMIs in (21). After solving the constrained optimization problem in (21), the optimal attenuation level ρ_{o }of genes 1, 2, 3, 4, 5 are 0.5723, 0.5328, 0.5489, 0.6980, 0.6964, for X_{1 }X_{2}, X_{3}, X_{4 }and X_{5}, respectively. This means that the noise attenuation levels of these genes can not exceed these values, respectively. The dynamic responses of different υ(t) are shown in Figure 3(g) and 3(h), from which we can calculate the attenuation levels of system simulation. The results are also consistent with the optimal attenuation levels calculated by the method we proposed.
Authors' contributions
BorSen Chen gave the topic and developed some results. YuChao Wang performed simulations and evaluated the results.
Acknowledgements
This study is supported by National Science Council of Taiwan under the contract NSC 9231112B007017.
References

Rao CV, Wolf DM, Arkin AP: Control, exploitation and tolerance of intracellular noise.
Nature 2002, 420:231237. PubMed Abstract  Publisher Full Text

White JA, Rubinstein JT, Kay AR: Channel noise in neurons.
Trends Neurosci 2000, 23:131137. PubMed Abstract  Publisher Full Text

Allen C, Stevens CF: An evaluation of causes for unreliability of synaptic transmission.
Proc Natl Acad Sci 1994, 91:1038010383. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

von Dassow G, Meir E, Munro EM, Odell GM: The segment polarity network is a robust developmental module.
Nature 2000, 406:188192. PubMed Abstract  Publisher Full Text

Houchmandzadeh B, Wieschaus E, Leibler S: Establishment of developmental precision and proportions in the early Drosophila embryo.
Nature 2002, 415:798802. PubMed Abstract  Publisher Full Text

van Oudenaarden A, Theriot JA: Cooperative symmetrybreaking by actin polymerization in a model for cell motility.
Nat Cell Biol 1999, 1:493499. PubMed Abstract  Publisher Full Text

Simon SM, Peskin CS, Oster GF: What drives the translocation of proteins?
Proc Natl Acad Sci 1992, 89:37703774. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Connell I, Agace W, Klemm P, Schembri M, Marild S, Svanborg C: Type 1 fimbrial expression enhances Escherichia coli virulence for the urinary tract.
Proc Natl Acad Sci 1996, 93:98279832. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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

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

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

Thattai M, van Oudenaarden A: Attenuation of noise in ultrasensitive signaling cascades.
Biophys J 2002, 82:29432950. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tegner J, Yeung MK, Hasty J, Collins JJ: Reverse engineering gene networks: integrating genetic perturbations with dynamical modeling.
Proc Natl Acad Sci 2003, 100:59445949. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Arkin AP: Signal processing by biochemical reaction networks. In SelfOrganized Biological Dynamics and Nonlinear Control: Toward Understanding Complexity, Chaos and Emergent Function in Living Systems. Edited by Walleczek J. London: Cambridge University Press; 2000:112144.

Zhang W, Chen BS, Tseng CS: Robust H_{∞ }filtering for nonlinear stochastic systems.
IEEE Trans Signal Processing 2005, 53:589598. Publisher Full Text

Paulsson J, Berg OG, Ehrenberg M: Stochastic focusing: fluctuationenhanced sensitivity of intracellular regulation.
Proc Natl Acad Sci 2000, 97:71487153. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Becskei A, Serrano L: Engineering stability in gene networks by autoregulation.
Nature 2000, 405:590593. PubMed Abstract  Publisher Full Text

Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli.
Nature 2000, 403:339342. PubMed Abstract  Publisher Full Text

Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: from theories to phenotypes.
Nat Rev Genet 2005, 6:451464. 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 2000, 97:46494653. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Arkin A, Ross J, McAdams HH: Stochastic kinetic analysis of developmental pathway bifurcation in phage λinfected Escherichia coli cells.
Genetics 1998, 149:16331648. PubMed Abstract  Publisher Full Text

Blake WJ, KAErn M, Cantor CR, Collins JJ: Noise in eukaryotic gene expression.
Nature 2003, 422:633637. PubMed Abstract  Publisher Full Text

Wolf DM, Arkin AP: Fifteen minutes of fim: control of type 1 pili expression in E. coli.
Omics 2002, 6:91114. PubMed Abstract  Publisher Full Text

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

Becskei A, Seraphin B, Serrano L: Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion.
EMBO J 2001, 20:25282535. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Isaacs FJ, Hasty J, Cantor CR, Collins JJ: Prediction and measurement of an autoregulatory genetic module.
Proc Natl Acad Sci 2003, 100:77147719. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gammaitoni L, Hanggi P, Jung P, Marchesoni F: Stochastic resonance.
Rev Mod phys 1998, 70:223287. Publisher Full Text

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

Barkai N, Leibler S: Circadian clocks limited by noise.
Nature 2000, 403:267268. PubMed Abstract  Publisher Full Text

Gonze D, Halloy J, Goldbeter A: Robustness of circadian rhythms with respect to molecular noise.
Proc Natl Acad Sci 2002, 99:673678. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Vilar JM, Kueh HY, Barkai N, Leibler S: Mechanisms of noiseresistance in genetic oscillators.
Proc Natl Acad Sci 2002, 99:59885992. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Barkai N, Leibler S: Robustness in simple biochemical networks.
Nature 1997, 387:913917. PubMed Abstract  Publisher Full Text

Chen BS, Wang YC, Wu WS, Li WH: A new measure of the robustness of biochemical networks.
Bioinformatics 2005, 21:26982705. PubMed Abstract  Publisher Full Text

Doyle J, Glover K, Khargonekar PP, Francis B: Statespace solutions to standard H_{2 }and H_{∞ }control problems.
IEEE Trans Automat Contr 1989, 34:831847. Publisher Full Text

Qu Z: Robust Control of Nonlinear Uncertain Systems. New York: John Wiley and Sons; 1998.

Chen BS, Zhang W: Stochastic H_{2}/H_{∞ }control with statedependent noise.
IEEE Trans Automat Contr 2004, 49:4557. Publisher Full Text

Zhang W, Chen BS: State feedback H_{∞ }control for a class of nonlinear stochastic systems.
SIAM J on Control and Optimization 2006, 44:19731991. Publisher Full Text

Boyd S, El Ghaoui L, Feron E, Balakrishnan V: Linear Matrix Inequalities in System and Control Theory. Philadelphia: SIAM; 1994.

Chen G, Chen Q, Hsu SH: Linear Stochastic Control Systems. New York: CRC Press; 1995.

Hlavacek WS, Savageau MA: Rules for coupled expression of regulator and effector genes in inducible circuits.
J Mol Biol 1996, 255:121139. PubMed Abstract  Publisher Full Text

Veflingstad SR, Almeida J, Voit EO: Priming nonlinear searches for pathway identification.
Theor Biol Med Model 2004, 1:8. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Kikuchi S, Tominaga D, Arita M, Takahashi K, Tomita M: Dynamic modeling of genetic networks using genetic algorithm and Ssystem.
Bioinformatics 2003, 19:643650. PubMed Abstract  Publisher Full Text