## Abstract

Cognitive function depends on an adaptive balance between flexible dynamics and integrative processes in distributed cortical networks. Patterns of zero-lag synchrony likely underpin numerous perceptual and cognitive functions. Synchronization fulfils integration by reducing entropy, while adaptive function mandates that a broad variety of stable states be readily accessible. Here, we elucidate two complementary influences on patterns of zero-lag synchrony that derive from basic properties of brain networks. First, mutually coupled pairs of neuronal subsystems—resonance pairs—promote stable zero-lag synchrony among the small motifs in which they are embedded, and whose effects can propagate along connected chains. Second, frustrated closed-loop motifs disrupt synchronous dynamics, enabling metastable configurations of zero-lag synchrony to coexist. We document these two complementary influences in small motifs and illustrate how these effects underpin stable versus metastable phase-synchronization patterns in prototypical modular networks and in large-scale cortical networks of the macaque (CoCoMac). We find that the variability of synchronization patterns depends on the inter-node time delay, increases with the network size and is maximized for intermediate coupling strengths. We hypothesize that the dialectic influences of resonance versus frustration may form a dynamic substrate for flexible neuronal integration, an essential platform across diverse cognitive processes.

## 1. Introduction

Understanding large-scale cortical dynamics and their relationship to the underlying anatomy is key to unlocking the computational principles of the brain. Despite cytological differences at the microscopic level, the repetitive design of mesoscopic motifs—circuits, columns and local hierarchies—is a stand-out feature of brain anatomy. By promising a way of scaling up in size, this repetitive principle is likely crucial to the success of the very large-scale computational modelling efforts currently underway. From this view, it is the larger scale structural connectivity and topological constraints of macroscopic anatomy that plays the crucial role in sculpting information processing in the cortex [1,2]. According to this modern connectionist approach, the anatomical network is fundamental because the functionally specialized and integrative roles of cortical regions come from their interactions with other brain regions, not only from their microscopic particularities. Functional network is a widely used approach to study one of the major current paradigms of neuroscience: the relationship between structure and dynamics [3]. While the physical structure remains roughly constant on the time scale of minutes to hours, the dynamic states of cortical regions evolve in a coordinated way to perform a multitude of tasks, which vary on far shorter time scales. Although some states of brain dynamics are very general and common, at least among individuals of a same species, others are very particular, and perhaps unique.

Spanning different neuronal scales, synchronization is an important dynamical feature that can increase the effectiveness of interactions between brain regions, for example by aligning neuronal spikes from different regions into critical windows for spike-time-dependent plasticity [4]. Indeed, an emerging paradigm considers synchronization as *the* key feature that modulates cortical interactions [5]: accordingly, its maximum effectiveness is achieved when two mutually interacting areas are phase coupled. The synchronization strength between regions is thus fundamental, and a functional network is composed of nodes that are linked by highly correlated pairs. The interactions between brain regions, in turn, are not static but rather erratic, depending on their dynamics and past activity [6–8]. Such a dynamic balance of integration and instability likely underlies the need for cortical function to switch according to the changing environmental needs [9–11]. Hence, variable and flexible dynamics, which depends on coupling delay [12] and the balance of excitation and inhibition [13], is a hallmark of adaptive cognitive behaviour.

Understanding dynamic functional networks is not only an important topic, but also a rather complex one. As with any problem endowed with complexity, a wise strategic approach is required. A standard way to reduce the level of complexity is to divide the problem into pieces. Often, however, it is not clear which is the best way to split the problem. Neuronal systems are clearly more than the sum of their parts: extrapolations thus demand careful analysis. The fingerprints of small motifs have already become a standard procedure to understand the structure of complex networks [14–19]. A recent, promising approach to understanding the relationship between structure/dynamics on complex networks consists of studying the structure/dynamics relationship in network motifs [20–23]. This rests on the premise that the dynamical features of network motifs plays a role in shaping the dynamics and synchronization of complex networks [21,22].

The connections between cortical regions in the brain, as derived from the primate cortical network, obey a few general principles that may inform the study of structure/dynamic analyses [24]: first, the brain topology is hierarchical with cortical regions occupying the highest positions in the hierarchy [25–27]. Second, spatial clustering leads to a modular architecture, in which nodes are strongly interconnected with regions within the same cluster (SC)—typically nearby in space—than with regions outside their cluster [28,29]. Third, some nodes can be classified as hubs because they have a degree larger than the average node degree of the network [17,30–32]. These principles have been linked to evolutionary and functional advantages as well as optimization arising from the physical and metabolic constraints under which brains operate on evolutionary time scales [33,34]. These general network properties can also shape the dynamics and synchronization of cortical networks [35–37]. However, a complete understanding of how the structural constraints of the anatomical connectivity influence dynamics is lacking, particularly regarding the mechanisms underlying the variability of synchronization.

