Synchronization in complex clustered networks

Abstract

Synchronization in complex networks has been an active area of research in recent years. While much effort has been devoted to networks with the small-world and scale-free topology, structurally they are often assumed to have a single, densely connected component. Recently it has also become apparent that many networks in social, biological, and technological systems are clustered, as characterized by a number (or a hierarchy) of sparsely linked clusters, each with dense and complex internal connections. Synchronization is fundamental to the dynamics and functions of complex clustered networks, but this problem has just begun to be addressed. This paper reviews some progress in this direction by focusing on the interplay between the clustered topology and network synchronizability. In particular, there are two parameters characterizing a clustered network: the intra-cluster and the inter-cluster link density. Our goal is to clarify the roles of these parameters in shaping network synchronizability. By using theoretical analysis and direct numerical simulations of oscillator networks, it is demonstrated that clustered networks with random inter-cluster links are more synchronizable, and synchronization can be optimized when inter-cluster and intra-cluster links match. The latter result has one counterintuitive implication: more links, if placed improperly, can actually lead to destruction of synchronization, even though such links tend to decrease the average network distance. It is hoped that this review will help attract attention to the fundamental problem of clustered structures/synchronization in network science.

Graphical abstract

Keywords

clustered networks / network analysis / synchronization / oscillators

1 Introduction

Recent years have witnessed a growing interest in the synchronizability of complex networks [1−21]. Earlier works [1−10] suggest that small-world [22], random [23] and scale-free [24] networks, due to their short network distances, are generally more synchronizable than regular networks. It has been found, however, that heterogeneous degree distributions typically seen in scale-free networks can inhibit their synchronizability [11], but add suitable weights to the network elements, i.e., the assignment of large weights to nodes with large degrees (the number of links) can enhance their chances to synchronize with each other [12−21]. Modifying local connecting structures, if done properly, could also change the synchronizability significantly [25−32]. Synchronizability of complex clustered networks has not been investigated until recently [33−36]. The purpose of this paper is to review recent progress in this new area of research in network science.

In general, a clustered network consists of a number of groups, in which nodes within each group are densely connected, but the linkages among the groups are sparse [37−45]. Clustered networks were first described by Zachary in 1977 as a model of social networks with group structures [46]. In particular, he examined the organization of martial-art clubs in a city and found that a number of clubs were actually originated from one root club, where the owners of those clubs had been former students in the root club. As a result, members within each club are close to each other, but interactions with members from a different club are much less likely. This type of organization can also be found commonly in the business world where restructuring and recombination are routine practices. The clustered structure can explain familiar social experiences such as quick identification of acquaintances. For example, once two people are introduced, they describe themselves in terms of their social characteristics (e.g., professions, places of work, and leisure activities, etc.). Next, each of them cites friends with social characteristics “close” to those of the other person. This is actually a second step in the process of introduction, but can be effectively seen as an attempt to find chains of acquaintances linking them. The success of this attempt depends on the clustered structure of the social network. Clustered networks have recently been systematically studied and analyzed in social science [37−39].

Besides in social science, clustered networks can arise in biological situations [40−42], as a key feature in a biological system is the tendency to form a clustered structure. For example, proteins with a common function are usually physically associated via stable protein-protein interactions to form larger macromolecular assemblies. These protein complexes (or clusters) are often linked by extended networks of weaker, transient protein-protein interactions to form interaction networks that integrate pathways mediating major cellular processes [40]. As a result, the protein-protein interaction network can naturally be viewed as an assembly of interconnected functional modules, or a clustered network. Furthermore, macroscopic tissues, a network of intercellular communication, typically exhibit clustered organizational structures in which organs largely consist of repeating, densely-connected, functional subunits such as the glomeruli in kidneys, liver lobules, and colon crypts. Interestingly, these organizational structures are greatly diminished in a cancer but still often persist as non-functional caricatures of the tissue of origin. The organizational structure of biological networks in organs represents an optimal strategy in the sense that synchronization must be maintained over a wide range of environments, e.g., an organ needs to be able to adapt itself to the environment and continue to function in the face of perturbations such as injury or infection. In addition, the strategy for optimal system dynamics within the cluster will probably be different from that for connecting the clusters into an organ. The clustered structure has also been identified in technological networks such as electronic circuits and computer networks [43−45].

A complex clustered network is typically small-world in that its average network distance is short. Moreover, its degree distribution can be made quite homogeneous. For the synchronization problem, an interesting question regarding the clustered structure is that for a given number of inter-cluster links, how would their distributions affect the synchronizability? By Use of linear stability analysis [47] and its generalization [48], the dependence of synchronizability on the number of clusters in a network has been studied [34], and it has been found that the network can become more synchronizable as the number of clusters is increased if the inter-cluster links are random. If those links of the clusters are mostly local or diametrical in a topological ring structure, the synchronizability would deteriorate continuously as more clusters appear in the network. Therefore, how links distribute among clusters has a significant influence on network synchronizability. A relevant question is that for a given distribution of the inter-cluster links, say random distribution, how does the number of links influence synchronizability? Given a complex network with a fixed (large) number of nodes, intuitively, its synchronizability can be improved by increasing the number of links, as a denser linkage makes the network more tightly coupled or, “smaller,” thereby facilitating synchronization. Our recent work [35, 36] on this problem has revealed a phenomenon that apparently contradicts this intuition. Namely, more links, which make the network smaller, do not necessarily lead to stronger synchronizability. There can be situations where more links can even suppress synchronization if placed improperly. In particular, it is found that the synchronizability of a clustered network is determined by the interplay between the inter-cluster and intra-cluster links in the network. Strong synchronizability requires that the numbers of the two types of links be approximately matched. In this case, increasing the number of links can indeed enhance synchronizability. However, if the matching is deteriorated, synchronization can be severely suppressed or even totally destroyed.

Our finding can have a potential impact on real network dynamics. In biology, synchronization is fundamental, on which many biological functions rely. Our result implies that, in order to achieve robust synchronization for a clustered biological network, the characteristics of links are more important than the number of them. Simply counting the number of links may not be enough to determine its synchronizability. Instead, links should be distinguished and classified to predict synchronization-based functions of the network. In technological applications, it is supposed that a large-scale, parallel computational task is to be accomplished by a computer network, for which synchronous timing is of paramount importance. Our result can provide clues as to how to design the network to achieve the best possible synchronization and consequently optimal computational efficiency.

In Section 2, a general linear-stability analysis is described for solving the synchronization problem in both continuous-time oscillator networks and discrete-time coupled-map networks. In Section 3, a theory is developed and numerical results are presented to demonstrate the effects of the distribution of inter-cluster links. In Section 4, the emphasis is placed on clustered networks with random inter-cluster links, and how the number of links affects the synchronizability of the oscillator network is examined. Two types of coupling schemes are studied in detail, theoretically and numerically. Extensive discussions of the main results and their broader implications are offered in Section 5.

2 Linear-stability analysis for synchronization

