Stochastic synchronization of genetic oscillator networks

Background The study of synchronization among genetic oscillators is essential for the understanding of the rhythmic phenomena of living organisms at both molecular and cellular levels. Genetic networks are intrinsically noisy due to natural random intra- and inter-cellular fluctuations. Therefore, it is important to study the effects of noise perturbation on the synchronous dynamics of genetic oscillators. From the synthetic biology viewpoint, it is also important to implement biological systems that minimizing the negative influence of the perturbations. Results In this paper, based on systems biology approach, we provide a general theoretical result on the synchronization of genetic oscillators with stochastic perturbations. By exploiting the specific properties of many genetic oscillator models, we provide an easy-verified sufficient condition for the stochastic synchronization of coupled genetic oscillators, based on the Lur'e system approach in control theory. A design principle for minimizing the influence of noise is also presented. To demonstrate the effectiveness of our theoretical results, a population of coupled repressillators is adopted as a numerical example. Conclusion In summary, we present an efficient theoretical method for analyzing the synchronization of genetic oscillator networks, which is helpful for understanding and testing the synchronization phenomena in biological organisms. Besides, the results are actually applicable to general oscillator networks.


Background
Elucidating the collective dynamics of coupled genetic oscillators not only is important for the understanding of the rhythmic phenomena of living organisms, but also has many potential applications in bioengineering areas. So far, many researchers have studied the synchronization in genetic networks from the aspects of experiment, numerical simulation and theoretical analysis. For instance, in [1], the authors experimentally investigated the synchronization of cellular clock in the suprachias-matic nucleus (SCN); in [2][3][4], the synchronization are studied in biological networks of identical genetic oscillators; and in [5][6][7], the synchronization for coupled nonidentical genetic oscillators is investigated. Gene regulation is an intrinsically noisy process, which is subject to intracellular and extracellular noise perturbations and environment fluctuations [8][9][10][11][12]14]. Such cellular noises will undoubtedly affect the dynamics of the networks both quantitatively and qualitatively. In [13], the authors numerically studied the cooperative behaviors of a multicell system with noise perturbations. But to our knowledge, the synchronization properties of stochastic genetic networks have not yet been theoretically studied. This paper aims to provide a theoretical result for the synchronization of coupled genetic oscillators with noise perturbations, based on control theory approach. We first provide a general theoretical result for the stochastic synchronization of coupled oscillators. After that, by taking the specific structure of many model genetic oscillators into account, we present a sufficient condition for the stochastic synchronization in terms of linear matrix inequalities (LMIs) [15], which are very easy to be verified numerically. To our knowledge, the synchronization of complex oscillator networks with noise perturbations, even not in the biological context, has not yet been fully studied. Recently, it was found that many biological networks are complex networks with small-world and scalefree properties [16,17]. Our method is also applicable to genetic oscillator networks with complex topology, directed and weighted couplings. To demonstrate the effectiveness of the theoretical results, we present a simulation example of coupled repressilators. Throughout this paper, matrix U ∈ R N × N is defined as an irreducible matrices with zero row sums, whose off-diagonal elements are all non-positive, and the other notations are defined in the Appendix A.