Here, we analyse the synchronization patterns as they arise in small motifs, prototypical modular structures and cortical networks. We aim to identify factors that enrich the landscape of synchronized cortical states. The paper is structured as follows: in §2 we describe mesoscopic dynamics as they arise from small motifs of coupled neural masses. In §3 we retrace recent work that highlights the stablizing role of reciprocally coupled pairs (which we call resonance pairs) on the motifs (small subnetworks of three nodes) in which they are embedded. This influence is contrasted with that of closed loops, which introduce frustration into the system. Frustrated motifs hence disrupt stability and synchrony, introducing metastability and variability. In §4–6 we study these two competing influences in small constructed modular networks. Subsequently in §7 we examine structural networks derived from tracing studies of the macaque brain, using the distribution of motifs to link the observed dynamics to those on our prototypical networks.

## 2. Neuronal network dynamics

We employed numerical methods to study neuronal dynamics on large-scale cortical networks. Each node was considered as a conductance-based neural mass model [36,38,39], described in detail in appendix A. This approach reduces the dimensionality of a cortical region in each node to three nonlinear differential equations, which govern the dynamics of the excitatory subpopulation, the inhibitory subpopulation and the fraction of open potassium channels. Analysis of such a reduced model is required to render the problem of large-scale cortical dynamics tractable. Distributed cortical regions are coupled together through long-range excitatory connections with a delay of 10 ms (unless otherwise stated) to account for finite axonal conduction speed [40]. In our mesoscale model, the dynamic variables start from random initial conditions. Exploring the phase space can lead to different basins of attraction and different transient dynamics, causing considerable trial-to-trial variability, consistent with experiments [41–43]. This variability is an intrinsic feature of many large-dimensional systems. Occasionally, such variability can also arise from transitions between metastable states. The following results should not be interpreted as deterministic outcomes of the model, but instead as likely statistical outcomes revealed through a large number of repetitions of the same numerical experiment. Deeper statistical features, in contrast, are remarkably persistent. For instance, nodes oscillate with two dominant frequencies, a slow approximately 10 Hz and a fast approximately 100 Hz time scale, and two mutually coupled nodes typically synchronize in anti-phase synchrony [23].

## 3. Resonance pairs in frustrated motifs promote variable synchronization patterns

We begin our analysis by focusing on three-node motifs, whose enumeration (M1, M2,…, M13) follows the seminal work of Sporns & Kötter [15]. As illustrated in the time traces for the common driving motif (M3) in the weak-coupling regime, robust synchronization between the driven nodes (1 and 3) does not occur (figure 1*a*), nor between directly coupled nodes (figure 1*b*). The lack of synchronization is also evident from the cross-correlation functions for one (figure 1*c*) and for an average of 40 trials (figure 1*d*). By contrast, strong zero-lag synchronization between these nodes (1 and 3) occurs in motifs (such as M6 and M9) that possess a reciprocal connection between at least one pair of nodes (figure 1*e*–*l*). A crucial feature of the dynamics of these synchronized motifs (M6 and M9) is that whereas the pair of edge nodes 1 and 3—which are indirectly connected via node 2—synchronize in-phase, the pairs of neighbouring nodes (1–2 and 2–3) synchronize in anti-phase at the slow rhythm. Even with a coupling delay of 10 ms, owing to the internal dynamics of the neural mass models, the phase difference between the neighbouring nodes amounts to about 60 ms. Hence, mutually connected nodes enhance the synchronization in these small motifs, an effect that we previously showed that can propagate along chains of connected nodes [23].

Importantly, this tendency of synchronized neighbouring nodes to oscillate in anti-phase cannot be satisfied for all motifs: for some configurations, for example when three nodes are mutually connected among themselves (motif M13), this tendency of synchronized neighbouring nodes to oscillate in anti-phase is frustrated. Typically in this motif (figure 1*m*–*p*), two out of the three pairs show anti-phase synchronization: when this happens, the third pair shows in-phase synchronization. For example, nodes 2 and 3 are in-phase for the case illustrated in figure 1*m*,*n*. Hence, this corresponds to an example of frustrated dynamics because although each pair shows a tendency to oscillate in anti-phase, only two out of the three pairs can accomplish this anti-phase configuration at any one time.