The approach taken to establish the result is to introduce nonlinear dynamics to each node in the network and then perform stability and eigenvalue analyses [48, 49].

2.1 Continuous-time oscillator networks

The goal is to establish synchronization conditions of clustered networks in a proper network-parameter space. Each oscillator, when isolated, is described by

\begin{array}{}\text{(1)}& \frac{\text{d}\mathit{x}}{\text{d}t}=\mathit{F}(\mathit{x})\end{array}

where x is a d-dimensional vector and F(x) is the velocity field. Without loss of generality a prototype oscillator model is chosen−the Rössler oscillator, for which x = [x, y, z]^T ([∗]^T denotes transpose), and

\begin{array}{}\text{(2)}& \mathit{F}(\mathit{x})=[-(y+z)\mathrm{,}x+ay\mathrm{,}b+z(x-c{)]}^{\text{T}}\end{array}

The parameters of the Rössler oscillator are chosen such that they exhibit chaotic oscillations. The network dynamics is described by

\begin{array}{}\text{(3)}& \frac{\text{d}{\mathit{x}}_i}{\text{d}t}=\mathit{F}({\mathit{x}}_i)-ϵ{\displaystyle \sum _{j=1}^N{G}_{ij}\mathit{H}({\mathit{x}}_j)}\end{array}

where H(x) = [x, 0, 0]^T is a linear coupling function, ϵ is a global coupling parameter, and G is the coupling matrix determined by the network topology. For theoretical convenience, G is assumed to satisfy the condition \begin{array}{l}{\displaystyle \sum _{j=1}^N}{G}_{ij}=0\hfill \end{array} for any i, where N is the network size; therefore, the system permits an exact synchronized solution: x¹ = x² = ··· = x^N = s, where ds/dt = F(s).

For the system described by Eq. (3), the variational equations governing the time evolution of the set of infinitesimal vectors \delta {\mathit{x}}_i(t)\equiv {\mathit{x}}_i(t)-\mathit{s}(t) are

\begin{array}{}\text{(4)}& \frac{\text{d}\text{δ}{\mathit{x}}_i}{\text{d}t}=\mathit{D}\mathit{F}(\mathit{s})\cdot \text{δ}{\mathit{x}}_i-ϵ{\displaystyle \sum _{j=1}^N{G}_{ij}\mathit{D}\mathit{H}(\mathit{s})\cdot \text{δ}{\mathit{x}}_j}\end{array}

where DF(s) and DH(s) are the Jacobian matrices of the corresponding vector functions evaluated at s(t). Diagonalizing the coupling matrix G yields a set of eigenvalues {λ_i} (i = 1, ···, N) and the corresponding normalized eigenvectors are denoted by e₁, e₂, ···, e_N. Generally, the eigenvalues are real and non-negatives [49] and thus can be sorted as 0 = λ₁ < λ₂ ≤ ··· ≤ λ_N. The smaller the ratio λ_N /λ₂, the stronger the synchronizability of the network (to be discussed below) [11−14, 17]. The transform \delta \mathit{y}=O^{-1}\cdot \delta \mathit{x}, where O is a matrix whose columns are the set of eigenvectors, leads to the following block-diagonally decoupled form of Eq. (4):

\begin{array}{}\text{(4)}& \frac{\text{d}{\mathit{y}}_i}{\text{d}t}=[\mathit{D}\mathit{F}(\mathit{s})-ϵ{\lambda }_i\mathit{D}\mathit{H}(\mathit{s})]\cdot \text{δ}{\mathit{y}}_i\end{array}

Letting K=ϵ{\lambda }_i(i=\mathrm{2,}\cdots \mathrm{,}N) be the normalized coupling parameter, it can be written as follows:

\begin{array}{}\text{(5)}& \frac{\text{d}\text{δ}\mathit{y}}{\text{d}t}=[\mathit{D}\mathit{F}(\mathit{s})-K\mathit{D}\mathit{H}(\mathit{s})]\cdot \text{δ}\mathit{y}\end{array}

The largest Lyapunov exponent from Eq. (5) is the master-stability function Ψ(K) [48]. If Ψ(K) is negative, a small disturbance from the synchronization state will be diminished exponentially, and the system is stable and can be synchronized; if Ψ(K) is positive, a small disturbance will be magnified and the system cannot be synchronized.

The function Ψ(K) is negative in an interval [K₁, K₂], where K₁ and K₂ solely depend on F(x) that governs the node dynamics, and the output function H(x). Thus, for K₁ < K < K₂, all the eigenvectors (eigenmodes) are transversely stable and the network can be synchronized, which gives the condition of the boundary of the synchronization region:

\begin{array}{}\text{(6)}& {\lambda }_2⩾\frac{K_1}{ϵ}\end{array}

\begin{array}{}\text{(7)}& {\lambda }_N⩽\frac{K_2}{ϵ}\end{array}

For the Rössler oscillators used in the simulation in Section 4(4.2), a = b = 0.2, c = 9 are chosen and the master stability function is shown in Fig. 1, where K₁ ≈ 0.2 and K₂ ≈ 4.62.

2.2 Discrete-time coupled-map network

The following general class of coupled-map networks is considered:

\begin{array}{}\text{(8)}& {\mathit{x}}_{m+1}^i=\mathit{f}({\mathit{x}}_m^i)-ϵ{\displaystyle \sum _j{G}_{ij}}\mathit{H}[\mathit{f}({\mathit{x}}_m^j)]\end{array}

where x_m+1 = f(x_m) is a d-dimensional map, ϵ is a global coupling parameter, G is the coupling matrix, and H is a coupling function. If the rows of the coupling matrix G have zero sum, Eq. (8) permits an exact synchronized solution: {\mathit{x}}_m^1={\mathit{x}}_m^2=\cdots ={\mathit{x}}_m^N={\mathit{s}}_m, where s_m+1 = f(s_m).

For the system described by Eq. (8), the variational equation governing the time evolution of the set of infinitesimal vectors \text{δ}{\mathit{x}}^i\equiv {\mathit{x}}^i-\mathit{s} is

\begin{array}{}\text{(9)}& \begin{array}{ll}\text{δ}{\mathit{x}}_{m+1}^i=\hfill & \mathit{D}\mathit{f}({\mathit{s}}_m)\cdot \text{δ}{\mathit{x}}_m^i-ϵ{\displaystyle \sum _j}{G}_{ij}\mathit{D}\mathit{H}[\mathit{f}({\mathit{s}}_m)]\hfill \\ \hfill & \cdot \mathit{D}\mathit{f}({\mathit{s}}_m)\cdot \text{δ}{\mathit{x}}_m^j\hfill \end{array}\end{array}

where Df and DH are the Jacobian matrices of the corresponding vector functions evaluated at s_m and f(s_m), respectively. Diagonalizing the coupling matrix G yields a set of eigenvalues {λ_i} (i = 1, ···, N), where 0 = λ₁ < λ₂ ≤ ··· ≤ λ_N. The transform \delta \mathit{y}=O^{-1}\cdot \text{δ}\mathit{x}, where O is a matrix whose columns are the set of eigenvectors, leads to the following block-diagonally decoupled form of Eq. (9):

