Microscopic interactions lead to mutual synchronization in a network of networks

Physics Letters A 375 (2011) 2809–2814

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Microscopic interactions lead to mutual synchronization in a network of networks Yao-Chen Hung Department of Physics, National Chung Cheng University, Chiayi 621, Taiwan

a r t i c l e

i n f o

Article history: Received 27 January 2011 Received in revised form 23 May 2011 Accepted 24 May 2011 Available online 22 June 2011 Communicated by C.R. Doering Keywords: Boolean networks Synchronization

a b s t r a c t This Letter proposes a stochastic coupling scheme to study the collective dynamics in a network that comprises random Boolean networks. Based on microscopic interactions, which are understood as the exchange of information among nodes, mutual synchronization can be achieved when the product of the assigning probability and influence probability exceeds a critical threshold. A mean field model is developed to approximate the dynamical behaviors of the original system. The effect of finite system size can be further mimicked by incorporating a noise term into the model. The dependence of the synchronization threshold on the degrees of connectivity and coupling configuration is analyzed. © 2011 Elsevier B.V. All rights reserved.

1. Introduction In past decades, synchronization phenomena have attracted considerable attention in a wide variety of disciplines, including physics [1–3], mathematics [4], biology [5], physiology [6,7], ecology [8], chemistry [9], and others [10]. Depending on the correlation among subsystems, synchronization can be classified into three major categories, i.e., complete synchronization (CS), phase synchronization (PS), and generalized synchronization (GS). CS means that the interacting subsystems adjust their states and finally converge to an identical trajectory [11–13]. PS means that the phases of the interacting subsystems are synchronized while the amplitudes still remain uncorrelated [14,15]. GS means the emergence of some functional relations between two interacting subsystems, i.e., xi (t ) = F[x j (t )], where xi (t ) and x j (t ) are the trajectories of the two interacting systems [16]. Other types of synchronization include lag synchronization [17], anticipation synchronization [18], and reduced order synchronization [19]. Nowadays, these definitions are regarded as different degrees of realization of a universal concept of synchronization [10]. Recently, interest has extended to synchronization in complex networks [20–22], with a focus on the emergence of collective behaviors under the influences of various topological structures [23,24]. Such studies are important because synchronization among all subsystems is required in biological systems when executing certain functions [25]. The emergence of clustering or partial synchronization, in which some of the subsystems synchronize but the others do not, is another particular feature represented in complex networks [26–30]. A good understanding of the transitions from turbulence, through partial synchronization, to mu-

E-mail address: [email protected]. 0375-9601/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2011.05.052

tual synchronization would be beneficial to understanding selforganization phenomena in nature [31], with potential implications for neuroscience, especially related to epilepsy [32]. Although synchronization in complex networks has been well examined, some essential problems remain to be investigated. For instance, previous studies have tended to focus on networks of dynamical units, which are usually formulated using differential equations [20,25] or iterated maps [27,30], and which are therefore coupled deterministically. Notwithstanding the fact that macroscopic descriptions can capture some of the dynamical features of the system, how coherent behaviors emerge from microscopic iterations remains to be understood. Exploring this problem is critical because dynamical processes in many real systems are inherently discrete and stochastic [33,34], and the stochasticity of microscopic processes may lead to novel behaviors [35]. Cellular automata (CA) are excellent models that mimic microscopic mechanisms that lead to collective behaviors in various systems [36]. However, because of their fully discrete nature, customary coupling methods cannot be employed to study synchronization in mutually interacting CA. Morelli and Zanette carried out pioneering work to solve the aforementioned problem [37]. They proposed a novel coupling mechanism, the stochastic coupling scheme, to explore the synchronization transition in a pair of coupled CA. The scheme has been successfully applied to other discrete systems, such as Kauffaman networks (KFs) [38] and random Boolean networks (RBNs) [39,40], which can abstractly represent the dynamical properties of genetic networks. The results thus obtained therefore provide useful information on the co-evolution of biological species with a horizontal exchange of genes. Along these lines, this Letter extends the stochastic coupling scheme to study the collective dynamics of more than two coupled RBNs, or a network of RBNs. This work is crucial owing to the ubiquity of multicellular interactions observed in nature [25].

2810

Y.-C. Hung / Physics Letters A 375 (2011) 2809–2814