Starting from random initial conditions in this frustrated motif (M13), one cannot predict which pair of nodes will synchronize in phase. Hence, the average zero-lag cross-correlation across trials between nodes 1 and 3 (figure 1*i*–*l*), as well as any other pair of nodes (figure 1*p*), is markedly diminished following the incorporation of the mutual connection between nodes 1 and 3. However, the selection of whichever two nodes have the maximum cross-correlation at zero-lag after each trial (and assigning them the labels 1 and 3 *a posteriori*) reveals that the maximum zero-lag synchronization in this frustrated motif can also be strong: we called the auxiliary motif M13 that has nodes 1 and 3 labelled *a posteriori* to maximize the zero-lag synchronization between them (figure 1*q*–*t*).

## 4. Synchronization mode selection by symmetry breaking

The preceding analyses illustrate the effect that the resonance pair has on promoting synchronization and, conversely, how frustration destabilizes the stable state in weakly coupled small motifs. We next studied the robustness of these results with respect to the coupling strength (figure 2*a*). Nodes 1 and 3 do not show in-phase synchrony for the common driving (M3) or frustrated (M13) motifs, but they synchronize even for very weak coupling for motifs which have resonance pairs (M6 and M9). Another motif that also exhibits zero-lag synchronization between nodes 1 and 3 is motif . This shows that for sufficiently strong coupling, one pair of nodes synchronizes in phase in the frustrated motif. Despite the unpredictable character of the particular pair of nodes that synchronize when parameters across all nodes are identical, a small mismatch (such as a 10% reduction in the input current over one neural mass model) breaks the symmetry, converting the distinct node to an apex. For example, a mismatch in node 2 favours phase synchronization between nodes 1 and 3 over the other pairs 1–2 and 2–3 (figure 2*b*).

## 5. Synchronization on modular networks

We now focus our analyses on canonical modular structures comprising two local communities connected through a single hub. In these prototypical modular networks, we set nodes to be mutually connected with every node within each module, and the hub is mutually connected to all other nodes. Such modular structures, even for the smallest cluster size of only two nodes per module (figure 3*a*), contain frustrated motifs within clusters. In addition, this type of modular structure shares another common feature with small motifs: the hub plays the role of the apex node in motif M9, and the clusters mimic the outer nodes (figure 3*a*).

To compare synchronization in modular structures against those in motifs, we focused on the cross-correlation between two nodes belonging either to the SC or to different clusters (DCs). For each trial, the pair of nodes was selected *a posteriori* as the two nodes exhibiting maximum zero-lag cross-correlation. This was then averaged over all trials. The maximum phase synchronization between nodes belonging to DCs resembles the phase synchronization between nodes 1 and 3 in motif M9, whereas the maximum phase synchronization between nodes belonging to the SC resembles motif (figure 3*b*,*c*). This confirms that the motif structure is indeed influential in shaping these modular dynamics.

In the weak-coupling regime, the maximum synchronization between DCs is remarkably similar to the synchronization between nodes 1 and 3 for motif M9 (figure 3*b*). However, for stronger coupling, when the synchronization within clusters grows, the synchronization between DCs decays. As a result, the synchronization between clusters is only stronger than local synchronization for very weak coupling (see black arrow in figure 3*b*).

We next tested the role of the intra-cluster connection strength in shaping the long-distance synchronization between clusters. This was achieved by multiplying the coupling strength of the intra-cluster connection *c* by a weight factor *w* ≤ 1. For this experiment, we used a modular structure with a larger cluster size, as shown in figure 4*a*. As illustrated in figure 4*b*,*c*, the synchronization strength increases for both SC and DC as the weight *w* is reduced. However, although our measures of synchronization improve for very low *w*, the variability of the synchrony pattern is reduced because the frustration is diminished. Hence, in the case of weak intra-cluster connectivity, the system exhibits a strong tendency to evolve to a state of global synchronization in which all cluster nodes are in phase with each other, and in anti-phase synchrony with the hub node.