\begin{array}{}\text{(10)}& \text{δ}{\mathit{y}}_{m+1}^i=\{\mathit{I}-ϵ{\lambda }_i\mathit{D}\mathit{H}[\mathit{f}({\mathit{s}}_m)]\}\cdot \mathit{D}\mathit{f}({\mathit{s}}_m)\cdot \text{δ}{\mathit{y}}_m^i\end{array}

Compared with Eq. (5), and letting K=ϵ{\lambda }_i, the largest Lyapunov exponent of Eq. (10) can be called its master stability function, denoted again by Ψ(K). Generally, there exist K₁ and K₂ such that Ψ(K) is negative in the interval [K₁, K₂], for which the synchronization solution of the coupled system is linearly stable. For a special class of coupled-map networks, K₁ and K₂ can be obtained explicitly as follows.

The system is stable if for any i that 2 ≤ i ≤ N. The following holds:

For a linear coupling function H, DH is a constant matrix. If the system is one dimensional, DH is simply a constant, say, γ. Equation (11) becomes

Recognizing that the second term in the above equation is the Lyapunov exponent µ of a single map, it can be obtained

which is

To gain insight, f(x) is set to be the one-dimensional logistic map f(x) = 1 − ax² and H(f) = f is chosen. Choosing \begin{array}{l}{\displaystyle \sum _{j\ne i}}{G}_{ij}=-1\hfill \end{array}, γ = 1 is obtained and Eq. (14) becomes [49]:

Because of the ordering of the eigenvalues, the above inequality will hold for all the is if it holds for i = 2 and i = N. Therefore, condition (15) can be further simplified as

\begin{array}{}\text{(16)}& {\lambda }_2>\frac{1}{ϵ}(1-{\text{e}}^{-\mu })\end{array}

Compared with Eqs. (6) and (7), it can be seen that for the coupled logistic-map network, K₁ = 1−e^−µ and K₂ = 1 + e^−µ. The boundary of the synchronization region in the phase diagram can be determined by setting {\lambda }_2=(1-{\text{e}}^{-\mu })/ϵ and {\lambda }_N=(1+{\text{e}}^{-\mu })/ϵ. In our simulation in Sec. IVC, a = 1.9 is used, and the corresponding Lyapunov exponent is µ = 0.55, so {\lambda }_2=0.423/ϵ, and {\lambda }_N=1.577/ϵ. If the coupling function H is nonlinear, DH[f(s_m)] will depend on the value of f(s_m) and generally Eq. (11) cannot be simplified further.

2.3 Physical understanding of synchronization boundaries

Now it can be seen that the synchronization conditions for both continuous-time oscillators and discrete-time coupled-map networks have the same form of Eqs. (6) and (7). Physically, the existence of the K₁ and K₂ boundaries can be understood by the roles of the coupling terms: (1) they serve to establish coherence among oscillators, and (2) they are effective perturbations to the dynamics of individual oscillators. Whether synchronization can occur depends on the interplay between these two factors. In particular, for small coupling, synchronization may not occur because of the fact that although the perturbing effect of the coupling terms is small, the amount of coherence provided by them is also small. For very large coupling, although the coupling terms can provide strong coherence, the effective perturbations are also large. As a large perturbation requires longer time for the system to reach an equilibrium state (e.g., synchronization), the system will have no time to respond to the perturbations, which consequently makes it unable to synchronize. Thus, synchronization may not occur if the coupling is too strong. In general, there exists a finite interval of the coupling parameter for which synchronization can occur [48]. These considerations are demonstrated in Fig. 1, the master stability function of the coupled Rössler system versus the generalized coupling parameter K. We see that synchronization can occur only when the coupling parameter K falls in the interval (K₁, K₂). Indeed, this behavior appears to be typical for a large class of coupled chaotic oscillators [48, 50].

The synchronizability of a complex network of oscillators for any linear coupling scheme can be inferred from Fig. 1. A given network can be characterized by the set of eigenvalues (λ₂ through λ_N) of the corresponding coupling matrix. For the fixed value of ϵ, the spread of the eigenvalues determines the range of possible variations in K. This suggests a general quantity that determines the synchronizability of a complex network, regardless of detailed oscillator models: the spread of the eigenvalues of the coupling matrix. In Fig. 1, it can be seen that in order to achieve synchronization, the spread of the eigenvalues must not be too large to fit the generalized coupling parameter K in the interval (K₁, K₂). That is, only if

there could exist values of ϵ of positive Lebesgue measure such that the synchronized state is linearly stable. Otherwise, the system is unstable for any ϵ value. For various chaotic oscillators, the value of β ranges from 5 to 100 [17]. The left-hand side of the inequality (18), the ratio λ_N /λ₂, depends only on the topology of interactions among oscillators. Hence, the impact of a particular coupling topology on the network’s ability to synchronize is represented by a single quantity λ_N/λ₂: the larger the ratio is, the more difficult it is to synchronize the oscillators, and vice versa [7].

3 Synchronizability transition phenomena in clustered networks

To address the synchronizability of a clustered network, it is insightful to explore the relationship between the eigenratio λ_N/λ₂ and the number of clusters. To construct an analyzable model of a clustered network, it is imagined that there is a network of N nodes grouped into M clusters located on a ring, each being connected to their nearest-neighboring clusters. The Laplacian matrix of the network is similar to that of a regular, “ring” network of M nodes. The eigenvector corresponding to the first non-zero eigenvalue of the ring network is

\begin{array}{}\text{(18)}& \sqrt{2/N}{\left[\sin \frac{2\text{π}j}{N}\right]}_{j=1}^N\end{array}

which can be considered as an envelope function of the eigenvectors in the clustered network, as shown in Fig. 2. Since the variance of the components in a cluster is small compared with the difference between the means of two consecutive clusters, components within a cluster are further approximated to obtain the same value and, as a result, the eigenvector of the network with M clusters is given as a piecewise continuous step function:

\begin{array}{}\text{(19)}& {\mathit{e}}_2^{\text{T}}=(\stackrel{N/M}{\overbrace{h_M^1\mathrm{,}\cdots }}\mathrm{,}\cdots \mathrm{,}\stackrel{N/M}{\overbrace{h_M^i\mathrm{,}\cdots \mathrm{,}{h}_M^i}}\mathrm{,}\cdots \mathrm{,}\stackrel{N/M}{\overbrace{\cdots \mathrm{,}{h}_M^M}})\end{array}

where all N/M components of the ith cluster have the same value \begin{array}{l}{h}_M^i=\sqrt{2/N}\sin (2\text{π}i/M)\hfill \end{array}. The first non-zero eigenvalue of the clustered network is

