Stochastic synchronization of genetic oscillator networks

Abstract

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 suprachiasmatic nucleus (SCN); in [2–4], the synchronization are studied in biological networks of identical genetic oscillators; and in [5–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–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 scale-free 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.

Results

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

$\frac{d{x}_i(t)}{dt}=F(x_i(t))+{\displaystyle \sum _{j=1}^N{G}_{ij}}D{x}_j(t)+v_i(t)n_i(t),i=1,\cdots ,N$

where F(·) defines the dynamics of individual oscillators, v _i (t) ∈ R^{n × 1}is called the noise intensity vector, belongs to L₂ [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:

$d{x}_i(t)=\left[F(x_i(t))+{\displaystyle \sum _{j=1}^N{G}_{ij}D{x}_j}(t)\right]dt+v_i(t)d{w}_i(t),i=1,\cdots ,N.\left(1\right)$

The work can be easily extended to the case that v _i (t) ∈ $R^{n\times l_i}$ and n _i (t) = [n_{i 1}(t), ⋯, $n_{i{l}_i}$(t)]^Tbe 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; $G_{ii}=-{\displaystyle {\sum }_{j=1,j\ne i}^N{G}_{ij}}$. Matrix G defines the coupling topology, direction, and the coupling strength of the network.

For network (1), a natural attempt is to study the mean-square 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:

Definition 1: For a given scalar γ > 0, the network (1) is said to be stochastically synchronous (under the combination matrix U) with disturbance attenuation γ if the network without disturbance (v _i = 0, ∀i) is asymptotically synchronous, and under the same initial conditions for all oscillators,

${\displaystyle \sum _{i<j}(-U_{ij}){\Vert x_i(t)-x_j(t)\Vert }_{E_2}^2<\gamma {\displaystyle \sum _i{\Vert v_i(t)\Vert }_2^2}},\left(2\right)$

for all nonzero v _i (t), where ${\Vert \cdot \Vert }_{E_2}^2=E\left({\displaystyle {\int }_0^{\infty }{|\cdot |}^{2}dt}\right)$. Here, by introducing the combination matrix U, we can flexibly select the form of the matrix to obtain different error combinations.

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(y₁ - y₂)^TP[F(y₁) - F(y₂) - T(y₁ - y₂)] + $\frac{\rho }{\gamma }$(y₁ - y₂)^T(y₁ - y₂) < 0, ∀y₁, y₂ ∈ Rⁿ(y₁ ≠ y₂),

(U ⊗ P)(G ⊗ D + I ⊗ T) + (G ⊗ D + I ⊗ T)^T(U ⊗ P) ≤ 0,     (3)

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:

$$\dot{y}$$

(1)

(t) = Ay(t) + B₁f₁(y(t)) + B₂f₂(y(t)),   (4)

where y(t) ∈ Rⁿrepresents the concentrations of proteins, RNAs and chemical complexes, A, B₁ and B₂ are matrices in R^{n × n}, f₁(y(t)) = [f₁₁(y₁(t)), ⋯, f_1n(y _n (t))]^Twith f_1j(y _j (t)) as a monotonic increasing function of the form $f_{1j}(y_j(t))={(y_j(t)/{\beta }_{1j})}^{H_{1j}}/[1+{(y_j(t)/{\beta }_{1j})}^{H_{1j}}]$, and f₂(y(t)) = [f₂₁(y₁(t)), ⋯ , f_2n(y _n (t))]^Twith f_2j(y _j (t)) as a monotonic decreasing function of the form $f_{2j}(y_j(t))=1/[1+{(y_j(t)/{\beta }_{2j})}^{H_{2j}}]$, where H_1jand H_2jare the Hill coefficients. Genetic oscillators of the form (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 $\dot{y}(t)=Ay(t)+{\displaystyle {\sum }_{j=1}^l{B}_j{f}_j}(y(t))$, 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 j th column of B_1,2 be zeros if f_1j,2j≡ 0. Since $f_{2j}(y_j(t))=\frac{1}{1+{(y_j(t)/{\beta }_{2j})}^{H_{2j}}}=1-\frac{{(y_j(t)/{\beta }_{2j})}^{H_{2j}}}{1+{(y_j(t)/{\beta }_{2j})}^{H_{2j}}}\equiv 1-g_j(y_j(t))$, and letting f(·) = f₁(·), we can rewrite (4) as follows:

$$\dot{y}$$

(2)

(t) = Ay(t) + B₁f(y(t)) - B₂g(y(t)) + B₂.    (5)

Obviously, f _i and g _i satisfy the sector conditions: $0\le \frac{f_i(a)-f_i(b)}{a-b}<k_{1i},0\le \frac{g_i(a)-g_i(b)}{a-b}\le k_{2i}$ or equivalently,

(f _i (a) - f _i (b)) [(f _i (a) - f _i (b)) - k_1i(a - b)] ≤ 0,   (6)

(g _i (a) - g _i (b)) [(g _i (a) - g _i (b)) - k_2i(a - b)] ≤ 0,

∀a, b ∈ R(a ≠ b);i = 1, ⋯, n,

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:

$d{x}_i(t)=[A{x}_i(t)+B_{1}f(x_i(t))-B_{2}g(x_i(t))+B_2+{\displaystyle {\sum }_{j=1}^N{G}_{ij}}D{x}_j(t)]dt+v_i(t)d{w}_i(t),i=1,\cdots ,N.\left(7\right)$

For this network, we have the following result:

Proposition 1: If there are matrices P > 0, Λ₁ = diag(λ₁₁, ⋯,λ_1n) > 0, Λ₂ = diag(λ₂₁, ⋯,λ_2n) > 0, Q ∈ R^{n × n}, U as defined above, and a positive real constant ρ such that the following matrix inequalities hold

$\begin{array}{l}\left[\begin{array}{ccc}(1,1)& P{B}_1+K_1{\Lambda }_1& -P{B}_2+K_2{\Lambda }_2\\ B_1^{T}P+K_1{\Lambda }_1& -2{\Lambda }_1& 0\\ -B_2^{T}P+K_2{\Lambda }_2& 0& -2{\Lambda }_2\end{array}\right]<0,\hfill \\ {(UG\otimes PD+U\otimes Q)}^T+(UG\otimes PD+U\otimes Q)\le 0,\hfill \\ U\otimes P\le \rho I,\hfill \end{array}\left(8\right)$

where (1,1) = PA + A^TP - Q - Q^T+ $\frac{\rho }{\gamma }$I, K₁ = diag(k₁₁, ⋯, k_1n), K₂ = diag(k₂₁, ⋯, k_2n). Then the network (7) is stochastic synchronization with disturbance attenuation γ.

Proposition 1 can be proved by replacing F in S₂ 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 cell-to-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.

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.

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 i th cell is denoted by S _i . Consequently, the mRNA and protein dynamics in the i th cell can be described by [5]:

$\begin{array}{c}\frac{d{a}_i}{dt}=-d_1{a}_i+\frac{\alpha }{\mu +C_i^m},\\ \frac{d{b}_i}{dt}=-d_2{b}_i+\frac{\alpha }{\mu +A_i^m},\\ \frac{d{c}_i}{dt}=-d_3{c}_i+\frac{\alpha }{\mu +B_i^m}+\frac{k{S}_i}{{\mu }_s+S_i},\\ \frac{d{A}_i}{dt}=-d_4{A}_i+{\beta }_1{a}_i,\\ \frac{d{B}_i}{dt}=-d_5{B}_i+{\beta }_2{b}_i,\\ \frac{d{C}_i}{dt}=-d_6{C}_i+{\beta }_3{c}_i,\\ \frac{d{S}_i}{dt}=-k_{s0}{S}_i+k_{s1}{A}_i+\eta (S_e-S_i),\end{array}\left(9\right)$

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]

$S_e=Q_0\overline{S}=\frac{Q_0}{N}{\displaystyle \sum _{j=1}^N{S}_j},$

where 0 <Q₀ < 1 is a constant. Thus the dynamics of S _i can be rewritten as

$\frac{d{S}_i}{dt}=-[k_{s0}+(1-Q_0)\eta ]S_i+k_{s1}{A}_i+\frac{\eta Q_0}{N}{\displaystyle \sum _{j=1}^N(S_j-S_i)}\left(10\right)$

Clearly, the individual model in (9) is of the form (5), in which f = [0, 0, 0, 0, 0, 0, S _i /(μ _s + S _i )]^T, B₁ is a 7 × 7 matrix with all zero entries except for B₁(3, 7) = k, g = [0, 0, 0, $A_i^m$/(μ + $A_i^m$), $B_i^m$/(μ +$B_i^m$), $C_i^m$/(μ + $C_i^m$), 0]^T, B₂ is a 7 × 7 matrix with all zero entries except for B₂(1,6) = B₂(2, 4) = B₂(3, 5) = α/μ, and all the other terms are in linear form. Obviously, the coupling term can also be written into the form defined previously.

The purpose of this example is to demonstrate the effectiveness and correctness of the theoretical result, instead of mimicking the real biological clock system. We consider a small size of network with N = 10 coupled oscillators. The parameters are set as m = 4, α = 1.8, d₁ = d₂ = d₃ = 0.4, μ = 1.3, k = 5, μ _s = 5, d₄ = d₅ = d₆ = 0.5, β₁ = β₂ = β₃ = 0.2, k_{s 0}= 0.016, k_{s 1}= 0.018, Q₀ = 0.8 and η = 0.4. Since Proposition 1 holds no matter what v _i (t) is and no matter where it is introduced, and the verification of Proposition 1 is independent of noise intensity v _i , for simplicity, we set v _i = 0.015 as a scalar for all i, and the noise term v _i n _i (t) is added to the first equation in (9), where n _i (t) is a scalar 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)|²) is still rather small because $v_i^2$ 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₁ 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

