Physics Reports 469(2008)93-153 Contents lists available at Science Direct Physics Reports LSEVIER ournalhomepagewww.elsevier.com/locate/physrep Synchronization in complex networks Alex Arenas a Albert Diaz-Guilera. Jurgen Kurths d, e, Yamir Moreno b. 4*, Changsong Zhou& Departament d'Enginyeria Informatica i Matematiques, Universitat Rovira i Virgili, 43007 Tarragona, Spair Institute for Biocomputation and Physics of Complex Systems(BIFI). University of zaragoza, Zaragoza 50009, Spain dEpartament de Fisica Fonamental, Universitat de barcelona, 08028 Barcelona, Spain d institute of physics Humboldt University, D-12489 Berlin, Newtonstrasse 15, Germany Potsdam Institute for Climate Impact Research, 14412 Potsdam, PF 601203, Germany f Department of Theoretical Physics, University of zaragoza, Zaragoza 50009, Spain epartment of physics and Centre for Nonlinear Studies, Hong Kong Baptist University, Kowloon Tong. Hong Kong china ARTICLE INFO ABSTRACT Article Synchronization processes in populations of locally interacting elements are the focus of ccepted 17 September 2008 intense research in physical, biological, chemical, technological and social systems. The Available online 24 September 2008 many efforts devoted to understanding synchronization phenomena in natural systems ake advantage of the recent theory of complex networks. In this review, we report the advances in the comprehension of synchronization phenomena when oscillating elements are constrained to interact in a complex network topology. We also take an overview of the new emergent features coming out from the interplay between the structure and the function of the underlying patterns of connections. Extensive numerical work as well as analytical approaches to the problem are presented. Finally, we review several application of synchronization in complex networks to different disciplines: biological systems and Complex networks euroscience, engineering and computer science, and economy and social sciences. o 2008 Elsevier B V. All rights reserved Contents 1. Introduction 2. Complex networks in a nutshell 9 3. Coupled phase oscillator models on complex networks 3. 1. Phase oscillators 3. 1.1. The Kuramoto model 3.1. 2. Kuramoto model on complex networks 1.4. Path towards synchronization in complex network 1.5. Kuramoto model on structured or modular network 3. 1.6. Synchronization by pacemakers plex networks 1. Master stability function formalism 112 lity and master stability function Corresponding author at: Institute for Biocomputation and Physics of Complex Systems(BIFI), University of Zaragoza, Zaragoza 50009, Spain. Tel +34 370-1573S-see front matter o 2008 Elsevier B V. All rights reserved. i:10.1016 physrep.200809002
Physics Reports 469 (2008) 93–153 Contents lists available at ScienceDirect Physics Reports journal homepage: www.elsevier.com/locate/physrep Synchronization in complex networks Alex Arenas a,b , Albert Díaz-Guilera c,b , Jurgen Kurths d,e , Yamir Moreno b,f,∗ , Changsong Zhoug a Departament d’Enginyeria Informàtica i Matemàtiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain b Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, Zaragoza 50009, Spain c Departament de Física Fonamental, Universitat de Barcelona, 08028 Barcelona, Spain d Institute of Physics, Humboldt University, D-12489 Berlin, Newtonstrasse 15, Germany e Potsdam Institute for Climate Impact Research, 14412 Potsdam, PF 601203, Germany f Department of Theoretical Physics, University of Zaragoza, Zaragoza 50009, Spain g Department of Physics and Centre for Nonlinear Studies, Hong Kong Baptist University, Kowloon Tong, Hong Kong, China a r t i c l e i n f o Article history: Accepted 17 September 2008 Available online 24 September 2008 editor: I. Procaccia PACS: 05.45.Xt 89.75.Fb 89.75.Hc Keywords: Synchronization Complex networks a b s t r a c t Synchronization processes in populations of locally interacting elements are the focus of intense research in physical, biological, chemical, technological and social systems. The many efforts devoted to understanding synchronization phenomena in natural systems now take advantage of the recent theory of complex networks. In this review, we report the advances in the comprehension of synchronization phenomena when oscillating elements are constrained to interact in a complex network topology. We also take an overview of the new emergent features coming out from the interplay between the structure and the function of the underlying patterns of connections. Extensive numerical work as well as analytical approaches to the problem are presented. Finally, we review several applications of synchronization in complex networks to different disciplines: biological systems and neuroscience, engineering and computer science, and economy and social sciences. © 2008 Elsevier B.V. All rights reserved. Contents 1. Introduction............................................................................................................................................................................................. 94 2. Complex networks in a nutshell ............................................................................................................................................................ 95 3. Coupled phase oscillator models on complex networks ...................................................................................................................... 96 3.1. Phase oscillators.......................................................................................................................................................................... 97 3.1.1. The Kuramoto model................................................................................................................................................... 97 3.1.2. Kuramoto model on complex networks..................................................................................................................... 98 3.1.3. Onset of synchronization in complex networks ........................................................................................................ 99 3.1.4. Path towards synchronization in complex networks................................................................................................ 103 3.1.5. Kuramoto model on structured or modular networks.............................................................................................. 105 3.1.6. Synchronization by pacemakers................................................................................................................................. 107 3.2. Pulse-coupled models ................................................................................................................................................................ 108 3.3. Coupled maps.............................................................................................................................................................................. 109 4. Stability of the synchronized state in complex networks .................................................................................................................... 112 4.1. Master stability function formalism.......................................................................................................................................... 112 4.1.1. Linear stability and master stability function............................................................................................................ 113 ∗ Corresponding author at: Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, Zaragoza 50009, Spain. Tel.: +34 976562212x223; fax: +34 976562215. E-mail address: yamir.moreno@gmail.com (Y. Moreno). 0370-1573/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.physrep.2008.09.002
A Arenas et aL / Physics Reports 469(2008)93-153 4.1.2. Measures of synchronizability.. 4.1.3. Synchronizability of 115 of networks 118 4. 1.5. Graph theoretical bo synchronizability 121 4.1.6. Synchronizability of weighted networks 4.1.7. Universal parameters controlling the synchronizability 4. 2. Design of synchronizable networks.. 4.2.1. Weighted couplings for enhancing synchronizability 4.2.2. Topological modification for enhancing synchronizability 2.3. Opti 4.3. Beyond the master stability function formalism 133 5.1. Biological systems and neuroscience 5.1.1. Genetic networks 5.1.2. Circadian rhythms 5.1.3. Ecology 5667 5.1.4. Neuronal networks 13 5.1.5. Cortical networks 5. 2. Computer science and engineering 5.2.1. Parallel/distributed computation .22.D 5.2.3. Consensus pi 5.2.4. Wireless co 5.2.5. Decentralized logistics 5.2.6. Power-grid 5.3. Social sciences and economy 3455 5.3.1. Opinion formation.. 5.3.2. Finance 533. World Trade Web Acknowledgments References 1. Introduction Synchronization, as an emerging phenomenon of a population of dynamically interacting units, has fascinated humans from ancestral times. Synchronization processes are ubiquitous in nature and play a very important role in many different contexts such as biology, ecology, climatology, sociology, technology, or even in arts [1, 2]. It is known that synchrony is rooted in human life from the metabolic processes in our cells to the highest cognitive tasks we perform as a group of individuals. For example, the effect of synchrony has been described in experiments of people communicating, or working gether with a background of shared, non-directive conversation, song or rhythm, or of groups of children interacting to an unconscious beat. In all cases the purpose of the common wave length or rhythm is to strengthen the group bond. The lack of such synchrony can indicate unconscious tension, when goals cannot be identified nor worked towards because the members are"out of sync"3]. Among the efforts for the scientific description of synchronization phenomena, there are several capital works that represented a breakthrough in our understanding of these phenomena. In 1665, the mathematician and physicist, inventor of the pendulum clock, C. Huygens, discovered an odd"kind of sympathy"in two pendulum clocks suspended side by side of each other. The pendulum clocks swung with exactly the same frequency and 180 out of phase; when the pendula were disturbed, the antiphase state was restored within half an hour and persisted indefinitely. huygens deduced that the crucial interaction for this effect came from"imperceptible movements"of the common frame supporting the two clocks. From that time on, the phenomenon became the focus of scientists. Synchronization involves, at least, two elements in interaction, and the behavior of a few interacting oscillators has been intensively studied in physics and mathematics literature. However, the phenomenon of synchronization of large populations is a different challenge and requires a different approach to be solved. We will focus our attention on this last challenge In the obituary of Arthur T Winfree, Strogatz [4 summarizes what can be considered the beginning of the modern quest to explain the synchronization of a population of interacting units: Wiener [5 posed a problem in his book Cybernetics: How is it that thousands of neurons or fireflies or crickets can suddenly fall into step with one another. all firing or flashing or chirping at the same time, without any leader or signal from the environment? wiener did not make significant mathematical progress on it, nor did anyone else until Winfree came along". Winfree[6] studied the nonlinear dynamics of a large population of weakly coupled limit-cycle oscillators with intrinsic frequencies that were distributed about so mean value, according to some prescribed probability distribution. The milestone here was to consider biological oscillators
94 A. Arenas et al. / Physics Reports 469 (2008) 93–153 4.1.2. Measures of synchronizability.................................................................................................................................... 114 4.1.3. Synchronizability of typical network models ............................................................................................................ 115 4.1.4. Synchronizability and structural characteristics of networks .................................................................................. 118 4.1.5. Graph theoretical bounds to synchronizability ......................................................................................................... 121 4.1.6. Synchronizability of weighted networks ................................................................................................................... 124 4.1.7. Universal parameters controlling the synchronizability .......................................................................................... 126 4.2. Design of synchronizable networks........................................................................................................................................... 128 4.2.1. Weighted couplings for enhancing synchronizability .............................................................................................. 128 4.2.2. Topological modification for enhancing synchronizability ...................................................................................... 131 4.2.3. Optimization of synchronizability.............................................................................................................................. 132 4.3. Beyond the master stability function formalism ...................................................................................................................... 133 5. Applications............................................................................................................................................................................................. 135 5.1. Biological systems and neuroscience ........................................................................................................................................ 135 5.1.1. Genetic networks......................................................................................................................................................... 135 5.1.2. Circadian rhythms ....................................................................................................................................................... 136 5.1.3. Ecology ......................................................................................................................................................................... 136 5.1.4. Neuronal networks ...................................................................................................................................................... 137 5.1.5. Cortical networks......................................................................................................................................................... 137 5.2. Computer science and engineering ........................................................................................................................................... 139 5.2.1. Parallel/distributed computation ............................................................................................................................... 140 5.2.2. Data mining.................................................................................................................................................................. 140 5.2.3. Consensus problems.................................................................................................................................................... 141 5.2.4. Wireless communication networks............................................................................................................................ 142 5.2.5. Decentralized logistics ................................................................................................................................................ 143 5.2.6. Power-grids.................................................................................................................................................................. 144 5.3. Social sciences and economy ..................................................................................................................................................... 145 5.3.1. Opinion formation ....................................................................................................................................................... 145 5.3.2. Finance ......................................................................................................................................................................... 146 5.3.3. World Trade Web ........................................................................................................................................................ 147 6. Perspectives............................................................................................................................................................................................. 147 7. Conclusions.............................................................................................................................................................................................. 148 Acknowledgments .................................................................................................................................................................................. 149 References................................................................................................................................................................................................ 149 1. Introduction Synchronization, as an emerging phenomenon of a population of dynamically interacting units, has fascinated humans from ancestral times. Synchronization processes are ubiquitous in nature and play a very important role in many different contexts such as biology, ecology, climatology, sociology, technology, or even in arts [1,2]. It is known that synchrony is rooted in human life from the metabolic processes in our cells to the highest cognitive tasks we perform as a group of individuals. For example, the effect of synchrony has been described in experiments of people communicating, or working together with a background of shared, non-directive conversation, song or rhythm, or of groups of children interacting to an unconscious beat. In all cases the purpose of the common wave length or rhythm is to strengthen the group bond. The lack of such synchrony can indicate unconscious tension, when goals cannot be identified nor worked towards because the members are ‘‘out of sync’’ [3]. Among the efforts for the scientific description of synchronization phenomena, there are several capital works that represented a breakthrough in our understanding of these phenomena. In 1665, the mathematician and physicist, inventor of the pendulum clock, C. Huygens, discovered an odd ‘‘kind of sympathy’’ in two pendulum clocks suspended side by side of each other. The pendulum clocks swung with exactly the same frequency and 180◦ out of phase; when the pendula were disturbed, the antiphase state was restored within half an hour and persisted indefinitely. Huygens deduced that the crucial interaction for this effect came from ‘‘imperceptible movements’’ of the common frame supporting the two clocks. From that time on, the phenomenon became the focus of scientists. Synchronization involves, at least, two elements in interaction, and the behavior of a few interacting oscillators has been intensively studied in physics and mathematics literature. However, the phenomenon of synchronization of large populations is a different challenge and requires a different approach to be solved. We will focus our attention on this last challenge. In the obituary of Arthur T. Winfree, Strogatz [4] summarizes what can be considered the beginning of the modern quest to explain the synchronization of a population of interacting units: ‘‘Wiener [5] posed a problem in his book Cybernetics: How is it that thousands of neurons or fireflies or crickets can suddenly fall into step with one another, all firing or flashing or chirping at the same time, without any leader or signal from the environment? Wiener did not make significant mathematical progress on it, nor did anyone else until Winfree came along’’. Winfree [6] studied the nonlinear dynamics of a large population of weakly coupled limit-cycle oscillators with intrinsic frequencies that were distributed about some mean value, according to some prescribed probability distribution. The milestone here was to consider biological oscillators
A. Arenas et aL Physics Reports 469(2008)93-153 Fig. 1. Small-world network construction from a regular lattice by rewiring links with a certain probability (randomness), as proposed by Watts and as phase oscillators, neglecting the amplitude. Working within the framework of a mean field model, winfree discovered that such a population of non-identical oscillators can exhibit a remarkable cooperative phenomenon. when the variance of the frequency distribution is large, the oscillators run incoherently, each one near its natural frequency. This behavior remains Then reducing the variance until a certain threshold is crossed. Below the threshold the oscillators begin to synchronize spontaneously(see[7). Note that the original winfree model was not solved analytically until recently 8. Although Winfree's approach proved to be successful in describing the emergence of spontaneous order in the system it was based on the premise that every oscillator feels the same pattern of interactions. However, this all-to-all connectivity between elements of a large population is difficult to conceive in the real world. When the number of elements is large enough, this pattern is incompatible with physical constraints as for example minimization of energy(or costs), and in local connectivity structure of the elements was missing(in fact, discarded )in these and subsequent approachs. particular general with the rare observation of long range interactions in systems formed by macroscopic elements. The In 1998. Watts and Strogatz presented a simple model of network structure, originally intended precisely to introduce the connectivity substrate in the problem of synchronization of cricket chirps, which show a high degree of coordination over long distances as though the insects were"invisibly"connected. Remarkably this work did not end in a new contribution to synchronization theory but as the seed for the modern theory of complex networks [9]. Starting with a regular lattice. they showed that adding a small number of random links reduces the distance between nodes drastically, see Fig. 1. This feature, known as small-world(SW)effect, had been first reported in an experiment conducted by Milgram [10 examining the average path length for social networks of people in the United States. Nowadays, the phenomenon has been detected in many other natural and artificial networks. the inherent complexity of the new model, from now on referred to as the Watts-Strogatz(WS) model, was in its mixed nature in between regular lattices and random graphs. Very soon, it turned out that the nature of many interaction patterns observed in scenarios as diverse as the Internet, the World-Wide Web, scientific collaboration networks and biological networks, was even more"complex"than the wS model. Most of them showed a heavy tailed distribution of connectivities with no characteristic scale. These networks have been since then called scale- free(SF)networks and the most connected nodes are called hubs. This novel structural complexity provoked an explosion of works, mainly from the physicists'community, since a completely new set of measures, models, and techniques, was needed to deal with these topological structures juring one decade we have witnessed the evolution of the field of complex networks, mainly from a static point of although some attempts to characterize the dynamical properties of complex networks have also been made. One of these dynamical implications, addressed since the very beginning of the subject, is the emergent phenomena of synchronization of a population of units with an oscillating behavior. The analysis of synchronization processes has benefited from the advance in the understanding of the topology of complex networks, but it has also contributed to the understanding of general emergent properties of networked systems. The main goal of this review is precisely to revise the research undertaken so far in order to understand how synchronization phenomena are affected by the topological substrate of interactions, in particular when this substrate is a complex network. The review is organized as follows. We first introduce the basic mathematical descriptors of complex networks that will be used henceforth. Next, we focus on the synchronization of populations of oscillators. Section IV is devoted to the analysis of the conditions for the stability of the fully synchronized state using the Master Stability Function(MSF)formalism Applications in different fields of science are presented afterwards and some perspectives provided Finally, the last section rounds off the review by giving our conclusions. 2. Complex networks in a nutshell There exist excellent reviews devoted to the structural characterization and evolution of complex networks [11-16]. Here we summarize the main features and standard measures used in complex networks. the goal is to provide the read with a brief overview of the subject as well as to introduce some notation that will be used throughout the review. The mathematical abstraction of a complex network is a graph g comprising a set of N nodes (or vertices)connected by a set of M links(or edges), being ki the degree(number of links )of node i. this graph is represented by the adjacency matrix A, with entries a i 1 if a directed link from j to i exists, and 0 otherwise. In the more general case of a weighted network, the graph is characterized by a matrix w, with entries wj, representing the strength(or weight)of the link fromj
A. Arenas et al. / Physics Reports 469 (2008) 93–153 95 Fig. 1. Small-world network construction from a regular lattice by rewiring links with a certain probability (randomness), as proposed by Watts and Strogatz [9]. as phase oscillators, neglecting the amplitude. Working within the framework of a mean field model, Winfree discovered that such a population of non-identical oscillators can exhibit a remarkable cooperative phenomenon. When the variance of the frequency distribution is large, the oscillators run incoherently, each one near its natural frequency. This behavior remains when reducing the variance until a certain threshold is crossed. Below the threshold the oscillators begin to synchronize spontaneously (see [7]). Note that the original Winfree model was not solved analytically until recently [8]. Although Winfree’s approach proved to be successful in describing the emergence of spontaneous order in the system, it was based on the premise that every oscillator feels the same pattern of interactions. However, this all-to-all connectivity between elements of a large population is difficult to conceive in the real world. When the number of elements is large enough, this pattern is incompatible with physical constraints as for example minimization of energy (or costs), and in general with the rare observation of long range interactions in systems formed by macroscopic elements. The particular local connectivity structure of the elements was missing (in fact, discarded) in these and subsequent approaches. In 1998, Watts and Strogatz presented a simple model of network structure, originally intended precisely to introduce the connectivity substrate in the problem of synchronization of cricket chirps, which show a high degree of coordination over long distances as though the insects were ‘‘invisibly’’ connected. Remarkably, this work did not end in a new contribution to synchronization theory but as the seed for the modern theory of complex networks [9]. Starting with a regular lattice, they showed that adding a small number of random links reduces the distance between nodes drastically, see Fig. 1. This feature, known as small-world (SW) effect, had been first reported in an experiment conducted by Milgram [10] examining the average path length for social networks of people in the United States. Nowadays, the phenomenon has been detected in many other natural and artificial networks. The inherent complexity of the new model, from now on referred to as the Watts–Strogatz (WS) model, was in its mixed nature in between regular lattices and random graphs. Very soon, it turned out that the nature of many interaction patterns observed in scenarios as diverse as the Internet, the World-Wide Web, scientific collaboration networks and biological networks, was even more ‘‘complex’’ than the WS model. Most of them showed a heavy tailed distribution of connectivities with no characteristic scale. These networks have been since then called scalefree (SF) networks and the most connected nodes are called hubs. This novel structural complexity provoked an explosion of works, mainly from the physicists’ community, since a completely new set of measures, models, and techniques, was needed to deal with these topological structures. During one decade we have witnessed the evolution of the field of complex networks, mainly from a static point of view, although some attempts to characterize the dynamical properties of complex networks have also been made. One of these dynamical implications, addressed since the very beginning of the subject, is the emergent phenomena of synchronization of a population of units with an oscillating behavior. The analysis of synchronization processes has benefited from the advance in the understanding of the topology of complex networks, but it has also contributed to the understanding of general emergent properties of networked systems. The main goal of this review is precisely to revise the research undertaken so far in order to understand how synchronization phenomena are affected by the topological substrate of interactions, in particular when this substrate is a complex network. The review is organized as follows. We first introduce the basic mathematical descriptors of complex networks that will be used henceforth. Next, we focus on the synchronization of populations of oscillators. Section IV is devoted to the analysis of the conditions for the stability of the fully synchronized state using the Master Stability Function (MSF) formalism. Applications in different fields of science are presented afterwards and some perspectives provided. Finally, the last section rounds off the review by giving our conclusions. 2. Complex networks in a nutshell There exist excellent reviews devoted to the structural characterization and evolution of complex networks [11–16]. Here we summarize the main features and standard measures used in complex networks. The goal is to provide the reader with a brief overview of the subject as well as to introduce some notation that will be used throughout the review. The mathematical abstraction of a complex network is a graph G comprising a set of N nodes (or vertices) connected by a set of M links (or edges), being ki the degree (number of links) of node i. This graph is represented by the adjacency matrix A, with entries aij = 1 if a directed link from j to i exists, and 0 otherwise. In the more general case of a weighted network, the graph is characterized by a matrix W, with entries wij, representing the strength (or weight) of the link from j
A Arenas et aL / Physics Reports 469(2008)93-153 to i. The investigation of the statistical properties of many natural and man-made complex networks revealed that, although representing very different systems, some categorization of them is possible. The most representative of these properties refers to the degree distribution P(k), that indicates the probability of a node to have a degree k. This fingerprint of complex networks has been taken for a long time to be its most differentiating factor. However, several other measures help to define the categorization. Examples are the average shortest path length (=(d;i where d;i is the length of the shortest path between node i and node j, and the clustering coefficient C that accounts for the fraction of actual triangles (three vertices forming a loop)over possible triangles in the graph The first classification of complex networks is related to the degree distribution P(k). The differentiation between homogeneous and heterogeneous networks AH is in general associated to the tail of the distribution. If it decays exponentially fast with the degree we refer to as homogeneous networks, the most representative example being the Erdos-Renyi(ER)random graph [17. On the contrary, when the tail is heavy one can say that the network is heterogeneous. In particular, SF networks are the class of networks whose distribution is a power-law, P(k)N kr, the Barabasi-Albert ( BA)model [18 being the paradigmatic model of this type of graph. This network is grown by a mechanism in which all incoming nodes are linked preferentially to the existing nodes. Note that the limiting case of lattices, or regular networks corresponds to a situation where all nodes have the same degree. e This categorization can be enriched by the behavior of e. For a lattice of dimension d containing N vertices, obviousy N/d For a random network, a rough estimate for e is also possible If the average number of nearest neighbors of a vertex is k then about k vertices of the network are at a distance e from the vertex or closer. Hence nN k and then e w In(N)/In(k), i.e. the average shortest-path length value is small even for very large networks. This smallness is usually referred to as the sw property. Associated to distances, there exist many measures that provide information about"centrality"of nodes. For instance, one can say that a node is central in terms of the relative distance to the rest of he network. One of the most frequently used centrality measures in the physics literature is the betweenness (load in some papers), that accounts for the number of shortest paths between any pair of nodes in the network that go through a given node or link The clustering coefficient C is also a discriminating property between different types of networks. It is usually calculated C G n; (1) where n; is the number of connections between nearest neighbors of node i, and k; is its degree. A large clustering coefficient implies many transitive connections and consequently redundant paths in the network, while a low C implies the opposite Finally, it is worth mentioning that many networks have a community structure, meaning that nodes are linked together densely connected groups between which connections are sparser. Finding the best partition of a network into communities is a very difficult problem. The most successful solutions, in terms of accuracy and computational cost [19 are those based on the optimization of a magnitude called modularity, proposed in [20], that precisely allows for the comparison of different partitionings of the network. the modularity of a given partition is, up to a multiplicative constant, the number of links falling within groups minus its expected number in an equivalent network with links placed at random. Given a network partitioned into communities, the mathematical definition of modularity is expressed in terms of the adjacency matrix qy nd the total number of links M=22iki as Q= ∑ is the community to which node i is assigned and the Kronecker delta function &c,c takes the value 1 if nodes re in the same community, and 0 otherwise. The larger the Q the more modular the network is. This property neralizations [278, 279] promises to be specially adequate to unveil structure-function relationships in complex 3. Coupled phase oscillator models on complex networks The need to understand synchronization, mainly in the context of biological neural networks, promoted the first studies of synchronization of coupled oscillators considering a network of interactions between them. In the late 80s, Strogatz and Mirollo [25 and later Niebur et al. 26] studied the collective synchronization of non-linear phase oscillators with random ntrinsic frequencies under a variety of coupling schemes in 2D lattices. beyond the differences with the actual conception of a complex network, the topologies studied in [26] can be thought of as a first approach to reveal how the complexity of the connectivity affects synchronization. the authors used a square lattice as a geometrical reference to construct three different connectivity schemes: four nearest neighbors, Gaussian connectivity truncated at 2o, and finally a random sparse connectivity. These results showed that random long-range connections lead to a more rapid and robust phase locking between oscillators than nearest-neighbor coupling or locally dense connection schemes. this observation is at the root of the recent findings about synchronization in complex networks of oscillators In the current section we review the results
96 A. Arenas et al. / Physics Reports 469 (2008) 93–153 to i. The investigation of the statistical properties of many natural and man-made complex networks revealed that, although representing very different systems, some categorization of them is possible. The most representative of these properties refers to the degree distribution P(k), that indicates the probability of a node to have a degree k. This fingerprint of complex networks has been taken for a long time to be its most differentiating factor. However, several other measures help to define the categorization. Examples are the average shortest path length ` = hdiji, where dij is the length of the shortest path between node i and node j, and the clustering coefficient C that accounts for the fraction of actual triangles (three vertices forming a loop) over possible triangles in the graph. The first classification of complex networks is related to the degree distribution P(k). The differentiation between homogeneous and heterogeneous networks AH is in general associated to the tail of the distribution. If it decays exponentially fast with the degree we refer to as homogeneous networks, the most representative example being the Erdös–Rényi (ER) random graph [17]. On the contrary, when the tail is heavy one can say that the network is heterogeneous. In particular, SF networks are the class of networks whose distribution is a power-law, P(k) ∼ k −γ , the Barabási–Albert (BA) model [18] being the paradigmatic model of this type of graph. This network is grown by a mechanism in which all incoming nodes are linked preferentially to the existing nodes. Note that the limiting case of lattices, or regular networks, corresponds to a situation where all nodes have the same degree. This categorization can be enriched by the behavior of `. For a lattice of dimension d containing N vertices, obviously, ` ∼ N 1/d . For a random network, a rough estimate for ` is also possible. If the average number of nearest neighbors of a vertex is k¯, then about k¯` vertices of the network are at a distance ` from the vertex or closer. Hence, N ∼ k¯` and then ` ∼ ln(N)/ ln(k¯), i.e. the average shortest-path length value is small even for very large networks. This smallness is usually referred to as the SW property. Associated to distances, there exist many measures that provide information about ‘‘centrality’’ of nodes. For instance, one can say that a node is central in terms of the relative distance to the rest of the network. One of the most frequently used centrality measures in the physics literature is the betweenness (load in some papers), that accounts for the number of shortest paths between any pair of nodes in the network that go through a given node or link. The clustering coefficient C is also a discriminating property between different types of networks. It is usually calculated as follows: C = 1 N XN i=1 Ci = 1 N XN i=1 ni ki(ki − 1)/2 , (1) where ni is the number of connections between nearest neighbors of node i, and ki is its degree. A large clustering coefficient implies many transitive connections and consequently redundant paths in the network, while a low C implies the opposite. Finally, it is worth mentioning that many networks have a community structure, meaning that nodes are linked together in densely connected groups between which connections are sparser. Finding the best partition of a network into communities is a very difficult problem. The most successful solutions, in terms of accuracy and computational cost [19], are those based on the optimization of a magnitude called modularity, proposed in [20], that precisely allows for the comparison of different partitionings of the network. The modularity of a given partition is, up to a multiplicative constant, the number of links falling within groups minus its expected number in an equivalent network with links placed at random. Given a network partitioned into communities, the mathematical definition of modularity is expressed in terms of the adjacency matrix aij and the total number of links M = 1 2 P i ki as Q = 1 2M X ij aij − kikj 2M δci ,cj (2) where ci is the community to which node i is assigned and the Kronecker delta function δci ,cj takes the value 1 if nodes i and j are in the same community, and 0 otherwise. The larger the Q the more modular the network is. This property and its generalizations [278,279] promises to be specially adequate to unveil structure–function relationships in complex networks [21–24]. 3. Coupled phase oscillator models on complex networks The need to understand synchronization, mainly in the context of biological neural networks, promoted the first studies of synchronization of coupled oscillators considering a network of interactions between them. In the late 80’s, Strogatz and Mirollo [25] and later Niebur et al. [26] studied the collective synchronization of non-linear phase oscillators with random intrinsic frequencies under a variety of coupling schemes in 2D lattices. Beyond the differences with the actual conception of a complex network, the topologies studied in [26] can be thought of as a first approach to reveal how the complexity of the connectivity affects synchronization. The authors used a square lattice as a geometrical reference to construct three different connectivity schemes: four nearest neighbors, Gaussian connectivity truncated at 2σ, and finally a random sparse connectivity. These results showed that random long-range connections lead to a more rapid and robust phase locking between oscillators than nearest-neighbor coupling or locally dense connection schemes. This observation is at the root of the recent findings about synchronization in complex networks of oscillators. In the current section we review the results
A Arenas et aL/Physics Reports 469(2008)93-153 obtained so far on three different kinds of oscillatory ensembles: limit cycle oscillators(Kuramoto), pulse-coupled models, and finally coupled map systems. We reserve for Section 4 those works that use the MSF formalism. Many other works whose major contribution is the understanding of synchronization phenomena in specific scenarios are discussed in the Appli ons 3. 1. Phase oscillators 3.1.1. The Kuramoto model The pioneering work by Winfree [6 spurred the field of collective synchronization and called for mathematical pproaches to tackle the problem. One of these approaches, as already stated, considers a system made up of a huge population of weakly-coupled nearly identical, interacting limit-cycle oscillators, where each oscillator exerts a phase dependent influence on the others and changes its rhythm according to a sensitivity function [ 27, 28 Even if these simplifications seem to be very crude, the phenomenology of the problem can be captured. Namely, the population of oscillators exhibits the dynamic analog to an equilibrium phase transition. When the natural frequencies of he oscillators are too diverse compared to the strength of the coupling, they are unable to synchronize and the system behaves incoherently. However, if the coupling is strong enough, all oscillators freeze into synchrony. The transition from one regime to the other takes place at a certain threshold At this point some elements lock their relative phase and a cluster of synchronized nodes develops. This constitutes the onset of synchronization. beyond this value, the population of oscillators is split into a partially synchronized state made up of oscillators locked in phase and a group of nodes whose atural frequencies are too different as to be part of the coherent cluster. Finally, after further increasing the coupling, more and more elements get entrained around the mean phase of the collective rhythm generated by the whole population and the system settles in the completely synchronized state. Kuramoto[29, 30] worked out a mathematically tractable model to describe this phenomenology. He recognized that the most suitable case for analytical treatment should be the mean field approach. He proposed an all-to-all purely sinusoidal coupling, and then the governing equations for each of the oscillators in the system are 选=a+∑sm(-B)(=1,…,N where the factor 1/N is incorporated to ensure a good behavior of the model in the thermodynamic limit, N-00, @i stands for the natural frequency of oscillator i, and K is the coupling constant. The frequencies i are distributed according to some function g(o), that is usually assumed to be unimodal and symmetric about its mean frequency 2. Admittedly, due to the rotational symmetry in the model, we can use a rotating frame and redefine O-O+52 for all i and set S2=0, so that the oi's denote deviations from the mean frequency The collective dynamics of the whole population is measured by the macroscopic complex order parameter, r()el9()-1 where the modulus<r(t)< 1 measures the phase coherence of the population and o (t) is the average phase. The values phase locked or move incoherently, respectively. Multiplying both parts of Eq (4)by n which all oscillators are either r y 1 and r z 0(where a stands for fluctuations of size O(N-1))describe the limits ci and equating imaginary parts gives rsin(φ一6 yielding 6=0+ Kr sin(中一)(i=1 Eq (6)states that each oscillator interacts with all the others only through the mean field quantities r and the first quantity provides a positive feedback loop to the system s collective rhythm: as r increases because the population becomes more coherent, the coupling between the oscillators is further strengthened and more of them can be recruited to take part in the coherent pack. Moreover, Eq (6)allows to calculate the critical coupling K and to characterize the order parameter limt-oo rr ( K)= r(k). Looking for steady solutions, one assumes that r(t)and (t)are constant. Next, without loss of generality, we can set =0, which leads to the equations of motion[29, 30 ei= o- kr sin(ei (i= N) The solutions of Eq (7)reveal two different types of long-term behavior when the coupling is larger than the critical value, K c. On the one hand, a group of oscillators for which oil s Kr are phase-locked at frequency s2 in the original frame according to
A. Arenas et al. / Physics Reports 469 (2008) 93–153 97 obtained so far on three different kinds of oscillatory ensembles: limit cycle oscillators (Kuramoto), pulse-coupled models, and finally coupled map systems. We reserve for Section 4 those works that use the MSF formalism. Many other works whose major contribution is the understanding of synchronization phenomena in specific scenarios are discussed in the Applications section. 3.1. Phase oscillators 3.1.1. The Kuramoto model The pioneering work by Winfree [6] spurred the field of collective synchronization and called for mathematical approaches to tackle the problem. One of these approaches, as already stated, considers a system made up of a huge population of weakly-coupled, nearly identical, interacting limit-cycle oscillators, where each oscillator exerts a phase dependent influence on the others and changes its rhythm according to a sensitivity function [27,28]. Even if these simplifications seem to be very crude, the phenomenology of the problem can be captured. Namely, the population of oscillators exhibits the dynamic analog to an equilibrium phase transition. When the natural frequencies of the oscillators are too diverse compared to the strength of the coupling, they are unable to synchronize and the system behaves incoherently. However, if the coupling is strong enough, all oscillators freeze into synchrony. The transition from one regime to the other takes place at a certain threshold. At this point some elements lock their relative phase and a cluster of synchronized nodes develops. This constitutes the onset of synchronization. Beyond this value, the population of oscillators is split into a partially synchronized state made up of oscillators locked in phase and a group of nodes whose natural frequencies are too different as to be part of the coherent cluster. Finally, after further increasing the coupling, more and more elements get entrained around the mean phase of the collective rhythm generated by the whole population and the system settles in the completely synchronized state. Kuramoto [29,30] worked out a mathematically tractable model to describe this phenomenology. He recognized that the most suitable case for analytical treatment should be the mean field approach. He proposed an all-to-all purely sinusoidal coupling, and then the governing equations for each of the oscillators in the system are: θ˙ i = ωi + K N XN j=1 sin(θj − θi) (i = 1, . . . , N), (3) where the factor 1/N is incorporated to ensure a good behavior of the model in the thermodynamic limit, N → ∞, ωi stands for the natural frequency of oscillator i, and K is the coupling constant. The frequencies ωi are distributed according to some function g(ω), that is usually assumed to be unimodal and symmetric about its mean frequency Ω. Admittedly, due to the rotational symmetry in the model, we can use a rotating frame and redefine ωi → ωi + Ω for all i and set Ω = 0, so that the ωi ’s denote deviations from the mean frequency. The collective dynamics of the whole population is measured by the macroscopic complex order parameter, r(t)e iφ(t) = 1 N XN j=1 e iθj (t) , (4) where the modulus 0 ≤ r(t) ≤ 1 measures the phase coherence of the population and φ(t) is the average phase. The values r ' 1 and r ' 0 (where ' stands for fluctuations of size O(N −1/2 )) describe the limits in which all oscillators are either phase locked or move incoherently, respectively. Multiplying both parts of Eq. (4) by e−iθi and equating imaginary parts gives r sin(φ − θi) = 1 N XN j=1 sin(θj − θi), (5) yielding θ˙ i = ωi + Kr sin(φ − θi) (i = 1, . . . , N). (6) Eq. (6) states that each oscillator interacts with all the others only through the mean field quantities r and φ. The first quantity provides a positive feedback loop to the system’s collective rhythm: as r increases because the population becomes more coherent, the coupling between the oscillators is further strengthened and more of them can be recruited to take part in the coherent pack. Moreover, Eq. (6) allows to calculate the critical coupling Kc and to characterize the order parameter limt→∞ rt(K) = r(K). Looking for steady solutions, one assumes that r(t) and φ(t) are constant. Next, without loss of generality, we can set φ = 0, which leads to the equations of motion [29,30] θ˙ i = ωi − Kr sin(θi) (i = 1, . . . , N). (7) The solutions of Eq.(7) reveal two different types of long-term behavior when the coupling is larger than the critical value, Kc . On the one hand, a group of oscillators for which |ωi | ≤ Kr are phase-locked at frequency Ω in the original frame according to
A Arenas et aL / Physics Reports 469(2008)93-153 the equation i= Kr sin(0i). On the other hand, the rest of the oscillators for which oil Kr holds, are drifting around the circle, sometimes accelerating and sometimes rotating at lower frequencies. Demanding some conditions for the stationary distribution of drifting oscillators with frequency ai and phases ei [27], a self-consistent equation for r can be derived as r=Kr/.(cose)g(@o)de where o= Kr sin (0). This equation admits a non-trivial solution, beyond which r >0. Eq. 8)is the Kuramoto mean field expression for the critical coupling at the onset of synchronization. Moreover, near the onset, the order parameter, r, obeys the usual square-root scaling law for mean field models, namely. (K-K) with B= 1/2. Numerical simulations of the model verified these results. The Kuramoto model(KM, from now on)approach to synchronization was a breakthrough for the understanding of synchronization in large populations of oscillators. Even in the simplest case of a mean field interaction, there are still unsolved problems that have resisted any analytical attempt. This is the case, e g. for finite populations of oscillators and some questions regarding global stability results 28 In what follows, we focus on another aspect of the models assumptions, namely that of the connection topology of real systems[14, 15], which usually do not show the all-to-all pattern of interconnections underneath the mean field approach. 3.1.2. Kuramoto model on complex networks To deal with the KM on complex topologies, it is necessary to reformulate eq ( 3)to include the connectivity =+∑qsi-)(=1,…N here oi is the coupling strength between pairs of connected oscillators and ai are the elements of the connectivity matrix. The original Kuramoto model is recovered by letting ai= 1, vi #j(all-to-all)and o k/N, vi, j The first problem when defining the KM in complex networks is how to state the interaction dynamics properly. In contrast with the mean field model, there are several ways to define how the connection topology enters in the governing uations of the dynamics. a good theory for Kuramoto oscillators in complex networks should be phenomenologically evant and provide formulas amenable to rigorous mathematical treatment. Therefore, such a theory should at least at the same time, should remain calculable in the thermodynamic limi k independently of the interaction dynamics, and For the original model, Eq (3), the coupling term on the right hand side of Eq (10 )is an intensive magnitude because the dependence on the size of the system cancels out. This independence on the number of oscillators n is achieved by choosing Ji= K/N. This prescription turns out to be essential for the analysis of the system in the thermodynamic limit N-0o in to become dependent on N. Therefore, in the thermodynamic limit, the coupling term tends to zero except for those nodes with a degree that scales with N. Note that the existence of such nodes is only possible in networks with power-law degree distributions[14,15, but this happens with a very small probability as k- ', with y > 2. In these cases, mean field solutions ndependent of N are recovered, with slight differences in the onset of synchronization of all-to-all and SF networks[31]. A second prescription consists in taking oi k/ki(where ki is the degree of node i)so that oi is a weighted interaction factor that also makes the right hand side of Eq. (10)intensive. This form has been used to solve the paradox o f heterogeneity [32 that states that the heterogeneity in the degree distribution, which often reduces the average distance between nodes, may suppress synchronization in networks of oscillators coupled symmetrically with uniform coupling strength. This result refers to the stability of the fully synchronized state, but not to the dependence of the order parameter on the coupling strength(where partially synchronized and unsynchronized states exist). Besides, the inclusion of weights in the interaction strongly affects the original KM dynamics in complex networks because it can impose a dynamic homogeneity that masks the real topological heterogeneity of the network. The prescription j= K/const, which may seem more appropriate, also causes some conceptual problems because the sum in the right hand side of Eq (10)could eventually diverge in the thermodynamic limit. The constant in the denominator could in principle be any quantity related to the topology such as the average connectivity of the graph. (k), or the maximum degree kmax. Its physical meaning is a re-scaling of the temporal scales involved in the dynamics. However, except for the case of oi =K/kmax, the other possible settings do not avoid the problems when N-o0. On the other hand, for a proper comparison of the results obtained for different complex topologies(e g SF or uniformly random), the global and local measures of coherence should be represented according to their respective time scales. Therefore, given two complex networks A and B with kmax = kA and kmax kB respectively, it follows that to make meaningful comparisons between
98 A. Arenas et al. / Physics Reports 469 (2008) 93–153 the equation ωi = Kr sin(θi). On the other hand, the rest of the oscillators for which |ωi | > Kr holds, are drifting around the circle, sometimes accelerating and sometimes rotating at lower frequencies. Demanding some conditions for the stationary distribution of drifting oscillators with frequency ωi and phases θi [27], a self-consistent equation for r can be derived as r = Kr Z π 2 − π 2 cos2 θ g(ω)dθ , where ω = Kr sin(θ ). This equation admits a non-trivial solution, Kc = 2 πg(0) (8) beyond which r > 0. Eq. (8) is the Kuramoto mean field expression for the critical coupling at the onset of synchronization. Moreover, near the onset, the order parameter, r, obeys the usual square-root scaling law for mean field models, namely, r ∼ (K − Kc ) β (9) with β = 1/2. Numerical simulations of the model verified these results. The Kuramoto model (KM, from now on) approach to synchronization was a breakthrough for the understanding of synchronization in large populations of oscillators. Even in the simplest case of a mean field interaction, there are still unsolved problems that have resisted any analytical attempt. This is the case, e.g. for finite populations of oscillators and some questions regarding global stability results [28]. In what follows, we focus on another aspect of the model’s assumptions, namely that of the connection topology of real systems [14,15], which usually do not show the all-to-all pattern of interconnections underneath the mean field approach. 3.1.2. Kuramoto model on complex networks To deal with the KM on complex topologies, it is necessary to reformulate Eq. (3) to include the connectivity θ˙ i = ωi + X j σijaij sin(θj − θi) (i = 1, . . . , N), (10) where σij is the coupling strength between pairs of connected oscillators and aij are the elements of the connectivity matrix. The original Kuramoto model is recovered by letting aij = 1, ∀i 6= j (all-to-all) and σij = K/N, ∀i, j. The first problem when defining the KM in complex networks is how to state the interaction dynamics properly. In contrast with the mean field model, there are several ways to define how the connection topology enters in the governing equations of the dynamics. A good theory for Kuramoto oscillators in complex networks should be phenomenologically relevant and provide formulas amenable to rigorous mathematical treatment. Therefore, such a theory should at least preserve the essential fact of treating the heterogeneity of the network independently of the interaction dynamics, and at the same time, should remain calculable in the thermodynamic limit. For the original model, Eq. (3), the coupling term on the right hand side of Eq. (10) is an intensive magnitude because the dependence on the size of the system cancels out. This independence on the number of oscillators N is achieved by choosing σij = K/N. This prescription turns out to be essential for the analysis of the system in the thermodynamic limit N → ∞ in the all-to-all case. However, choosing σij = K/N for the governing equations of the KM in a complex network makes them to become dependent on N. Therefore, in the thermodynamic limit, the coupling term tends to zero except for those nodes with a degree that scales with N. Note that the existence of such nodes is only possible in networks with power-law degree distributions [14,15], but this happens with a very small probability as k −γ , with γ > 2. In these cases, mean field solutions independent of N are recovered, with slight differences in the onset of synchronization of all-to-all and SF networks [31]. A second prescription consists in taking σij = K/ki (where ki is the degree of node i) so that σij is a weighted interaction factor that also makes the right hand side of Eq. (10) intensive. This form has been used to solve the paradox of heterogeneity [32] that states that the heterogeneity in the degree distribution, which often reduces the average distance between nodes, may suppress synchronization in networks of oscillators coupled symmetrically with uniform coupling strength. This result refers to the stability of the fully synchronized state, but not to the dependence of the order parameter on the coupling strength (where partially synchronized and unsynchronized states exist). Besides, the inclusion of weights in the interaction strongly affects the original KM dynamics in complex networks because it can impose a dynamic homogeneity that masks the real topological heterogeneity of the network. The prescription σij = K/const, which may seem more appropriate, also causes some conceptual problems because the sum in the right hand side of Eq. (10) could eventually diverge in the thermodynamic limit. The constant in the denominator could in principle be any quantity related to the topology, such as the average connectivity of the graph,hki, or the maximum degree kmax. Its physical meaning is a re-scaling of the temporal scales involved in the dynamics. However, except for the case of σij = K/kmax, the other possible settings do not avoid the problems when N → ∞. On the other hand, for a proper comparison of the results obtained for different complex topologies (e.g. SF or uniformly random), the global and local measures of coherence should be represented according to their respective time scales. Therefore, given two complex networks A and B with kmax = kA and kmax = kB respectively, it follows that to make meaningful comparisons between
A Arenas et aL/Physics Reports 469(2008)93-153 8 0.2 口N=104 ◇N=5x10 Fig. 2. Order parameter r(Eq(4))as a function of a for several BA networks of different sizes. Finite size scaling analysis shows that the onset of synchronization takes place at a critical value a =0.05(1). The inset is a zoom around ae From 35] observables, the equations of motion(Eq (10))should refer to the same time scales, i.e. d= ka/ka ke/kB= o. with this formulation in mind, Eq (10)reduces to B=m+∑q(-B)(=,……N however, that the same order parameter, Eq (4), is often used to describe the coherence of the synchronized state ZAaB ndently of the specific topology of the network. this allows us to study the dynamics of Eq (11on different topologies pare the results, and properly inspect the interplay between topology and dynamics in what concerns synchronization As we shall see, there are also several ways to define the order parameter that characterizes the global dynamics of the system, some of which were introduced to allow for analytical treatments at the onset of synchronization. We advance. 3.1.3. Onset of synchronization in complex networks Studies on synchronization in complex topologies where each node is considered to be a Kuramoto oscillator, were first reported for WS networks [33, 34 and Ba graphs 35, 36]. These works are mainly numerical explorations of the onset of synchronization, their main goal being the characterization of the critical coupling beyond which groups of nodes beating coherently first appear. In [34, the authors considered oscillators with intrinsic frequencies distributed according to a aussian distribution with unit variance arranged in a wS network with varying rewiring probability, p, and explored how the order parameter, Eq. (4), changes upon addition of long-range links. Moreover, they assumed a normalized coupling der variatio.K/(k), where (k)is the average degree of the graph. Numerical integration of the equations of motion(10) The results confirm that networks obtained from a regular ring by just rewiring a tiny fraction of links(p> 0)can be synchronized with a finite K. Moreover, in contrast with the arguments provided in [34 we notice that their results had been obtained for a fixed average degree and thus the Kuramoto's critical coupling cannot be recovered by simply taking p-1, which produces a random ER graph with a fixed minimum connectivity. This limit is recovered by letting(k)increase. Actually, numerical simulations of the same model in 33]showed that the Kuramoto limit is approached when the average connectivity grows. In[35] the same problem in BA networks is considered. The natural frequencies and the initial values of e; were randomly drawn from a uniform distribution in the interval (-1/2, 1/2)and(-, T), respectively. The global dynamics of the system, Eq (11). turns out to be qualitatively the same as for the original KM as shown in Fig. 2, where the dependence of the order parameter Eq (4)with o is shown for several system sizes. The existence of a critical point for the KM on SF networks came as a surprise. Admittedly, this is one of the few cases in which a dynamical process shows a critical behavior when the substrate is described by a power-law connectivity distribution with an exponent y 0, the order parameter reaches a stationary value as N oo(though still with O(1//N) fluctuations). Therefore, plots of r versus N allow us to
A. Arenas et al. / Physics Reports 469 (2008) 93–153 99 Fig. 2. Order parameter r (Eq. (4)) as a function of σ for several BA networks of different sizes. Finite size scaling analysis shows that the onset of synchronization takes place at a critical value σc = 0.05(1). The inset is a zoom around σc . From [35]. observables, the equations of motion (Eq. (10)) should refer to the same time scales, i.e. σij = KA/kA = KB/kB = σ. With this formulation in mind, Eq. (10) reduces to θ˙ i = ωi + σ X j aij sin(θj − θi) (i = 1, . . . , N), (11) independently of the specific topology of the network. This allows us to study the dynamics of Eq.(11) on different topologies, compare the results, and properly inspect the interplay between topology and dynamics in what concerns synchronization. As we shall see, there are also several ways to define the order parameter that characterizes the global dynamics of the system, some of which were introduced to allow for analytical treatments at the onset of synchronization. We advance, however, that the same order parameter, Eq. (4), is often used to describe the coherence of the synchronized state. 3.1.3. Onset of synchronization in complex networks Studies on synchronization in complex topologies where each node is considered to be a Kuramoto oscillator, were first reported for WS networks [33,34] and BA graphs [35,36]. These works are mainly numerical explorations of the onset of synchronization, their main goal being the characterization of the critical coupling beyond which groups of nodes beating coherently first appear. In [34], the authors considered oscillators with intrinsic frequencies distributed according to a Gaussian distribution with unit variance arranged in a WS network with varying rewiring probability, p, and explored how the order parameter, Eq. (4), changes upon addition of long-range links. Moreover, they assumed a normalized coupling strength σij = K/hki, where hki is the average degree of the graph. Numerical integration of the equations of motion (10) under variation of p shows that collective synchronization emerges even for very small values of the rewiring probability. The results confirm that networks obtained from a regular ring by just rewiring a tiny fraction of links (p & 0) can be synchronized with a finite K. Moreover, in contrast with the arguments provided in [34], we notice that their results had been obtained for a fixed average degree and thus the Kuramoto’s critical coupling cannot be recovered by simply taking p → 1, which produces a random ER graph with a fixed minimum connectivity. This limit is recovered by letting hkiincrease. Actually, numerical simulations of the same model in [33] showed that the Kuramoto limit is approached when the average connectivity grows. In [35] the same problem in BA networks is considered. The natural frequencies and the initial values of θi were randomly drawn from a uniform distribution in the interval(−1/2, 1/2) and (−π , π ), respectively. The global dynamics of the system, Eq. (11), turns out to be qualitatively the same as for the original KM as shown in Fig. 2, where the dependence of the order parameter Eq. (4) with σ is shown for several system sizes. The existence of a critical point for the KM on SF networks came as a surprise. Admittedly, this is one of the few cases in which a dynamical process shows a critical behavior when the substrate is described by a power-law connectivity distribution with an exponent γ ≤ 3 [14,15,37]. In principle it could be a finite size effect, but it turned out from numerical simulations that this was not the case. To determine the exact value of σc , one can make use of standard finite-size scaling analysis. At least two complementary strategies have been reported. The first one allows bounding the critical point and is computationally more expensive. Consider a network of size N, for which no synchronization is attained below σc , where r(t) decays to a small residual value of size O(1/ √ N). Then, the critical point may be found by examining the N-dependence of r(σ , N). In the sub-critical regime (σ σc , the order parameter reaches a stationary value as N → ∞ (though still with O(1/ √ N) fluctuations). Therefore, plots of r versus N allow us to
A Arenas et aL / Physics Reports 469(2008)93-153 locate the critical point o. Alternatively, a more accurate approach can be adopted Assume the scaling form for the order parameter[38: r=N (N(o-ac) where f(x) is a universal scaling function bounded as x ->+oo and a and v are two critical exponents to be determined. ofo for g s oc, the value of the functionf is independent of N, the estimation of oc can be done by plotting@ras a function a by-product, the method also allows us to calculate the two scaling exponents a and v, the latter can be obtained from the In[(dr/do )lo ]=(v-a)In n+const, once a is computed Following these scaling procedures, it was estimated a value for the critical coupling strength o 0.05(1)[35, 39, 40] Moreover, r M(o-o) when approaching the critical point from above with B= 0.46(2)indicating that the square-root behavior typical of the mean field version of the model(B= 1/2)seems to hold as well for BA networks. Before turning our attention to some theoretical attempts to tackle the onset of synchronization, it is worth briefly summarizing other numerical results that have explored how the critical coupling depends on other topological features of the underlying SF graph. Recent results have shed light on the influence of the topology of the local interactions on the route to, and the onset of, synchronization. In particular, the authors in[41-43 explored the Kuramoto dynamics on networks in which the degree distribution is kept fixed, while the clustering coefficient(C)and the average path length (e)of the graph change. The results suggest that the onset of synchronization is mainly determined by C, namely, networks with a high clustering coefficient promote synchronization at lower values of the coupling strength On the other hand, when the coupling is increased beyond the critical point, the effect of e dominates over C and the phase diagram is smoothed out(a sort of stretching), delaying the appearance of the fully synchronized state as the average shortest path length In a series of recent papers [31, 44-48 the onset of synchronization in large networks of coupled oscillators has been analyzed from a theoretical point of view under certain(sometimes strong) assumptions. Despite these efforts no exact analytical results for the KM on general complex networks are available up to now. Moreover, the analytical approaches ponent y s 3, the critical coupling va anishes as n-o. in contrast to numerical simulations on BA model networks. It appears that the strong heterogeneity of real networks and the finite average connectivity strongly hampers analytical solutions of the model Following [31]. consider the system in Eq (11), with a symmetric adjacency matrix i ji. Defining a local order parameter r as re=∑ae) (14) re(.)t stands for a time average, a new global order parameter to measure the macroscopic coherence is readily (15) Now, rewriting Eq (11)as a function of ri, yields, 6=0一 or sin(一中)一ah(t) (16) In Eq (16);(t)= Im(e- L, ai ((ei)t-e)) depends on time and contains time fluctuations. Assuming the terms in the previous sum to be statistically independent hi (t)is expected to be proportional to ki above the transition, where T ki. Therefore, except very close to the critical point, and assuming that the number of connections of each node is large enough(ki > 1 as to be able to neglect the time fluctuations entering h, i. e, hi <ri). the equation describing the dynamics of node i can be reduced to (17) I The reader can find the extension of the forthcoming formalism to directed networks in(44). This obviously restricts the range of real networks to which the approximation can be applied
100 A. Arenas et al. / Physics Reports 469 (2008) 93–153 locate the critical point σc . Alternatively, a more accurate approach can be adopted. Assume the scaling form for the order parameter [38]: r = N −α f(N ν (σ − σc )), (12) where f(x) is a universal scaling function bounded as x → ±∞ and α and ν are two critical exponents to be determined. Since at σ = σc , the value of the function f is independent of N, the estimation of σc can be done by plotting N α r as a function of σ for various sizes N and then finding the value of α that gives a well-defined crossing point, the critical coupling σc . As a by-product, the method also allows us to calculate the two scaling exponents α and ν, the latter can be obtained from the equality ln[(dr/dσ )|σc ] = (ν − α)ln N + const, (13) once α is computed. Following these scaling procedures, it was estimated a value for the critical coupling strength σc = 0.05(1) [35,39,40]. Moreover, r ∼ (σ − σc ) β when approaching the critical point from above with β = 0.46(2) indicating that the square-root behavior typical of the mean field version of the model (β = 1/2) seems to hold as well for BA networks. Before turning our attention to some theoretical attempts to tackle the onset of synchronization, it is worth briefly summarizing other numerical results that have explored how the critical coupling depends on other topological features of the underlying SF graph. Recent results have shed light on the influence of the topology of the local interactions on the route to, and the onset of, synchronization. In particular, the authors in [41–43] explored the Kuramoto dynamics on networks in which the degree distribution is kept fixed, while the clustering coefficient (C) and the average path length (`) of the graph change. The results suggest that the onset of synchronization is mainly determined by C, namely, networks with a high clustering coefficient promote synchronization at lower values of the coupling strength. On the other hand, when the coupling is increased beyond the critical point, the effect of ` dominates over C and the phase diagram is smoothed out (a sort of stretching), delaying the appearance of the fully synchronized state as the average shortest path length increases. In a series of recent papers [31,44–48], the onset of synchronization in large networks of coupled oscillators has been analyzed from a theoretical point of view under certain (sometimes strong) assumptions. Despite these efforts no exact analytical results for the KM on general complex networks are available up to now. Moreover, the analytical approaches predict that for uncorrelated SF networks with an exponent γ ≤ 3, the critical coupling vanishes as N → ∞, in contrast to numerical simulations on BA model networks. It appears that the strong heterogeneity of real networks and the finite average connectivity strongly hampers analytical solutions of the model. Following [31], consider the system in Eq. (11), with a symmetric1 adjacency matrix aij = aji. Defining a local order parameter ri as rie iφi = XN j=1 aijhe iθjit, (14) where h· · ·it stands for a time average, a new global order parameter to measure the macroscopic coherence is readily introduced as r = PN i=1 ri PN i=1 ki . (15) Now, rewriting Eq. (11) as a function of ri , yields, θ˙ i = ωi − σri sin(θi − φi) − σ hi(t). (16) In Eq. (16), hi(t) = Im{e −iθi PN j=1 aij(he iθjit − e iθj)} depends on time and contains time fluctuations. Assuming the terms in the previous sum to be statistically independent, hi(t) is expected to be proportional to √ ki above the transition, where ri ∼ ki . Therefore, except very close to the critical point, and assuming that the number of connections of each node is large enough2 (ki 1 as to be able to neglect the time fluctuations entering hi , i.e., hi ri), the equation describing the dynamics of node i can be reduced to θ˙ i = ωi − σri sin(θi − φi). (17) 1 The reader can find the extension of the forthcoming formalism to directed networks in [44]. 2 This obviously restricts the range of real networks to which the approximation can be applied
A. Arenas et aL Physics Reports 469(2008)93-153 Next, we look for stationary solutions of Eq (17), i.e. sin (0i -i)=o/or In particular, oscillators whose intrinsic frequency satisfies oil s or become locked. Then, as in the Kuramoto mean field model, there are two contributions(though in this case to the local order parameter), one from locked and the other from drifting oscillators such that (9-中) To move one step further, some assumptions are needed. Consider a graph such that the average degree of nearest neighbors is high(i. e if the neighbors of node i are well-connected ) Then it is reasonable to assume that these nodes are not affected by the intrinsic frequency of i. This is equivalent to assuming solutions (ri, i) that are, in a statistical sense, independent of the natural frequency oi. With this assumption, the second summand in Eq (18)can be neglected. Taking into account that the distribution g(o) is symmetric and centered at 32=0. after some algebra one is left with [31 n=∑0(-1-( 19) The critical coupling o is given by the solution of Eq (19)that yields the smallest o. It can be argued that it is obtained when cos(j -i)=1 in Eq (19), thus which is the main equation of the time average approximation(recall that time fluctuations have been neglected). Note however, that to obtain the critical coupling, one has to know the adjacency matrix as well as the particular values of o; for all i and then solve eq. (20)numerically for the ril. Finally the global order parameter defined in Eq (15)can be computed Even if the underlying graph satisfies the other aforementioned topological constraints, it seems unrealistic to require knowledge of the oils. A further approach, referred to as the frequency distribution approximation can be adopted. According to the assumption that ki >>1 for all i, or equivalently that the number of connections per node is large(a dense graph), one can also consider that the natural frequencies of the neighbors of node i follows the distribution g(o) Then, Eq (20)can be rewritten avoiding the dependence on the particular realization of oil to yield, n-2aC-(m)-业 (21) with x=o/(ori). This equation allows us to readily determine the order parameter r as a function of the network topology ai), the frequency distribution (g(o))and the control parameter(o ) On the other hand, Eq- (20)still does not provide approximating(xar])Ng(O)which is valid for small, but nonzero, values of r. Namely, when x uce explicit expressions for the order parameter and the critical coupling strength. To this end, one introduces a first-order r where K= 2/(g(O)is Kuramoto's critical coupling. Moreover, as the smallest value of o corresponds to oc, it follows that the critical coupling is related to both Ke and the largest eigenvalue amax of the adjacency matrix, yielding Eq (22)states that in complex networks, synchronization is first attained at a value of the coupling strength that inversely depends ong(O)and on the largest eigenvalue Amax of the adjacency matrix. Note that this equation also recovers Kuramoto's result when ai = 1, vi j. since Amax = N-1. It is worth stressing that although this method allows us to calculate oc analytically, it fails to explain why for uncorrelated random SF networks with y 3 and in the thermodynamic limit N-oo, the critical value remains finite. This disagreement comes from the fact that in these SF networks, Amax is proportional to the cutoff of the degree distribution, kmax which in turn scales with the system size. Putting the two dependencies together, one obtainsλmax~kax~Nx→∞asN→∞, thus predicting o=0 in the thermodynamic limit, in contrast to finite size scaling analysis for the critical coupling via numerical solution of the equations of motion. Note, however, that the difference may be due to the use of distinct order parameters. Moreover, even in the case of SF networks with y>3, max still diverges when we take the thermodynamic limit, so that o -0 as well. As we shall see soon, this is not the case when other approaches are adopted, at least for y>3
A. Arenas et al. / Physics Reports 469 (2008) 93–153 101 Next, we look for stationary solutions of Eq. (17), i.e. sin(θi−φi) = ωi/σri . In particular, oscillators whose intrinsic frequency satisfies |ωi | ≤ σri become locked. Then, as in the Kuramoto mean field model, there are two contributions (though in this case to the local order parameter), one from locked and the other from drifting oscillators such that ri = XN j=1 aijhe i(θj−φi ) it = X |ωj |≤σrj aije i(θj−φi ) + X |ωj |>σrj aijhe i(θj−φi ) it . (18) To move one step further, some assumptions are needed. Consider a graph such that the average degree of nearest neighbors is high (i.e. if the neighbors of node i are well-connected). Then it is reasonable to assume that these nodes are not affected by the intrinsic frequency of i. This is equivalent to assuming solutions (ri, φi) that are, in a statistical sense, independent of the natural frequency ωi . With this assumption, the second summand in Eq. (18) can be neglected. Taking into account that the distribution g(ω) is symmetric and centered at Ω = 0, after some algebra one is left with [31] ri = X |ωj |≤σrj aij cos(φj − φi) s 1 − ωj σrj 2 . (19) The critical coupling σc is given by the solution of Eq. (19) that yields the smallest σ. It can be argued that it is obtained when cos(φj − φi) = 1 in Eq. (19), thus ri = X |ωj |≤σrj aijs 1 − ωj σrj 2 , (20) which is the main equation of the time average approximation (recall that time fluctuations have been neglected). Note, however, that to obtain the critical coupling, one has to know the adjacency matrix as well as the particular values of ωi for all i and then solve Eq. (20) numerically for the {ri}. Finally, the global order parameter defined in Eq. (15) can be computed from ri . Even if the underlying graph satisfies the other aforementioned topological constraints, it seems unrealistic to require knowledge of the {ωi}’s. A further approach, referred to as the frequency distribution approximation can be adopted. According to the assumption that ki 1 for all i, or equivalently, that the number of connections per node is large (a dense graph), one can also consider that the natural frequencies of the neighbors of node i follows the distribution g(ω). Then, Eq. (20) can be rewritten avoiding the dependence on the particular realization of {ωi} to yield, ri = X j aij Z σrj −σrj g(ω)s 1 − ω σrj 2 dω = σ X j aijrj Z 1 −1 g(xσrj) p 1 − x 2dx, (21) with x = ω/(σrj). This equation allows us to readily determine the order parameter r as a function of the network topology (aij), the frequency distribution (g(ω)) and the control parameter (σ). On the other hand, Eq. (20) still does not provide explicit expressions for the order parameter and the critical coupling strength. To this end, one introduces a first-order approximation g(xσrj) ≈ g(0) which is valid for small, but nonzero, values of r. Namely, when rj → 0 + r 0 i = σ Kc X j aijr 0 j , where Kc = 2/(πg(0)) is Kuramoto’s critical coupling. Moreover, as the smallest value of σ corresponds to σc , it follows that the critical coupling is related to both Kc and the largest eigenvalue λmax of the adjacency matrix, yielding σc = Kc λmax . (22) Eq. (22) states that in complex networks, synchronization is first attained at a value of the coupling strength that inversely depends on g(0) and on the largest eigenvalue λmax of the adjacency matrix. Note that this equation also recovers Kuramoto’s result when aij = 1, ∀i 6= j, since λmax = N − 1. It is worth stressing that although this method allows us to calculate σc analytically, it fails to explain why for uncorrelated random SF networks with γ ≤ 3 and in the thermodynamic limit N → ∞, the critical value remains finite. This disagreement comes from the fact that in these SF networks, λmax is proportional to the cutoff of the degree distribution, kmax which in turn scales with the system size. Putting the two dependencies together, one obtains λmax ∼ k 1 2 max ∼ N 1 2(γ −1) → ∞as N → ∞, thus predicting σc = 0 in the thermodynamic limit, in contrast to finite size scaling analysis for the critical coupling via numerical solution of the equations of motion. Note, however, that the difference may be due to the use of distinct order parameters. Moreover, even in the case of SF networks with γ > 3, λmax still diverges when we take the thermodynamic limit, so that σc → 0 as well. As we shall see soon, this is not the case when other approaches are adopted, at least for γ > 3
A Arenas et aL / Physics Reports 469(2008)93-153 It is possible to go beyond the latter approximation and to determine the behavior of r near the critical point. In 31] a perturbative approach to higher orders of Eq. (21)is developed, which is valid for relatively homogeneous degree distributions(y>5). They showed that for(o/oc)-1-0+ k where n1=-Tg(O)Kc/16 and n (24) here u is the normalized eigenvector of the adjacency matrix corresponding to Amax and(u)=2i u4 The analytical insights discussed so far can also be reformulated in terms of a mean field approximation [31, 46-48] for complex networks. This approach(valid for large enough( k))considers that every oscillator is influenced by the local field created in its neighborhood, so that T is proportional to the degree of the nodes ki, i.e i ki. Assuming this is the case and introducing the order parameter r through k∠a(e吗) after summing over i and substituting ri rki in Eq (21)we obtain [31 k g(xr)√1-x2d The above relation, Eq (26), was independently derived in [46 who first studied analytically the problem of synchronization in complex networks, though using a different approach. Taking the continuum limit, Eq (26)becomes kP(k)dk=o/kP(k)dk/g(xo rk)v1-x2dx which for r-0+ verifies /k=o/koo一 og(O)T KP(k)dk, which leads to the condition for the onset of synchronization(r>0)as og(o)T/ P(dk> kP(k)dk, that is (0)(k2)-(k2 The mean field result, Eq (29), gives as a surprising result that the critical coupling o in complex networks is nothing else but the one corresponding to the all-to-all topology ke re-scaled by the ratio between the first two moments of the degree coupling strongly depends on the topology of the underlying graph. In particular, ac ->0 when the second momenta e distribution, regardless of the many differences between the patterns of interconnections. Precisely, it states that the critic distribution( k)diverges, which is the case for SF networks with y s 3. Note, that in contrast with the result in Eq (22). for y >3, the coupling strength does not vanish in the thermodynamic limit On the other hand, if the mean degree is kept fixed and the heterogeneity of the graph is increased by decreasing y, the onset of synchronization occurs at smaller values of oc. Interestingly enough, the dependence gathered in Eq (29)has the same functional form for the critical points of other dynamical processes such as percolation and epidemic spreading processes [14, 15, 37]. while this result is still under numerical scrutiny, it would imply that the critical properties of many dynamical processes on complex networks are essentially determined by the topology of the graph, no matter whether the dynamics is nonlinear or not. The corroboration of this last claim will be of extreme importance in physics, probably changing many preconceived ideas about the nature dynamical phenomena. The approach holds if the fourth of the degree distribution, ()=/i P(k)kdk remains finite when N00
102 A. Arenas et al. / Physics Reports 469 (2008) 93–153 It is possible to go beyond the latter approximation and to determine the behavior of r near the critical point. In [31] a perturbative approach to higher orders of Eq. (21) is developed, which is valid for relatively homogeneous degree distributions (γ > 5).3 They showed that for (σ /σc ) − 1 ∼ 0 + r 2 = η η1K 2 c σ σc − 1 σ σc 3 , (23) where η1 = −πg 00(0)Kc/16 and η = hui 2λ 2 max Nhki 2 hu 4 i , (24) where u is the normalized eigenvector of the adjacency matrix corresponding to λmax and hu 4 i = PN j u 4 j /N. The analytical insights discussed so far can also be reformulated in terms of a mean field approximation [31,46–48] for complex networks. This approach (valid for large enough hki) considers that every oscillator is influenced by the local field created in its neighborhood, so that ri is proportional to the degree of the nodes ki , i.e., ri ∼ ki . Assuming this is the case and introducing the order parameter r through r = ri ki = 1 ki XN j=1 aijhe iθjit , (25) after summing over i and substituting ri = rki in Eq. (21) we obtain [31] XN j kj = σ XN j k 2 j Z 1 −1 g(xσrkj) p 1 − x 2dx. (26) The above relation, Eq.(26), was independently derived in [46], who first studied analytically the problem of synchronization in complex networks, though using a different approach. Taking the continuum limit, Eq. (26) becomes Z kP(k)dk = σ Z k 2 P(k)dk Z 1 −1 g(xσrk) p 1 − x 2dx, (27) which for r → 0 + verifies Z kP(k)dk = σ Z k 2 P(k)dk Z 1 −1 g(0) p 1 − x 2dx = σg(0)π 2 Z k 2 P(k)dk, (28) which leads to the condition for the onset of synchronization (r > 0) as σg(0)π 2 Z k 2 P(k)dk > Z kP(k)dk, that is, σc = 2 πg(0) hki hk 2 i = Kc hki hk 2 i . (29) The mean field result, Eq. (29), gives as a surprising result that the critical coupling σc in complex networks is nothing else but the one corresponding to the all-to-all topology Kc re-scaled by the ratio between the first two moments of the degree distribution, regardless of the many differences between the patterns of interconnections. Precisely, it states that the critical coupling strongly depends on the topology of the underlying graph. In particular, σc → 0 when the second moment of the distribution hk 2 i diverges, which is the case for SF networks with γ ≤ 3. Note, that in contrast with the result in Eq. (22), for γ > 3, the coupling strength does not vanish in the thermodynamic limit. On the other hand, if the mean degree is kept fixed and the heterogeneity of the graph is increased by decreasing γ , the onset of synchronization occurs at smaller values of σc . Interestingly enough, the dependence gathered in Eq. (29) has the same functional form for the critical points of other dynamical processes such as percolation and epidemic spreading processes [14,15,37]. While this result is still under numerical scrutiny, it would imply that the critical properties of many dynamical processes on complex networks are essentially determined by the topology of the graph, no matter whether the dynamics is nonlinear or not. The corroboration of this last claim will be of extreme importance in physics, probably changing many preconceived ideas about the nature of dynamical phenomena. 3 The approach holds if the fourth moment of the degree distribution, hk 4 i = R ∞ 1 P(k)k 4dk remains finite when N → ∞