\begin{array}{}\text{(20)}& {\lambda }_2={\mathit{e}}_2^{\text{T}}\mathit{G}{\mathit{e}}_2\end{array}

where G is the Laplacian matrix of the network. Inserting Eq. (19) into Eq. (20), and assuming there are only two nearest connections per cluster, the following equation can be obtained

\begin{array}{}\text{(21)}& {\lambda }_2=\frac{2}{N}{\displaystyle \sum _{i=1}^M\sin \frac{i}{{\theta }_0}\left(2\sin \frac{i}{{\theta }_0}-\sin \frac{i-1}{{\theta }_0}-\sin \frac{i+1}{{\theta }_0}\right)}\end{array}

where θ₀ = M/2π. For M\gg 1, we have

\begin{array}{}\text{(22)}& {\lambda }_2\approx \frac{4\text{π}}{NM}{\displaystyle {\int }_0^{2\text{π}}}\sin \theta (-{\nabla }^2\sin \theta )\text{d}\theta =\frac{4{\text{π}}^2}{NM}\end{array}

Although in general no approximation to the largest eigenvalue can be obtained in a similar way, approximations for it are known in some special cases. For example, when the cluster is a regular “ring” network with 2k nearest connections, the largest eigenvalue of the cluster is λ_N ≈ (2k + 1)(1 + 2/3π) [51]. Since λ_N depends on the number of nearest neighbors only, the same λ_N may be used as the largest eigenvalue of the network. Therefore, the ratio is obtained

\begin{array}{}\text{(23)}& \frac{{\lambda }_N}{{\lambda }_2}\approx \frac{(2k+1)(1+2/3\text{π})NM}{4{\text{π}}^2}\end{array}

which indicates that the eigenratio increases linearly with M, as shown in Fig. 3(a). That is, only with local connections between clusters, it is more difficult to achieve synchronization as the number of clusters increases.

Now a more general clustered network model that allows long range interactions can be considered. That is, two arbitrary clusters can make a connection with a probability p(l), where l is the distance between the clusters. It is reasonable to assume that p(l) has an exponential dependence on l,

\begin{array}{}\text{(24)}& p(l)=\frac{{\text{e}}^{-\alpha l}}{N_l}\end{array}

where α is a parameter and N_l is a proper normalization constant. Inserting this p(l) into Eq. (20) and again replacing the sum with an integral, it can be obtained

\begin{array}{}\text{(25)}& {\lambda }_2\approx \frac{4{\text{π}}^2{K}_0}{NM}\frac{{\displaystyle {\int }_1^{M/2}}{l}^2{\text{e}}^{-\alpha l}\text{d}l}{{\displaystyle {\int }_1^{M/2}}{\text{e}}^{-\alpha l}\text{d}l}\end{array}

where the factor K₀ comes from the fact that each cluster can have multiple (2K₀) connections. As an example, for the case of a “globally” connected clustered network, α → 0 and K₀ ∼ M/2 can be obtained. The integral can then be simplified as M ² and λ₂ ∼ M². Thus, the eigenratio decreases as M increases:

\begin{array}{}\text{(26)}& \frac{{\lambda }_N}{{\lambda }_2}\sim M^{-2}\end{array}

as shown in Fig. 3(b). This is opposite to the result in Fig. 3(a), where no long-range links among clusters are allowed. It is noted that for α → ∞and K₀ = 1, Eq. (25) becomes Eq. (22).

An interesting consequence of Eq. (25) is that a transition phenomenon may occur in the synchronizability of clustered networks. For \alpha \gtrsim 0, connections between clusters are decentralized and the integral in Eq. (25) is proportional to M ², and, as a result, the eigenratio decreases as M increases. However, for \alpha \gg 0, connections between clusters are centralized among nearest neighboring clusters and the integral becomes independent of M for sufficiently large M\gg {\alpha }^{-1} and the eigenratio eventually increases as M increases. A similar transition phenomenon is also expected when K₀ is regarded as a function of M (i.e., K₀ ∼ M^η). For example, it is supposed that α → 0 and the total number of connections between clusters is fixed. Then K₀ is inversely proportional to M; as a result, the eigenratio remains constant. Distinct synchronization behavior can arise depending on the value of η.

Our considerations so far have been limited to the case where the parameter α is positive. From Eq. (24), it can be seen that a relatively large, positive value of α stipulates smaller probabilities for long-range links among clusters as compared with those for short-range links. Thus, as α is increased from zero, long-range links become increasingly improbable, reducing the network synchronizability. An interesting question is what happens if α is negative. Intuitively, it is expected that for a negative value of α, long-range links among clusters can be much more probable than short-range links; as a result, network synchronizability should improve as α is decreased from zero. But is this really the case?

Insight has been provided in Ref. [34], where a ring clustered-network model is employed and different networks are generated for visualization with several different values of α, as shown in Fig. 4(a)−(c). For Fig. 4(a), the value of α is positive so that short-range links among the clusters are favored. In this case, the average network distance can be large. For α = 0 [Fig. 4(b)], short-range and long-range links are equally probable, making the connections among the clusters small-world like with the small average network distance. For α < 0, long-range links are favored but, the most favorable links are those that are squarely across the ring configuration. For example, for a circular ring, the inter-cluster links close to the diameter of the circle are the most favorable ones, as shown in Fig. 4(c). This, in fact, makes the average network distance large.

In Ref. [11], it is shown that scale-free networks are generally more difficult to be synchronized, despite their smaller average network distances as compared with small-world networks. Intuitively, this is mainly because of the highly heterogeneous degree distribution in scale-free networks that stipulates the existence of a small subset of nodes with extraordinarily larger numbers of links as compared with most nodes in networks. It is speculated in Ref. [11] that communication can be blocked at these nodes, significantly reducing the synchronizability of the whole network as compared with networks with more homogeneous degree distributions, such as small-world networks or random networks. However, for networks with similar characteristics, either homogeneous or heterogeneous, the average network distance (or diameter) is the determining factor for synchronizability [3]. These considerations suggest that, quite counterintuitively, as α becomes negative from a positive value, the network synchronizability is expected to increase and reach its maximum for α = 0, and then to decrease as α is decreased from zero.

For the ring clustered-network configuration, estimates of the average network distance can be obtained for some limiting cases. In particular, for α → ∞, there is a strong tendency for a cluster to connect only to its nearest neighboring clusters. For such a configuration with a large number (M) of clusters, the average network distance is d ∼ M/4. For α = 0, the probabilities for a cluster to connect to other clusters are equal, so the inter-cluster links appear random. In this case, the average network distance is d ∼ ln M, as for random networks [23]. In the limiting case where α → −∞, there is a tendency for a cluster to connect to diametrically opposite clusters as seen in Fig. 4(c), making most links nearly pass through the center of the ring configuration. In this case, the average network distance is d∼M/8. For reasonably large values of M, it can be obtained

suggesting that ring clustered networks with α near zero are most synchronizable.