We numerically demonstrate that mutual synchronization can be achieved using the modified coupling scheme. To understand the underlying mechanism and predict the dynamical behaviors, a mean-field model that is based on the concepts of statistical probability is formulated. Transitions from the unsynchronized state to the synchronized state are analyzed under the conditions of unfrozen coupling configuration and frozen coupling configuration. The effects of finite size and degree of connectivity on the synchronization threshold are also studied. This Letter is organized as follows. The next section describes the system of interest and the observed synchronization phenomena. In Section 3, a mean-field model that consists of density maps is constructed to describe the dynamics of the original system. Moreover, a refined model by including a noise term into the mean-field model, is proposed to mimic the effect of finite system size. Section 4 exhaustively analyzes the mutual synchronization that is revealed by the system. Finally, a brief conclusion is drawn and suggestions for further work are given. 2. Model description Consider M coupled RBNs that each consist of N nodes. Each node cni , where the superscript i = 1, 2, . . . , M is the network index and the subscript n = 1, 2, . . . , N is the node index, is described by only two states, 0 or 1. Suppose the evolution of all cells c˜ is governed by two operators,





ˆ c˜ (t ) , c˜ (t + 1) = Sˆ ◦ R

(1)

in which c˜ is the state vector of the studied system. The operator ˆ models intra-network interactions between nodes, while Sˆ repR resents inter-network interactions based on a stochastic coupling scheme. In this work, the states of all nodes are updated simultaneously. Within an individual RBN, a node cni is randomly assigned a set of neighbors which dominate the state of cni based upon the

ˆ The number of neighbor nodes (or namely deBoolean rule R. i gree) of cn is denoted by kni . At each time step, the neighbors of cni are re-determined randomly, but the value of kni keeps fixed dur-

ˆ stands for the independent evolution ing iteration. The operator R of each RBN. Here we adopt a totalistic rule, the generalized ECA rule 126, as the working rule for its potential implication in biological evolution and rich dynamics (especially chaos) [41–43]. The generalized ECA rule 126 reads: if the state of cni and its neighbors

ˆ otherwise cni = 1. are all 1 or 0, cni = 0 after the operation of R; Other totalistic rules that reveal biological implications and chaotic motions, such as the rule 22 [44], can also be applied to study the emerged collective dynamics. The effects of different evolving rules on mutual synchronization will be reported elsewhere. The operator Sˆ governs the process of information exchange between nodes in adjacent RBNs [38]. Here the algorithm is generalized to couple more than two RBNs. At each time step, the ith RBN is randomly assigned a set of neighboring RBNs Q i containing ˆ each node cni within qi neighboring RBNs. Under the operation of S, j

the ith RBN is randomly assigned a driver cm from its neighboring RBNs ( j ∈ Q i ) with the assigning probability ρ . There is an equal j

chance of choosing cm from any of the q i RBNs in the set Q i . On the other hand, the node cni is non-assigned with a complementary probability 1 − ρ and its state will not change. When cni is as-

signed a driver cm , the operator Sˆ compares cni with cm . If cni = cm , j

j

the state of cni is not modified. If cni = j cm

j cm ,

j

influence probability. Consequently, the stochastic coupling scheme for the node cni can be expressed as [39]

Sˆ ni =

⎧ ci ⎪ ⎪ ⎪ ni ⎨

if cni is non-assigned,

⎪ cni ⎪ ⎪ ⎩ j

with probability 1 − g , if cni = cm ,

j

j

if cni is assigned a driver cm and cni = cm ,

cn

j

(2)

j

with probability g , if cni = cm .

cm

Note that the coupling relations between RBNs are not bidirectional, i.e., j ∈ Q i does not indicate i ∈ Q j simultaneously. At each time step, the neighboring RBNs of the ith RBN is re-determined from all RBNs with equal probability, which is referred to as unfrozen coupling configuration [45]. For simplicity, the neighboring set Q i may contain the ith RBN itself, and the possibility of choosing the same neighbor to a RBN several times is not ruled out. j Similarly, the information exchange process of cni and cm is also j

unidirectional (that cni is assigned a driver cm does not imply that

j cm

cni ).