Setting the weight ratio *w* = 1 and comparing the maximum zero-lag synchronization in modular structures with DC sizes (grey curve of figure 3*b* against red curve of figure 4*b*), we find that the synchronization inside a cluster is stronger for larger cluster sizes. Having the maximum synchronization inside a cluster that grows with the size is an expected result because for larger clusters, the total coupling is stronger by construction (nodes inside a cluster are all-to-all connected), and the number of candidate pairs to synchronize also grows. On the other hand, an unexpected result appears when comparing the synchronization between clusters for DC sizes (blue curve of figure 3*b* against red curve (*w* = 1, dots) of figure 4*c*): in this case, the larger cluster gives rise to a weaker synchronization between clusters. These results indicate that for *w* = 1, there is a competition between SC and DC synchronization, in which the presence of SC synchronization prevents the synchronization between DCs.

## 6. Variability of synchronization patterns

Even frustrated motifs show a small set of stable patterns of synchronization. Owing to the larger number of nodes and the presence of several frustrated motifs within each cluster, modular structures exhibit a larger number of stable patterns of synchronization. By performing multiple trials with random initial conditions, it is possible to estimate the variability of phase synchrony by counting the number of different outcomes that are found in a quasi-stationary regime.

To measure the variability of phase synchrony, we dichotomized pairs as either synchronized or not synchronized for each trial of 2 s (discarding a transient period of 500 ms) depending on whether the cross-correlation coefficient was greater than a threshold *θ*. Figure 5*a* illustrates different patterns of zero-lag synchronization obtained from the network studied in figure 4*a* for different coupling strengths, in which white (black) stripes represent synchronized (unsynchronized) pairs. We then estimated the variability by counting the number of different patterns and dividing by 500, the total number of trials.

The analysis of the phase synchronization—defined as the maximum cross-correlation at any lag (figure 5*b*)—reveals a strong dependence of the synchronization variability on the threshold *θ*. The appearance of a maximum synchronization variability at intermediary coupling strength is evident for sufficiently high threshold. Next, we restricted the synchronization to zero-lag phase synchronization, such that a synchronized pair is identified only if the zero-lag cross-correlation coefficient was above *θ* (figure 5*c*). Hence, we see that the zero-lag synchronization variability is also optimized for an intermediary coupling strength: the synchronization is too sparse for very weak coupling and conversely too stable for strong coupling. Moreover, as the number of nodes, resonance pairs and frustrated motifs grow with size, the maximal variability also grows with the cluster size (figure 5*d*).

## 7. Macaque cortical networks

### (a) Statistical features of macaque structural networks

We first re-visit some basic network features, such as modularity and motif distribution, as expressed in cortical networks of the macaque. These will guide our structure–function analyses. In particular, we study three binary CoCoMac networks: the visual cortex (with *N* = 30 nodes; figure 6*a*), the visual and sensorimotor areas (with *N* = 47 nodes; figure 6*b*) and the larger macaque cortical network (composed of *N* = 71 nodes; figure 6*c*). These are collated cortical networks of the macaque (CoCoMac [44,45]), available at the brain connectivity toolbox [18].

A modular decomposition of these networks [18,46,47] suggests that these three networks have two, four and four modules, respectively. Analyses of the frequency of structural counts of the 13 configurations of three-node motifs (figure 6*d*) for the macaque networks reveal substantial differences with respect to the presence of inter-modular connectivity (figure 6*e*–*g*). Motifs M12 and M13 are typically composed entirely by intra-modular connections (blue bars, bottom). By contrast, motifs M4, M6 and M9 appear more often connecting inter-modular nodes, either by exclusively inter-modular connections (red bars, top) or by a combination of inter- and intra-modular connections (white bars, middle). The participation index (see appendix A), which quantifies the inter- to intra-modular connectivity ratio of each node [17,18,48], shows that many nodes with low degree exhibit strictly intra-modular connections (null participation index; lower left-hand corner of figure 6*h*). Conversely, high-degree nodes typically interconnect modules (upper right-hand corner). Between these extremes, there exists a tremendous variety, such that nodes of intermediate degree may connect within and/or between modules. Interestingly, network motifs also show some degree-dependent statistical patterns in these macaque networks. High-degree nodes, for example, tend to occupy the apex position of motif M9 (figure 6*i*) [17]. This feature suggests that nodes in the apex position of this bi-directionally connected three-node chain may play an influential role on the dynamics of distant segregated areas. Additionally, excluding nodes with very low degree, the prototypical frustrated motif M13 appears to be expressed roughly independently of the node degree (figure 6*j*). This indicates that although this frustrated motif appears more often inside clusters, frustration can be a general feature affecting both hubs as well as peripheral nodes.