For the ring clustered network model, the eigenratio and the average network distance versus α are shown in Fig. 5. The network consists of M clusters, each containing five nodes that are connected in a one-to-all manner. For comparison, the link ratio (the ratio between the number of links among the clusters and the total number of links in the network) is fixed at a small value (0.01 for Fig. 5). In Fig. 5(a), the eigenratio versus α is shown for three values of M. It is seen that for relatively large value of M, the eigenratio exhibits the expected behavior in that it decreases as α is decreased from a positive value to zero, and it increases as α becomes negative from zero. Figure 5(b) shows the behavior of the average network distance, which is consistent with that exhibited by the eigenratio, as expected.

Similar transition behavior persists in internally scale-free clustered networks, hierarchical clustered networks and Zachary networks [34].

4 Quantitative analysis of synchronization in clustered networks

The preceding section discusses how the distribution of inter-cluster links affects the synchronizability of the network. Here, it will be examined, for a given distribution of links, how the number of links affects synchronizability. The emphasis will be placed on a random clustered network model for developing a theory and providing numerical support that theoretical predictions are typical for clustered networks.

To be concrete, the following clustered network model is considered: N nodes are classified into M groups, where each group has n = N/M nodes. In a group, a pair of nodes are connected with probability p_s, and nodes of different groups are connected with probability p_l. This forms a clustered random network. Typically, the number of inter-connections is typically far less than that of intra-connections. As a result, the parameter region of small p_l values is more relevant.

Theoretical analysis should be developed for synchronization, which yields the stability regions for synchronization in the two-dimensional parameter space defined by the probabilities of the two types of links. The analytical predictions are verified by direct numerical simulations of corresponding dynamical networks. The following coupling scheme is considered: for any i (1 ≤ i ≤ N), G_ii = 1, G_ij = −1/k_i if there is a link between node i and j, and G_ij = 0 otherwise, where k_i is the degree of node i. The coupling matrix G is not symmetric since G_ij = −1/k_i while G_ji = −1/k_j. The Gerschgorin theorem stipulates that all the eigenvalues be located within a disc centered at 1 with radius 1, thus λ_N ≤ 2. One of the synchronization conditions, , can usually be satisfied. Thus the synchronizability of the system is determined by λ₂. In the following section, a theoretical formula will be derived to understand the dependence of λ₂ on p_l and p_s for small values of p_l, the typical parameter regime for realistic clustered networks.

4.1 Dependence of λ₂ on p_l and p_s

For a clustered network, the components of the eigenvector e₂ have approximately the same value within any cluster, while they can be quite different among different clusters, as illustrated in Fig. 6. It is written that {\mathit{e}}_2\approx {[{\tilde{e}}_1\mathrm{,}\cdots \mathrm{,}{\tilde{e}}_1\mathrm{,}{\tilde{e}}_2\mathrm{,}\cdots \mathrm{,}{\tilde{e}}_2\mathrm{,}\cdots \mathrm{,}{\tilde{e}}_M\mathrm{,}\cdots \mathrm{,}{\tilde{e}}_M]}^{\text{T}}, and for each I, 1 ≤ I ≤ M, there are n {\tilde{e}}_I’\text{s} in e₂. By definition, G · e₂ = λ₂e₂ and e₂ · e₂ = 1, {\lambda }_2={\mathit{e}}_2^{\text{T}}\cdot \mathit{G}\cdot {\mathit{e}}_2={\displaystyle \sum _{i\mathrm{,}j=1}^N{e}_{2i}{G}_{ij}{e}_{2j}} is obtained, where e_2i is the ith component of e₂. Expanding the summation in j gives

\begin{array}{}\text{(27)}& \begin{array}{ll}{\lambda }_2=\hfill & {\displaystyle \sum _{i=1}^N}{e}_{2i}\{G_{i1}{\tilde{e}}_1+G_{i2}{\tilde{e}}_1+\cdots +G_{in}{\tilde{e}}_1\hfill \\ \hfill & +G_{in+1}{\tilde{e}}_2+\cdots +G_{iN}{\tilde{e}}_M\}\hfill \end{array}\end{array}

It is recalled that G_ii = 1; and if i and j belong to the same cluster, G_ij equals −1/k_i with probability p_s and 0 with probability 1−p_s; while if i and j belong to different clusters, G_ij equals−1/k_i with probability p_l and 0 with probability 1−p_l, where k_i is the degree of node i. Thus, the following equation can be obtained:

\begin{array}{}\text{(27)}& \begin{array}{ll}{\lambda }_2=\hfill & {\displaystyle \sum _{i=1}^N}{e}_{2i}(-n\frac{p_l}{k_i}{\tilde{e}}_1-n\frac{p_l}{k_i}{\tilde{e}}_2+\cdots \hfill \\ \hfill & +{\tilde{e}}_I-n\frac{p_s}{k_i}{\tilde{e}}_I+\cdots -n\frac{p_l}{k_i}{\tilde{e}}_M)\hfill \end{array}\end{array}

where {\tilde{e}}_I is the eigenvector component value corresponding to the cluster that contains node i. When 1−np_s/k_i = (N−n)p_l/k_i is noted, the following equation can be obtained

\begin{array}{}\text{(27)}& \begin{array}{ll}{\lambda }_2\hfill & ={\displaystyle \sum _{i=1}^N}{e}_{2i}\left[(N-n)\frac{p_l}{k_i}{\tilde{e}}_I-n\frac{p_l}{k_i}{\displaystyle \sum _{J\ne I}^M}{\tilde{e}}_J\right]\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^N}{e}_{2i}\left(N\frac{p_l}{k_i}{\tilde{e}}_I-n\frac{p_l}{k_i}{\displaystyle \sum _{J=1}^M}{\tilde{e}}_J\right)\hfill \end{array}\end{array}

For the clustered random network models, the degree distribution has a narrow peak centered at k = np_s + (N−n)p_l, thus k_i ≈ k. The summation over i can be carried out as follows:

\begin{array}{}\text{(27)}& \begin{array}{ll}{\lambda }_2\hfill & \approx {\displaystyle \sum _{I=1}^{M}n{\tilde{e}}_I\left(N\frac{p_l}{k}{\tilde{e}}_I-n\frac{p_l}{k}{\displaystyle \sum _{J=1}^M{\tilde{e}}_J}\right)}\hfill \\ \hfill & =N\frac{p_l}{k}{\displaystyle \sum _{I=1}^{M}n{\tilde{e}}_I^2}-{\left(n{\displaystyle \sum _{J=1}^M{\tilde{e}}_J}\right)}^2\frac{p_l}{k}\hfill \end{array}\end{array}

Since \begin{array}{l}{\displaystyle \sum _{I=1}^M}n{\tilde{e}}_I^2\approx {\displaystyle \sum _{i=1}^N}{e}_{2i}^2=1\hfill \end{array}, and \begin{array}{l}n{\displaystyle \sum _{J=1}^M}{\tilde{e}}_J={\displaystyle \sum _{i=1}^N}{e}_{2i}\hfill \end{array}, it is obtained