is assigned a driver The inter-network couplings between two nodes are also shuffled at each time step. The topological structure of the studied model is a random network which consists of random networks. The inter-network coupling considered here is the microscopic interaction which relies on the probabilistic assimilating process between two nodes (i.e., the state of cni is set to coincide with its randomly assigned driver j

cn with a probability). Since the coupling scheme is strictly different from deterministic coupling, an interesting question therefore arises: can the coherent behaviors, especially the mutual synchronization, be achieved based on the aforementioned coupling scheme? To answer this question,  N we define the density of the ith RBN at time t as p i (t ) = n=1 cni (t )/ N (0  p i (t )  1). A suitable approach to study synchronization is to compute the standard mean square deviation of p i (t ) [46]



M

1  2 σ (t ) =  p i (t ) − p¯ (t ) , M

i =1

(3)

M

1 i known as the synchronization error. Here p¯ (t ) = M i =1 p (t ) is the average density at time t. All RBNs are synchronized whenever limt →∞ σ (t ) ≈ 0. Parameters M = 10, N = 106 , and the assigning probability ρ = 0.9 are chosen. The initial state of each node cni (0) within the ith RBN is assigned 1 with a biased probability p i (0); otherwise cni (0) = 0. The value of p i (0) is determined by an independent random number distributed homogeneously in the interval (0, 1), so that the p i (t ) of each RBN evolves from different initial values. To simplify the analysis, all ki s are set as an identical constant, ki = k = 12, and the number of neighboring RBNs q i = q = 3. Fig. 1 illustrates σ (t ) as a function of t for various values of influence probability g. In the condition of g = 0.00, all RBNs evolve independently and the errors fluctuate around a constant, as shown in Fig. 1(a). With the increase of g, the mean values of errors decrease correspondingly (Figs. 1(b) and 1(c)), and finally approach to zero when g = 0.75 (Fig. 1(d)). Similar phenomena can be observed for other values of k which ensure the chaotic motions of RBNs. These results indicate that the proposed coupling scheme is efficient enough to overcome the chaos revealed by RBNs to induce the mutual synchronization. It is important to explicitly determine under what conditions the synchronization can occur. To this aim, the mean error, σ = t0 +T 1 t =t 0 σ (t ), is evaluated as a function of coupling strength ε ≡ T

ρ g.1 Here t 0 is the length of the discarded transients and T = 105

the state of cni is set

to coincide with with a probability g, or keeps invariant with a complementary probability 1 − g. Here g is referred to as the

1 The definition of coupling strength section.

ε used here will be clarified in the next

Y.-C. Hung / Physics Letters A 375 (2011) 2809–2814

2811

At the second sub-step, the stochastic coupling operator Sˆ is applied to all RBNs. To facilitate the formulation, the notations f i and f j are used in place of the density function f i ( p i (t ), ki ) and ˆ the f j ( p j (t ), k j ). Within an arbitrary RBN, after being acted by S, number of nodes still staying in 1 can be expressed by

N 1i →1 = N (1 − ρ ) f i +

Nρ qi



fi f j + fi 1− f

j



 (1 − g ) . (6)

j∈ Q i

The first term in the right-hand side of Eq. (6) corresponds to the number of unconnected nodes. The second term indicates in the j j respective conditions of (cni , cm ) = (1, 1) and (cni , cm ) = (1, 0) the

ˆ Based on a similar alnumber of cni still staying invariant after S. gorithm, we may get the number of nodes changing their states from 1 to 0, Fig. 1. The evolution of synchronization error σ (t ) in the condition of (a) g = 0.00, (b) g = 0.35, (c) g = 0.65, and (d) g = 0.75. When the influence probability g is sufficient large, the value of σ (t ) approaches to zero indicating the achievement of mutual synchronization among all RBNs. System parameters M = 10, N = 106 , ρ = 0.9, k = 12, and q = 3 are used.

N 1i →0 =

Nρ qi

j∈ Q



fi 1− f

j

  g .

(7)

i

Similarly,





N 0i →0 = N (1 − ρ ) 1 − f i +



+ 1− f

i



N ρ  qi



1− fi 1− f

j



j∈ Q i

 f (1 − g ) , j

(8)

and

N 0i →1 =

  N ρ  1 − f i f jg . i q i

(9)

j∈ Q

