Hub-driven remote synchronization in brain networks

The phenomenon of “remote synchronization” (RS), first observed in a star network of oscillators, involves synchronization of unconnected peripheral nodes through a hub that maintains independent dynamics. In the RS regime the central hub was thought to serve as a passive gate for information transfer between nodes. Here, we investigate the physical origin of this phenomenon. Surprisingly, we find that a hub node can drive remote synchronization of peripheral oscillators even in the presence of a repulsive mean field, thus actively governing network dynamics while remaining asynchronous. We study this novel phenomenon in complex networks endowed with multiple hub-nodes, a ubiquitous feature of many real-world systems, including brain connectivity networks. We show that a change in the natural frequency of a single hub can alone reshape synchronization patterns across the entire network, and switch from direct to remote synchronization, or to hub-driven desynchronization. Hub-driven RS may provide a mechanism to account for the role of structural hubs in the organization of brain functional connectivity networks.

network 20 , we show that a degree-dependent distribution of natural frequencies results in complex patterns of remote synchronization, and that a shift in a single hub's frequency can induce dramatic changes in synchronization patterns.
Finally, we study the conditions for this phenomenon to emerge in brain connectivity networks. Specifically, we leverage recent electrophysiological 21 and structural connectivity data 22,23 to model the dynamics of spontaneous activity in the macaque brain, and demonstrate a potential role for hub-driven remote synchronization in shaping patterns of coherent activity, sometimes referred to as functional connectivity.