Theoretical results
Since we know very little about how the cellular noises act on the genetic networks, a simple way to incorporate random effects is to assume that certain noises randomly perturb the genetic networks in an additive manner. We consider the following networks of N coupled genetic oscillators with random noise perturbations where F(·) defines the dynamics of individual oscillators, v i (t) ∈ R n × 1 is called the noise intensity vector, belongs to L 2 [0, ∞). As we will see in the following analysis, the results hold no matter what v i (t) is and no matter where it is introduced, so we do not explicitly express the form of v i (t) here. v i (t) can also be a function of the variables (if so, some minor modifications are needed in the following). n i (t) is a scalar zero mean Gaussian white noise process. Recall that the time derivative of a Wiener process is a white noise process [19], hence we can define dw i (t) = n i (t)dt, where w i (t) is a scalar Wiener process. Thus, the above equation can be rewritten as the following stochastic differential equation form: The work can be easily extended to the case that v i (t) ∈ and n i (t) = [n i1 (t), ʜ, (t)] T be an l i -dimensional mutually independent zero mean Gaussian white noise process. D ∈ R n × n defines the coupling between two genetic oscillators. G = (G ij ) N × N is the coupling matrix of the network. If there is a link from oscillator j to oscillator i (j ≠ i), then G ij equals to a positive constant denoting the coupling strength of this link; otherwise, G ij = 0; . Matrix G defines the coupling topology, direction, and the coupling strength of the network.
For network (1), a natural attempt is to study the meansquare asymptotic synchronization. But existing experimental results show that usually the genetic oscillators can not achieve mean-square synchronization (see, e.g. experimental results in [1] and Appendix B for a theoretical discussion). Analogue to the stochastic stability with disturbance attenuation (see, e.g. [20]), we give a less restrictive (but more realistic) definition of the stochastic synchronization as follows: By using the techniques described in the Appendix C, we know that if there exist matrices P > 0, T ∈ R n × n and U, and a scalar ρ > 0, such that the following conditions are satisfied, S 2 ≡ 2(y 1 -y 2 ) T P[F(y 1 ) -F(y 2 ) -T(y 1 -y 2 )] + (y 1 -y 2 ) T (y 1 -y 2 ) < 0, ∀y 1 , y 2 ∈ R n (y 1 ≠ y 2 ), U P ≤ ρI, then, the network (1) will achieve stochastic synchronization with disturbance attenuation γ.
The above condition (3) is a general result for the stochastic synchronization of coupled oscillators. But we do not have general efficient method for verifying the first inequality in (3) for arbitrary F due to its nonlinearity. Next we consider a special structure of genetic oscillators to obtain an easy-verified result.
Genetic oscillators are biochemically dynamical networks, which can usually be modelled as nonlinear dynamical systems. Mathematically many genetic oscillators can be expressed in the form of multiple additive terms, and the terms are monotonic functions of each variable, which particularly, are of linear, Michaelis-Menten and Hill forms. In our previous papers [7,23], we have taken such special structure properties of gene networks into account, and have shown that these genetic oscillators can be transformed into Lur'e form and can be further analyzed by using Lur'e system method in control theory [22]. In this paper, we also consider such special structure.
To make our paper more understandable and self-contained, we will first introduce the approach briefly, and after that we will analyze the stochastic genetic oscillator networks theoretically. We consider the following general form of genetic oscillator: (t) = Ay(t) + B 1 f 1 (y(t)) + B 2 f 2 (y(t)), (4) where y(t) ∈ R n represents the concentrations of proteins, RNAs and chemical complexes, A, B 1 and B 2 are matrices in R n × n , f 1 (y(t)) = [f 11 (y 1 (t)), ʜ, f 1n (y n (t))] T with f 1j (y j (t)) as a monotonic increasing function of the form , and  (4) is by no mean peculiar. Many well-known genetic system models can be represented in this form, such as the Goodwin model [24], the repressilator [25], the toggle switch [26], and the circadian oscillators [27]. Undoubtedly, this work can be easily generalized to the case of , where there are more than two nonlinear terms in each equation of the genetic oscillator. From the synthetic biology viewpoint, genetic oscillators with only linear, Michaelis-Menten and Hill terms can also be implemented experimentally.
To avoid confusion, we let the jth column of B 1,2 be zeros if f 1j,2j ≡ 0. Since and letting f(·) = f 1 (·), we can rewrite (4) as follows: Obviously, f i and g i satisfy the sector conditions: or equivalently, Recall that a Lur'e system is a linear dynamic system, feedback interconnected to a static nonlinearity that satisfies a sector condition [22]. Hence, the genetic oscillator (5) can be seen as a Lur'e system, which can be investigated by using the fruitful Lur'e system approach in control theory.
By substituting the individual genetic oscillator dynamics (5) for F in the network (1), we obtain the following network of N coupled genetic oscillators: For this network, we have the following result: , Then the network (7) is stochastic synchronization with disturbance attenuation γ.
Proposition 1 can be proved by replacing F in S 2 of (3) by the dynamics of (5), and using the sector conditions (6). The details are given in Appendix D. If we choose U beforehand, the matrix inequalities in (8) are all LMIs, which are very easy to be verified numerically [15]. For some special G and D, we can further simplify the verification process [21,7].

