Remote synchronization reveals network symmetries and functional modules

We study a Kuramoto model in which the oscillators are associated with the nodes of a complex network and the interactions include a phase frustration, thus preventing full synchronization. The system organizes into a regime of remote synchronization where pairs of nodes with the same network symmetry are fully synchronized, despite their distance on the graph. We provide analytical arguments to explain this result and we show how the frustration parameter affects the distribution of phases. An application to brain networks suggests that anatomical symmetry plays a role in neural synchronization by determining correlated functional modules across distant locations.

We study a Kuramoto model in which the oscillators are associated with the nodes of a complex network and the interactions include a phase frustration, thus preventing full synchronization. The system organizes into a regime of remote synchronization where pairs of nodes with the same network symmetry are fully synchronized, despite their distance on the graph. We provide analytical arguments to explain this result and we show how the frustration parameter affects the distribution of phases. An application to brain networks suggests that anatomical symmetry plays a role in neural synchronization by determining correlated functional modules across distant locations. Synchronization of coupled dynamical units is a ubiquitous phenomenon in nature [1]. Remarkable examples include phase locking in laser arrays, rhythms of flashing fireflies, wave propagation in the heart, and also normal and abnormal correlations in the activity of different regions of the human brain [2][3][4][5]. In 1975 Y. Kuramoto proposed a simple microscopic model to study collective behaviors in large populations of interacting elements [6]. In its original formulation the Kuramoto model describes each unit of the system as an oscillator which continuously readjusts its frequency in order to minimize the difference between its phase and the phase of all the other oscillators. This model has shown very successful in understanding the spontaneous emergence of synchronization and, over the years, many variations have been considered [7][8][9]. Recently, the Kuramoto model has been also extended to sets of oscillators coupled through complex networks [2,10,11], and it has been found that the topology of the interaction network has a fundamental role in the emergence and stability of synchronized states [12,13]. In particular, the presence of communities -groups of tightly connected nodes-has a relevant impact on the path to synchronization [14][15][16][17][18], and units that are close to each other on the network, or belong to the same module or community [19], have a higher chance to exhibit similar dynamics. This implies that nodes in the same structural module share similar functions, which is a belief often supported by empirical findings [3,20]. However, various examples are found in nature where functional similarity is instead associated with morphological symmetry. In these cases, units with similar roles, which could potentially swap their position without altering the overall functioning of the system, appear in remote locations of the network. Some examples include cortical areas in brains [21], symmetric organs in plants and vertebrates [22,23], and even atoms in complex molecules [24]. Therefore, identifying the sets of symmetric units of a complex system might be helpful to understand its organization. Finding the global symmetries in a graph, i.e., constructing its automorphism group, is a classical problem in graph theory. However, it is still unknown if this problem is polynomial or NP-complete [25,26], even if there exist polynomial-time algorithms for graphs with bounded maximum degree [27]. Recent works have focused instead on defining and detecting local symmetries in complex networks [28,29]. Nevertheless, the interplay between the structural symmetries of a network and the dynamics of processes occurring over the network has been studied only marginally [30][31][32], or for specific small network motifs [33][34][35].
In this Letter we show that network symmetries play a central role in the synchronization of a system. We consider networks of identical Kuramoto oscillators, in which a phase frustration parameter forces connected nodes to maintain a finite phase difference, thus hindering the attainment of full synchronization. We prove that the configuration of phases at the synchronized state reflects the symmetries of the underlying coupling network. In particular, two nodes with the same symmetry have identical phases, i.e., are fully synchronized, despite the distance between the two nodes on the graph. Such a remote synchronization behavior is here induced by the network symmetries and not by an initial ad hoc choice of different natural frequencies [30].
Let us consider N identical oscillators associated to the nodes of a connected graph G(N , L), with N = |N | nodes and K = |L| links. Each node i is characterized, at time t, by a phase θ i (t) whose time evolution is governed by the equatioṅ Here ω is the natural frequency, identical for all the oscillators, and A ≡ {a ij } is the adjacency matrix of the arXiv:1211.5390v3 [nlin.AO] 27 Apr 2013 coupling graph. The model has two control parameters: λ > 0 accounting for the strength of the interaction, and α, the phase frustration parameter ranging in [0, π/2]. When α = 0, the model reduces to a network of identical Kuramoto oscillators. In this case, the fully synchronized state is globally stable for a set of initial conditions having finite measure [6,9] and the transient dynamics closely reflects the structure of the graph, so that nodes belonging to the same structural module evolve similarly in time [14]. However, the synchronized state can coexist with other nontrivial attractors, e.g., uniformly twisted waves, especially if the coupling topology is regular and sparse (see Ref. [12] for a discussion about the size of the sync basin). Instead, if the oscillators are not identical the frequency distribution tends to separate their phases and, as a result, there is a transition from an incoherent state (with order parameter r = 1 N N j=1 e iθj equal to 0) to a synchronized one (r = 0) at a critical value λ c of the coupling strength.
The introduction of a phase frustration α = 0 forces directly connected oscillators to maintain a constant phase difference [36]. In particular, we found that for a wide range of α > 0 the dynamics in Eqs. (1) reaches a stationary state in which the oscillators at two symmetric nodes have exactly the same phase, and this phase differs from the phases of nodes with different symmetries. Let us first illustrate this behavior and the effect of α on the three graphs G a , G b and G c shown in Fig. 1. In the three topmost panels of Fig. 2 we report the results of the numerical integration of Eqs. (1) on the graph G a for three different values of α. We find that, after a transient, the system settles into a stationary state in which,  at any time t, the phases are grouped into four different trajectories: θ 1 (t), θ 2 (t) = θ 3 (t), θ 4 (t) = θ 7 (t) and θ 5 (t) = θ 6 (t). In general, by increasing α up to a certain value α c we better separate the four trajectories.
The four clusters of nodes obtained for α < α c are identified by a color code in Fig. 1. We notice that each cluster groups together all the nodes with the same symmetry. In this way two distant nodes of the graph, e.g. node 4 and node 7, are fully synchronized even if the other nodes in the paths connecting them have different phases. In this respect, what we observe is a remote synchronization [30]. We have found similar results for the linear chain and for the Bethe lattice (see nodes with the same colors in G b and G c in Fig. 1).
Notice that if the system reaches a synchronized state and α is small enough, Eq. (1) can be linearized by replacing the sinus with its argument. We obtaiṅ where k i = j a ij is the degree of node i, L ij are the entries of the Laplacian matrix of the graph L ≡ D − A, and D is a diagonal matrix such that D ii = k i . Without loss of generality, we can set λ = 1, ω = 0. If the system is synchronized thenθ i = Ω, ∀i, so that the phases must satisfy the equations where k = N −1 i k i is the average degree of the network. This corresponds to a synchronization frequencẏ θ i = Ω = −α k ∀i. In a connected graph the Laplacian matrix has one null eigenvalue and the system of Eqs. (3) is singular. Consequently, at each time t we can solve the system by computing the phase difference between each node and a given node chosen as reference.
. This is in agreement with the results of the simulations: the phases are clustered into four groups, with nodes with the same symmetry having the same phase, and nodes with different symmetries being separated by a phase lag that depends on α as in the relations found above. An analogous analytical expression can be derived for a finite chain (graph G b in Fig.1), for which we Consequently, two nodes symmetrically placed with respect to node 0 will have identical phases.
We now provide a general argument to explain why the synchronization of Eqs. (1) is related to graph symmetries. A graph G(N , L) has a symmetry if and only if it is possible to find a bijection π : N → N which preserves the adjacency relation of G, i.e., which is an automorphism for G. Formally, this means that there exists a permutation matrix P = P (π) such that P AP −1 = A. If P corresponds to an automorphism of G then P commutes with A, i.e., P A = AP , and P AP −1 performs a relabeling of the nodes of the original graph which preserves the adjacency matrix [37]. In general a graph can admit more than one automorphism. For instance, graph G a in Fig.1 has at least three nontrivial bijections which preserve the adjacency matrix, namely π 1 : (1, 2, 3, 4, 5, 6, 7) → (1, 3, 2, 4, 5, 6, 7) π 2 : (1, 2, 3, 4, 5, 6, 7) → (1, 2, 3, 7, 6, 5, 4) π 3 : (1, 2, 3, 4, 5, 6, 7) → (1, 3, 2, 7, 6, 5, 4) Node 2 and node 3 are symmetric because we can relabel the nodes of G a (e.g. by means of either π 1 or π 3 ) so that node 2 is mapped into node 3 and vice versa, and the adjacency matrix of G a is left unchanged. Similarly, for the pairs {4, 7} and {5, 6}, there are two different relabelings which preserve adjacency relations, i.e., π 2 and π 3 . In terms of symmetries, the graph G has the following four different classes of nodes: . Now, if a permutation of the nodes is an automorphism of G, then P LP −1 = P DP −1 − P AP −1 = D − A = L, i.e., the associated permutation matrix P also commutes with the Laplacian matrix of the graph. By left-multiplying both sides of Eq. (3) by P , we get P Lθ = αP [ k 1 − k]. Since P L = LP (P commutes with L) and P k = k (symmetric nodes have the same degree) then we have Combining Eq. (3) and Eq. (4) we finally obtain the linear system which is singular, i.e., has one free variable. Again, it can be solved by leaving free one of the N variables θ i , setting φ j = θ j − θ i and considering the new systemLP φ =Lφ. The matrixP is obtained from P by removing the row and the column corresponding to node i. If P does not permute node i with another node, thenP is still a permutation matrix. Similarly,L is the reduced Laplacian, i.e., the matrix obtained from the Laplacian by deleting the i−th row and the i−th column. By left-multiplying byL −1 , which is not singular, we obtaiñ SinceP φ is a permutation of the phases of symmetric nodes, Eq.(6) implies that the phases of symmetric nodes will be equal at any time, whereas by solving Eq. (4) we can get the values of the corresponding phases. This argument is valid for small values of α, since the lineariza- but as shown in Fig. 2(a)-(c) we observe the formation of the same perfectly synchronized clusters of symmetric nodes for a wide range of α. However, when α becomes larger than a certain value α c , the assumptionθ i = Ω, ∀i does not hold any more and the global synchronized state loses stability. By looking at Fig. 2(d)-(e) we notice that for α > α c , with α c 1.05 for the graph G a , the value of r steadily decreases while the dispersion of phases increases, until it reaches the expected valueσ θ 1.39 for a system of seven incoherent oscillators (see Fig. 2(d) and Appendix). Moreover, for α > α c the maximal Lyapunov exponent of the system λ max becomes positive and the system enters a chaotic regime (see Fig. 2(e)). Interestingly, the results reported in Fig. 3(a)-(d) confirm that in this regime the coherence of symmetric nodes, measured by the pairwise order parameter r 2 , is higher than expected for incoherent oscillators (refer to Appendix for additional details). Figure 3(e) shows that for α > α c the system exhibits metastable, partially synchronized states, in which pairs of symmetric nodes alternates intervals of perfect synchronization with intervals of complete incoherence. We point out that in this regime chimera states could potentially occur [38][39][40][41] and could even coexist with remote synchronization for α < α c . Qualitatively similar results are obtained for different coupling topologies, but the actual value of α c seems to depend on the structure of the coupling network in a nontrivial way.
Application to the brain.-As an example, we investigate here the role of symmetry in the human brain by considering anatomical and functional brain connectivity graphs defined on the same set of N = 90 cortical areas (see details in Appendix). We have first constructed a graph of anatomical brain connectivity as obtained from DW-MRI data [42], where links represent axonal fibers, and we used this graph as a backbone network to integrate Eqs. (1). We identified candidate pairs of anatomically symmetric areas by means of agglomerative clustering, i.e., grouping together nodes having close phases at the stationary state (full dendrogram and details are provided in Appendix). Then, we considered the graph of functional brain connectivity, in which links represent statistically significant correlations between the BOLD fMRI time-series of cortical areas (see details in Appendix). Figure 4 illustrates the results for α = 0.5 (we obtained qualitatively similar results in a wide range of α). Consider nodes 57 and 74, corresponding respectively to the green and blue areas in panel (a). Not only the two areas are spatially separated, but there is no edge connecting the two corresponding nodes in the anatomical connectivity network. However, the two nodes are detected as a candidate symmetric pair since at the stationary state of the Kuramoto dynamics in Eq. (1) the oscillators associated to these two nodes have very close phases (see dendrogram in Appendix). As shown in Fig. 4(b), also the BOLD fMRI signals corresponding to nodes 57 and 74 also are strongly synchronized. We obtain remarkably different results when we consider node 74 and node 76. These nodes correspond to two spatially adjacent areas of the brain (the red and blue regions in Fig. 4(a) and are directly connected in the anatomical connectivity network. However, at the stationary state of Eq. (1) the phase difference of the oscillators associated to node 74 and 76 is quite large. Interestingly, in this case the fMRI time-series associated to these nodes are much less similar to each other [see the two bottom trajectories reported in Fig. 4(b)].
To quantify this effect, we plot in Fig. 4c the average functional correlation Z between the fMRI activity of pairs of brain areas as a function of the phase differences ∆θ between the phases of the corresponding oscillators, obtained from the dynamics of Eq. (1) on the anatomical connectivity network. The fact that Z decreases with ∆θ suggests that structural symmetry plays an important role in determining human brain functions. In fact, the functional activities of anatomically symmetric areas can be strongly correlated, even if the areas are distant in space. These results suggest that the study of anatomical symmetries in neural systems might provide meaningful insights about the functional organization of distant neural assemblies during diverse cognitive or pathological states [21]. Applied to other connectivity networks as a method to spot potential network symmetries, our study could provide new insights on the interplay be-tween structure and dynamics in complex systems.

APPENDIX
Expected order parameter for two incoherent oscillators. -In Fig.3 of the main text we showed that when α approaches π/2 the order parameter r 2 for a pair of symmetric oscillators remains higher than the expected valuer 2 for a pair of incoherent oscillators. We report here the derivation ofr 2 for a pair of oscillators that are not synchronized. We assume that the phases θ 1 and θ 2 of the two oscillators are drawn uniformly in [0, 2π]. Thanks to the rotational symmetry, we can set one of the phases equal to 0, e.g. θ 1 = 0, while drawing the other uniformly in [0, 2π]. We notice that the value of r 2 for the generic pair of phases (θ 1 = 0, θ 2 ) reads: Consequently, the expected value of r 2 is obtained as the average over all the possible choices of θ 2 , namely: Phase dispersion and numerical estimate for incoherent oscillators. -In Fig.2d of the main text we reported, as a function of α, the dispersion of the phases σ θ for a system of seven oscillators coupled through graph G a in Fig.1a. We noticed that the dispersion approaches the value 1.39 when α → π/2. Given a set of N unitary vectors having phases {θ 1 , θ 2 , . . . , θ N }, consider the average vector: having polar coordinates (r, ψ). The dispersion of the N phases of the set around the phase ψ of the average vector is defined as: where the difference (θ j − ψ) is computed mod(2π) and takes in [0, π]. An estimate of the expected phase dispersion for a set of N incoherent oscillators can be obtained as the average over M independent realizations of σ θ computed on a set of N phases uniformly drawn in [0, 2π]. For the system of N = 7 oscillators coupled through graph G a we averaged σ θ over M = 10 7 samples, obtaining the estimateσ θ = 1.394 ± 0.248.
Computation of the maximum Lyapunov exponent λ max . -The computation of the maximum Lyapunov exponent for the system coupled through graph G a of Fig. (1) in the main text was performed using equally separated values of α between 0 and 1.57 π/2 (∆α = 0.01). Then, for each value of α, we considered 500 different initial configurations of the phases of the seven oscillators. For each initial condition, we let evolve the trajectory θ(t) according to Eq. (1) in the main text, until it reached the attractor (or the stationary state, for α π/2). To properly take into account the rotational symmetry of the system, we studied the evolution of θ(t) in Cartesian coordinates, i.e., by looking at the set of 2N variables x = {cos θ 1 , cos θ 2 , . . . , cos θ N , sin θ 1 , sin θ 2 , . . . , sin θ N }. We considered a perturbation of x of magnitude ε = 10 −4 , i.e., a trajectoryx such that d(x(t),x(t)) = ε. Here d(·, ·) denotes the Euclidean distance in R 2N . Then, we integrated both trajectories for one integration step h (using a standard fourth-order Runge-Kutta integration scheme), we measured the distance d 1 = d(x(t+h),x(t+ h)) and we computed λ = log (d 1 /ε). The quantity λ is a one-step approximation of the largest Lyapunov exponent of the system [43]. Then, we realigned the perturbed trajectoryx so that the distance between x(t + h) and the realigned perturbed trajectoryx(t + h) was equal to ε in the same direction of d 1 , and we iterated the procedure. The value of λ max for a set of initial conditions was obtained by averaging the values of λ computed at each iteration over 10 4 subsequent integration steps.
Brain data acquisition and pre-processing. -The anatomical connectivity network is based on the connectivity matrix obtained by Diffusion Magnetic Resonance Imaging (DW-MRI) data from 20 healthy participants, as described in [44]. The elements of this matrix represent the probabilities of connection between the 90 anatomical regions of interest (N = 90 nodes in the network) of the Tzourio-Mazoyer brain atlas [45]. These probabilities are proportional to the density of fibers between different areas, so each element of the matrix represents an approximation of the connection strength between the corresponding pair of brain regions.
The functional brain connectivity was extracted from BOLD fMRI resting state recordings obtained as described in [46]. All fMRI data sets (segments of 5 minutes recorded from 15 healthy subjects) were co-registered to the anatomical data set and normalized to the standard MNI (Montreal Neurological Institute) template image, to allow comparisons between subjects. As for DW-MRI data, normalized and corrected functional scans were sub-sampled to the anatomical labeled template of the human brain [45]. Regional time series were estimated for each individual by averaging the fMRI time series over all voxels in each region (data were not spatially smoothed before regional parcellation). To eliminate low frequency noise (e.g. slow scanner drifts) and higher frequency artifacts from cardiac and respiratory oscillations, time-series were digitally filtered with a finite impulse response (FIR) filter with zero-phase distortion (bandwidth 0.01 − 0.1 Hz) as in [46].
Functional synchrony. -A functional link between two time series x i (t) and x j (t) (normalized to zero mean and unit variance) was defined by means of the linear cross-correlation coefficient computed as r ij = x i (t)x j (t) , where · denotes the temporal average. For the sake of simplicity, we only considered here correlations at lag zero. To determine the probability that correlation values are significantly higher than what is expected from independent time series, r ij (0) values (denoted r ij ) were firstly transformed by the Fisher's Z transform Z ij = 0.5 ln 1 + r ij 1 − r ij (S-3) Under the hypothesis of independence, Z ij has a normal distribution with expected value 0 and variance 1/(df − 3), where df is the effective number of degrees of freedom [47][48][49]. If the time series consist of independent measurements, df simply equals the sample size, N . Nevertheless, autocorrelated time series do not meet the assumption of independence required by the standard significance test, yielding a greater Type I error [47][48][49]. In presence of auto-correlated time series df must be corrected by the following approximation: where r xx (τ ) is the autocorrelation of signal x at lag τ .
To estimate a threshold for statistically significant correlations, a correction for multiple testing was used. The False Discovery Rate (FDR) method was applied to each matrix of Z ij values [50]. With this approach, the threshold of significance Z th was set such that the expected fraction of false positives is restricted to q ≤ 0.05.
Clustering of phase values. -To identify brain areas that could be related by a topological symmetry, we used the anatomical connectivity obtained from the DW-MRI data as the connectivity matrix (N = 90 nodes) in Eq. (1) of the main text. A standard hierarchical agglomerative clustering algorithm was then used to identify nodes with similar phases [51]. The resulting dendrogram is depicted in Fig. S-1.