The above equations also satisfy the conservation condition, N 1i →1 + N 1i →0 + N 0i →0 + N 0i →1 = N. The MF model then can be constructed based on these equations. The probability of finding a node in state 1 at time t + 1 is given by Fig. 2. (Color online.) The average synchronization errors σ versus coupling strength ε of the coupled RBNs for N = 106 (red squares), N = 104 (blue triangles), and N = 5 × 103 (green circles). Other system parameters are identical to those used in Fig. 1. The dashed line labels the results predicted by the mean-field model. The inset plots the deviation of σ between the original coupled RBNs and the mean-field model for different values of N.

is the averaged time steps. The results for N = 106 have been shown in Fig. 2 marked by (red) squares, which reveal a typical scenario of transition from unsynchronized state to mutually synchronized state. For the present parameters set, ε > 0.60 ensures the achievement of synchronization. 3. Mean-field model To understand the underlying mechanism of synchronization, a mean-field model (MF model) is constructed to capture the dynamical properties of the coupled RBNs. Let N 1i (t ) and N 0i (t ) be the number of nodes in the state 1 and 0 respectively within the ith RBN at time t. The two variables satisfy the conservation condition, N 1i (t ) + N 0i (t ) = N. Now the MF model can be formulated progressively. At the first sub-step, all RBNs are acted by the free ˆ and then N i and N i turn out to be [39,40] evolving operator R, 1 0





N 1i = N f i p i (t ), ki ,







N 0i = N 1 − f i p i (t ), ki ,

(4)

where f i ( p i (t ), ki ) is the density map determined by the generalized ECA rule 126 [41],

 i

 i

ki +1 i f p i (t ), k = 1 − p i (t )k +1 − 1 − p i (t ) .

(5)

p i (t + 1) =

N 1i →1 + N 0i →1 N

(10)

.

Inserting Eqs. (6) and (9) into Eq. (10), we can get

p i (t + 1) = (1 − ρ g ) f i +

ρg

Defining coupling strength formulated as



qi

j∈ Q

f j.

(11)

i

ε ≡ ρ g, the MF model finally can be 

p i (t + 1) = (1 − ε ) f i p i (t ), ki +

ε qi

f

j





p j (t ), k j ,

(12)

j∈ Q i

where i = 1, 2, . . . , M, Q i is the set of neighboring RBNs of the ith RBN, and q i is the number of RBNs within Q i . Regardless of the stochastic nature of the original system, the MF model is represented in a familiar form of deterministically coupled maps: each map in the random graph diffuses a certain fraction of density to its neighboring maps and receives a comparable fraction from others [27,30]. Using the parameters set identical to aforementioned results, i.e., M = 10, ki = k = 12, and qi = q = 3, we plot mean errors σ versus ε in Fig. 2 (the black dashed line). As one can observe, the values of σ predicted by MF model coincide well with the results from the original model, excluding those near the synchronized threshold εc .2 Similar phenomena can be observed for other parameters sets q and k which ensure

2 In this Letter, synchronization threshold εc is determined numerically by meas of the MF model. When ε exceeds εc , synchronization error σ (t ) will converge to 0 asymptotically. Detailed analysis is included in the next section.

2812

Y.-C. Hung / Physics Letters A 375 (2011) 2809–2814

the deviations from the deterministic descriptions that are caused by the finite system size. Importantly, the stochastic model used here is based on a very simple assumption: the introduced noise is a biased white noise. However, further analysis demonstrates that the probability distribution function of the intrinsic noise that is measured from the original model is neither a uniform function with bias nor a Gaussian function, but has a complex form. Some preliminary calculations indicate that the measured noise is weakly autocorrelated and its variance √ is inversely proportional to the square root of the system size N. Detailed studies that involve the refinement of the proposed stochastic MF model are still in progress, and these will be published in the future. Fig. 3. (Color online.) The average synchronization errors σ versus coupling strength ε of the stochastic mean-field model for noise intensity η = 10−4 (red solid line), η = 10−3 (blue dash-dot), and η = 10−2 (green dotted line). Other system parameters are identical to those used in Fig. 1. The dashed line labels the results of the deterministic mean-field model (η = 0). The inset plots the deviation of σ between the stochastic mean-field model and the mean-field model for different values of noise intensity η .