$J(t)\le E\{{\displaystyle {\int }_0^t[\alpha LV}(x(s))+{\displaystyle {\sum }_{i<j}(-U_{ij}){(x_i(s)-x_j(s))}^T(x_i(s)-x_j(s))-\gamma {\displaystyle {\sum }_i{v}_i^T(s)v_i(s)]ds}}\}+E(V(x(0)).$

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.

Figure 2

Simulation results of the coupled repressilators with the same initial values. (a) The evolution dynamics of the mRNA concentrations of tetR (a _i ) of all the genetic oscillators. (b) Zooming in the range t ∈ [600, 700] of (a). (c) The evolution of the synchronization error of a _i - a₁ for i = 2, ⋯ ,10.

Figure 3

Coupled genetic oscillators with different initial conditions: (a) The evolution dynamics of the mRNA concentrations of tetR (a _i ) of all the genetic oscillators. (b) Zooming in the range t ∈ [600, 700] of (a). (c) The evolution of the synchronization error of a _i - a₁ for i = 2, ⋯ ,10.

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.

Figure 4

Coupled genetic oscillators without noise perturbation: (a) The evolution dynamics of the mRNA concentrations of tetR (a _i ) of all the genetic oscillators. (b) Zooming in the range t ∈ [600, 700] of (a). (c) The evolution of the synchronization error of a _i - a₁ for i = 2, ⋯ ,10.

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 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 synchronization 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 case such that f _i , g _i , the component of f(y(t)), g(y(t)), are functions of y(t), not only of y _i (t), and f and g can also be of non-Hill form, provided that $0\le \frac{f_i(a)-f_i(b)}{c_{1i}^T(a-b)}\le k_{1i},0\le \frac{g_i(a)-g_i(b)}{c_{2i}^T(a-b)}\le k_{2i}$, ∀a, b ∈ Rⁿ(a ≠ b), i = 1, ⋯ , n, where c_1i, c_2i∈ Rⁿare real vectors. (II) Biologically, the genetic oscillators are usually nonidentical. We can consider genetic networks with both parametric mismatches and stochastic perturbations in similar ways as those presented in this paper and [7]. (III) There are significant time delays in the gene regulation, due to the slow processes of transcription, translation and translocation. Our result can be easily extended to the case that there are delays both in the coupling and the individual genetic oscillators.

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 $\dot{x}$(t) = f(x) + g(x)ξ(t), the well-known Euler-Maruyama scheme is most frequently used, which is also used in this paper. In this scheme, the numerical trajectory is generated by x_n+1= x _n + hf(x _n ) + $\sqrt{h}$g(x _n )η _n , where h is the time step and η _n is a discrete time Gaussian white noise with <η _n >= 0 and <η _n η _m > = δ _nm . For more details, see e.g. [18].

Appendices

A. Notations

Throughout this paper, A^Tdenotes the transpose of a square matrix A. The notation M > (<) 0 is used to define a real symmetric positive definite (negative definite) matrix. R^mdenotes the m-dimensional Euclidean space; and R^{n × m}denotes the set of all n × m real matrices. In this paper, if not explicitly stated, matrices are assumed to have compatible dimensions. E(·) denotes the expectation operator; L₂[0, ∞) is the space of square-integrable vector functions over [0, ∞); |·| stands for the Euclidean vector norm, and ||·||₂ stands for the usual L₂[0, ∞) norm. The Kronecker product A ⊗ B of an n × m matrix A and a p × q matrix B is the np × mq matrix defined as

$A\otimes B=\left[\begin{array}{ccc}{A}_{11}B& \cdots & A_{1m}B\\ ⋮& \ddots & ⋮\\ A_{n1}B& \cdots & A_{nm}B\end{array}\right].$

For a general stochastic systems

dx = f(t, x(t))dt + g(x(t))dw(t),

the diffusion operator L acting on Y(t, x(t)) is defined by

LY(t, x(t)) = Y _t (t, x(t)) + Y _x (t, x(t))f(t, x(t)) + $\frac{1}{2}$trace[g(x(t))g^T(x(t))Y _xx (t, x(t))].

B. Mean-square synchronization

For network (1), a natural attempt is to study the mean-square 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 $E{\left|x_i(t)-x_j(t)\right|}_{t>0}^2$ <ε for |x _i (0) - x _j (0)| <δ(ε), ∀i, j. If in addition, lim_t→∞E|x _i (t) - x _j (t)|² = 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 $x(t)={[x_1^T(t),\cdots ,x_n^T(t)]}^T\in R^{Nn\times 1}$. According to [21], this Lyapunov function is equivalent to V(x(t)) = ∑_{i <j}(-U _ij )(x _i (t) - x _j (t))^TP(x _i (t) - x _j (t)).

By Itô's formula [19], we obtain the following stochastic differential along (1)

dV(x(t)) = LV(x(t))dt + 2x^T(t)(U ⊗ P)v(t)dw(t)

where v(t) = diag(v₁, ⋯, v _N ) ∈ R^{Nn × N}, L is the diffusion operator, and

LV(x(t)) = 2∑_{i <j}(-U _ij )(x _i - x _j )^TP[F(x _i ) - F(x _j ) - T(x _i - x _j )]

+ 2x^T(t)(U ⊗ P)(G ⊗ D + I ⊗ T)x(t)

+ trace(v(t)v^T(t)(U ⊗ P))

We discuss two special cases of the stochastic terms:

1. The genetic oscillators are perturbed by the same noise, which can occur in the situation that genetic oscillators communicate via a common environment. In this case, v _i d _wi are the same for all i. We let [$v={[v_i^T,\cdots ,v_N^T]}^T$] and dw = dw _i Since U is a matrix with zero row sums and v _i is the same for all i, it is easy to show that the last term of LV is zero. Thus if the following conditions hold, we will have E[dV(x(t))] = E[LV(x(t))dt] < 0.

(y₁ - y₂)^TP[F(y₁) - F(y₂) - T(y₁ - y₂)] < 0,

∀ y₁, y₂ ∈ Rⁿ(y₁ ≠ y₂)

(U ⊗ P)(G ⊗ D + I ⊗ T) + (G ⊗ D + I ⊗ T)^T(U ⊗ P) ≤ 0.

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

$v_i^T{v}_i\le {\left[{\displaystyle \sum _j{G}_{ij}D{x}_j(t)}\right]}^T{H}_i\left[{\displaystyle \sum _j{G}_{ij}D{x}_j(t)}\right],H_i\ge 0.$

Defining v = diag(v₁, ⋯, v _N ) ∈ R^{Nn × N}, w(t) = [w₁(t), ⋯, w _N (t)]^T, and H = diag(H₁, ⋯, H _N ), and assuming U ⊗ P ≤ ρI, we have

trace(v(t)v^T(t)(U ⊗ P)) ≤ λ _max (U ⊗ P)trace(v(t)v^T(t))

≤ ρ∑ _i $v_i^T$(t)v _i (t)

≤ ρx^T(t) (G ⊗ D)^TH(G ⊗ D)x(t).

So, the conditions for the mean-square asymptotically synchronization of the network (1) in this case are

(y₁ - y₂)^TP[F(y₁) - F(y₂) - T(y₁ - y₂)] < 0,

∀y₁, y₂ ∈ Rⁿ(y₁ ≠ y₂)

(U ⊗ P)(G ⊗ D + I ⊗ T) + (G ⊗ D + I ⊗ T)^T(U ⊗ P) + ρ(G ⊗ D)^TH(G ⊗ D) ≤ 0,

U ⊗ P ≤ ρI.

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.

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)|² = 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₁ = ⋯ = 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)) = $\left({\displaystyle {\int }_0^{t}LV}(x(s))ds\right)$. For γ > 0, we define