\begin{array}{}\text{(28)}& {\lambda }_2\approx \frac{N{p}_l}{n{p}_s+(N-n)p_l}-{\left({\displaystyle \sum _{i=1}^N{e}_{2i}}\right)}^2\frac{p_l}{k}\end{array}

The normalized eigenvector e₁ of λ₁ corresponds to the synchronized state, so its components have constant values: {\mathit{e}}_1=[1/\sqrt{N}\mathrm{,}\cdots \mathrm{,1}/\sqrt{N}{]}^{\text{T}}. If G is symmetric, eigenvectors associated with different eigenvalues are orthogonal: {\mathit{e}}_i\cdot {\mathit{e}}_j={\text{δ}}_{ij}, where {\text{δ}}_{ij}=1 for i = j and 0 else. When i = 1 and j = 2, {\displaystyle \sum _{l=1}^N{e}_{2l}=0} is obtained. If the coupling matrix G is slightly asymmetric, \begin{array}{l}{\displaystyle \sum _{i=1}^N}{e}_{2i}\hfill \end{array} is nonzero but small, and the second term in Eq. (28) can be omitted, thus obtaining

\begin{array}{}\text{(29)}& {\lambda }_2\approx \frac{N{p}_l}{n{p}_s+(N-n)p_l}\end{array}

For fixed p_l and large p_s, λ₂ decreases as p_s increases, indicating that the network becomes more difficult to be synchronized. This is an abnormal behavior in the network synchronizability, which will be verified numerically. Furthermore, since λ₂ depends only on the ratio of p_l/p_s in Eq. (29), the synchronization-desynchronization boundaries in the (p_s, p_l) parameter plane should consist of straight-line segments.

The above analysis can be extended to more general clustered networks, i.e., those with different cluster sizes or heterogeneous degree distributions in each cluster, by replacing n with n_I−the size of the Ith cluster−for each I, and using the degree distribution P_I(k) in the summation over 1/k. In this case, p_s and p_l can be regarded as effective parameters, and may vary in different clusters. A formula similar to Eq. (29) can be obtained, because even in such a case, the contribution of the second term in Eq. (28) to λ₂ is small. Thus, the abnormal synchronization phenomenon is due to the clustered network structure; it does not depend on the details of dynamics.

4.2 Numerical support with coupled Rössler networks

Depending on the initial conditions and the network realization, the Rössler system may have desynchronization bursts [52, 53]. It is thus necessary to characterize the network synchronizability statistically. P_syn is defined as the probability that the fluctuation width of the system W (t) is smaller than a small number δ (chosen somewhat arbitrarily) at all times during a long observational period T₀ in the steady state, say, from T₁ to T₁ + T₀, where W(t)=\langle |x(t)-\langle x(t)\rangle |\rangle, and \langle \cdot \rangle denotes the average over the nodes of the network. If δ is small enough, the system can be deemed as being synchronized in the period T₀, so P_syn is the probability of synchronization of the system in the period T₀, with P_syn = 1 if the network for the given parameters can synchronize. Practically, P_syn can be calculated by the ensemble average, i.e. the ratio of the number of synchronized cases over the number of all random network realizations. Since P_syn can change drastically from 0 to 1 in a small region in the parameter space, it is possible to define the boundary between the synchronizable region and the unsynchronizable region as follows: for a fixed p_s, the boundary value p_lb is such that the quantity

\begin{array}{}\text{(29)}& \parallel \nabla P_{syn}(p_s\mathrm{,}{p}_l)\parallel \equiv {\sqrt{{\left(\frac{\partial P_{syn}}{\partial p_s}\right)}^2+{\left(\frac{\partial P_{syn}}{\partial p_l}\right)}^2}|}_{(p_s\mathrm{,}{p}_l)}\end{array}

is maximized at (p_s, p_lb). Figure 7 shows the synchronization boundary in the parameter space (p_s, p_l) from both numerical calculation and theoretical prediction of Eqs. (6) and (7). It can be seen that both results agree with each other. If the number of inter-cluster connections is fixed, say, p_l = 0.2, as the number of intra-cluster links exceeds a certain value (about 0.78), the system becomes desynchronized. When p_s is small, the number of the inter-cluster connections and the number of the intra-cluster connections are approximately matched, and the networks are synchronizable. As p_s becomes larger, the matching condition deteriorates and the network loses its synchronizability, even though its average distance becomes smaller. That is, too many intra-cluster links tend to destroy the global synchronization. The same phenomenon persists for different parameter values. This is precisely the abnormal synchronization phenomenon predicted by theory.

4.3 Numerical support with coupled logistic-map networks

For the coupled logistic-map network, if the system is synchronizable, starting from a random initial condition, it will approach the synchronization state. In the simulation, synchronization is defined as , where \langle \cdot \rangle denotes the average over the network. The average time T required for the system to become synchronized can be conveniently used to characterize the ability of the system to synchronize. If the system is unsynchronizable, the time T is infinite. Figure 8 shows the behavior of T in the two-dimensional parameter space (p_l, p_s) for networks with two clusters (a) and ten clusters (b). This gives the synchronizable region (grey regions in Fig. 8) in the parameter space in which the system is able to synchronize within a certain time, and the unsynchronizable region (white regions in Fig. 8). The shape of the figure depends on the coupling strength ε and the contour lines of λ₂ and λ_N. For the two-cluster network, if ε = 1, the shape appears to be symmetric, while if ε < 1, the boundary is asymmetric. Figure 8(a) demonstrates that for a given p_l (e.g., 0.2) as p_s is increased from 0.2, synchronization time T is also increased, and at a certain point (about 0.75 in this case), the system becomes unsynchronizable. The same phenomenon persists for different networks and dynamical parameters. Again, when the number of inter-cluster links is fixed, too many intra-cluster links violate the matching condition and thus tend to destroy the global synchronization.

One remark concerning the physical meaning of the result, as exemplified by Figs. 7 and 8, is in order. Consider two clustered networks where (A) the two types of links are approximately matched and (B) there is a substantial mismatch. Theory predicts that network A is more synchronizable than network B. This statement is meaningful in a probabilistic sense, as whether or not a specific system may achieve synchronization is also determined by many other factors such as the choice of the initial condition, possible existence of multiple synchronized states, and noise, etc. Our result means that, under the influence of these random factors, there is a higher probability for network A to be synchronized than for network B.

5 Discussion

In summary, some recent results about the synchronizability of complex networks with a cluster structure have been reviewed. The first result is that random, long-range couplings among clusters can enhance the synchronizability, while connections among nodes within individual clusters have little impact on the network’s ability to synchronize. In terms of the relationship between the synchronizability and the number of clusters in the network, an interesting transition phenomenon is uncovered, where the network synchronizability exhibits a different behavior depending on parameter that controls the probability of random, long-range links among the clusters. In particular, when these links are less probable, the synchronizability tends to deteriorate as the number of clusters is increased; the opposite occurs when the links are more probable. There has been some theoretical understanding of this phenomenon based on the analysis of a class of simplified networks with clusters distributed according to a ring topology. These findings imply that, in the context of social networks, a viable strategy to achieve synchronization is to devote resources to establishing and enhancing connections among distant communities.