the chaotic motion of individual map. The inset in Fig. 2 clearly draws the deviation of σ between the MF model and the original RBNs, σδ ≡ |σMF − σRBNs |, near εc . The deviation would be further magnified if N (the number of nodes within a RBN) is further decreased. Examples are labelled by blue triangles (N = 104 ) and green circles (N = 5 × 103 ) in Fig. 2. In fact, the predictions of the MF model are in good agreement on the original system only in the limit N → ∞. When the number of interacting elements N is finite or presented in a small value, the MF model fails to well capture the dynamical properties of the original system, as shown in Fig. 2. From the perspective of statistical physics, the inherent stochasticity and the discreteness of the dynamics in finite systems results in the existence of intrinsic noise, which causes deviations from the deterministic descriptions [47]. Similar phenomena have been observed in biological systems, and especially in the processes of gene expression in living cells [48]. Accordingly, an additional noise term is introduced into the deterministic density map to mimic the effect of finite size,





ki +1

f i p i (t ), ki = 1 − p i (t )k +1 − 1 − p i (t ) i

+ ηξ i (t ),

(13)

where η is the noise intensity, ξ (t ) ∈ [−1, 1] is a uniform white noise which satisfies ξ i (t ), ξ i (t  ) = δt ,t  and ξ i (t ), ξ j (t  ) = δi , j . Note that ξ i (t ) is chosen with the constraint 0  f i ( p i (t ), ki )  1. Thus if a chosen ξ i (t ) exceeds the reasonable bounds, it would be discarded and re-chosen. Such a noise term is referred to as the biased noise [49]. Applying the stochastic map to Eq. (12), we can get a stochastic MF model and then calculate σ versus ε under the influence of different noise intensities η , as shown in Fig. 3. In this diagram, the (black) dashed line, (red) solid line, (blue) dashdot, and (green) dotted line represent the results of noise intensity η = 0, η = 10−4 , η = 10−3 , and η = 10−2 respectively. The inset shows σδ , the deviation of σ between the stochastic MF model and the deterministic MF model, versus ε near the synchronization threshold εc . Compared with Fig. 2, the stochastic MF model truly reproduces the dynamical properties observed in the coupled RBNs with finite number of nodes. Our results indicate that the stochastic MF model with η = 10−4 approaches the original system with N = 106 while η = 10−3 is comparable with N = 104 . Since systems that consist of a finite number of interacting elements are ubiquitous in nature [34], understanding the effect of finite system size is crucial in dealing with real systems. As displayed in Fig. 2, a finite value of N may increase the threshold of transition to the synchronous state or desynchronize the system. Although the deterministic MF model is no longer completely accurate, the refined model with a noise term can accurately describe i

4. Synchronization analysis In this section, the synchronous behaviors revealed by the system will be further analyzed based on the MF model. Setting ki = k = 12 and qi = q = 3, we depict mean errors σ as a function of ε for various values of the number of RBNs M in Fig. 4. The (black) circles, (red) triangles, and (blue) squares respectively denote the numerical results for M = 10, M = 100, and M = 1000. We clearly see that the mutual synchronization can still be achieved even for the network with a larger number of RBNs. The inset in Fig. 4 draws log |σ | against ε for M = 1000 in detail. It displays that the value of σ undergoes a sharp cut-off near εc , which indicates the transition to synchronization in the studied system is of first order [46]. The aforementioned first-order transition scenario persists in the networks with different values of M (the largest number we have investigated is M = 214 ) and other values of k which ensure chaotic motions of individual RBN. Another interesting phenomenon exhibited in Fig. 4 is the dependence of critical coupling strength εc on the number of RBNs M. In order to characterize this relationship more clearly, εc is calculated as a function of M ranging from 23 to 214 , as illustrated in Fig. 5. Empty circles in the figure are averages over 10 to 100 independent initial conditions and coupling configurations depending on the required time for simulations. The step distance in determining εc for each realization is ε = 10−4 . The figure shows that a higher synchronization threshold εc is required when the value of M increases. The solid line is the fitting curve of the numerically obtained data. The fitting function reads,

εc = a −

b

(1 + cM )1/d

,

(14)

