Structure-Function Discrepancy: Inhomogeneity and Delays in Synchronized Neural Networks

The discrepancy between structural and functional connectivity in neural systems forms the challenge in understanding general brain functioning. To pinpoint a mapping between structure and function, we investigated the effects of (in)homogeneity in coupling structure and delays on synchronization behavior in networks of oscillatory neural masses by deriving the phase dynamics of these generic networks. For homogeneous delays, the structural coupling matrix is largely preserved in the coupling between phases, resulting in clustered stationary phase distributions. Accordingly, we found only a small number of synchronized groups in the network. Distributed delays, by contrast, introduce inhomogeneity in the phase coupling so that clustered stationary phase distributions no longer exist. The effect of distributed delays mimicked that of structural inhomogeneity. Hence, we argue that phase (de-)synchronization patterns caused by inhomogeneous coupling cannot be distinguished from those caused by distributed delays, at least not by the naked eye. The here-derived analytical expression for the effective coupling between phases as a function of structural coupling constitutes a direct relationship between structural and functional connectivity. Structural connectivity constrains synchronizability that may be modified by the delay distribution. This explains why structural and functional connectivity bear much resemblance albeit not a one-to-one correspondence. We illustrate this in the context of resting-state activity, using the anatomical connectivity structure reported by Hagmann and others.


Introduction
Much of the current focus in the empirical study of large-scale neuronal networks has been on their intrinsic activity and the degree to which the coherent patterns of this intrinsic activity reflect anatomy. The use of fMRI and diffusion spectrum imaging has allowed for a comprehensive evaluation of the structurefunction map of resting-state networks (RSNs). In fMRI the spatial patterns of spontaneous changes in blood oxygenation leveldependent signals seem to reflect the generating neural architecture of RSNs. Despite the very slow changes of these signals, Biswal and co-workers [1] defined RSNs as networks of brain areas that exhibit temporally coherent activity in the absence of identifiable externally imposed or measurable events. More recently, RSNs penetrated the field of encephalography [2,3]. For M/EEG, locally synchronized neural activity is considered to yield macroscopic oscillations that provide a basis for defining functional brain networks [4]. In most studies, structural connectivity is considered a good predictor of functional connectivity [5,6]: Structural connectivity agrees with the anatomical connections between network nodes and functional connectivity covers the statistical relationship of nodal activity.
The predictive value of structure for function found support in recent modeling work using full brain systems with realistic anatomy, which demonstrated the structural dependency of functional network configurations [7]. There, functional connectivity has been estimated between all nodes over several hundred seconds of simulated time yielding the pattern of functional connectivity over this time window that largely reproduced the structural connectivity. At smaller time windows, however, shorter-living patterns of functional connectivity emerged that had not been predicted by anatomy. To understand this discrepancy we investigated effects of time delays vis-à-vis effects of structural inhomogeneity on synchronization patterns of neuronal networks.
Delays are inherent in neuronal networks due to finite conduction velocities [8] and synaptic transmission [9]. Ignoring delays may be a valid starting point for mathematical analysis but when doing so one runs the risk of loosing biological plausibility. However, incorporating delays in oscillatory networks does come with immense challenges. Already for low-dimensional oscillatory systems (or for high-dimensional ones with strong symmetry) the presence of delays is known to change the dynamical repertoire significantly [10,11]. Yeung and Strogatz showed for very large networks how time delays can alter synchronization properties, even if the structure is isotropic and homogeneous [12]; see also [13,14]. Numerical assessments revealed similar results for biologically motivated and hence more inhomogeneous connectivities. Delays seem to be crucial in establishing the spatiotemporally organized fluctuations typically observed in resting state brain recordings [15][16][17]. In the present study we sought to tackle this issue and separated the effect of time delays from that of inhomogeneous connectivity by studying networks consisting of distinct neural masses. Neural mass models offer a low-dimensional description of the dynamics of a large neuronal population and exist in a variety of forms [18]. We chose for Freeman's seminal model [19][20][21], since it covers the dynamics of mean membrane potential changes that relate closely to encephalographic signals. A network of such entities may constitute RSNs if we regard the neural masses to be representative of individual brain areas.
Throughout the paper we describe functional connectivity by means of phase synchronization whose dynamics can be estimated in voltage-based and firing-rate models using a combination of rotating wave and slowly varying amplitude approximations, or in brief averaging, see [22,23]. In the Methods section this combination of approximations is briefly summarized for Freeman neural mass models in the oscillatory regime. Central outcome measure is thus the phase dynamics of the individual nodes in the network or, to be more precise, the density of the nodes' phases as a function of time, often also referred as time-dependent population distributions. We note that we applied this approach before to instantaneously coupled Wilson-Cowan firing rate models [24] (see also [25]) but, as said, we here chose for the Freeman model for an easier comparison with M/EEG studies. For coupled Freeman models we could analytically determine the corresponding stationary distributions even in the presence of delays and inhomogeneous coupling between neural masses. We could not only prove the existence of these solutions, but we were also able to determine the loss of stability of the desynchronized state as soon as the overall coupling strength exceeded a critical value. More complicated scenarios including biological plausible anatomical adjacencies were treated numerically to illustrate the non-trivial relationship between structural and functional connectivity.