Remote synchronization in a star network of Kuramoto oscillators. We adopt a Kuramoto-
is degree of the node i. In case of star-like network ( Fig. 1) we set φ = ϕ 0 to be a hub and denote ε = A a i0 and α δ = i 0 for i = 1… N, ε = B a j 0 and β δ = j 0 for j = 1 … N. Thus, for identical phase oscillators described in ref. 25  where φ denotes the phase of the hub (or leader) and ϕ k the phases of the leaf oscillators. In this case, a synchronous solution can include constant phase differences between oscillators. System (2) was analytically analyzed in ref. 25 by virtue of the Watanabe-Strogatz ansatz 26,27 , a variable transformation akin to the Mobius transform (see section Methods for details). This method provides a description of the dynamics of a system like (2) through a complex global variable, namely the order parameter 28 . As shown in ref. 25, the dynamics of this system includes hysteretic transitions between asynchronous and synchronous states. Depending on model's parameters, system (2) can have three different stable solutions.
(i). In the region of relatively small absolute values of frequency mismatch ω ω − 0 , the system (2) has one stable synchronous solution that is stationary in the rotating reference frame. For this solution, the phases are identical, while the phase of the leader φ = Φ − ∆Φ, where ∆Φ = const. Hence, the synchronous solution has non-zero but constant phase shift between the leader and the leaves. (ii). With increasing frequency mismatch, an asynchronous regime emerges. Stability of the asynchronous solution depends on the sign of the expression α β ω ω + − sign(sin( ))( ) 0 . For positive values, the asynchronous solution is stable, while for negative values it becomes unstable. (iii). When the frequency mismatch is too large to lock the phases of the leader and the leaves, the synchronous solution (i) turns into the RS regime, whereby all the leaf oscillators have same frequency and phase, while the leader's frequency and phase are different. In ref. 25 this is called "synchronous limit cycle" solution. We note that this regime corresponds to the definition of RS according to Bergner et al. 14 . A specific condition for the RS regime to be stable is:   Figure 2a demonstrates the RS solution of system (2) for the model parameters satisfying expression (3). This figure shows that, under these conditions, the difference between the synchronized leaves and the hub's phases builds up in time. The Watanabe-Strogatz approach is applicable only if the natural frequencies of the peripheral nodes are identical; however the RS regime can be also observed in the inhomogeneous case, as illustrated on Fig. 2b.
From the brief description reported above and expression (3) it follows that the case of zero overall phase shift (α β = − ) is special. Indeed, stability analysis 25 shows that in this case the asynchronous fixed point is neutrally stable. This implies that the RS regime is neutrally stable as well, thus making it difficult, if not impossible, to detect it numerically in a pure Kuramoto model with zero phase shift in the coupling term. This is the reason why the RS regime was not observed with pure Kuramoto oscillators in ref. 14 and 18, where the authors concluded that RS can only be observed in the presence of an additional degree of freedom, like amplitude in the Stuart-Landau equations. This additional degree of freedom was thought to be necessary for the appearance of RS regime, as it enables a hidden transfer of information through the amplitudes of the oscillators. However, the appearance of RS regime can arise from direct action of the hub on the leaves, as discussed below.
In summary, in this section we have presented the main results of ref. 25 and conclude that the regime therein described for Kuramoto oscillators corresponds to remote synchronization as defined by Bergner et al. 14 , in sharp contrast with their conclusion that remote synchrony cannot occur in pure phase oscillators' systems 14, 18 . Hubs actively drive synchronization. In order to show that a leader (or hub) can have a direct synchronizing or desynchronizing effect on peripheral nodes depending on the model's parameters, we consider a star network where the leaf oscillators are additionally subjected to a Kuramoto-Sakaguchi mean field. This system was introduced but not analyzed in ref. 25 in the form: System (4) can be obtained from (1) by setting ε = C a ij and γ δ = ij for = … i j N , 1 . Below we demonstrate that hub nodes can induce non-trivial regimes in the presence of attractive or repulsive mean fields. Specifically, we find: (i) asynchronous solutions in the presence of an attractive ( γ > cos 0) mean field. Figure 3a shows that in this regime the leader desynchronizes leaf oscillators in the bounded manifold of the positive frequency mismatch ω ω − 0 . (ii) Synchronous solutions in the presence of a repulsive ( γ < cos 0) mean field. Figure 3b shows that the leader synchronizes the leaves in the bounded region of the negative frequency mismatch ω ω − 0 . As in the conventional RS regime, the leader maintains free dynamics, asynchronous with respect to the leaves. The above examples show that the leader can have a synchronizing action even in the presence of a desynchronizing mean field or, vice versa, a desynchronizing one in the presence of a synchronizing mean field. This is inconsistent with the notion that the hub serves the sole purpose of enabling information transfer between the peripheral nodes, and suggests that it directly drives the dynamics of the leaves.
By virtue of the Watanabe-Strogatz transformation, we performed a stability analysis of the synchronous solutions (see section Methods for details). The synchronous solution is stable if λ < 0 s (5), with: . From expression (5) it follows that even for relatively large mean-field strength the hub counteracts the effects of the field when its frequency mismatch falls in a certain range (Fig. 4). For attractive mean-fields, active hub-driven desynchronization is observed for positive frequency mismatches (ω ω − > 0 0 , with ω 0 indicating the natural frequency of the hub). Vice versa, negative mismatches can drive synchronization in the presence of repulsive mean-field.
Hub-driven remote synchronization in a prototypical complex network. In order to investigate hub-driven remote synchrony in more complex systems, we implement the Kuramoto-Sakaguchi model (1) on the Karate club network 20 , a widely studied social network. This network represents a middle ground between a star-like toy model and more complex real world networks, and possesses clearly identifiable hub nodes.
As shown in the previous section, if the sine of the cumulative phase shift is positive, the frequencies of the hubs must be sufficiently larger than the frequencies of the peripheral oscillators (leaves) to observe the RS regime. By way of example, we adopt the distribution of frequencies suggested in ref. 9, where the frequency of the oscillator is proportional to the degree of the node ω = k i i (by rescaling time the coefficient of proportionality can be absorbed in the coupling strength ε). This distribution of frequencies has no physical meaning for this particular network, and was taken solely to demonstrate the possibility of a RS regime in a small complex network with hub nodes.
As a measure of synchrony between nodes, we take the time-averaged order parameter called synchronization index: where ⋅ t denotes an average over large period of time. Note that this order parameter is not sensitive to constant phase shifts, and takes large values when the average frequencies of the two oscillators are similar. We consider nodes i and j to be synchronized if > . r 0 75 ij , a value that corresponds to a minimum in the histogram of node-wise synchronization indices. Other thresholds in the range 0.7-0.9 produce qualitatively and quantitatively similar results.
In order to study how the frequencies of the hubs affect synchronization of the network a multiplier ω x is introduced to selectively vary the frequencies of the two hub nodes (33 and 34). Figure 5 shows the number of synchronized clusters of nodes N c and the size of the largest cluster S c as a function of ω x .
For each value of ω x , several numerical simulations were performed with different, random initial conditions. With increasing ω x the size of the largest cluster follows a bell shaped curve, with a maximum corresponding to ω x = 0.4. Conversely, the number of clusters shows a complementary behavior, albeit less pronounced.
Patterns of synchronization on the Karate-club network for different values of ω x are shown in Fig. 6, where grey lines denote structural links, and red and green dashed lines connect nodes that are directly or remotely synchronized, respectively. In the case of non-zero phase shifts δ, we note that the highly connected hubs 33 and 34 generate a cluster of remotely synchronized nodes while remaining asynchronous.
At ω x = 0.7 remote synchronization starts to emerge, comprising small clusters of nodes (Fig. 6a).
Above ω x = 0.7, the frequencies of the hubs are large enough to remotely synchronize their leaf nodes and the size of the remotely synchronized cluster increases sharply (Fig. 6b). We observed a similar behavior also for other hub nodes, like node 1 (data not shown), thus suggesting that the mechanism is not idiosyncratic to the particular configuration of node 33 and 34, but is generally applicable to hubs.
In the case of zero phase shift δ = 0, occasional pairwise correlations may appear (e.g. between nodes 15 and 16 in Fig. 6c) depending on the initial conditions, but extended clusters of remotely synchronized oscillators do not emerge.
In summary, hub-driven remote synchrony is not limited to toy star-like networks, but can be found in more complex networks endowed with hubs, like the Karate club network. In this model, we found that synchronization patterns strongly depend on the hubs' natural frequencies.  Fig. 3a and b, respectively. Hub-driven remote synchronization in brain connectivity networks. The interplay between structural connections and synchronization of dynamical processes on networks is reminiscent of the concept of functional connectivity in the realm of neuroscience. In this context, functional connectivity is defined in terms of correlations or coherence between oscillatory behaviors observed, e.g., in the electrical or hemodynamic spontaneous activity of the brain.
The question we intend to address in the following is whether the conditions for remote synchronization exist in brain networks and may play a role in determining patterns of correlated activity as observed, e.g., in EEG or MEG neuroimaging experiments in the brain under resting conditions. The key ingredients for this phenomenon to emerge are: (i) the presence of hubs within brain networks of structural connectivity; (ii) a phase shift between remote nodes resulting from a delay in the interaction terms; (iii) a difference in the oscillatory frequency of hubs with respect to their peripheral nodes.
A number of studies (see refs 29 and 30 for recent reviews) based on Diffusion Tensor Magnetic Resonance Imaging have shown that certain brain regions are characterized by high degree and "betweenness", i.e. they act as gateways of many shortest paths connecting pairs of nodes. These connector hubs play an important role in the integration of the network and in ensuring efficient transfer of information across the graph. This has been observed in humans 31 as well as in other species, including non-human primates and rodents 32 . Ex-vivo studies using anterograde or retrograde tracers corroborate this evidence in experimental laboratory animals like the macaque. Hence, condition (i) appears to be fulfilled.
Finite signal propagation speed in axons generates distance-dependent time delays in the interaction terms in brain networks [33][34][35] . In a model of coupled phase oscillators, these delays can be represented as phase shifts proportional to internodal distance 33,36 , as in condition (ii).
Finally, condition (iii) has been recently addressed in a meta-analysis of electrophysiology experiments in the macaque, showing anatomical dependence of spontaneous oscillations of populations of neurons in a number of brain areas as measured by invasive electrophysiology 21 . Comparison with degree distribution in macaque 32 shows that the fast nodes from 21 are also structural hubs.
In the light of this evidence, we have modeled the synchronization phenomenon in the macaque brain using experimental data and empirically determined parameters from the literature.
For our simulations, we adopted the structural connectivity graph described in 22,23 , where connectivity data was obtained by retrograde tracer injections in 29 areas of the macaque cerebral cortex. Extrinsic fraction of labeled neurons (the ratio between the number of labeled neurons in the source area over the total number of labeled cortical neurons extrinsic to the injected area) for each pathway determines the weight of the connection between areas. The 29 by 29 connectivity matrix was thresholded by percolation analysis of the giant component [37][38][39][40] and binarized 41 .
Simulations were performed based on this connectivity graph by adopting the model (1) without normalization of the coupling strengths by node degree. Following 25,33,42 we assumed that the phase shifts in the coupling terms are proportional to the distances between nodes. Internodal distances were taken from ref. 23 as the length of the shortest trajectory interconnecting areas via the white matter, approximating the axonal distance.
Murray et al. 21 collected measurements of timescales of intrinsic fluctuations in spiking activity in different areas of the macaque brain. In ref. 43, the selection of regions was extended using modeling of the macaque neocortex. Specifically, Chaudhuri et al. 43 calculated autocorrelation functions of activity in response to white-noise input to all areas. Time constants of the decay of autocorrelation were calculated for all nodes included in our model. The dominant time constants in various areas of the network were extracted by fitting single or double exponentials to the autocorrelation. In case of double exponentials, the timescales of the two components were calculated as the weighted average of the two time constants. For the frequencies in our model we took inverse timescales.
The results of our simulations are shown on Fig. 7. In the structural connectivity network, a node in the prefrontal cortex denoted as 10 plays the role of structural hub that connects areas in the frontal and temporal cortices. Using the frequency distribution calculated from ref. 43, we obtain a wide pattern of synchronization (red dashed lines), including all frontal and temporal cortices (Fig. 7a). When we switch the frequency of node 10 to that of the fast component associated with this node 43 , coherent synchronization among those areas is preserved, despite the fact that node 10 becomes asynchronous with the rest of the cluster, consistent with the definition of remote synchronization (green dashed lines in Fig. 7b). When we further increase the frequency of the oscillator associated with node 10, frontal and temporal areas become functionally disconnected, and form two separate clusters of synchronized nodes. Hence, a switch between the regime of normal synchronization (Fig. 7a), remote synchronization (Fig. 7b) and asynchrony (Fig. 7c) can be driven by a single parameter in the model, namely the hub's frequency.
This finding is particularly interesting as it may indicate a new and potentially important mechanism in the emergence of functional connectivity patterns in the brain. Here, we have shown that the appearance of clusters of synchronized areas can be driven by structural hub regions that do not appear to be functionally connected to these areas. This is consistent with the observation that structural and functional hubs in the brain do not necessarily coincide 44,45 . Importantly, changes in the frequency of spontaneous fluctuations in a hub region can dynamically reconfigure patterns of functional connectivity in the brain, switching between the regimes of direct and remote synchronization, or actively desynchronizing functional modules even when couplings between nodes are attractive. Changes in spontaneous fluctuations may occur in response to external stimuli, or perhaps changes in brain state, resulting in a reconfiguration of the synchronization pattern.
While the factors enabling hub-driven RS appear to be present in the brain, a role of this phenomenon in shaping functional connectivity patterns remains to be proven. However, this hypothesis can be experimentally tested. By way of example, it may be envisaged that optogenetic technology 46 , whereby photosensitive proteins are expressed in specific populations of neurons, could provide a means to drive the oscillatory behavior of selected hub regions while measuring patterns of synchronized activity across the brain.