with the coefficients a = 0.5615, b = 0.3525, c = 6.6362, and d = 2.8820. In the limit M → ∞, the second term in the right-hand side of Eq. (14) approaches to zero and thus the corresponding critical value approaches ε ∗  a = 0.5615. The inset plots the deviation between εc and ε ∗ against M. The fitting line has slop −0.3451. The number of neighboring RBNs q has a great impact on the collective behaviors of the system. In the cases of q = 1 and q = 2, no mutually complete synchronization is emerged for all possible ε owing to the sparse wiring structure (treelike topology). For q  3, there exist critical values εc which separate the unsynchronized state (US) and the mutually synchronized state (MS). The synchronization transitions appearing at these critical points are also of first order. Fig. 6 illustrates the phase diagram in ε –q plane. The solid circles separating US and MS are synchronization thresholds obtained in the condition of unfrozen coupling configuration. Obviously, the value decreases in response to the ascending order of q, and asymptotically converges to the limitation ε g = 0.436 (horizontal dashed line), the threshold in the condition of global coupling configuration. Here ε g is determined analytically by means of the diffusive synchronization stability matrix (DSSM) [50]. To

Y.-C. Hung / Physics Letters A 375 (2011) 2809–2814

Fig. 4. (Color online.) The average synchronization errors σ versus coupling strength ε of the mean-field model for the number of RBNs M = 10 (black circles), M = 100 (red triangles), and M = 1000 (blue squares). Other system parameters are identical to those used in Fig. 1. The inset depicts the transition from unsynchronized state to synchronized state near the synchronization threshold (in the case of M = 1000). The sharp variation of σ ensures the appearance of first order transition.

Fig. 5. Dependence of the averaged synchronization threshold εc on the network size M. Each empty circle is averaged over 10–100 different realizations with independent initial conditions and coupling structures. The vertical axis in the inset indicates the deviation between εc and ε ∗ . The fitting function reads |εc − ε ∗ | ∝ M Γ , where Γ = −0.3451. System parameters k = 12, and q = 3 are used.

2813

Fig. 7. Dependence of the departure between εc and ε g on the value of q, the number of neighboring RBNs. Empty circles, squares, and triangles draw the results of unfrozen coupling configuration in the conditions of k = 12, k = 9, and k = 15, respectively. Solid circles, squares, and triangles represent the results of frozen coupling configuration under k = 12, k = 9, and k = 15. The fitting slops are γUF ≈ −1.5 (unfrozen configuration) and γ F ≈ −0.82 (frozen configuration).

Indeed, both εc and ε g depend on the maximal Lyapunov exponent λ of the density map, which is a function of k [39,44]. It is therefore interesting (and important) to explore how the parameter k influences the aforementioned scaling behavior. Here we focus on the values of k ensuring the chaotic motion of maps (λ > 0). Numerical results of two examples, k = 9 (empty squares) and k = 15 (empty triangles), are drawn in Fig. 7. As one can see, the scaling coefficient γUF ≈ −1.5 persists even for the adoption of other values of k. We have tested more than 8 different values of k and gotten an identical value of γUF . Such a result implies these systems with various k share an universal converging route toward the global coupling limitation. We also consider the frozen coupling configuration, i.e., the neighboring RBNs of the ith RBN are randomly chosen in the beginning and fixed during evolution. The connecting graph has been ensured to be fully connected. In other words, it cannot be separated into two or more than two disconnected parts. In this case, only q  4 ensures the achievement of MS and a higher εc is required comparing to unfrozen coupling configuration, as shown in Fig. 6 (solid circles k = 12, solid squares k = 9, solid triangles k = 15). With the increase of q, the value of εc also approaches to the limitation of ε g but with a smaller slope (empty triangles in Fig. 7). The fitting functions read,

εcF − ε g ∝ qγ F ,

(16)

where the superscript F indicates the frozen coupling configuration and the scaling coefficient γ F ≈ −0.82. Similarly, the systems with different k share an identical value of γ F . 5. Conclusion

Fig. 6. Phase diagram in ε –q plane for the deterministic MF model defined in Eq. (12). The symbols US and MS correspond to the unsynchronized state and the mutual synchronization state. The solid circles and the empty triangles separating US and MS are obtained in the conditions of unfrozen coupling configuration and frozen coupling configuration, respectively. The horizontal dashed line labels the synchronization threshold of the global coupling architecture.

clarify the converging process, the departure of εc from ε g is calculated as a function of q and the results are depicted in Fig. 7 (empty circles). Numerically, we have

εcUF − ε g ∝ qγUF ,