An example
To demonstrate the effectiveness of our theoretical results, we consider a population of N coupled biological clocks, and the individual genetic oscillator is the repressilator [25]. The repressilator is a network of three genes, the products of which inhibit the transcription of each other in a cyclic way (10). Specifically, the gene lacI expresses protein LacI, which inhibits transcription of the gene tetR.
The protein product TetR, inhibits transcription of the gene cI, the protein product CI of which in turn inhibits expression of lacI, thus forming a negative feedback cycle.
The quorum-sensing system is used for the coupling purpose, which was described in [5]. The system achieves cellto-cell communication through a mechanism that makes use of two proteins, the first one of which (LuxI), under the control of the repressilator protein LacI, synthesizes a small molecule known as an autoinducer (AI), which can diffuse freely through the cell membrane. When a second protein (LuxR) binds to this molecule, the resulting complex activates the transcription lacI, as shown in Fig. 1. The noise perturbations in the model can arise both intracellularly, due to the intrinsically noisy property of the gene regulation process, and extracellularly, due to environment fluctuations.
To model the system, we use a i , b i , c i , and A i , B i , C i to represent the dimensionless concentrations of the genes tetR, cI, lacI and their product proteins TetR, CI, LacI, respectively. As in [5], assuming equal lifetimes of the TetR and Luxl proteins, their dynamics are identical, and hence we can use the same variable to describe both protein concentrations. The concentration of AI inside the ith cell is denoted by S i . Consequently, the mRNA and protein dynamics in the ith cell can be described by [5]: where m is the Hill coefficient, S e denotes the extracellular AI concentration, and the meaning of the other parameters are standard in genetic network models. We assume that the release of the AI is fast with respect to the timescale of the oscillators and becomes approximately homogeneous to establish an average AI level outside the cells.
In the quasi-steady-state approximation, the extracellular AI concentration can be approximated by [5] where 0 <Q 0 < 1 is a constant. Thus the dynamics of S i can be rewritten as Clearly, the individual model in (9)  Gaussian white noise process. According to Proposition 1 (by letting U = -G, and using MATLAB LMI Toolbox), we know that the above all-to-all coupled network can achieve stochastic synchronization with disturbance attenuation γ = 6. Although γ is a large value, it is easy to show from (2) that the time average of E(∑ i ∑ j |x i (t) - x j (t)| 2 ) is still rather small because is very small. We omit the computational details here. In Fig. 2(a) &2(b), when starting from the same initial values, we plot the time evolution of the mRNA concentrations of tetR (a i ) of all the oscillators, which behaviors are similar to the experimental results (see, e.g. [1]). Fig. 2(c) shows the synchronization error of a i -a 1 for i = 2, ʜ, 10.
In Definition 1, it requires that all the genetic oscillators have the same initial conditions, so that V(x(0)) = 0. If the genetic oscillators have different initial conditions, V(x(0)) ≠ 0, and thus (12) in the Appendix C is replaced by Since in genetic networks, the variables usually represent the concentrations of mRNAs, proteins and chemical complexes, which are of (not so large) limited values, and so is V(x(0)). For a long time scale, the last term of the above inequality is usually much smaller than the absolute value of the first term in the right-hand side, and thus the last term can be ignored roughly. In Fig. 3, we show the same computations as those in Fig. 2 except that the oscillators are with different initial values (randomly in the range (0, 1)). After a period of evolution, the network behaviors are similar to those in Fig. 2, which verifies our above argument. In other words, rigorously, according to Definition 1, we need that all the oscillators have the same initial conditions, but practically, for oscillators with different initial conditions, we can obtain almost the same results.
For the purpose of comparison, in Fig. 4 we show the simulation results of a network without noise perturbations. As we can conclude from Figs. 2, 3, 4, the networks with noise perturbation, though can't achieve perfect synchronization, can indeed achieve synchronization with small error fluctuation, and the network behaviors are similar to those of networks without noise perturbations.

Synthesis
In addition to providing a sufficient condition for the stochastic synchronization, Proposition 1 can also be used for designing genetic oscillator networks, which is a v i Schematic representation of the coupled repressilator network Figure 1 Schematic representation of the coupled repressilator network. In the left big circle, detailed regulation and coupling mechanism are presented. The repressilator module is located at the left of the vertical dotted line, and the coupling module appears at the right. 1 cell byproduct of the main results. From the synthetic biology viewpoint, to minimize the influence of the noises (on the synchronization), we can design genetic oscillator networks according to the following rule: min γ, such that the LMIs (8) hold, (11) which is obviously from the above theoretical result. This is similar to the H ∞ synthesis problem in control theory.