### (b) Synchronization patterns on macaque cortex

We finally extend our analyses of the variability of synchronization patterns to these macaque cortical networks. Consistent with the dynamics on our prototypical modular networks (figure 5*b*), the variability of network patterns of zero-lag synchronization for the visual cortical network also exhibits a maximum for a moderate coupling strength (figure 7*a*). While our results have so far been restricted to a delay coupling of 10 ms, it is important to note that the coupling delay also has a strong effect in the synchronization of the network and its variability (figure 7*b*): maximum variability appears for short (10 ms or less) and very long (more than 40 ms) delays, with a minimum for intermediate values (15–20 ms). Also consistent with our prototypical networks, the dynamics on the larger cortical networks exhibit a richer variability of patterns of zero-lag synchronization (figure 7*c*).

## 8. Discussion

Adaptive brain function rests upon the ability of cortical networks to balance the competing constraints of variability and stability. Using the dynamics on network motifs, we investigated the influence of structural connectivity on the diversity of synchronization patterns. We identified resonance and frustration as two major principles in sculpting the variability of the functional networks: reciprocal pairs enhance zero-lag synchronization via resonance-induced synchronization; Frustration decomposes stable synchronized dynamics into multiple competing stable solutions. Together, the two principles appear to give rise to metastable dynamics in systems of different sizes, from motifs to cortical networks.

We propose that in larger networks, these dyadic influences reflect the contrasting deployment of these motifs with respects to their modular organization. Assuming that modules contain larger clustering, motif M13 represents a key connectivity pattern within modules. An important contribution of this frustrated motif is to enrich the dynamical landscape of synchronized cortical states. This property of frustrated systems to enrich and diversify the dynamics is general and occurs in other systems [49–53]. However, to the best of our knowledge, this is the first time it has been related to large-scale cortical dynamics.

The over-representation of motifs M9 and M6 is a consistent and characteristic property of cortical networks [15,17,23]. Taking into account that apex regions typically connect distant regions, in agreement with previous work [54–56], our results suggest that one functional role for these hubs is to promote the synchronization between segregated regions belonging to DCs. According to this proposition, these connections would enrich the variability of synchronization between clusters, which in our model occurs more sparsely than inside clusters.

Changes in brain connectivity have been associated with pathological brain states, such as schizophrenia [57] and multiple sclerosis [58]. Our results show that two crucial features of the structural connectivity that shape the cortical dynamics are the number of resonant pairs and frustrated motifs. Future work could investigate whether these two structural elements (resonant pairs and frustrated motifs) are altered in neurological and psychiatric illnesses characterized by network dysplasia [59].