$J(t)=E{\{{\displaystyle {\int }_0^t[{\displaystyle \sum _{i<j}(-U_{ij})(x_i(s)}}-x_j(s))}^T(x_i(s)-x_j(s))-\gamma {\displaystyle \sum _i{v}_i^T(s)v_i(s)]ds\}}\left(12\right)$

Then, it is easy to show that for α > 0,

$\begin{array}{c}J(t)\le E{\{{\displaystyle {\int }_0^t[\alpha LV(x(s))+{\displaystyle {\sum }_{i<j}(-U_{ij})(x_i(s)-x_j(s)}})}^T\\ \cdot (x_i(s)-x_j(s))-\gamma {\displaystyle {\sum }_i{v}_i^T}(s)v_i(s)]ds\}\\ \equiv E\{{\displaystyle {\int }_0^t{S}_1}(s)ds\}.\end{array}\left(13\right)$

Assuming U ⊗ P ≤ ρI, and letting α = γ/ρ, we have

$\begin{array}{ccc}{S}_1& \equiv & \frac{\gamma }{\rho }\{2{\displaystyle {\sum }_{i<j}(-U_{ij}){(x_i(t)-x_j(t))}^T}\\ \cdot P[F(x_i(t))-F(x_j(t))-T(x_i(t)-x_j(t))]\\ +2x^T(t)(U\otimes P)(G\otimes D+I\otimes T)x(t)\\ +\text{trace}(v(t)v^T(t)(U\otimes P))\\ +\frac{\rho }{\gamma }{{\displaystyle {\sum }_{i<j}(-U_{ij})(x_i(t)-x_j(t))}}^T(x_i(t)-x_j(t))\\ -\rho {\displaystyle {\sum }_i{v}_i^T}(t)v_i(t)\}\\ \le & \frac{\gamma }{\rho }\{{\displaystyle {\sum }_{i<j}(-U_{ij})[2{(x_i(t)-x_j(t))}^T}\\ \cdot P(F(x_i)-F(x_j)-T(x_i-x_j))\\ +\frac{\rho }{\gamma }{(x_i(t)-x_j(t))}^T(x_i(t)-x_j(t))]\\ +2x^T(t)(U\otimes P)(G\otimes D+I\otimes T)x(t)\}.\end{array}$

