Kuramoto model simulation of neural hubs and dynamic synchrony in the human cerebral connectome
BMC Neuroscience volume 16, Article number: 54 (2015)
The topological structure of the wiring of the mammalian brain cortex plays an important role in shaping the functional dynamics of large-scale neural activity. Due to their central embedding in the network, high degree hub regions and their connections (often referred to as the ‘rich club’) have been hypothesized to facilitate intermodular neural communication and global integration of information by means of synchronization. Here, we examined the theoretical role of anatomical hubs and their wiring in brain dynamics. The Kuramoto model was used to simulate interaction of cortical brain areas by means of coupled phase oscillators—with anatomical connections between regions derived from diffusion weighted imaging and module assignment of brain regions based on empirically determined resting-state data.
Our findings show that synchrony among hub nodes was higher than any module’s intramodular synchrony (p < 10−4, for cortical coupling strengths, λ, in the range 0.02 < λ < 0.05), suggesting that hub nodes lead the functional modules in the process of synchronization. Furthermore, suppressing structural connectivity among hub nodes resulted in an elevated modular state (p < 4.1 × l0−3, 0.015 < λ < 0.04), indicating that hub-to-hub connections are critical in intermodular synchronization. Finally, perturbing the oscillatory behavior of hub nodes prevented functional modules from synchronizing, implying that synchronization of functional modules is dependent on the hub nodes’ behavior.
Our results converge on anatomical hubs having a leading role in intermodular synchronization and integration in the human brain.
The human brain has to continuously process a large variety of information relating to vision, hearing, motor function, and numerous associative processes. The topological organization of the brain’s wiring pattern has an important role in enabling complex functionality, evidenced by studies examining structural wiring architecture of mammalian and nematode neural systems that revealed several network attributes of an efficient processing and communication structure [1, 2]. Neural systems have been shown to exhibit synchrony patterns associated with communication between segregated anatomical neural components in terms of neuronal discharges, beta and gamma band fluctuations and resting-state blood oxygenation level dependent (BOLD) signals of brain regions [3–9]. Synchronized activity of specialized groups of neurons has been proposed to provide a neural basis for binding and integrating information [10–12]. Moreover, brain disorders such as schizophrenia, epilepsy, autism and Alzheimer and Parkinson disease have been associated with abnormal neural synchronization .
These observations of (dys)synchrony called for the formulation of mathematical models to provide insight into the neural dynamics and its relation to the complex architecture of the underlying structural pathways [13, 14].
Embracing network science as a general mathematical framework to understand topological organizational features of neural wiring, studies have revealed high levels of clustering and community structure, suggesting the existence of functional sub-systems, as well as short communication pathways and the existence of high-degree hub nodes that take a central role in the overall network organization [3, 15–18]. Furthermore, recent research has shown that a highly interconnected collection of high-degree nodes, the ‘rich club’, is embedded across many functional domains and may thus serve as an anatomical infrastructure for global neural communication between otherwise segregated systems [1, 2, 16, 19, 20].
Using the Kuramoto model  to simulate brain interactions through synchrony on the basis of structural connectivity, functional implications of the organization of brain connectivity can be examined. These include for example how cross-domain integration can be influenced by the underlying wiring anatomy of the brain. Simulation studies employing the Kuramoto model have investigated synchronization patterns in the cortical brain networks of the cat , the macaque  and the human [23–25], showing correspondence with resting-state functional data supportive of the applicability of the model [24, 26]. Here, extending to previous work on random or targeted attack in structural brain networks [22, 27, 28], we studied the synchronization effects of suppressed structural connectivity and perturbed intrinsic oscillatory behavior of brain regions in the human cortical network. Simulating these phenomena we shed light on the role of anatomical hubs in functional dynamics of the brain network, which may ultimately contribute to the understanding of effects of structural brain pathology on brain function.
The cohort consisted of the 40 healthy subjects (mean age ± standard deviation 27 ± 6.9; male/female: 27/13) previously reported on by Van den Heuvel et al. and used as a replication set in . Each subject gave informed written consent according to the Declaration of Helsinki and the study was approved by the medical ethics committee for research into humans of the University Medical Center Utrecht. Scanning sessions lasted 45 min and were performed on a 3 Tesla Philips Achieva Clinical scanner. For each subject, diffusion weighted imaging data (acquisition parameters: SENSE-p = 3; two sets of 30 different weighted directions, and 2 × b = 0 images; repetition time (TR)/echo time (TE) = 7,035/68 ms, 2 × 2 × 2 mm, 75 slices covering whole brain, b weighting of 1000 s/mm2, second set with reversed k-space readout) and an anatomical T1 image for anatomical reference (3D fast field echo (FFE) using parallel imaging; TR/TE = 10 ms/4.6 ms; field of view (FOV) = 240 × 240 mm, 200 slices covering whole brain, 0.75 × 0.75 × 0.75 mm) were acquired.
The structural network from  was recomputed at a 219 × 219 network resolution using a subdivision of the Desikan–Killiany atlas . Each area of the brain was designated as a node in the network, and edges represent anatomical connections between these nodes. For each individual, if reconstructed white matter fiber(s) were found between two nodes i and j, then a 1, indicating the presence of a connection, was placed in cell c(i,j) of the 219 × 219 connectivity matrix, M. If no fibers were found connecting regions i and j, then these nodes were taken to be non-connected and the corresponding entry in the weighted connectivity matrix M was set to zero. From the individual connectivity matrices a binary group matrix was constructed having a 1 for connections present in at least 40 % of the subjects and a 0 otherwise .
Cortical network, rich club, and functional modules
A template of resting-state networks (RSNs) [32–35], distinguishing between subnetworks of regions that show correlated activity during rest based on resting-state functional MRI (fMRI) data from a group of healthy subjects, was used to define 11 commonly reported RSN networks (Additional file 1: Figure S1) . Van den Heuvel and Sporns 2013 used Independent Component Analysis (ICA) decomposition on the resting-state fMRI data to define this RSN template. The 11 functional network maps were realigned to the group-averaged normalized T1 images of the subjects and resampled to the 219 parcellation scheme of the diffusion data. As performed in the Van den Heuvel and Sporns paper , each node was assigned a single RSN by means of a majority vote based on its voxels’ participation across the 11 RSNs. This resulted in each network node participating in a single functional module, with the functional modules varying in size from 11 to 31 nodes. Following the line of thought of previous studies [32, 33], these 11 RSNs were taken as a description of 11 functional domains (i.e. modules) of the human brain.
Hub nodes were selected on the basis of degree of connectivity. Consistent with the size and layout of a neural rich club as reported in literature [2, 3, 36], the thirty-nine nodes with the highest number of connections (17.8 % of all nodes) were designated as members of the rich club. In concurrence with previous findings, the anatomical rich club included the insula, precuneus, superior frontal and superior parietal regions and participated in all of the 11 functional modules [16, 20].
In order to verify our findings on the human cortical network, a reconstruction of the macaque cortical network was used based on collated data from tract-tracing experiments. To this end the open source CoCoMac neuroinformatics database  was queried for white matter axonal projections linking cortical regions of the Felleman and Van Essen 91 atlas . The Felleman and Van Essen 91 atlas describes a single hemisphere and includes a parcellation of the cerebral cortex into 78 non-overlapping regions. For all node pairs a 1 was included in the symmetric structural connectivity (SC) matrix if the query returned at least 1 report and if of these reports at least 40 % were positive reports; otherwise a 0 was included in the matrix. The resulting SC matrix had a 27.6 % density.
Since no predefined functional modules were available for the macaque, modules were created based on structural connectivity data by maximizing the number of within-group edges . Modules were selected so that each module contained a minimum of 8 nodes out of the total of 78 nodes.
The Kuramoto model
The Kuramoto model simulates dynamic behavior of a set of coupled oscillators. Applying the Kuramoto model to the group-averaged anatomical brain network, each node was assigned an oscillator with a fixed, random internal angular frequency and an initial random phase. The model is defined as,
where θ i (t) and ω i are the phase and internal angular frequency of oscillator i respectively. The cortical coupling strength that is applied to the edges is denoted by λ. W ij is the symmetrical group matrix containing all connections between cortical nodes (Fig. 1) and N is the total number of nodes. The set of differential equations (Eq. 1) was solved numerically using a Runge–Kutta solver. During every run of the model, each node was assigned a random initial phase, θ, uniformly distributed between [−π, π] and a random internal frequency, ω, uniformly distributed between [0, 1]. Simulations were carried out for T = 700 with a transient time of τ = 300, keeping these parameters the same as in previous work . Explorative longer runs of T = 5000 did not change the nature of our findings. Output data from the last 400 time points was then used in the analyses.
Synchronization metrics and computational details
Two order parameters r and r link were used to describe global dynamic coherence.Footnote 1 The first order parameter, r, is based on a complex number z defined as,
where r(t) is the modulus of z(t). The variable r(t) is a measure of the phase synchrony among the population of N oscillators. Φ(t) is the average combined phase of all oscillators. Averaging r(t) over time yields the first order parameter, r.
In addition to the average phase coherence, r, the average edgewise synchrony also provides information about the state of the network [3, 40]. The edgewise synchrony between nodes i and j, C ij , is defined by
together forming C, the simulated functional connectivity matrix. Average edgewise synchronization levels are taken together to yield rlink,
where N is the total number of nodes of the network. The phase coherence, r, and the fraction of synchronized node pairs, r link, together describe the global dynamics of the system.
Going beyond the examination of global dynamics, synchronization dynamics were also studied on a modular level, analyzing nodes together on the basis of their functional module assignment. In order to assess the likelihood of particular nodes to be in synchrony, node pairs need to be classified as either synchronized or non-synchronized. To this end the r link parameter, ranging from 0 (total incoherence) to 1 (perfect synchrony), was used to define the proportion of nodes that are in synchrony. Assuming a threshold level of synchronization, C ij , exists above which the nodes’ oscillatory behavior converge to full synchrony, a theoretically equivalent system with equal order parameter r link, was constructed composed of N(N-1) × r link perfectly coherent and N(N-1) × (1-r link) completely incoherent node pairs, described by a binary synchronization matrix F ij :
The number of edges present in F now matches the number of synchronized connections predicted by r link. The node pairs showing the highest synchronization levels are thus categorized as synchronized, and those with lower synchrony as unsynchronized. In order to study synchronization levels independent of initial conditions, n = 103 trials were run for each cortical coupling strength λ. The F ij values were averaged across trials into r ij such that,
where r ij reflects the probability of regions i and j to be synchronized and where l is the trial number. Based on the probability of synchronization between any two areas, modular synchronization, r αβ , is calculated as,
where i describes the nodes in module α, and j are nodes in module β. N α and N β are the number of nodes in module α and β respectively. In the case of intramodular synchrony (α = β), Eq. 7 was slightly adapted to Eq. 8 so that synchrony within a node (i.e. r ii ) was eliminated, to avoid module size to have an effect on intramodular synchrony.
The mean field representation is an alternative representation of the Kuramoto model, describing each oscillator’s dynamics by a centroid with a mean radius and a mean phase:
where R k is the radius and ϕ k the phase of the centroid representing node k and D k is the degree of node k (i.e. the number of other nodes it is connected to). The centroid’s radius R k captures the phase coherence among the nodes that node k is connected with and, as can be seen from Eq. 9, it denotes the strength by which the mean field influences the frequency of node k. Rewriting Eq. 10 as follows yields each module’s contribution to the mean field and therefore one obtains a measure of how strongly each of the modules influence the frequency of node k:
where D k,α is the number of connections between node k and module α, R k,α stands for the influence module α has on the frequency of node k, and ϕ k,α is the mean phase of module α nodes that are connected to node k. A similar expression can be used to describe the influence of hub nodes’ on the frequency of node k:
If node k is connected to only one node of module α, D k,α = 1, it follows that R k,α = 1 which cannot be interpreted as a measure of phase coherence. Therefore analyses were based on instances of D k,α > 1.
As a measure of the impact of module nodes or hub nodes on nodes in another module, values of R k,α were averaged over all nodes k belonging to the collection of nodes of interest, e.g.:
similarly, for hubs:
where R β,α (or R β,hubs ) denotes the average influence of module α nodes (or hub nodes) on the frequencies of module β nodes.
To assess the effects of suppressed connectivity on synchrony, selected edges were computationally removed from the adjacency matrix . Hub connectivity was suppressed by decoupling all hub nodes with respect to each other. Since this type of hub connectivity suppression alters the overall connectivity density (which may have general consequences for network synchronization), an equal number of random edges was removed for reference purposes in the test condition. For trials with suppressed connectivity, the structural adjacency network, W, was adjusted by,
where K is the subset of nodes selected to be modulated. It is important to note that in the modulation conditions all hub nodes remained connected to the rest of the network as normal.
Nodal frequency perturbation
In the trials where oscillator frequency was tracked, the data was grouped by module and a Fourier transform was applied in order to find the most powerful frequency within the modules. In the model trials, in which a set of nodes (i.e. hubs, random nodes or all nodes belonging to one of the modules) had their internal frequencies perturbed, the nodes were first assigned a random initial phase and internal frequency (as described above). Then, nodes within the selected set had their angular frequency increased by one.
To test the variation among dominant frequencies of the functional modules, a Bartlett’s test was performed on the unperturbed modules at the cortical coupling factor immediately prior to when the perturbed set of nodes became synchronized with the rest of the modules. For this purpose, whole brain synchronization was defined as the point at which lowest modular frequency was within 5 % of the highest modular frequency.
Assessing global synchronization in the anatomical brain network with the Kuramoto model revealed a modular state for low cortical coupling strengths, λ, that quickly transitioned into whole brain synchrony at increased cortical coupling. The level of global synchronization was measured by the Kuramoto model’s order parameters r and r link, denoting the extent to which nodes are in phase and the share of in-phase nodes respectively, where r (or r link) close to zero indicates global incoherence and r (or r link) near one reflects a state of global synchrony.
Figure 2a shows how levels of global synchronization progress with respect to cortical coupling strength. The order parameter r has been shown to follow a square-root function of the coupling strength λ , therefore the shape of the curve results from a square-root curve smoothed out by the heterogeneity in the intrinsic frequencies and by the network topology of the oscillators. A critical regime was observed between λ = 0.02 and λ = 0.04, indicated in blue in Fig. 2a, during which the network abruptly switches from brain regions being synchronized locally to a state of global synchrony, consistent with previous work on the cat cerebral cortex . Figure 2b shows three panels of synchronization levels, or simulated functional connectivity, illustrating the modular state (small clusters of synchronized nodes), transition state (widespread synchrony arising) and global synchrony state (nearly all nodes in synchrony), evaluated by the Kuramoto model at cortical coupling factors λ = 0.02, λ = 0.03 and λ = 0.04, respectively.
Overlaying the reconstructed structural brain network, described by a binary 219 × 219 matrix derived from diffusion weighted imaging (DWI) data of 40 healthy subjects, with 11 commonly reported functional modules (see “Methods”), revealed that functional modules were structurally more densely connected as compared to the overall density of the SC matrix (Fig. 1, p < 10−3, permutation test, 103 random selections of nodes). As a result, the highly internally connected functional modules produced higher levels of synchrony in our Kuramoto simulations, peaking at λ = 0.02. Synchrony was found to be significantly higher among node pairs that were directly connected by a structural connection compared with structurally unconnected node pairs (factor 3.9, p < 10−4, permutation test, 104 random selections of unconnected node pairs), reflecting the apparent overlap between the lower triangular part of the matrix of Fig. 1, the SC matrix, and the upper triangular part containing the simulated synchronization levels. The high structural density of the functionally defined modules and their high levels of synchronization are thus indicative of a positive structure–function relationship, in accordance with recent literature [43–45].
Hub selection and participation
The anatomical brain network showed a clear rich club organization [2, 20, 46], with nodes of high degree, k, or hub nodes, being more densely interconnected than expected by chance (11 < k < 44). Next, using k = 26, we selected the 39 nodes (18 % of all nodes) having the largest number of connections (degree), as the set of anatomical hubs [2, 3, 36]. These hub nodes were found to be distributed throughout the brain network and to participate in all 11 functional modules (Fig. 1), in accordance with the central embedding of hubs in the anatomical network and their suggested ability to facilitate communication and integration between anatomically segregated brain regions [3, 16].
Intramodular synchrony, reflecting the share of in-phase nodes within a functional module, was observed to evolve as an s-curve with respect to coupling strength in all 11 functional modules (Fig. 3a) again with a critical regime between λ = 0.02 and λ = 0.04 showing a steep increase in intramodular synchrony similar to global synchrony (Fig. 2). In Fig. 3b the synchronization among the hub nodes (intra-hub synchrony) is contrasted with the intramodular synchrony levels observed in each of the 11 functional modules. While similar in shape, the curve corresponding to the hub nodes is observed to be shifted towards lower cortical coupling strengths with respect to the graphs of the functional modules. In other words, intra-hub synchrony is higher than any module’s intramodular synchrony (p < 10−4, 104 random permutations of connection labels ‘intra-hub’ and ‘intramodular’, 0.02 < λ < 0.05), implying that the rich club structure, spread out across the brain network and involved in all functional modules, requires less cortical coupling in order to reach synchrony than the similarly dense or even denser functional modules.
Furthermore, nodes participating in the same functional module were found to be more closely synchronized to each other than nodes from different modules, reflected by a 1.5 times higher synchrony between nodes within a module than between nodes in different modules (p < 10−3, 103 random permutations of connection labels ‘intramodular’ and ‘intermodular’, evaluated at the onset of the critical regime, λ = 0.02). Examining hubs showed this class of nodes to have a significantly higher synchronization among themselves than the nodes belonging to any functional module or than an equally sized random set of nodes (p < 10−4, 104 random permutations of node labels ‘hub’ and ‘modular’, 0.005 < λ < 0.075), indicating that global hub nodes are leading the functional modules in establishing local modular synchrony. In addition, synchrony between a functional module and the set of hub nodes was found to be significantly higher than intermodular synchrony (p < l0−4, 104 random permutations of connection labels ‘hub-module’ and ‘intermodular’) for 53 out of all 55 combinations of functional modules (11 × 10/2) (Additional file 2: Figure S2), suggesting the structure of the network supports intermodular synchronization through the strongly tied hub nodes.
The mean field representation of the Kuramoto model and the break down of synchrony contributions across the functional modules and the hub nodes Eqs. 13 and 14, allowed for the comparison between the contributions of modules and hub nodes towards global synchrony. Figure 3c shows the synchrony contributions received by the Default Mode Network from each of the functional modules and the hub nodes (for the other modules see Additional file 3: Figure S3). In the critical regime between λ = 0.02 and λ = 0.04, where the hub nodes become more synchronized among each other, the hubs exert stronger influence on the modules and take on a dominant role in exerting influence on the rest of the network, evidenced by the steepest increase in impact on functional modules compared to the change in influences exerted by the modules (Additional file 3: Figure S3, paired t-test between a module’s and the hub nodes’ impact on the set of 11 modules, p < 0.02). Moreover, at the end of the critical regime, λ = 0.04, the hub nodes’ influence on each of the modules is greater than that of 9 out of 11 modules (paired t-test, p < 0.04) and not significantly different from the two remaining modules (Visual and Bilateral Parietal, paired t-test, p > 0.85).
Hub connectivity suppression
Next, we examined the role of hub nodes in establishing global synchrony. Targeted suppression of structural connectivity was inflicted by computationally removing edges between hub nodes. For comparison, in a separate simulation, an equal number of random edges were removed. In simulations with suppressed connectivity (either between hubs or between random nodes), stronger cortical coupling was required for synchronization to compensate for the loss of connectivity, while the general shape of global synchronization progression towards a state of global synchrony remained the same (Additional file 4: Figure S4). However, simulations in which hub-to-hub connectivity was suppressed showed a significantly higher intramodular to global synchrony ratio than randomly suppressed connectivity simulations (Fig. 4) (p < 4.1 × l0−3, 104 random permutations of hub-to-hub and random connectivity suppressions, 0.015 < λ < 0.04). This increased intramodular to global synchrony ratio as seen in the hub suppressed network indicates that connections spanning between hub nodes (so-called ‘rich club connections’ ) play an important role in the integration of modules needed to reach global synchrony. Also in comparison with the unsuppressed, normal condition, suppression of hubs resulted in increased modularity (Additional file 5: Figure S5), again illustrating that without interconnected hub nodes, functional modules become more separated and the emergence of global synchrony is hampered.
Apart from manipulating connectivity through edges as presented in the previous analysis, synchronization dynamics in the Kuramoto model can also be perturbed by offsetting the initial internal frequencies of a set of nodes. As lower levels of synchronization are a trivial consequence of these frequency perturbations, we tracked nodal frequencies along the path towards global synchrony to examine how different sets of nodes affected the synchronization of the entire system. Figure 5a shows that with the frequency of the hub nodes perturbed, the functional modules remained out of synchrony with respect to each other until the hub nodes stabilized at a frequency near the whole brain frequency corresponding to global synchrony. In contrast, in all other simulations where a different set of nodes was perturbed, the frequencies of the unperturbed modules converged to a common global frequency, reflecting global synchronization, even though the perturbed set of nodes were still operating at a markedly higher frequency. For example, when a random set of nodes (equal in number to the hub nodes) had their internal frequencies perturbed, the remaining unperturbed portion of the system was found to synchronize without the perturbed set of nodes (Fig. 5b). Similarly, perturbing the nodes of any of the 11 functional modules did not keep the other modules from synchronizing among one another (see Fig. 5c for an example in which the Default Mode Network module was perturbed). Nearly identical findings for the remaining 10 functional modules are provided in Additional file 6: Figure S6.
Testing the variation of frequencies among modules in the case of hub node perturbation, compared to if either a functional module or a random set of nodes was perturbed, showed, as depicted in Fig. 5, that hub node perturbation caused the highest frequency variation among the functional modules (p < 6× 10−9, Bartlett’s test). This effect of relative global asynchrony with perturbed hub nodes remained present when a smaller rich club comprising 22 hub nodes, closer to the average functional module size of 17. 9 ± 5.1 nodes, was examined (Additional file 7: Figure S7). Taken together, these frequency perturbation findings suggest that the hub nodes are essential components that cannot operate at frequencies isolated from the rest of the network in order for global synchronization to occur.
A reconstruction of the macaque connectome, based on an extraction from the CoCoMac database  containing anatomical tract-tracing data (see Methods), was used to verify synchronization effects observed in the human cortical network. Analogous to the main analysis of the human cortical network, identical Kuramoto simulations were run on the SC data describing the macaque connectome, of which the results are summarized below. Findings in the macaque were in close agreement with those in the human:
Global synchrony Global synchronization in the macaque brain was observed to increase with cortical coupling strength and revealed a modular state and a state of global synchrony separated by a critical regime (Additional file 8: Figure S8).
Network connectivity Synchrony between structurally connected nodes was higher than synchrony between nodes not directly connected; the ratio between these two values peaked early in the critical regime at λ = 0.015, where synchrony in directly structurally connected nodes was 3.07 times higher (p < 10−5).
Hub selection and participation Rich club organization was found for k > 14 to k = 56. The 15 nodes (19 %) having degree k > 38 were selected as the set of hub nodes. These hubs were again found to be distributed across the cortical network, participating in all functional modules.
Modular synchronization The synchronization among hub nodes and the intramodular synchronization showed the characteristic s-curve progression for increasing cortical coupling strength, with the hub nodes synchronizing at lower cortical coupling than the modules. Furthermore, synchrony between a module and the hub nodes during the modular state (λ = 0.01) was found to be significantly higher than intramodular synchrony for all 6 out of 6 modules and higher than intermodular synchrony (p < 7.2 × 10−4) for all 15 combinations of modules (6 × 5/2).
Hubs showed a steeper increase in their influence on the frequencies of modules during the critical regime than 4 out of 6 modules (paired t-test, p < 0.03) and at the end of the critical regime, the hubs’ impact on each of the modules is greater than that of all 6 modules (paired t-test, p < 0.04).
Hub connectivity suppression Simulations in which hub connectivity was suppressed showed increased modularity (higher intramodular to global synchrony ratio, p < 10−4, 0.01 < λ < 0.03) compared with the removal of an equal sized random set of connections (Fig. 6a).
Hub perturbation Perturbation of the hub nodes revealed that the remaining unperturbed nodes of the network were not able to synchronize before the hub nodes had converged to the frequency of global synchrony (Fig. 6b). In contrast, for a perturbed random set of nodes or a perturbed module, the unperturbed part of the network was still able to synchronize (Fig. 6c, d). Significantly higher variation was observed for the hub nodes perturbed simulations than if a random set of nodes or a module was perturbed (p < 7 × 10−5).
Our simulations of synchrony patterns based on a structural reconstruction of the human cortical network produce converging evidence that hub nodes act as spatially distributed but functionally central integrators of neural information between otherwise segregated functional domains. First, hub nodes showed high levels of synchronization among themselves and were found to be distributed across all functional modules. Second, suppression of connectivity among hub nodes caused the network to show stronger modularity in synchronization patterns. Third, perturbation of hub nodes’ intrinsic oscillatory behavior prevented the functional modules from synchronizing.
We verified these findings in the macaque, using a connectome reconstruction based on collated tract-tracing data.
Functional modules based on resting-state fMRI data displayed high structural density as well as strong intramodular synchrony, providing simulation based indications of a positive structure–function relationship in empirical data [29, 44, 47].
Previous modeling work revealed that the wiring of cortical networks gives rise to the formation of functional modules [14, 48–50]. Our simulations over a range of coupling strengths showed a critical regime [3, 25] separating a modular state from the state of global synchrony. Although caution is needed to directly interpret the simulated cortical coupling strengths of cortico-cortical connections as a mechanism to allow a cortical network to switch between a modular state and global synchrony, our results do support the notion of topological organization of the cortical network to enable theoretical synchrony both at the modular and at the global level . The state of global synchrony itself is not to be directly interpreted as biologically meaningful, being merely the attractor state of the Kuramoto model. It is the progression towards global synchrony that reveals how synchronization—and therewith hypothesized binding of information [10, 12]—among segregated brain regions may dynamically occur.
Hub nodes were found to be distributed across the cortical network and to participate in all functional modules, consistent with previous observations [16, 36, 44–46], which are characteristics well-suited to a facilitating role in intermodular communication and integration contributing to global synchrony . In further support of the importance of hub nodes in global synchrony, our simulation results show synchronization among hub nodes to be stronger than the intramodular synchronization of the functional modules, suggesting that the network topology is such that hub nodes form a closely tied and functionally linked ensemble leading the functional modules in establishing synchrony. Anatomically densely intra-connected functional modules showed the highest levels of intramodular synchrony. Hub nodes stood out due to their strong mutual synchronization while having a structural density near the average density of the functional modules, indicating strong synchrony among hub nodes was not purely a consequence of structural density.
Our observation of strong synchrony between hub nodes aligns well with the results of two other Kuramoto simulation studies, one on the macaque cortex, reporting hubs to synchronize at shorter timescales compared to other nodes  and another on the cat cortex reporting anatomical hubs to be leading the transition to whole brain synchrony . Moreover, a recent study of our group in which we used an alternative model for the simulation of brain dynamics (an Ising spin glass model ) showed the tendency of hub nodes to collectively be in an ‘activated state’ and to enrich the overall functional dynamics of the system .
Suppressing connectivity, irrespective of the targeted edges in the network, resulted in a reduction of global and intramodular synchrony, directly related to the sparser structural input. The ratio between intramodular and global synchrony provided more insight into the relative effect of connectivity suppression, with targeting the connections linking hub nodes leading to an increased intramodular to global synchrony ratio. Such an increase in modularity was not observed when a random set of connections was removed (Additional file 5: Figure S5), indicating connections between hub nodes to be particularly important for intermodular communication and integration.
Further evidence of the importance of hub nodes in shaping intermodular functional dynamics was demonstrated by perturbing the intrinsic oscillatory behavior of the hub nodes. While offsetting the internal frequencies of the hub nodes did not have a permanent effect on the network dynamics other than a trivial change in whole brain synchrony state frequency, the functional modules remained desynchronized until the hub nodes approached a matching synchronization frequency (Fig. 5a). Perturbing a functional module or an equal number of random nodes, by contrast, did not prevent the unperturbed modules from synchronizing with each other, leaving the perturbed set of nodes isolated (Fig. 5b, c).
An important consideration resulting from the interpretation of r link as a fraction of synchronized nodes is that the cut-off value of C ij , above which node pairs are considered to be in synchrony, may vary. As a consequence, node pairs with a relatively high level of synchronization, falling just below the cut-off value, can be classified as completely incoherent using this approach. Interpreting r link as the fraction of synchronized nodes is most accurate where there is a clear divide between node pairs in their levels of synchronization. Furthermore, it should be noted that a binary, undirected adjacency matrix was used as structural input to our simulations, resulting in each edge to contribute equally in the model drawing all nodes to a global frequency converging to the average of all internal node frequencies. Future studies could focus on modified Kuramoto model implementations, e.g. including connection weights, for enhanced biological plausibility  to further elucidate the observed patterns of modular and global synchronization. Another limitation arises from the crudeness of the suppression applied to the structural input. Reduced white matter integrity is linked to various psychiatric and neurological disorders, but at the scale of this cortical network, consisting of 219 nodes, disease effects including the loss of entire anatomical connections are improbable. It would therefore be of interest to explore the effects of connectivity suppression in a weighted approach, in which strengths of individual connections can be adjusted.
In the macaque verification analyses ‘functional’ modules were derived from anatomical connectivity. For consistency with our module definition in the human, macaque resting-state fMRI acquisitions would have been ideal. However, with structural and functional connectivity and modularity strongly overlapping both in the human [43, 54] and in the macaque [47, 55], the structurally defined modules are a valid approximation for macaque functional modules.
The synchronization patterns examined in this study provide insight into how anatomical hub nodes may form an infrastructure facilitating communication between functional modules enabling global binding and integration of information.
The Kuramoto model’s ability to simulate structurally dependent functional synchrony could be further utilized to simulate neurological and neuropsychiatric diseases that affect the anatomical connectivity of the brain, and could potentially aid in understanding the relationships between disease structural abnormalities and changes in functional brain dynamics and connectivity.
r link is equivalent to r* link in , in which the Kuramoto model was applied to the cat brain.
Kaiser M, Varier S. Evolution and development of brain networks: from Caenorhabditis elegans to Homo sapiens. Netw Comput Neural Syst. 2011;22:143–7.
Van den Heuvel MP, Sporns O. Rich-club organization of the human connectome. J Neurosci. 2011;31:15775–86.
Gómez-Gardeñes J, Zamora-López G, Moreno Y, Arenas A. From modular to centralized organization of synchronization in functional areas of the cat cerebral cortex. PLoS One. 2010;5:e12313.
Kreiter AK, Singer W. Stimulus-dependent synchronization of neuronal responses in the visual cortex of the awake macaque monkey. J Neurosci. 1996;76:2381–96.
Lachaux JP, Rodriguez E, Martinerie J, Varela FJ. Measuring phase synchrony in brain signals. Hum Brain Mapp. 1999;8:194–208.
Engel AK, Fries P, Singer W. Dynamic predictions: oscillations and synchrony in top-down processing. Nat Rev Neurosci. 2001;2:704–16.
Uhlhaas PJ, Singer W. Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron. 2006;52:155–68.
Uhlhaas PJ, Singer W. Abnormal neural oscillations and synchrony in schizophrenia. Nat Rev Neurosci. 2010;11:100–13.
Ponce-Alvarez A, Deco G, Hagmann P, Romani GL, Mantini D, Corbetta M. Resting-State Temporal Synchronization Networks Emerge from Connectivity Topology and Heterogeneity. PLoS Comput Biol. 2015;11:e1004100.
Tononi G, Edelman GM, Sporns O. Complexity and coherency: integrating information in the brain. Trends Cogn Sci. 1998;2(12):474–484.
Sporns O, Gally JA, Reeke GN, Edelman GM. Reentrant signaling among simulated neuronal groups leads to coherency in their oscillatory activity. Proc Natl Acad Sci USA. 1989;86:7265–9.
Singer W, Engel AK, Kreiter AK, Munk MH, Neuenschwander S, Roelfsema PR. Neuronal assemblies: necessity, signature and detectability. Trends Cogn Sci. 1997;1:252–61.
Honey CJ, Kötter R, Breakspear M, Sporns O. Network structure of cerebral cortex shapes functional connectivity on multiple time scales. Proc Natl Acad Sci USA. 2007;104:10240–5.
Cabral J, Hugues E, Sporns O, Deco G. Role of local network oscillations in resting-state functional connectivity. Neuroimage. 2011;57:130–9.
Kaiser M, Hilgetag CC. Nonoptimal component placement, but short processing paths, due to long-distance projections in neural systems. PLoS Comput Biol. 2006;2:e95.
Van den Heuvel MP, Sporns O. An anatomical substrate for integration among functional networks in human cortex. J Neurosci. 2013;33:14489–500.
Towlson EK, Vértes PE, Ahnert SE, Schafer WR, Bullmore ET. The rich club of the C. elegans neuronal connectome. J Neurosci. 2013;33:6380–7.
Shanahan M. The brain’s connective core and its role in animal cognition. Philos Trans R Soc Lond B Biol Sci. 2012;367:2704–14.
Zamora-López G, Zhou C, Kurths J. Graph analysis of cortical networks reveals complex anatomical communication substrate. Chaos. 2009;19:015117.
Harriger L, van den Heuvel MP, Sporns O. Rich club organization of macaque cerebral cortex and its role in network communication. PLoS One. 2012;7:e46497.
Kuramoto Y. Chemical Oscillations, Waves, and Turbulance. Berlin: Springer-Verlag; 1984.
Honey CJ, Sporns O. Dynamical consequences of lesions in cortical networks. Hum Brain Mapp. 2008;29:802–9.
Kitzbichler MG, Smith ML, Christensen SR, Bullmore E. Broadband criticality of human brain network synchronization. PLoS Comput Biol. 2009;5:e1000314.
Cabral J, Kringelbach ML, Deco G. Exploring the network dynamics underlying brain activity during rest. Prog Neurobiol. 2014;114:102–31.
Villegas P, Moretti P, Muñoz MA. Frustrated hierarchical synchronization and emergent complexity in the human connectome network. Sci Rep. 2014;4:5990.
Vuksanović V, Hövel P. Functional connectivity of distant cortical regions: role of remote synchronization and symmetry in interactions. Neuroimage. 2014;97:1–8.
Alstott J, Breakspear M, Hagmann P, Cammoun L, Sporns O. Modeling the impact of lesions in the human brain. PLoS Comput Biol. 2009;5:e1000408.
Cabral J, Hugues E, Kringelbach ML, Deco G. Modeling the outcome of structural disconnection on resting-state functional connectivity. Neuroimage. 2012;62:1342–53.
Van den Heuvel MP, Kahn RS, Goñi J, Sporns O. High-cost, high-capacity backbone for global brain communication. Proc Natl Acad Sci USA. 2012;109:11372–7.
Cammoun L, Gigandet X, Meskaldji D, Thiran JP, Sporns O, Do KQ, Maeder P, Meuli R, Hagmann P. Mapping the human connectome at multiple scales with diffusion spectrum MRI. J Neurosci Methods. 2012;203:386–97.
de Reus MA, Van den Heuvel MP. Estimating false positives and negatives in brain networks. Neuroimage. 2013;70:402–9.
Beckmann CF, DeLuca M, Devlin JT, Smith SM. Investigations into resting-state connectivity using independent component analysis. Philos Trans R Soc Lond B Biol Sci. 2005;360:1001–13.
Damoiseaux JS, Rombouts SARB, Barkhof F, Scheltens P, Stam CJ, Smith SM, Beckmann CF. Consistent resting-state networks across healthy subjects. Proc Natl Acad Sci USA. 2006;103:13848–53.
Greicius MD, Krasnow B, Reiss AL, Menon V. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc Natl Acad Sci USA. 2003;100:253–8.
Van den Heuvel M, Mandl R, Pol HH. Normalized cut group clustering of resting-state fMRI data. PLoS One. 2008;3(4):e2001. doi:10.1371/journal.pone.0002001.
Van den Heuvel MP, Sporns O, Collin G, Scheewe T, Mandl RCW, Cahn W, Goni J, Hulshoff Pol HE, Kahn RS. Abnormal rich club organization and functional brain dynamics in schizophrenia. JAMA. Psychiatry. 2013;70(8):783–92.
Stephan KE, Kamper L, Bozkurt A, Burns GA, Young MP, Kötter R. Advanced database methodology for the Collation of Connectivity data on the Macaque brain (CoCoMac). Philos Trans R Soc Lond B Biol Sci. 2001;356:1159–86.
Felleman DJ, Van Essen DC. Distributed hierarchical processing in the primate cerebral cortex. Cereb Cortex. 1991;1:1–47.
Rubinov M, Sporns O. Complex network measures of brain connectivity: uses and interpretations. Neuroimage. 2010;52:1059–69.
Brede M. Locals vs. global synchronization in networks of non-identical Kuramoto oscillators. Eur Phys J B. 2008;62:87–94.
Breakspear M, Heitmann S, Daffertshofer A. Generative models of cortical oscillations: neurobiological implications of the kuramoto model. Front Hum Neurosci. 2010;4(November):190.
de Reus MA, van den Heuvel MP. Simulated rich club lesioning in brain networks: a scaffold for communication and integration? Front Hum Neurosci. 2014;8(3):1–5.
Van den Heuvel MP, Mandl RCW, Kahn RS. Hulshoff Pol HE: Functionally linked resting-state networks reflect the underlying structural connectivity architecture of the human brain. Hum Brain Mapp. 2009;30:3127–41.
Honey CJ, Sporns O, Cammoun L, Gigandet X, Thiran JP, Meuli R, Hagmann P. Predicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci USA. 2009;106:1–6.
Honey CJ, Thivierge J-P, Sporns O. Can structure predict function in the human brain? Neuroimage. 2010;52:766–76.
Dereus MA, Vandenheuvel MP. Rich club organization and intermodule communication in the cat connectome. J Neurosci. 2013;33:12929–39.
Miranda-Dominguez O, Mills BD, Grayson D, Woodall A, Grant KA, Kroenke CD, Fair DA. Bridging the gap between the human and macaque connectome: a quantitative comparison of global interspecies structure-function relationships and network topology. J Neurosci. 2014;34:5552–63.
Nicosia V, Valencia M, Chavez M, Díaz-Guilera A, Latora V. Remote synchronization reveals network symmetries and functional modules. Phys Rev Lett. 2013;110:174102.
Senden M, Deco G, de Reus MA, Goebel R, van den Heuvel MP. Rich Club organization supports a diverse set of functional network configurations. Neuroimage. 2014;96:174–82.
Zhou C, Zemanová L, Zamora G, Hilgetag C, Kurths J. Hierarchical organization unveiled by functional connectivity in complex brain networks. Phys Rev Lett. 2006;97:238103.
Arenas A, Díaz-Guilera A, Pérez-Vicente C. Synchronization reveals topological scales in complex networks. Phys Rev Lett. 2006;96:114102.
Dehaene S, Kerszberg M, Changeux JP. A neuronal model of a global workspace in effortful cognitive tasks. Proc Natl Acad Sci USA. 1998;95:14529–34.
Deco G, Senden M, Jirsa V. How anatomy shapes dynamics: a semi-analytical study of the brain at rest by a simple spin model. Front Comput Neurosci. 2012;6(September):68.
Van den Heuvel MP, Sporns O. Network hubs in the human brain. Trends Cogn Sci. 2013;17:683–96.
Scholtens LH, Schmidt R, de Reus MA, van den Heuvel MP. Linking macroscale graph analytical organization to microscale neuroarchitectonics in the macaque connectome. J Neurosci. 2014;34(36):12192–205.
RS, KJRL, and MPvdH conceived of the study design. RS, KJRL, MAdR and MPvdH were involved in study implementation. RS and KJRL performed the analyses. LHvdB contributed to the interpretation of data. All authors contributed to drafting the first manuscript and read and approved the final manuscript.
KJRL was sponsored by a Fulbright Grant (http://www.fulbright.nl). MPvdH was supported by a VENI (#451-12-001) grant of the Netherlands Organization for Scientific Research (NWO) and by a Fellowship of the Brain Center Rudolf Magnus. LHvdB received funding from the Netherlands Organization for Scientific Research VICI Grant, the ALS Foundation Netherlands, the Prinses Beatrix Fonds and from the European Community’s Health Seventh Framework Program (FP7/2007-2013) under grant agreement no. 259867. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Compliance with ethical guidelines
Competing interests The authors declare that they have no competing interests.
Ruben Schmidt and Karl J. R. LaFleur contributed equally to this work
Figure S1. Functional modules derived from resting-state fMRI data. This figure from Van den Heuvel and Sporns, 2013, shows the location of the functional modules based on independent component analysis (ICA) of resting-state fMRI data.
Figure S2. Inter- and intramodular synchrony. Each bar plot corresponds to the inter- and intramodular synchrony of each of the 11 functional resting-state modules and the hub nodes. Shown synchronization levels were evaluated at the onset of the critical regime with a cortical coupling factor λ = 0.02. The hub nodes showed particularly high synchrony compared with the functional modules even though they were distributed across the modules. In some cases, synchrony between the hub nodes and a module was even higher than the module’s intramodular synchrony (see plot 1, 5, 6 and 8). This strong level of synchrony suggests the hub nodes’ importance for synchrony across the network.
Figure S3. Influences on oscillation frequencies of the modules and hub nodes. For each of the 11 modules (1 plot per module), the influences on its frequencies of the modules and hub nodes (12 lines per plot corresponding to the modules and the hub nodes) are shown. The hub nodes become dominant in the process of global synchronization during the critical regime.
Figure S4. Hub versus random connectivity suppression. When edges were removed from the adjacency matrix, attaining whole brain synchrony required a higher cortical coupling factor. Removing both random and edges between hub nodes produces this result.
Figure S5. Modularity increased with hub connectivity suppressed. The ratio between intramodular synchrony and whole brain synchrony is shown for normal, hub connectivity suppressed, and random edge suppressed states. During the critical regime, the random edge suppressed state was observed to be almost identical to the normal state. In contrast, the hub connectivity suppressed state shows a significantly increased intramodular synchrony relative to whole brain synchrony, reflecting stronger modularity in the hub connectivity suppressed state.
Figure S6. Modular frequency tracking during perturbation. Next to the Default Mode module (Figure 5), the effects of perturbation were examined for the remaining 10 functional modules as well. In each instance, the non-perturbed modules synchronized before the perturbed module joined in whole brain synchrony, consistent with the Default Mode perturbation and distinctly different from the hub node perturbation shown in figure 5.
Figure S7. Small set of hub nodes perturbed. Perturbation of a smaller set of hub nodes also prevented the functional modules from synchronizing until the hub nodes’ frequency joined in whole brain synchrony.
Figure S8. Global synchrony progression in the macaque. The order parameters r and r link for the macaque model network progressed in a manner similar to the human network (Figure 2), displaying a modular and a whole brain synchrony state separated by a critical regime.
About this article
Cite this article
Schmidt, R., LaFleur, K.J.R., de Reus, M.A. et al. Kuramoto model simulation of neural hubs and dynamic synchrony in the human cerebral connectome. BMC Neurosci 16, 54 (2015). https://doi.org/10.1186/s12868-015-0193-z