(15)

where the superscript UF indicates the unfrozen coupling configuration and γUF ≈ −1.5.

Mutual synchronization in complex networks has recently inspired a considerable number of investigations. However, to our knowledge, less attention has been paid to microscopic couplings that lead to globally mutual synchronization. In this Letter, transitions to synchronization in a random network of random Boolean networks (RBNs) is studied by introducing a microscopically stochastic coupling scheme. The mutually synchronous stat is numerically demonstrated able to be realized when the coupling strength exceeds a threshold value. The proposed mean-field model, a random network of density maps, further provides an analytical model for predicting the dynamics of the original system. The results thus obtained have been carefully checked for various parameters sets. Recently, the effect of noise on synchronization phenomena has been studied based on the macroscopically coupled Langevin

2814

Y.-C. Hung / Physics Letters A 375 (2011) 2809–2814

equations [4]. Noise may increase the threshold of transition to a synchronous state, and the synchronization threshold strongly depends on the nature of the noise. Here, no specific type of fluctuation is introduced; rather, in the system presented herein, the noise that results from the finite number of interacting elements is inherent and thus more realistic. The system herein therefore can be potentially applied to explore the destructive effect of intrinsic noise on synchronization. From the perspective of nonlinear dynamics, some systems may be able to utilize noise to improve the efficiency of functional performance [51,52]. Our further works include an exploration of constructive noise in the proposed microscopically interacting systems. Acknowledgements This work was supported by the National Science Council of Republic of China (Taiwan) under the Grant No. NSC 99-2112-M194-001. Ted Knoy is appreciated for his editorial assistance. References [1] K. Otsuka, R. Kawai, S.-L. Hwong, J.-Y. Ko, J.-L. Chern, Phys. Rev. Lett. 84 (2000) 3049. [2] K. Otsuka, T. Ohtomo, A. Yoshioka, J.-Y. Ko, Chaos 12 (2002) 678. [3] J. Hassel, L. Grönberg, P. Helistö, H. Seppä, Appl. Phys. Lett. 89 (2006) 072503. [4] B.C. Bag, K.G. Petrosyan, C.-K. Hu, Phys. Rev. E 76 (2007) 056210. [5] S. Rajesh, S. Sinha, S. Sinha, Phys. Rev. E 75 (2007) 011906. [6] C. Schäfer, M.G. Rosenblum, J. Kurths, H.-H. Abel, Nature (London) 392 (1998) 239. [7] M.-C. Wu, C.-K. Hu, Phys. Rev. E 73 (2006) 051917. [8] R.E. Amritkar, G. Rangarajan, Phys. Rev. Lett. 96 (2006) 258102. [9] I.Z. Kiss, Q. Lv, J.L. Hudson, Phys. Rev. E 71 (2005) 035201(R). [10] A.S. Pikovsky, M.G. Rosenblum, J. Kurths, Synchronization: A Univeral Concept in Vonlinear Sciences, Cambridge University Press, Cambridge, 2001. [11] L.M. Pecora, T.L. Carroll, Phys. Rev. Lett. 64 (1990) 821. [12] K. Kaneko, Phys. Rev. Lett. 65 (1990) 1391. [13] M.-C. Ho, Y.-C. Hung, Phys. Lett. A 301 (2002) 424. [14] A. Pikovsky, M. Zaks, M. Rosenblum, G. Osipov, J. Kurths, Chaos 7 (1997) 680; A. Pikovsky, M. Zaks, M. Rosenblum, G. Osipov, J. Kurths, Phys. Rev. Lett. 79 (1997) 47. [15] M.-C. Ho, Y.-C. Hung, I.-M. Jiang, Phys. Lett. A 324 (2004) 450. [16] U. Parlitz, L. Junge, L. Kocarev, Phys. Rev. Lett. 79 (1997) 3158.