Discussion and Conclusion
In the first part of this paper we have summarized the analytical results obtained in ref. 25 for the Kuramoto-Sakaguchi model on star networks. Specifically, we have explicitly shown that "remote synchronization" (RS), in Bergner's definition 14 , can be observed in systems of Kuramoto phase oscillators if a non-zero phase shift is introduced in the coupling terms. For pure star networks, the RS regime is stable when a sufficiently large frequency mismatch is imposed between leaves and the hub. Importantly, we have shown that the hub can exert a synchronizing action even in the presence of additional repulsive mean field acting on the peripheral nodes. Hence, the hub does not simply enable transfer of information between the leaves, but actively drives synchronization while remaining asynchronous with the rest of the network. Conversely, the hub can actively desynchronize the leaves in the presence of an attractive mean field. We have analytically derived the conditions whereby these phenomena can occur in a star network by a stability analysis of the synchronous solutions.
In the second part we have explored the role of this mechanism in the synchronization of an exemplary complex network. As an example, we have chosen the Karate-club network, a prototypical social network in which a few highly connected individuals play the role of hubs. When a degree-dependent distribution of natural frequencies is introduced in the model, remote synchrony emerges and plays a substantial role in the formation of clusters of synchronized nodes within the network.
Finally, we have explored the potential role of this mechanism in brain connectivity networks. Indeed, the prerequisites for hub-driven remote synchronization appear to exist in the brain, including the presence of structural hubs, delays in the couplings between nodes, and region-dependent frequency of spontaneous fluctuations. We have leveraged connectomic data from the macaque brain and recent electrophysiological measurements to derive a dynamical model of synchronization in the macaque cortex. Our simulations show that a change in the intrinsic frequency of a hub can dramatically reshape synchronization patterns, shifting from direct to remote synchronization, and to a hub-driven desynchronization regime. This experimentally testable hypothesis may explain the mismatch between structural and functional hubs sometimes observed in brain connectivity networks.