Results
We considered a set of N coupled neural masses whose mean membrane potentials V k follow the dynamics In this expression a k and b k represent mean rise and decay times of neural responses in population k, q k stands for an external input, and S : ½ denotes a sigmoidal activity function covering the effects of pulsecoupled neurons in populations l~1 . . . N [19,26]. Corresponding mean activities V l arrive at population k after yet arbitrary delays t kl . The structural connectivity matrix C kl served to introduce both excitatory and inhibitory connections in the thus asymmetric coupling; see Fig 1A. We first considered the case in which a large degree of homogeneity was present in C kl to define a 'baseline'. Subsequently we introduced inhomogeneity to mimic, e.g., the sparse connectivity presumably underlying RSNs. Two seminal coupling schemes are sketched in Fig 1. Excitatory and inhibitory populations were always properly balanced to stabilize oscillatory behavior [27,28]. This translates to the condition that at least one pair of the eigenvalues of the linearized system around the fixed points V (0) k ,U (0) k must be imaginary with positive real part. From the Methods section it follows that the considered coupling schemes did satisfy this condition.
For the convolution we chose an exponentially decaying kernel W to represent the dynamics at the synaptic junction. In the particular case of infinitesimal memory, i.e. for c??, we found the corresponding phase dynamics by transforming the neural mass dynamics to polar coordinates around an unstable focus, i.e. V k ?V (0) k zR k cos Vtzw k ð Þ ; R k denotes the amplitude of oscillation and w k its phase corresponding to the central frequency V. Subsequently, we averaged the dynamics over one period 2p=V. We assumed that the characteristic times of amplitude and phase kl . In both cases excitatory/inhibitory-pairs have a symmetric 'internal' coupling with strength +K 1 but the coupling between these pairs differs. In A the between-pair coupling is homogeneous at strength K 2 whilst in B the between-pair coupling may differ across pairs; i.e. it is randomly chosen from a certain distribution U with the constraint that inhibitory units map (on average) with negative and excitatory with positive coupling strength. See text for more details and Figs 7 and 8A for the more complicated coupling schemes employed below. doi:10.1371/journal.pcbi.1003736.g001