In addition to changes in structural connectivity, many pathological brain states are also associated with alterations in neuronal synchrony, which may be either reduced (e.g. autism [60] and schizophrenia [61,62]) or increased in levels of synchrony (e.g. epileptic seizures [63] and Parkinson's disease [64]). A balanced amount of synchronization is an essential feature of healthy neuronal dynamics [9–11]. Consistent with experiments and simulations in networks of spiking neurons [13], our numerical analyses of large-scale oscillatory cortical (modular and macaque) networks show that the maximum variability of synchronization patterns occurs for moderate coupling: weaker coupling reduces the variability because fewer pairs can synchronize and the set of possible outcomes is restricted; stronger coupling reduces the variability because the networks often become globally synchronized, and a smaller set of stable states is favoured over others. Such changes might plausibly underlie the reduction in consciousness associated with anaesthesia [65]. This finding of maximal variability for moderate coupling reinforces the importance of homeostatic regulatory mechanisms to maintain the brain close to its optimal region of the parameter space [11]. Our results are also consistent with prior analyses of the role of delays [12,66] and network size in shaping the cortical dynamics.

While resonance-induced synchronization has been shown to hold for synchronized pairs, regardless of their phase relation [23], frustration requires the synchronization between pairs of nodes to be in anti-phase. Anti-phase synchronization is a common phase relationship between separated cortical regions found in a variety of experimental sets [13,67] and models [23,54,68]. The phase relation of the synchronization depends on many factors such as the balance between excitation and inhibition within a cortical region [13], and the coupling strength between regions [68]. More importantly, however, anti-phase becomes the dominant regime when inter-region synapses involve a substantial conduction delay [68]. For this reason, we argue that our model, which favours anti-phase synchronization, is a suitable model for the purpose of understanding the dynamics in large-scale cortical networks, in which long-distance connections are invariably associated with conduction latencies [40]. Extension of our conclusions to other systems should fundamentally depend on whether anti-phase synchronization, the condition for frustration to occur, is the stable solution for two coupled nodes.

We have used a simple way of estimating the complexity of the landscape of cortical activity states. The advantage of this measure is that it is clear and intuitive to count the number of different possible dynamic outcomes. Future work could focus on more sophisticated measures of dynamic variability, based on the system entropy [69], phase reconfiguration [70] and the temporal dwelling in the various stable solutions [71]. Links to empirical measures of dynamic functional connectivity should also be forged [8]. In addition, while we here forge a link between motifs and modules when bridging scales, future work could focus more closely on these related larger scales (of modules) as well as the shaping influence of the largest structural backbone in the primate cortex, namely the rich club [72].

## Acknowledgements

We are grateful to Olaf Sporns, Claudio Mirasso, James Roberts, Andrew Zalesky and Chris Honey for valuable discussions.

## Appendix A

**(a) Neural mass models**

The conductance-based neural mass model we used represents a mesoscale cortical region [38], which is derived from the biophysical Morris–Lecar model [73]. After adaptation to connect areas [74] with synaptic interactions [75], it reached its most recent formulation, which is suitable for modelling whole brain activity in large-scale networks [23,36,39]. Each cortical region was described by the neural mass model of spontaneous cortical dynamics
A 1
A 2
A 3
where *V* is the mean membrane potential of the excitatory pyramidal neurons, *Z* is the mean membrane potential of the inhibitory interneurons and *W* is the average number of open potassium ion channels. The neural-activation functions exhibit a sigmoidal-saturating growth with *V* that governs the fraction of open channels *m*_{ion}
A 4
Finally, and are the average neuronal firing rate of excitatory and inhibitory subpopulations of region *i*. Assuming Gaussian distributions, they are described by the following sigmoidal activation functions:
A 5
A 6

In the above equations, the following set of parameters was used [75]: *g*_{Ca}(=1.1), *r*_{NMDA}(=0.25), *a _{ee}*(=0.4),

*V*

_{Ca}(=1),

*g*

_{Na}(=6.7),

*V*

_{Na}(=0.53),

*g*

_{K}(=2),

*V*

_{K}(=−0.7),

*g*

_{L}(=0.5),

*V*

_{L}(=−0.5),

*a*(=2),

_{ie}*a*(=1),

_{ne}*I*

_{0}(=0.3),

*b*(=0.1),

*a*(=0.4),

_{ni}*a*(=2),

_{ei}*T*

_{Ca}(=−0.01),

*T*

_{Na}(=0.3),

*T*

_{K}(=0),

*δ*

_{Ca}(=0.15),

*δ*

_{Na}(=0.15),

*δ*

_{K}(=0.3),

*ϕ*(=0.7),

*τ*(=1),

_{W}*Q*

_{V}_{max}(=1),

*V*(=0),

_{T}*δ*(=0.65),

_{V}*Q*

_{Z}_{max}(=1),

*Z*(=0) and

_{T}*δ*(=0.65). The above equations also include other varying parameters: the coupling strength between cortical regions

_{Z}*c*, the synaptic delay between cortical regions

*τ*and

*j*(=1,…,

*N*), which are the afferent regions of region

*i*. Consistent with the average values reported in macaques [40], we typically used a constant delay of 10 ms. Although axonal delays between two given regions are typically distributed and specific roles have been attributed to non-homogeneous delays (e.g. to shape the EEG power spectrum [76]), as a first approximation here we consider only homogeneous delays. Small delay mismatches have been shown not to affect the synchronization in neuronal motifs qualitatively [23], but the functions of delay diversity and distributed delays remain to be investigated in large-scale network models. The model was simulated in Matlab (Math Works) using the function

*dde23*.

**(b) Participation index**

The participation index *P* quantifies how connected nodes are with regions belonging to other modules [17,18,48]. It is defined as
A 7
where the summation *m* runs over all identified modules *N _{M}*;

*k*and

_{i}*κ*stand for the degree of node

_{im}*i*and the number of edges from node

*i*to nodes within module

*m*, respectively.

## Footnotes

One contribution of 12 to a Theme Issue ‘Complex network theory and the brain’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.