The second result concerns theoretical and numerical evidence that the optimal synchronization of complex clustered networks can be achieved by matching the probabilities of inter-cluster and intra-cluster links. That is, at a global level, the network has the strongest synchronizability when these probabilities are approximately equal. Overwhelmingly, strong intra-cluster connection can counterintuitively weaken the network synchronizability. This phenomenon persists for another typical coupling scheme, i.e., for any i (1 ≤ i ≤ N), G_ii = k_i, G_ij = −1 if there is a link between node i and j, and G_ij = 0 otherwise. A new set of analysis and numerical justification has been provided in Ref. [36]. While the network model used to arrive at this result is somewhat idealized, it can be argued that similar phenomena should persist in more general clustered networks. In real systems with a clustered structure, if global synchronization is the best performance of the system, special attention needs to be devoted to distinguishing the inter-cluster and intra-cluster connections as in this case, where a proper distribution of the links is more efficient than adding links blindly.

For biological networks, such as metabolic networks and protein-protein interaction networks, certain nodes may have many more links than others, which forms a hierarchical cluster structure [58]. This indicates a power-law distribution of the degree k: P(k) ∼ k^−γ, and the network is scale-free. Therefore, it is interesting to study clustered scale-free networks, in which each cluster contains a scale-free subnetwork. The synchronizability of such clustered networks has been studied. In particular, for each cluster, the subnetwork is generated via the preferential attachment rule [24]. Initially, there is a fully connected small subset of size m₀, then a new node is added with m links, and the probability that a previous node i is connected to this new node is proportional to its current degree k_i. New nodes are continuously added until a prescribed network size n is reached. In our simulation, m₀ = 2m + 1 is set so that the average degree of this network is 2m. Given M such scale-free subnetworks, each pair of nodes in different clusters are connected with probability p_l. For this model, p_l controls the number of inter-clustered links, and m controls the number of intra-clustered links. Numerical simulations have been carried out, and it is found that the patterns for the eigenvalues λ_N and λ₂ are essentially the same as that for the clustered network where each cluster contains a random subnetwork. This indicates that optimization of synchronization by matching different types of links is a general phenomenon, regardless of the detailed topology in each individual cluster.

Some networks, e.g., intercellular communication networks, may have a locally regular structure. For example, a tissue network can be defined such that the nodes are cells and the links are the interactions between cells, i.e., the transmission of signal molecules. These interactions mainly occur between adjacent cells, which form a locally regular linkage structure. In addition, larger diffusing growth factors provide long-range links. A clustered network with each cluster having a regular backbone is thus a plausible model for biological tissue organization. The synchronizability of such clustered regular networks has been studied. First, for each cluster, a one-dimensional regular lattice with the periodic boundary condition is generated, i.e. each node is connected with 2m of its neighbors. For example, node i connects with nodes i−m, i−m + 1, ···, i−1, i + 1, ···, i + m. Then each pair of nodes in different clusters are connected with probability p_l. Thus p_l and m are the control parameters for the number of inter- and intra-clustered links, respectively. Numerical simulations show that for large m and small p_l, which is typical for clustered networks, optimization of synchronization can also be achieved by constraining m such that the number of inter- and intra- links are approximately matched. For intermediate m values, an interesting synchronization phenomenon is uncovered. That is, for a locally regular clustered networks the synchronizability exhibits an alternating, highly non-monotonic behavior as a function of the intra-cluster link density. In fact, there are distinct regions of the density for which the network synchronizability is maximized, but there are also parameter regions in between for which the synchronizability is diminished [59].

A basic assumption in the existing works is that all the clusters in a network are on the equal footing in the sense that their sizes are identical and the interactions between any pair of clusters are symmetrical. In realistic applications the distribution of the cluster size can be highly uneven. For example, in a clustered network with a hierarchical structure, the size of a cluster can in general depend on the particular hierarchy to which it belongs. More importantly, the interactions between clusters in different hierarchies can be highly asymmetrical. For instance, the coupling from a cluster at the top of the hierarchy to a cluster in a lower hierarchy can be much stronger than the other way around. An asymmetrically interacting network can in general be regarded as the superposition of a symmetrically coupled network and a directed network, both being weighted. A weighted, directed network is actually a gradient network [60, 61], a class of networks for which the interactions or couplings among nodes are governed by some gradient field on the network. For a complex gradient network, a key parameter is the strength of the gradient field (the extent of the directness of links), denoted by g. A central issue is how the network synchronizability depends on g. As g is increased, the interactions among various clusters in the network become more directed. From a dynamical-system point of view, uni-directionally coupled systems often possess strong synchronizability [62]. Thus, intuitively, we expect to observe enhancement of the network synchronizability with the increase of g. The question is whether there exists an optimal value of g for which the network synchronizability can be maximized. This is in fact the problem of optimizing synchronization in clustered gradient networks, and our recent findings [63] suggest an affirmative answer to the question. In particular, we are able to obtain solid analytic insights into a key quantity that determines the network synchronizability. The theoretical formulas are verified by both numerical eigenvalue analysis and direct simulation of oscillatory dynamics on the network. The existence of an optimal state for gradient clustered networks to achieve synchronization may have broad implications for fundamental issues such as the evolution of biological networks and for practical applications such as the design of efficient computer networks.

After all, the general observation is that the synchronizability of the clustered networks is mainly determined by the underlying clustered structure. Insofar as there is a clustered structure, details such as the link topology within each cluster, node dynamics and parameters, etc. do not appear to have a significant influence on the synchronization of the coupled oscillator networks supported by the clustered backbone. A practical usage is that, even if the details about the dynamics of a biological system are not available, insofar as the underlying network has a clustered structure, it is possible to make predictions about synchronization of the network. Such insights may begin to provide, for example, first principles for the organizational dynamics of normal and abnormal (i.e. cancer) tissue which currently remain largely unknown.

The clustered topology has also been identified in technological networks such as certain electronic circuit networks and computer networks [43−45]. For a computer network, the main functions include executing sophisticated codes to carry out extensive computations. Suppose a large-scale, parallel computational task is to be accomplished by the network, for which synchronous timing is of paramount importance. Our results can provide useful clues as to how to design the network to achieve the best possible synchronization and consequently optimal computational efficiency.

References

[1] Lago-Fernandez L. F., Huerta R., Corbacho F., and Siguenza J. A., Phys. Rev. Lett., 2000, 84(12): 2758

[2] Gade P. M. and Hu C.-K., Phys. Rev. E, 2000, 62(5): 6409

[3] Wang X. F. and Chen G., Int. J. of Bifur. Chaos, 2002, 12: 187