Discussion and conclusion
In this paper, we presented a general theoretical method for analyzing the stochastic synchronization of coupled genetic oscillators based on systems biology approach. By taking the specific structure of genetic systems into account, a sufficient condition for the stochastic synchro-nization was derived based on LMI formalism, which can be easily verified numerically. Although the method and results are presented for genetic oscillator networks, it is also applicable to other dynamical systems. In coupled genetic oscillator networks, since there is a maximal activity of fully active promoters, it is more realistic to consider a Michaelis-Menten form of the coupling terms. As argued in [7], our theoretical method is also applicable to this case. To make the theoretical method more understandable and to avoid unnecessarily complicated notation, we discussed only on some simplified forms of the genetic oscillators, but more general cases regarding this topic can be studied in a similar way. For example: (I) The genetic oscillator model (5) can be generalized to more general Simulation results of the coupled repressilators with the same initial values As we know, noises can play both beneficial and harmful roles (for synchronization) in biological systems. For the latter case, the noise is a kind of perturbation, and it is interesting to study the robustness of the synchronization with respect to noise. In this paper, we addressed this question. For the former case, in [13,14], the authors studied the mechanisms of noise-induced synchronization.

Methods
To simulate the stocahstic differential equaiton (t) = f(x) + g(x)ξ(t), the well-known Euler-Maruyama scheme is most frequently used, which is also used in this paper. In 0 0

B. Mean-square synchronization
For network (1), a natural attempt is to study the meansquare asymptotic synchronization. Analogue to the definition of mean-square stability [19], we can define the mean-square synchronization as follows: Definition A1: The network (1) is said to be mean-square synchronous if for every ε > 0, there is a δ(ε) > 0, such that <ε for |x i (0) -x j (0)| <δ(ε), ∀i, j. If in addition, lim t→∞ E|x i (t) -x j (t)| 2 = 0 for all initial conditions, the network is said to be mean square asymptotically synchronous.
In analyzing the synchronization of the network (1), we use the Lyapunov function V(x(t)) = x T (t)(U P)x(t) [21], where is the Kronecker product, and . According to [21], this Lyapunov function is equivalent to V(x(t)) = ∑ i <j (- By Itô's formula [19], we obtain the following stochastic differential along (1) where v(t) = diag(v 1 , ʜ, v N ) ∈ R Nn × N , L is the diffusion operator, and We discuss two special cases of the stochastic terms: (y 1 -y 2 ) T P[F(y 1 ) -F(y 2 ) -T(y 1 -y 2 )] < 0, ∀ y 1 , y 2 ∈ R n (y 1 ≠ y 2 ) Hence, if there are matrices P > 0, T ∈ R n × n and U, such that the above conditions hold, the network (1) will achieve mean-square asymptotically synchronization. In this case, roughly speaking, the noise will not affect the synchronous state (since they are common for all oscillators), but it will affect the individual oscillator dynamics.
2. The noise intensity matrix v i is a function of ∑ j G ij D xj , which means that if there is no coupling from oscillator j to i, then j does not have contribution to the perturbation of oscillator i. We further assume that v i can be estimated by So, the conditions for the mean-square asymptotically synchronization of the network (1) in this case are (y 1 -y 2 ) T P[F(y 1 ) -F(y 2 ) -T(y 1 -y 2 )] < 0, If we consider genetic oscillators of the form of (5), the conditions for the mean-square asymptotically synchronization can be analyzed by the same method as that in the following Appendix D. 1 2 From Definition Al, we know that the definition of the mean-square asymptotically synchronization is rather restrictive, which requires that lim t→∞ E|x i (t) -x j (t)| 2 = 0, ∀i, j. If it is neither of the above two cases, that is neither V i d wi are the same for all i, nor v i (∀i) reduce to zero when x 1 = ʜ = x n , the network is hardly to achieve mean-square asymptotically synchronization. Experimental results also show that usually the genetic oscillators can not achieve mean-square synchronization (see for example [1]). So, we argue that the study of mean-square synchronization is unrealistic (and therefore meaningless) in genetic networks.
In Ref. [28], the authors studied the mean-square asymptotic synchronization of two master-slave coupled Chua's circuits. They assume that the noise intensity depends on the difference of the states of the two systems, which is also somewhat unrealistic.

C. Analysis of the general synchronization condition
To obtain the general synchronization condition (3) of the network (1), we also use the Lyapunov function V(x(t)) = x T (t)(U P)x(t). By Itôs formula [19], we obtain the stochastic differential dV(x(t))= LV(x(t))dt + 2x T (t)(U P)v(t)dw(t). According to Definition 1, we assume that the oscillators have the same initial conditions, thus we can derive E(V(x(t)) = . For γ > 0, we define Then, it is easy to show that for α > 0, Assuming U P ≤ ρI, and letting α = γ/ρ, we have If E(S 1 ) < 0, we will have J(t) < 0, and thus, (2) follows immediately from (3).