Author Summary
Separating the time scale of oscillations from that of the phase dynamics allowed for reducing a network of coupled neural mass models to a system of phase oscillators. We studied the dynamics of networks of phases and their synchronization characteristics as being seminal for functional neural networks. We put particular focus on effects of time delays in the coupling on the network dynamics and contrasted that to effects due to altered structural connectivity. Does neuroanatomical structure prescribe all the macroscopic activity patterns that we observe through electrophysiological brain recordings? We found that heterogeneity in structural coupling and distributed delays have equivalent effects on the shape of phase distributions, i.e., on functional connectivity. The contribution of changes in structural connectivity to network synchronization can therefore not readily be distinguished from that of distributed delays. Interestingly, the emergence of phase clusters in networks requires a subtle interplay between coupling and delays, which may form a window into disentangling structural effects from those induced by delay distributions. Therefore, when investigating neural network behavior, both structural connectivity and delay distribution should be addressed.
dynamics are significantly larger than the period 2p=V, i.e. amplitude and phase dynamics are slow compared to the oscillation. We also assumed the time delays to be of the same order of magnitude or smaller than the period. As a result the delays t kl reduced to mere phase shifts D kl and the phase dynamics became In the dynamics (3) the frequencies obeyed the form and the phase shifts read D kl~p 2 {Vt kl . As said, delays could be transformed to phase shifts due to a time scale separation in the system, i.e. the phase dynamics was set sufficiently slow compared to the oscillation and the delays t kl were up to the same order of magnitude as one period of oscillation (cf. Fig 2). At first glance the dynamics (3) seemed to largely resemble the Kuramoto network of phase oscillators that is well known for its bifurcation scheme from desynchronized to synchronized states. The latter, i.e. the fully synchronized state, would imply large-scale -if not whole-brain -synchronization, which, apart for very pathological cases, is not observed experimentally. Nonetheless we considered linking our study to Kuramoto's profound work to be very insightful, as understanding the Kuramoto network is essential for understanding synchronizability in our more general framework. A closer look revealed that, although similar, (3) differs from the Kuramoto network in two important aspects. First, in the phase dynamics the coupling is given by D kl , which in general is not entirely homogeneous as in Kuramoto's case. This expression for D kl agrees with the previously derived phase dynamics of coupled Wilson-Cowan oscillators in that the amplitude relation R l =R k between nodes l and k affects the corresponding (relative) phase [24]. Second, the finite delay t kl yields a non-trivial phase shift D kl~p 2 {Vt kl . This phase shift might alter the spectrum of phase synchronization entirely; see, e.g., [29,30] for the case. Instead of studying the mere collection of phases fw 1 , . . . ,w N g we investigated the dynamics of their probability density P~P(w 1 , . . . ,w N ,t), because it forms a proper measure of synchronization. To explain: Synchronization around a certain phase value h manifests itself as a peak in the probability density around that value, i.e. a phase cluster around h. For our set of Freeman neural masses we found that the phase distribution follows the dynamics This is the continuity equation, which equals the zero-noise limit of the corresponding Fokker-Planck equation when presuming ergodicity and very large networks (N??). Despite the presence of delays t kl we here succeeded to specify stationary solutions P stat w 1 , . . . ,w N ð Þ . Crucial in finding these solutions was the fact that they can always be written as a mere sum of distinct clustered states, i.e. they always obey the form In order to illustrate stationary solutions of (4) we first specified the coupling matrix C kl to be either homogeneous or inhomogeneous. In both cases we labeled excitatory and inhibitory populations by E and I , respectively, and defined the corresponding sets . As mentioned above we guaranteed throughout analysis that all excitatory populations had an inhibitory counterpart generating and stabilizing oscillatory activity at around the central frequency V [27,28]. For the sake of simplicity we ignored self-coupling of neuronal populations. The stationary phase distribution could thus be written as with within-pair coupling K 1~1 , between-pair coupling strength K 2~0 :2, and distributed delays t kl ; see text for details. The presence of between-pair coupling K 2 and even heterogeneity in coupling, C (hom) kl ?C (inh) kl , did not yield qualitatively different results.

Homogeneous coupling
In the homogeneous case we employed the coupling scheme sketched in Fig 1A. Excitatory and inhibitory populations were fully connected (apart from self-coupling) with coupling values K 1 and K 2 discriminating within-pair and between-pair coupling. For the numerical assessment we always fixed the within-pair coupling to K 1~1 . In more detail, we chose the overall homogeneous coupling matrix as with sub-population connectivities In the absence of delays, i.e. for t kl~0 (D kl~p 2 ) and sufficiently strong coupling K 2 we found robust distributions with M~2 phase clusters, one containing all the excitatory populations and one all the inhibitory ones:  As for strong coupling we found synchronized solutions P stat,2 in the weak coupling case K 2~0 :2. In addition solutions P stat,4 were also present. The number of clusters did depend on the delay value t. That is, altering the delay from t~0 to t~0:2 caused a switch in stability between these two solutions. Since for t~0:05 the two clusters within an E=I -group only showed a small phase difference, we conjecture that t~0:05 is close to the critical value of this bifurcation parameter. doi:10.1371/journal.pcbi.1003736.g004 An example of this solution is illustrated in Fig 3A; the remaining panels in that figure refer to cases of non-vanishing delay that will be summarized below. We note that due to symmetry the homogeneous case with t kl~0 can be readily transformed, proving its resemblance with the Kuramoto network [31]. The somewhat lengthy analytic derivations are given in the Methods section.
Next to the homogeneously synchronized state P stat,2 we found a solution with M~4 phase clusters given by Again we refer to the Methods for the analytical treatment with respect to the existence of this stationary solution; see Figs 3A and 4A for the corresponding numerical assessments. The specific form of the dynamics (3) and (4) already suggested that t kl~0 is just a special case of t kl~t =0. We therefore expected the existence of the solutions above not to be affected by introducing homogeneous, finite delays t kl~t =0, even if this appeared somewhat counterintuitive. In fact, numerics confirmed this expectation as displayed in Figs 3B, 3C, 4B and 4C. The introduction of distributed delays t kl instead of a single t changed results profoundly. To exemplify this, we distributed delays t kl by drawing them at random from a uniform distribution over a certain interval. Recall that according to the transformation from (1) to (3), a distribution of delays t kl generally implies an equivalent distribution of phase shifts D kl . If D kl differed for all populations k and l, the stationary solution P stat (w 1 , . . . ,w N ) of the continuity equation (4) required the presence of many distinct phase clusters. We could prove the existence of that set and,  although the generic solution appeared similar to the homogeneous delay case t kl~t , it did contain M~N centroid values h k instead of the small number M~2 or M~4 shown above. We depict examples of phase distributions for several parameter settings in Fig 5. Interestingly, the heterogeneity in t kl , or equivalently in D kl , agreed with weakening the between-population coupling K 2 in that both cause a profound widening of the phase distributions; compare Fig 5A with 5C. That is, for a network with homogeneous structural connectivity, it is not the presence of delays per se that hinders synchronization but rather the distribution of delays (or the lack of coupling strength).

Inhomogeneous coupling
According to (3) both distributed phase shifts (or delays) and heterogeneous coupling may in principle result in inhomogeneity of phase coupling. In other words, distributed delays t kl and structural heterogeneity may yield inhomogeneity in functional connectivity. We therefore expected a heterogeneous coupling matrix with homogeneous delays t kl~t to be accompanied by desynchronization equivalent to the case of homogeneous coupling and distributed delays. To verify this, we used the inhomogeneous coupling sketched in Fig 1B, which can be given more formally as EE ,r kl ). C: Overall synchronizability Sr kl T. Between-pair coupling strength K 2 appears on the horizontal axis; while blue, red, and black lines correspond to St kl T~0,0:05, and 0:1 respectively. Note that an increase in delay distribution width increased synchronizability due to increased St kl T, but the accompanying higher t kl variance decreased structure-function correspondence as predicted by (3). D-F: Corresponding spatial synchronization matrices r kl for delay distributions t kl~0 , t kl *U(0,0:1) and t kl *U(0,0:2) respectively (K 2~0 :1). Color coding is displayed in the most right panel; only excitatory nodes are shown. doi:10.1371/journal.pcbi.1003736.g008 with sub-population connectivities denoting that the off-(sub)diagonal entries were randomly drawn from a uniform The numerical simulations depicted in Fig 6 confirmed our hypothesis. For C (inh) kl and t kl~t~0 :05 we observed a widening of the phase distribution similar to that shown in Figs 5B and 5C where coupling was established by C (hom) kl but with delay distributions t kl *U(0,0:2) and t kl *U(0,0:05) respectively. Increasing coupling strength reduced the width of the phase distribution comparable to the switch from Fig 5C to 5A. Functional connectivity thus seems to result from an interplay between structural connectivity C kl and delay structure t kl . Therefore both should be taken into account when studying functional connectivity in neuronal networks.

Anatomical coupling
The coupling matrices C kl considered so far were admittedly quite academic. However, these seminal examples did provide important insights that | as we will show here | generalize to more complicated and biologically plausible cases. We performed simulations using the coupling scheme displayed in Fig 7. The matrix C (N) kl had the same structure as C (hom) kl but now the blocks were given by C (N) EE~S C, C (N) II~0 , and C (N) IE~{ C (N) EI~I K1 . The acronym SC stands for 'structural connectivity' that here refers to a neuroanatomical connection matrix as can be derived using DTI/DSI imaging [32,33] and I K1 is the identity matrix with K 1 along the diagonal. To be precise, we used a binary form of the Hagmann connection matrix; see [24] for specifics of preprocessing.
In line with earlier studies we quantified functional connectivity in terms of phase uniformity or phase locking value of the pair-wise relative phases, i.e. r kl~D Se i(w l (t){w k (t)) TD. Using this synchronization measure, simulation results can be best summarized in the form of functional connectivity matrices constructed from r kl values for all available pairs. The effects of delay structure t kl on these functional connectivity matrices are depicted in Figs 8D-F with the underlying structural connectivity given in Fig 8A. The functional connectivity matrix appeared rather sensitive for parameter values, as increasing coupling strength from K 2~0 :1, which was the value used in Fig 8D-F, up to K 2~0 :12 resulted in a fully synchronized network as can be seen in Fig 8C. The sudden synchronization is reminiscent of the phase transition towards the fully synchronized state at the critical coupling strength K c in the Kuramoto network.
In a nutshell, from (4) we could deduce the mechanism responsible for the general finding that structural and functional connectivity are positively correlated [5,6]; see also Fig 8B. Our results clearly show that delay distribution affects both the spatial distribution of functional connectivity (Figs 8B; 8D-F) and the overall level of synchronization in the network (Fig 8C). The increase of overall synchronization is caused by a decreased phase shift D kl by which the phase dynamics (3) converges towards the Kuramoto model, i.e. the delay induces a change in stability of the (partially) synchronized state.

Discussion
We investigated the effect of time delays in the coupling between neural mass dynamics, where we consider an oscillatory regime, established by creating pairs of excitatory/inhibitory neural masses. Although we employed a specific neural mass model, we do consider our results generic because the mappings C kl ?D kl and t kl ?D kl are largely independent of the generating nodal dynamics, presuming that the time scales in the system are sufficiently separated; cf. Methods.
By using this oscillatory dynamics to describe activity in certain brain areas, our approach links directly to the ongoing dispute about changes of functional connectivity in resting state networks (RSNs). There is growing evidence from experimental research that spontaneous brain activity during rest is highly structured into characteristic RSN patterns [1,3,34,35]. These activity patterns seem not to be the result of structural connectivity alone [5,36], but to reflect a non-trivial interplay between the neuroanatomical structure and dynamics [37]. The distribution of time delays involved in this dynamics may have an important role in shaping patterns of activity per se and neuronal synchronization in particular [15][16][17]38].
Key to our analysis was the reduction of a neural mass network to a system of phase oscillators summarized in (3). Several previous studies struggled with computational complexity when trying to unravel effects of delays vis-à-vis coupling on network dynamics [15][16][17]38]. By contrast, our analytic reduction 'readily' allowed for disentangling the contributions of both structural connectivity C kl and delays t kl to the phase dynamics (3). Delays t kl entered the phase dynamics as phase shifts D kl , given a proper time scale separation of oscillatory and phase dynamics. Furthermore, we found that heterogeneity in delays yields effects equivalent to those of heterogeneity in structural connectivitiy. That is, connectivity and delay effects cannot be easily distinguished when solely looking at functional connectivity patterns.
The decrease of V as a function of delay, as depicted in Fig 2, agreed with the analytical findings of [13] as well as with our small-delay approximations outlined in the Methods. However, when further comparing the current results with the literature, one has to realize that some fundamental differences exist between general phase oscillator networks and our dynamics (3). One of those differences is the finite dimensionality of the system (3). We assumed every excitatory/inhibitory neural mass pair to represent a single brain area, by which the dimension of the system under study may be fairly low. On the contrary, most analytical work on phase oscillator networks considered the limit N?? [14,[39][40][41] rendering one-to-one comparisons all but trivial. This can already be appreciated by the rather dramatic finite-size effects in Kuramoto networks [42][43][44][45]. Moreover, our structural connectivity C (hom) kl is rather atypical due to its strong +-asymmetry. Usually, the connectivity structure in similar phase oscillator networks comprises either fully homogeneous coupling or the entries C kl are distributed according to some unimodal distribution [42,46]. An exception is [47], who investigated repulsive coupling, which is similar to the excitatory/inhibitory connections in C (hom) kl . That study reported the presence of two anti-phase clusters reminiscent of the separate E/I -groups observed here.
The time scale separation in (3) and the resulting simplification of delays t kl as phase shifts D kl also hinders direct comparison with studies on delays in the Kuramoto network. This may indeed explain our seemingly contradicting results. For example, we did not observe the emergence of multi-stability mediated by specific (t,K)-combinations as reported in [12,48]. In those studies, regions in the (t,K) phase plane were found in which synchronization was entirely absent. This is clearly not the case in our study. We did find synchronized solutions irrespective of delays. This is not trivial, because for t kl~0 the phase dynamics equation (3) may be reduced to a cosine-variant of the traditional Kuramoto model, which is known for its inability to display synchronized behavior [49]. By exploiting the pairing of the E=I -groups and the +-asymmetry in C (hom) kl we could map our averaged neural mass network to a conventional, fully homogeneous Kuramoto network via a mere transformation of variables. Hence our system can display synchronized behavior even for vanishing t values. This is consistent with [50] who did not find any qualitative effects of a phase shift a on the stability of the Kuramoto network.
Apart from choosing random values as entries of the C kl -matrix, inhomogeneous coupling might also stem from creating distinctively different sub-populations in the network. It can then be studied by modifying within-versus between-network interactions. Particularly interesting in this respect is the occurrence of clustering in the network. When delays are not incorporated one needs either structural inhomogeneity [51] or higher-order Fourier harmonics, in combination with an appropriate phase shift [30], to achieve clustering. The phenomenon of clustering is important in the light of the study of RSNs, i.e. the strong spatiotemporal organization observed in brain activity during resting state conditions [1,3,34,35].
We briefly consider a simple, low-dimensional example of three isolated excitatory nodes. We define a cluster as a number of (excitatory) masses attaining the same centroid phase value h k . By denoting the excitatory and inhibitory centroid values as h HE and h HI , respectively and assuming t k,kzN=2~0 , from the phase derivation (3) and its corresponding continuity equation (4)  . In fact, these forms already hint at the interference between coupling and delays and its effect on synchronization structure. First, all terms on the right-hand side must have equal magnitude requiring specific combinations of D kl and D kl . However, both D kl and D kl are constrained by biology: D kl by the neuroanatomical coupling SC as part of C kl ; and D kl by the spatial structure of the brain, as delays are proportional to distance between masses k and l due to finite conduction velocity. Second, because the left-hand side does not vanish, D kl must have some lower bound. If C kl %1, then sin D kl cannot compensate this and the equality cannot be satisfied. Because C kl determines D kl , there must be some minimal coupling strength between nodes for synchronization to emerge. This explains the positive correlation between structural and functional connectivity, see, e.g., [6,17,38]. It also shows the intricate interplay between structure and delays in establishing synchronization structure. Interestingly, D kl may be regarded as the effective coupling matrix that is typically encountered in dynamic causal modeling approaches [52]. The fact that D kl is directly determined by C kl also explains the finding that models using the structural connectivity as a prior do show more evidence than models using other priors [53]. That is, models that have structural connectivity as a starting point, perform better in terms of data explanation. The sparsity of D kl induced by C kl may yield coexisting synchronized and desynchronized groups within the network, which are often labelled chimera states in the study of phase oscillator systems. It has been found that they crucially depend on the combination of coupling strengths and phase shifts [45,54,55] (or delays [56]), confirming that there has to be a specific matching of coupling and delays for synchronization to occur.
Against the background of the aforementioned t kl -dependence of functional connectivity and the functional-structural connectivity correlation for a biological plausible network, we numerically investigated this by performing simulations of (1) with structural connectivity given by the anatomical connectivity matrix C (N) kl reported by Hagmann and co-workers [32]. Functional connectivity was quantified as pair-wise phase uniformity, i.e. the phase locking value. Our numerical assessment is summarized in Fig 8. It clearly revealed off-diagonal patches with synchronization between nodes that are not coupled (contradicting what has been sketched above). The topology of the Hagmann et al. network shares similarities with the Watts and Strogatz' small-world network [57], i.e. both have a relatively large clustering coefficient with a small average path length. This kind of topology is often believed to be generic in biological neural networks like our brain [5] and enhances synchronizability compared to random networks [57][58][59]. The presence of sparsely connected clusters establishes synchronization between nodes that are only indirectly coupled via their clusters. This causes 'blurring' of the structural connectivity matrix: The functional connectivity matrix is less sparse than the structural one [60]. Although this 'blurring' is similar to the effects attributed to volume conduction [61], in this case it is solely due to network topology.
Next to such clustering phenomena, we can make even more general predictions about the effect of delays in this network. Structural and functional connectivity are most prominently correlated for homogeneous delays, since D kl~D yields an interaction term in (3) that is merely a scaled version of C (N) kl . Hence, the resulting spatial synchronization distribution largely resembles D kl and thus C (N) kl , presuming that the overall coupling strength is not excessively large. This effect can be seen in Fig 8B, where we depicted the Pearson correlation coefficient between the lower triangular parts of C (N) EE and the functional connectivity matrix r kl . Increasing the width of the delay distribution results in a decrease in structure-function correspondence. The positive correlations are consistent with the finding that the pattern of resting state activity is spanned by the eigenmodes of the underlying connectivity matrix in [62]. This is not as trivial as it may seem because the node dynamics in this study were noisedriven fluctuations around a stable fixed point and therefore entirely different from the self-sustained oscillations considered in the current study.
Widening the delay distribution also had another effect: It increased its average value and consequently the mean phase shift SD kl T~S p 2 {Vt kl T tended to vanish. Therefore, the interaction term D kl sin (w l {w k zD kl ) became more similar to an 'ordinary' sine-term, which is known for its capacity to enhance synchronizability [63]. We illustrate this effect in Fig 8C, where overall synchrony SrT is shown as a function of coupling strength K 2 for different average delay-values St kl T. A similar phenomenon has been reported for a system of coupled Hindmarsh-Rose neurons, where a stable synchronized region appears to exist despite the presence of a (constant) delay [64]. Throughout this study we assumed the amplitudes to be constant. The relation between the envelope dynamics of M/ EEG and fMRI-BOLD signals [3] suggests that considering the temporal change of the amplitude may be very important for unravelling the spatio-temporal structure of resting state brain activity. Given that we focused on phase synchronization together with the slow time scale on which the BOLD dynamics evolve (v0:1 Hz, [1]), we believe that the assumption of constant amplitude is justified here. Investigating this assumption in depth, however, is beyond the scope of the current study.
We summarize that the dynamics of a system of coupled Freeman neural masses (1) can be captured by the averaged phase dynamics (3), in which the role of the structural connectivity C kl and delay distribution t kl become explicit. By this, one can identify the relative contributions of structure and delay to phase synchronization, i.e. to the functional connectivity of the neural network. Heterogeneity in structural coupling and distributed delays have equivalent effects on the observed phase distributions. Overall, this supports the notion that structure and delay are both crucial determinants of network behavior and should therefore be taken into account in unison whenever modeling realistic neural networks [37]. Our examples on clustering detailed how the intricate interplay between coupling and delays determines the form and spatial distribution of clustering in these networks. Pinpointing the explicit contributions of t kl and C kl in the phase dynamics (3) enabled us to understand their roles in establishing synchronization structure and why functional and structural connectivity are so closely correlated. This implies that the observed temporal changes in synchronization structure in resting state and task conditions can be modulated through either t kl or amplitudes R k .

Methods
To analyze the gross membrane voltage of a neural mass we first wrote its dynamics (1) as a two-dimensional system using the auxiliary variable U k~_ V V k . For the sake of simplicity we further rewrote the convolution integral in (1) by means of . Then the dynamics (1) reads In the following we discuss this dynamics in its oscillatory regime after transforming the system into polar coordinates to derive the corresponding phase dynamics. That transform, however, was not applied to the original state variables V k ,U k ½ but to the deviations dV k ,dU k ½ around the unstable fixed points in the sense of Taylor and obtained With this expansion the dynamics of the deviations dV k ,dU k could be approximated as A closer look at this linear system revealed that the fixed points were indeed unstable nodes, provided that a proper balance between excitatory and inhibitory masses was present. That is, the eigenvalues of the linear system (8) were complexvalued with positive real-parts -explicit expressions for the eigenvalues (as a function of delay) can be found below. Next, we transformed system (8) into polar coordinates by means of ½dV k ,dU k ? R k cos Vtzw k ð Þ ,{VR k sin Vtzw k ð Þ ½ so that the phases could be cast into the generic dynamics Importantly, in the current case we could assume that the phase dynamics (9) in general contains two (or more) distinct time scales: the rate of change given by the oscillation defined via the frequency V and the time scale(s) of the amplitudes R k and the phases w k ; the latter are much slower than 2p=V and can hence be separated. More formally we used D _ R R k =R k D%V and D _ w w k =w k D%V : In first order approximation we could thus consider R k and w k to be constant during one period of oscillation, T~2p=V. This approach is referred to as a combination of a rotating wave approximation and a slowly varying amplitude approximation [65]. It enabled us to average the dynamics over the interval 0,T ½ Þ; see also [23]. As will be shown in the following, this averaging procedure decoupled amplitude and phase dynamics, which ultimately resulted in the dynamics (3). Below we will provide a more formal discussion regarding these approximations.

Averaging-towards the phase oscillator model
To average the dynamics (9) we defined the averaging operator as Sf s ð ÞT~Ð T 0 f (s)ds. We first substituted (8) into (9) and used The convolution integral on the right-hand side of the second equation of (8) required more attention. Recall the definition of W in (2) When multiplied by dV k (t), this yielded the two averaged trigonometric expressions After substituting this in (9) we obtained the phase dynamics where we defined the constants A kl and B kl as A kl~c 2 sin (Vt kl )zVc cos (Vt kl ) By this procedure we omitted all fast oscillating terms as they averaged out (cf. rotating wave approximation).
As mentioned in the Results section, we focused on the case in which our convolution kernel W did not contain any memory. That is, we considered the limit c??. In this limit only the first terms in the numerators of A kl and B kl remained non-zero and we could cast (11) in the form (3) using Note that in this form, the delays t kl only appeared in the phase shift D kl . Last but not least we simplified expression (3) by exploiting the homogeneity of C (hom) kl to explicitly formulate the D kl matrix multiplication. In particular for equal delays, i.e. for t kl~t , this led to a greatly simplified form of (4) that we summarized in the Results section and will be discussed in more detailed below.

Fixed points and amplitudes
The neural masses do not oscillate around the origin but around the fixed points presuming the fixed points exist. Therefore we were free to choose t kl~0 , such that under the limit c?? the coupling term in (6) reduced to After inserting the form of C (hom) kl we explicitly found That is, in the case of the homogeneous coupling the fixed points of the excitatory masses are equal and the same holds for the inhibitory masses.
The coupling D kl also depended on the amplitudes R k of the neural masses -see also [24]. Accounting for the high degree of homogeneity in the system, we assumed the amplitudes to be equal for equal types of neural masses, i.e. R k~RE , k [ E and R k~RI , k [ I . Furthermore we randomized the parameters a k~a z", b k~b z" by introducing " as a mean-centered random variable. Whenever appropriate we chose " sufficiently small to restrict discussion to the mean values a~1=N P N k~1 a k and b~1=N

Homogeneous coupling -existence of solutions
Since the phase oscillator system (3) can be cast in Kuramoto form, fully synchronized solutions may be stable despite the presence of equal delays t kl~t . But how about solutions other than the fully synchronized ones? In what follows we discuss existence and form of partially synchronized solutions of (3) for general delays t kl . We concentrated on homogeneous coupling and varied the distribution of t kl . In the homogeneous case we found the dynamics of the k-th node's phase distribution P k to be where the subscript 0~E when k [ E and 0~I when k [ I . Note ð12Þ that equation (12) may differ for every k. Homogeneity of C (hom) kl enable us to express D kl explicitly and to define the following constants Sufficient for the existence of a stationary solution is the case in which the drift coefficient in the dynamics of the probability density vanishes, here the bracketed term on the right-hand side of (12). From the dynamics (12) it readily follows that the phase distribution obeys the form (5), i.e. P~P stat (w 1 , . . . ,w N ), containing N different centroid phase values h k .
Equal delays. t kl~t . As said, for homogeneous coupling the sums in (12) could be evaluated in the form of the constants in (13). Furthermore in the case of equal delays t kl~t (i.e. D kl~D ) we only had to account for two distinct populations each with just a single type of density dynamics, by which the system (12) of N equations could be reduced to merely two distinct ones: To show the existence of P stat,2 (w E ,w I ) consisting of two phase clusters, one for the excitatory units and one for the inhibitory units, i.e.
we followed a constructive approach and substituted P stat,2 (w E ,w I ) into (14). Vanishing of the drift coefficient required when abbreviating y~h E {h I . Using DC IE zC EI DwDC EE zC II D, which followed from (13) irrespective of the value of K 2 , we could conclude that a solution y satisfying (15) exists. Note that only y appeared in (15) and not the centroid values h 0 , which allows for the mapping to the conventional Kuramoto model as will be discussed below.
We could readily generalize this line of reasoning to an arbitrary number of clusters after defining the stationary phase probability density containing M clusters as We defined generalized phase differences y ab,kl~ha,k {h b,l where a,b [ fE,Ig and k,l [ f1, . . . ,Mg, by which we obtained the set of constraining equations as In contrast to (15) we here considered general delays t kl , but note that the form of (16) did not change for t kl~t . As the values of the constants C 0 still fulfilled the inequality DC IE zC EI DwDC EE zC II D we again could conclude that this system can be solved. Thus, even for equal delays t kl~t (or D kl~D ) the stationary phase distribution may contain any number of clusters MƒN. Of course, this existence proof does not imply that all these solutions will be found in reality since we have not yet addressed their stability; cf. Figs 4A and 3A. Note that apart from the fully synchronized state, we here restricted our stability analysis to numerical evaluations.
The distributed delay case (16) does not require D kl~D . Hence, for general D kl the stationary distribution P stat,M will contain Stability for t kl~t -a fling with the Kuramoto network. Before proving the existence of distinct solutions in the case of arbitrary delays t kl , we first briefly discuss the case of constant delay because it readily links the model (3) to the Kuramoto network of coupled phase oscillators [63]. For t kl~t~0 , the interaction term in (3) reduces to a mere cosine term because of D kl~p 2 . We note that in Kuramoto's traditional form this cosine term leads to a system that is unable to show synchronized solutions [49]. Our numerical simulations, however, revealed synchronized solutions; see, e.g., Figs 3 and 4. In what follows we show that due to the fact that the centroid phase values h k are irrelevant one can distill the traditional Kuramoto model. It is the coupling matrix C (hom) kl that made the difference with Kuramoto's network as we assumed a (balanced) combination of excitatory and inhibitory units. In the given form illustrated in Fig 1A, the matrix C (hom) kl is strictly speaking not homogeneous but contains an asymmetry. This, however, did allow for mapping our model to the Kuramoto network by adapting an approach of Frank and co-workers [31]. In essence, we applied a change in variables yielding a fully homogeneous coupling matrix. For the sake of legibility we set R E~RI , S' Then D kl simplified to D kl~+ ab V where the '+'-sign discriminates rows that correspond to the excitatory/inhibitory neural masses E=I , respectively. By this simplification the key feature of the system (9), the asymmetry in the coupling matrix, remained untouched. We further transformed the system into a rotating frame at frequency v~a b{V 2 2V by defining [65,66] ð14Þ After abbreviating the constants ab In line with Kuramoto's analysis we introduced a mean (or cluster) phase h and amplitude r in terms of Here we dropped the subscript k because the equation applies for all nodes. This dynamics agrees entirely with the Kuramoto model, apart from the fact that we here considered a k~a and b k~b by which all oscillators' natural frequencies v k~v agree. Due to this correspondence to the Kuramoto model, we could conclude that the system with equal delays t kl~t and 'homogeneous' coupling matrix C (hom) kl can generate synchronized states provided the overall coupling strength ab V K 2 S' V (0) Â Ã is properly chosen.

Numerical assessments
Both distributed delays t kl and heterogeneous coupling called for numerical assessments, particularly when it came to the stability of solutions. We performed numerical simulations of the coupled neural masses (1) using a conventional Euler-forward integration scheme with time step Dt~10 {3 s over a time span of 10 s. We verified the appropriateness of this simple implementation against a more elaborate predictor/corrector integrator [67], which revealed little to no difference but demanded far more numerical resources. The simulated network consisted of 500 nodes (250 E=I pairs) with a k and b k being randomly drawn from uniform distributions (a k [ ½25,75 and b k [ ½175,225) to mimic distributed natural frequencies per E=I -pair of nodes. Although randomly drawn, these sets were fixed across simulations trials. Initial conditions where chosen randomly and did vary between trials; V k (t~0) [ ½{100,100,k [ E, V k (t~0) [ ½{10,60,k [ I and similarly for _ V V k (t~0). The external input q k was set to q k~2 0, k [ E, q k~0 , k [ I . Coupling between masses was achieved by using the sigmoidal activation function S½V that was given as The fraction in this equation may be interpreted as the cumulative distribution function of the normal distribution N V {V th ,s 2 À Á of the firing thresholds V th across the population, whereas the constant c is just a scaling factor [68]. In the simulations we used the following values for the activation function parameters: c~250, V th~1 5 and s~1.
In order to compare the numerics with our analytical results, we determined the phase values w k from the simulated potentials V k (t) by means of the Hilbert phase. To this end, we first determined oscillation frequency V as the lowest frequency with a coinciding peak in the power spectra for all nodes. Voltage traces V k (t) were band-pass filtered using a 1-st order bi-directional Butterworth filter in the frequency band ½0:8V,1:2V. For each sample in the interval t [ ½7,8, phase values w k (t) were then calculated as the angle of the analytical signal. By restricting analysis to that interval we avoided transient behavior as well as possible filter artifacts. The sodetermined w k contained the frequency component Vt, which we first subtracted to obtain Q k . Then, we opted to compute phase distributions over the time interval T and over successive trials. The (circular) mean phase, however, differed from trial to trial because of the randomly chosen initial conditions (see above). Hence for every trial we shifted phases by the mean phase of the excitatory population prior to concatenating trials. By this, the phase distributions of the excitatory phases always became centered around zero and the inhibitory phases were considered relative values. The mean phase per trial was given as where arctan 0 denotes the quadrant-corrected inverse tangents. The distributions displayed in the figures are the phase distributions obtained from 100 simulated trials with different initial conditions. As said, parameter values were identical across trials. For the simulations involving C (N) kl connectivity, we used a binary form of the 66 areas parcelated Hagmann et al. matrix [24,32] as C (N) EE block; see Fig 7 and 8A for the coupling scheme and C (N) EE , respectively. Functional connectivity was quantified as phase coherence given by r kl~D Se i(w k (t){w l (t)) TD. We performed one hundred simulation runs of 10 s for each t kl distribution with different initial conditions and a k ,b k and averaged the r kl matrices over these runs. This was done to avoid high r kl values due to common oscillation frequency alone. For each run, data in the interval t [ ½4,8 was used to determine r kl . The overall coupling strengths were set to K 1~1 and K 2~0 :1 for the r kl matrices displayed in Fig 8. Structure-function correspondence r C (N) EE ,r kl was quantified as the Pearson correlation coefficient between the lower triangular parts of both matrices to avoid spurious correlation values due to common terms along the diagonal.

Eigenvalues of the linear system
To estimate the oscillatory regime of the system of coupled Freeman neural masses (1) we considered the linearized dynamics (8) for c? ? which in general reads Assuming all delays t kl to be (very) small we expanded dV l t{t kl ð Þ in the sense of Taylor and approximated up to the linear order in t kl : For the sake of legibility we here considered a single isolated E=Ipair with C EE~CII~0 and C IE~{ C EI~K1 . We also assumed equal delays, i.e. t kl~t . Then, we found the resulting linear dynamics as _ u ũ u u~Aũ u where we abbreviatedũ ud V E ,dU E ,dV I ,dU I ð Þ T and The matrix A came with eigenvalues l j where we abbreviated y~ab A necessary condition for the existence of a stable limit cycle, and hence for the system (1) to display oscillatory behavior, is that the fixed point V (0) E ,U (0) E ,V (0) I ,V (0) I is unstable. This means that for at least one of the conjugate pairs l j , <(l j )w0 must hold, which was indeed the case irrespective of t. The corresponding =(l j ) then provided a rough estimate for the frequency V as a function of t, as shown in Fig 2 ( Further, for the real-part to be positive the coupling constant K 1 had to be sufficiently large and the intrinsic damping a and/or b sufficiently small but finite, because of < l 1,2 f g t~0 w0 u K 1 w azb where we note that < l 3,4 f gv0.

Distinct time scales
For the separation of time scales underlying all the major approximations in the current study we considered the case in which two distinct time scales are present in the system of coupled neural masses: the oscillation described by the (mean) frequency V and its corresponding period 2p=V -from here-on referred to as fast time scale -as well as a slower time scale, on which the dynamics of w k evolve. In what follows we will verify the expression for _ w w k in (9) and show how the separation of time scales enabled us to determine the role of t kl in the convolution W t{t kl ð ÞÃS V l t ð Þ ½ . As is conventional in multiple-scaling approaches, we set the time t as 'fast' time and the 'slow' time as g~"t with 0ƒ"%1. For the sake of legibility we here adopted the dot-notation for temporal derivatives d=dt and further abbreviated partial derivatives as L L( : )~L ( : ) .
We denoted the deviation of the voltage from its fixed-point as dV k (t)~R k cos (Vtzw k ) where we assumed that w evolved on the slow time g, i.e. w k~wk (g). Note that an equivalent approach can be adopted for the amplitude dynamics, i.e. R k~Rk (g), which is referred to as the slowly varying amplitude approximation. As we were primarily interested in the w k dynamics, we regarded amplitude R k as constant on both time scales; see Discussion section. By this we could readily apply the chain rule and obtained d _ V V k~Lt dV k z"L g dV k :~{ Vz"L g w k À Á R k sin Vtzw k ð ÞdU k where the last equality follows from (8). For the derivative d _ U U k we found d _ U U k~{ Vz"L g w k À Á 2 R k cos (Vtzw k ){" 2 L 2 g w k R k sin Vtzw k ð Þ: Next we considered the expression (9) for _ w w k , where we emphasize that _ w w k could here be identified as L g w k , since w k evolves only on the slow time scale. To anticipate: (9) discarded all terms evolving on the fast time scale, i.e. all O(1) terms in favor of O(") terms and higher. To show this, recall that the right-hand side of (9) read : In words, only expressions evolving on the slow time scale were retained, i.e. only O(") and O(" 2 )-order terms. When focusing on the slow time sale O(") and discarding the even slower time scale O(" 2 ), we could conclude that, up to a constant, _ w w k is given by (9); note that we here applied the so-called two-timing method; see, e.g., [23,69,70].
As said, we used the constancy of w k (on the fast time scale t) to evaluate the convolution term W (t{t kl ) Ã S V (0) l h i . We exploited the description in two time scales to justify the transformation of the delay t kl into phase shifts Vt kl . We explicitly evaluated the integral Ð t{t kl { ? e cs cos (Vszw l )ds to show that (10) is its O(1) result. For the sake of readability we dropped the explicit time dependence of w l whenever possible.
First, by integrating by parts twice we obtained