[4] Wang X. F. and Chen G., IEEE Trans. on Circ. Sys., Part I, 2002, 49: 54

[5] Hong H., Choi M. Y., and Kim B. J., Phys. Rev. E, 2002, 65(2): 026139

[6] Jalan S. and Amritkar R. E., Phys. Rev. Lett., 2003, 90(1): 014101

[7] Barahona M. and Pecora L. M., Phys. Rev. Lett., 2002, 89: 054101

[8] Lv J., Yu X., Chen G., and Cheng D., IEEE Trans. on Circ. Sys., Part I, 2004, 51: 787

[9] Fan J. and Wang X. F., Physica A, 2005, 349: 443

[10] Gong B., Yang L., and Yang K., Phys. Rev. E, 2005, 72: 037101

[11] Nishikawa T., Motter A. E., Lai Y.-C., and Hoppensteadt F. C., Phys. Rev. Lett., 2003, 91: 014101

[12] Motter A. E., Zhou C., and Kurths J., Europhys. Lett., 2005, 69(3): 334

[13] Motter A. E., Zhou C., and Kurths J., Phys. Rev. E, 2005, 71: 016116

[14] Zhou C., Motter A. E., and Kurths J., Phys. Rev. Lett., 2006, 96: 034101

[15] Zhou C. and Kurths J., Phys. Rev. Lett., 2006, 96: 164102

[16] Hwang D.-U., Chavez M., Amann A., and Boccaletti S., Phys. Rev. Lett., 2005, 94: 138701

[17] Chavez M., Hwang D.-U., Amann A., Hentschel H. G. E., and Boccaletti S., Phys. Rev. Lett., 2005, 94: 218701

[18] Huang D., Phys. Rev. E, 2006, 74: 046208

[19] Chavez M., Hwang D.-U., Martinerie J., and Boccaletti S., Phys. Rev. E, 2006, 74: 066107

[20] Wang X., Lai Y.-C., and Lai C.-H., Phys. Rev. E, 2007, 75: 056205

[21] Nishikawa T. and Motter A. E., Phys. Rev. E, 2006, 73: 065106

[22] Watts D. J. and Strogatz S. H., Nature, 1998, 393: 440

[23] Erdös P. and Rényi A., Publ. Math. Inst. Hung. Acad. Sci., 1960, 5: 17; Bollobaás B., Random Graphs, London: Academic, 1985

[24] Barabási A.-L. and Albert R., Science, 1999, 286: 509

[25] Zhao M., Zhou T., Wang B.-H., and Wang W.-X., Phys. Rev. E, 2005, 72: 057102

[26] Yin C.-Y., Wang W.-X., Chen G., and Wang B.-H., Phys. Rev. E, 2006, 74: 047102

[27] Atay F. M. and Bıyıko˘glu T., Phys. Rev. E, 2005, 72: 016217

[28] Donetti L., Hurtado I., and Mu˜noz M. A., Phys. Rev. Lett., 2005, 95: 188701

[29] Atay F. M., Bıyıko˘glu T., and Jost J., IEEE Trans. Circuits Syst. I: Regular Papers, 2006, 53: 92

[30] Wang X. F. and Chen G., J. Systems Science and Complexity, 2003, 16: 1

[31] Boccaletti S., Hwang D.-U., Chavez M., Amann A., Kurths J., and Pecora L. M., Phys. Rev. E, 2006, 74: 016102

[32] Restrepo J. G., Ott E., and Hunt B. R., Phys. Rev. Lett., 2006, 97: 094102

[33] Oh E., Rho K., Hong H., and Kahng B., Phys. Rev. E, 2005, 72: 047101

[34] Park K., Lai Y.-C., Gupte S., and Kim J.-W., Chaos, 2006, 6: 015105

[35] Huang L., Park K., Lai Y.-C., Yang L., and Yang K., Phys. Rev. Lett., 2006, 97: 164101

[36] Huang L., Lai Y.-C., Park K., and Gatenby R. A., submitted

[37] Watts D. J., Dodds S., and Newman M. E. J., Science, 2002, 296: 1302

[38] Girvan M. and Newman M. E. J., Proc. Natl. Acad. Sci. USA, 2002, 99: 7821

[39] Motter A. E., Nishikawa T., and Lai Y.-C., Phys. Rev. E, 2003, 68: 036105

[40] Spirin V. and Mirny L. A., Proc. Natl. Acad. Sci. USA, 2003, 100: 12123

[41] Ravasz E., Somera A. L., Mongru D. A., Oltvai Z., and Barabási A.-L., Science, 2002, 297: 1551

[42] Palla G., Derényi I., Farkas I., and Vicsek T., Nature, 2005, 435: 814

[43] Milo R., Shen-Orr S., Itzkovitz S., Kashtan N., Chklovskii D., and Alon U., Science, 2002, 298(5594): 824

[44] Vázquez A., Pastor-Satorras R., and Vespignani A., Phys. Rev. E, 2002, 65(6): 066130

[45] Eriksen K. A., Simonsen I., Maslov S., and Sneppen K., Phys. Rev. Lett., 2003, 90: 148701

[46] Zachary W. W., J. Anthropol. Res., 1977, 33:452

[47] Fujisaka H. and Yamada T., Prog. Theor. Phys., 1983, 69: 32

[48] Pecora L. M. and Carroll T. L., Phys. Rev. Lett., 1998, 80: 2109

[49] Jost J. and Joy M. P., Phys. Rev. E, 2002, 65: 016201

[50] Fink K. S., Johnson G., Carroll T., Mar D., and Pecora L., Phys. Rev. E, 2000, 61: 5080

[51] Monasson R., Eur. Phys. J. B, 1999, 12: 555

[52] Restrepo J. G., Ott E., and Hunt B. R., Phys. Rev. Lett., 2004, 93: 114101

[53] Restrepo J. G., Ott E., and Hunt B. R., Phys. Rev. E, 2004, 69: 66215

[54] Wigner E. P., Ann. Math., 1955, 62: 548

[55] Wigner E. P., Ann. Math., 1957, 65: 203

[56] Mehta M. L., Random Matrices, 2nd ed., New York: Academic, 1991

[57] Farkas I. J., Derényi I., Barabási A.-L., and Vicsek T., Phys. Rev. E, 2001, 64: 026704

[58] Barabási A.-L. and Oltvai Z. N., Nature Reviews—Genetics, 2004, 5: 101

[59] Huang L., Lai Y.-C., and Gatenby R. A., submitted

[60] Toroczkai Z. and Bassler K.E., Nature, 2004, 428: 716

[61] Park K., Lai Y.-C., Zhao L., and Ye N., Phys. Rev. E, 2005, 71: 065105

[62] Kocarev L. and Parlitz U., Phys. Rev. Lett., 1996, 76: 1816

[63] Wang X., Huang L., Lai Y.-C., and Lai C.-H., submitted

Rights & permissions

Higher Education Press and Springer-Verlag 2007