Methods
The Watanabe-Strogatz approach. The Watanabe-Strogatz (WS) theory 26,27 can be applied to any system of the form k i k where f(t) is an arbitrary real and F(t) is an arbitrary complex functions. Note that a global coupling can enter one or both functions. The dynamics of the general system (7) can be characterized by one global complex variable = z z t ( ) and one real global variable Ψ = Ψ t ( ) and N constants of motion ψ k (of which only N − 3 are independent) by virtue of WS variable transformation Equations for z(t) and Ψ(t) are obtained by substituting (8) into the system (7). 2 Where z(0) and Ψ(0) together with the constants ψ k are determined by initial conditions of original variables ϕ (0) k . Stability analysis. Stability analysis of the system (4) is performed with the help of WS theory presented above. Following 25 we perform variable transformation to the phase differences ϕ ∆ k H t e Im( ( )) Im( ) Im( ( ) ),  (11) and (7) gives As shown in refs 25 and 28, in the thermodynamic limit → ∞ N and for a uniform distribution of constants of motion ψ the order parameter j is equal to z t ( ). In ref. 47 it was demonstrated that, in case of small perturbations, initially non-uniform distributions of constants tend toward the vicinity of the uniform one. Therefore, in this case From the expressions (13) for the global fields, it is apparent that Ψ(t) does not enter the equation for z(t) in (9). Thus the equation for z(t) (14) describes the dynamics of the system (11) and consequently of the original system (4). where − corresponds to γ ∆ − > x qsin 1 and + to ∆ γ − <− x qsin 1. Taking into account that out of the two steady solutions only ΔΦ s1 is of interest and returning to the original parameters, we obtain the expression (5) for λ s .