If E(S₁) < 0, we will have J(t) < 0, and thus, (2) follows immediately from (3).

D. Proof of proposition 1

Proposition 1 can be proved by replacing F in S₂ of (3) by the dynamics of (5), and using the sector conditions (6). We have

$\begin{array}{ccc}{S}_2& \equiv & 2{(y_1(t)-y_2(t))}^{T}P[(A-T)(y_1(t)-y_2(t))\\ +B_1(f(y_1(t))-f(y_2(t)))-B_2(g(y_1(t))-g(y_2(t)))]\\ +\frac{\rho }{\gamma }{(y_1(t)-y_2(t))}^T(y_1(t)-y_2(t))\\ \le & 2{(y_1(t)-y_2(t))}^{T}P(A-T)(y_1(t)-y_2(t))\\ +\frac{\rho }{\gamma }{(y_1(t)-y_2(t))}^T(y_1(t)-y_2(t))\\ +2{(y_1(t)-y_2(t))}^{T}P{B}_1(f(y_1(t))-f(y_2(t)))\\ -2{(y_1(t)-y_2(t))}^{T}P{B}_2(g(y_1(t))-g(y_2(t)))\\ -2{\displaystyle {\sum }_{l=1}^n{\lambda }_{1l}}(f_l(y_{1l}(t))-f_l(y_{2l}(t)))\\ \cdot [(f_l(y_{1l}(t))-f_l(y_{2l}(t)))-k_1(y_{1l}(t)-y_{2l}(t))]\\ -2{\displaystyle {\sum }_{l=1}^n{\lambda }_{2l}(g_l(y_{1l}(t))-g_l(y_{2l}(t)))}\\ \cdot [(g_l(y_{1l}(t))-g_l(y_{2l}(t)))-k_2(y_{1l}(t)-y_{2l}(t))].\end{array}$

By letting Q = PT and denoting the first matrix in (8) by M₁, we have S₂ ≤ ξ(t)M₁ξ(t) < 0 for all y₁, y₂ except for y₁ = y₂, where ξ(t) = [(y₁(t) - y₂(t))^T, (f(y₁(t) - f(y₂(t)))^T, (g(y₁(t) - g(y₂(t)))^T]^T∈ R^{3n × 1}. So, the first condition in (3) is satisfied. Substituting Q = PT, the second inequality in (3) is equivalent to the second inequality in (8). Thus, Proposition 1 is proved.