- Research article
- Open access
- Published:

# Stochastic synchronization of genetic oscillator networks

*BMC Systems Biology*
**volume 1**, Article number: 6 (2007)

## 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*_{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:

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*)]^{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; {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} ≡ 2(*y*_{1} - *y*_{2})^{T}*P*[*F*(*y*_{1}) - *F*(*y*_{2}) - *T*(*y*_{1} - *y*_{2})] + \frac{\rho}{\gamma}(*y*_{1} - *y*_{2})^{T}(*y*_{1} - *y*_{2}) < 0, ∀*y*_{1}, *y*_{2} ∈ *R*^{n}(*y*_{1} ≠ *y*_{2}),

(*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:

(*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 {f}_{1j}({y}_{j}(t))={({y}_{j}(t)/{\beta}_{1j})}^{{H}_{1j}}/[1+{({y}_{j}(t)/{\beta}_{1j})}^{{H}_{1j}}], and *f*_{2}(*y*(*t*)) = [*f*_{21}(*y*_{1}(*t*)), ⋯ , *f*_{2n}(*y*_{
n
}(*t*))]^{T}with *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*_{1j}and *H*_{2j}are 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*_{1}(·), we can rewrite (4) as follows:

(*t*) = *Ay*(*t*) + *B*_{1}*f*(*y*(*t*)) - *B*_{2}*g*(*y*(*t*)) + *B*_{2}. (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, Λ_{1} = diag(λ_{11}, ⋯,λ_{1n}) > 0, Λ_{2} = diag(λ_{21}, ⋯,λ_{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*^{T}*P* - *Q* - *Q*^{T}+ \frac{\rho}{\gamma}*I*, *K*_{1} = diag(*k*_{11}, ⋯, *k*_{1n}), *K*_{2} = diag(*k*_{21}, ⋯, *k*_{2n}). 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 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.

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*_{0} < 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*_{1} is a 7 × 7 matrix with all zero entries except for B_{1}(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*_{2} is a 7 × 7 matrix with all zero entries except for *B*_{2}(1,6) = *B*_{2}(2, 4) = *B*_{2}(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*_{1} = *d*_{2} = *d*_{3} = 0.4, *μ* = 1.3, *k* = 5, *μ*_{
s
}= 5, *d*_{4} = *d*_{5} = *d*_{6} = 0.5, *β*_{1} = *β*_{2} = *β*_{3} = 0.2, *k*_{s 0}= 0.016, *k*_{s 1}= 0.018, *Q*_{0} = 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*)|^{2}) 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*_{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

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.

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 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*^{n}(*a* ≠ *b*), *i* = 1, ⋯ , *n*, where *c*_{1i}, *c*_{2i}∈ *R*^{n}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*^{T}denotes the transpose of a square matrix *A*. The notation *M* > (<) 0 is used to define a real symmetric positive definite (negative definite) matrix. *R*^{m}denotes 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*_{2}[0, ∞) is the space of square-integrable vector functions over [0, ∞); |·| stands for the Euclidean vector norm, and ||·||_{2} stands for the usual *L*_{2}[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\\ \vdots & \ddots & \vdots \\ {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*)|^{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 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*))^{T}*P*(*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* + 2*x*^{T}(*t*)(*U* ⊗ *P*)*v*(*t*)*dw*(*t*)

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

*LV*(*x*(*t*)) = 2∑_{i <j}(-*U*_{
ij
})(*x*_{
i
}- *x*_{
j
})^{T}*P*[*F*(*x*_{
i
}) - *F*(*x*_{
j
}) - *T*(*x*_{
i
}- *x*_{
j
})]

+ 2*x*^{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*_{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})

(*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*_{1}, ⋯, *v*_{
N
}) ∈ *R*^{Nn × N}, *w*(*t*) = [*w*_{1}(*t*), ⋯, *w*_{
N
}(*t*)]^{T}, and *H* = diag(*H*_{1}, ⋯, *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*)^{T}*H*(*G* ⊗ *D*)*x*(*t*).

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,

∀*y*_{1}, *y*_{2} ∈ *R*^{n}(*y*_{1} ≠ *y*_{2})

(*U* ⊗ *P*)(*G* ⊗ *D* + *I* ⊗ *T*) + (*G* ⊗ *D* + *I* ⊗ *T*)^{T}(*U* ⊗ *P*) + *ρ*(*G* ⊗ *D*)^{T}*H*(*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*)|^{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* + 2*x*^{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))]\\ +2{x}^{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))]\\ +2{x}^{T}(t)(U\otimes P)(G\otimes D+I\otimes T)x(t)\}.\end{array}

If **E**(*S*_{1}) < 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*_{2} 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*_{1}, we have *S*_{2} ≤ *ξ*(*t*)*M*_{1}*ξ*(*t*) < 0 for all *y*_{1}, *y*_{2} except for *y*_{1} = *y*_{2}, where *ξ*(*t*) = [(*y*_{1}(*t*) - *y*_{2}(*t*))^{T}, (*f*(*y*_{1}(*t*) - *f*(*y*_{2}(*t*)))^{T}, (*g*(*y*_{1}(*t*) - *g*(*y*_{2}(*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.

## References

Yamaguchi S, Isejima H, Matsuo T, Okura R, Yagita K, Kobayashi M, Okamura H: Synchronization of cellular clocks in the suprachiasmatic nucleus. Science. 2003, 302: 1408-1412.

McMillen D, Kopell N, Hasty J, Collins JJ: Synchronizing genetic relaxation oscillators by intercell signalling. Proc Natl Acad Sci USA. 2002, 99: 679-684.

Wang R, Chen L: Synchronizing genetic oscillators by signaling molecules. J Biol Rhythms. 2005, 20: 257-269.

Kuznetsov A, Kaern M, Kopell N: Synchronization in a population of hysteresis-based genetic oscillators. SIAM J Appl Math. 2004, 65: 392-425. 10.1137/S0036139903436029.

Garcia-Ojalvo J, Elowitz MB, Strogatz SH: Modelling a synthetic multicellular clock: Repressilators coupled by quorum sensing. Proc Natl Acad Sci USA. 2004, 101: 10955-10960.

Gonze D, Bernard S, Waltermann C, Kramer A, Herzel H: Spontaneuous synchronization of coupled circadian oscillators. Biophys J. 2005, 89: 120-129.

Li C, Chen L, Aihara K: Synchronization of coupled nonidentical genetic oscillators. Phys Biol. 2006, 3: 37-44.

Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Science. 2002, 297: 1183-1186.

Paulsson J: Summing up the noise in gene networks. Nature. 2004, 427: 415-418.

Raser JM, O'Shea EK: Noise in gene expression: Origins, consequences, and control. Science. 2005, 309: 2010-2013.

Blake WJ, Kaern M, Cantor CR, Collins JJ: Noise in eukaryotic gene expression. Nature. 2003, 422: 633-637.

Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: From theories to pheno-types. Nature Reviews Genetics. 2005, 6: 451-464.

Chen L, Wang R, Zhou T, Aihara K: Noise-induced cooperative behavior in a multicell system. Bioinfor. 2005, 22: 2722-2729. 10.1093/bioinformatics/bti392.

Li C, Chen L, Aihara K: Transient resetting: A novel mechanism for synchrony and its biological examples. PLoS Comp Biol. 2006, 2: e103-10.1371/journal.pcbi.0020103.

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

Tong AH, Lesage G, Bader GD, Ding H, Xu H, Xin X, Young J, Berriz GF, Brost RL, Chang M, Chen Y, Cheng X, Chua G, Friesen H, Goldberg DS, Haynes J, Humphries C, He G, Hussein S, Ke L, Krogan N, Li Z, Levinson JN, Lu H, Menard P, Munyana C, Parsons AB, Ryan O, Tonikian R, Roberts T, Sdicu AM, Shapiro J, Sheikh B, Suter B, Wong SL, Zhang LV, Zhu H, Burd CG, Munro S, Sander C, Rine J, Greenblatt J, Peter M, Bretscher A, Bell G, Roth FP, Brown GW, Andrews B, Bussey H, Boone C: Global mapping of the yeast genetic interaction network. Science. 2004, 303: 808-813.

Barabási AL, Oltvai ZN: Network biology: Understanding the cell's functional organization. Nature Reviews Genetics. 2004, 5: 101-114.

Kloeden PE, Platen E, Schurz H: Numerical solution of SDE through computer experiments. 1994, Springer-Verlag: Berlin

Arnold L: Stochastic Differential Equations: Theory and Applications. 1974, John Wiley & Sons: London

Xu S, Chen T: Robust

*H*_{∞}control for the uncertain stochastic systems with state delay. IEEE Trans Auto Contr. 2002, 47: 2089-2094. 10.1109/TAC.2002.805670.Wu CW: Synchronization in Coupled Chaotic Circuits and Systems. 2002, World Scientific: Singapore

Vidyasagar M: Nonlinear Systems Analysis. 1993, Englewood Cliffs, NJ: Prentice-Hall

Li C, Chen L, Aihara K: Stability of genetic networks with SUM regulatory logic: Lur'e system and LMI approach. IEEE Trans. Circuits and Systems-I. 2006, 53: 2451-2458. 10.1109/TCSI.2006.883882.

Goodwin BC: Oscillatory behavior in enzymatic control processes. Adv Enzyme Regul. 1965, 3: 425-438.

Elowitz MB, Leibler S: A synthetic oscillatory network of transcriptional regulators. Nature. 2000, 403: 335-338.

Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in

*Escherichia coli*. Nature. 2000, 403: 339-342.Goldbeter A: A model for circadian oscillations in the

*Drosophila*period protein (PER). Proc R Soc Lond B: Biol Sci. 1995, 261: 319-324. 10.1098/rspb.1995.0153.Lin W, He Y: Complete synchronization of the noise-perturbed Chua's circuits. Chaos. 2005, 15: 023705-10.1063/1.1938627.

## Acknowledgements

This research was supported by Grant-in-Aid for Scientific Research on Priority Areas 17022012 from MEXT of Japan, the Fok Ying Tung Education Foundation under Grant 101064, the Program for New Century Excellent Talents in University, the Distinguished Youth Foundation of Sichuan Province, and the National Natural Science Foundation of China (NSFC) under grant 60502009.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### Authors' contributions

CL gave the topic, developed theoretical results, designed the numerical experiments, analyzed the data, contributed materials/analysis tools. CL, LC and KA wrote the paper.

## Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

## Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## About this article

### Cite this article

Li, C., Chen, L. & Aihara, K. Stochastic synchronization of genetic oscillator networks.
*BMC Syst Biol* **1**, 6 (2007). https://doi.org/10.1186/1752-0509-1-6

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/1752-0509-1-6