[17] M.G. Rosenblum, A.S. Pikovsky, J. Kurths, Phys. Rev. Lett. 78 (1997) 4193. [18] H.U. Voss, Phys. Rev. E 61 (2000) 5115; H.U. Voss, Phys. Rev. E 64 (2001) 039904(E); H.U. Voss, Phys. Rev. Lett. 87 (2001) 014102. [19] M.-C. Ho, Y.-C. Hung, Z.-Y. Liua, I.-M. Jiang, Phys. Lett. A 348 (2006) 251. [20] S.H. Strogatz, Nature (London) 410 (2001) 268. [21] T. Nishikawa, A.E. Motter, Y.-C. Lai, F.C. Hoppensteadt, Phys. Rev. Lett. 91 (2003) 014101. [22] A.E. Motter, M.A. Matias, J. Kurths, E. Ott, Physica D 224 (2006) vii. [23] P.M. Gade, C.-K. Hu, Phys. Rev. E 62 (2000) 6409. [24] M. Barahona, L.M. Pecora, Phys. Rev. Lett. 89 (2002) 054101. [25] J. Garcia-Ojalvo, M.B. Elowitz, S.H. Strogatz, Proc. Natl. Acad. Sci. 101 (2004) 10955. [26] W.-X. Qin, G.-R. Chen, Physica D 197 (2004) 375. [27] S. Jalan, R.E. Amritkar, Phys. Rev. Lett. 90 (2003) 014101. [28] S. Jalan, R.E. Amritkar, C.-K. Hu, Phys. Rev. E 72 (2005) 016211; S. Jalan, R.E. Amritkar, C.-K. Hu, Phys. Rev. E 72 (2005) 016212. [29] R.E. Amritkar, C.-K. Hu, Chaos 16 (2006) 015117. [30] Y.-C. Hung, Y.-T. Huang, M.-C. Ho, C.-K. Hu, Phys. Rev. E 77 (2008) 016202. [31] S. Camazine, J.L. Deneubourg, N.R. Franks, J. Sneyd, G. Theraulaz, E. Bonabeau, Self-Organization in Biological Systems, Princeton University Press, Princeton, 2001. [32] R. Ferri, C.J. Stam, B. Lanuzza, F.I. Cosentino, M. Elia, S.A. Musumeci, G. Pennisi, Clin. Neurophysiol. 115 (2004) 1202. [33] B. Blasius, A. Huppert, L. Stone, Nature 399 (1999) 354. [34] M. Kaern, T.C. Elston, W.J. Blake, J.J. Collins, Nature Rev. Genet. 6 (2005) 451. [35] L.M. Bishop, H. Qian, Biophys. J. 98 (2010) 1. [36] S. Wolfram, A New Kind of Science, Wolfram Media, Champaign, 2002. [37] L.G. Morelli, D.H. Zanette, Phys. Rev. E 58 (1998) R8. [38] L.G. Morelli, D.H. Zanette, Phys. Rev. E 63 (2001) 036204. [39] Y.-C. Hung, M.-C. Ho, J.-S. Lih, I.-M. Jiang, Phys. Lett. A 356 (2006) 35. [40] M.-C. Ho, Y.-C. Hung, I.-M. Jiang, Phys. Lett. A 344 (2005) 36. [41] M. Andrecut, M.K. Ali, Int. J. Mod. Phys. B 15 (2001) 17. [42] M.T. Matache, J. Heidel, Phys. Rev. E 69 (2004) 056214; M.T. Matache, J. Heidel, Phys. Rev. E 71 (2005) 026232. [43] G.L. Becka, M.T. Matache, Physica A 387 (2008) 4947. [44] J.-Y. Ko, Y.-C. Hung, M.-C. Ho, I.-M. Jiang, Chaos Solitons Fractals 36 (2008) 934. [45] H. Chaté, P. Manneville, Chaos 2 (3) (1992) 307. [46] P.G. Lind, J.A.C. Gallas, H.J. Herrmann, Phys. Rev. E 70 (2004) 056207. [47] T. Galla, Phys. Rev. Lett. 103 (2009) 198702. [48] P.S. Swain, M.B. Elowitz, E.D. Siggia, Proc. Natl. Acad. Sci. 99 (2002) 12795. [49] A. Maritan, J.R. Banavar, Phys. Rev. Lett. 72 (1994) 1451. ´ [50] A. Stefanski, J. Wojewoda, T. Kapitaniak, S. Yanchuk, Phys. Rev. E 70 (2004) 026217. [51] J. Wang, J. Zhang, Z. Yuan, T. Zhou, BMC Systems Biology 1 (2007) 50. [52] Y.-C. Hung, C.-K. Hu, Comput. Phys. Commun. 182 (2011) 249.