Phase-lags in large scale brain synchronization: Methodological considerations and in-silico analysis

Architecture of phase relationships among neural oscillations is central for their functional significance but has remained theoretically poorly understood. We use phenomenological model of delay-coupled oscillators with increasing degree of topological complexity to identify underlying principles by which the spatio-temporal structure of the brain governs the phase lags between oscillatory activity at distant regions. Phase relations and their regions of stability are derived and numerically confirmed for two oscillators and for networks with randomly distributed or clustered bimodal delays, as a first approximation for the brain structural connectivity. Besides in-phase, clustered delays can induce anti-phase synchronization for certain frequencies, while the sign of the lags is determined by the natural frequencies and by the inhomogeneous network interactions. For in-phase synchronization faster oscillators always phase lead, while stronger connected nodes lag behind the weaker during frequency depression, which consistently arises for in-silico results. If nodes are in anti-phase regime, then a distance π is added to the in-phase trends. The statistics of the phases is calculated from the phase locking values (PLV), as in many empirical studies, and we scrutinize the method’s impact. The choice of surrogates do not affects the mean of the observed phase lags, but higher significance levels that are generated by some surrogates, cause decreased variance and might fail to detect the generally weaker coherence of the interhemispheric links. These links are also affected by the non-stationary and intermittent synchronization, which causes multimodal phase lags that can be misleading if averaged. Taken together, the results describe quantitatively the impact of the spatio-temporal connectivity of the brain to the synchronization patterns between brain regions, and to uncover mechanisms through which the spatio-temporal structure of the brain renders phases to be distributed around 0 and π. Trial registration: South African Clinical Trials Register: http://www.sanctr.gov.za/SAClinicalbrnbspTrials/tabid/169/Default.aspx, then link to respiratory tract then link to tuberculosis, pulmonary; and TASK Applied Sciences Clinical Trials, AP-TB-201-16 (ALOPEXX): https://task.org.za/clinical-trials/.


Introduction
Many processes in nature are oscillatory, from heart beats and birds flapping their wings, to firing of neurons [1] and brain rhythms [2]. Oscillators are rarely isolated and they interact when coexisting in the same environment, thus synchronizing by adjusting their rhythms [3]. Synchronization, or consistent phase relationships, of distant regions of the brain has been detected by a variety of measures and may be a key mechanism for the regulation of cortical processing and communication [4,5]. Advances of non-invasive structural brain imaging [6,7] have made feasible large-scale network modeling approaches using biologically realistic connectivity, defined by the connection topology and delays, the so-called connectome, which is a crucial determinant of the network behavior [8][9][10][11][12][13]. The Kuramoto model (KM) [14] as a paradigm for the emergent group dynamics of coupled oscillatory subsystems [15,16] is well suited for assessing how the connectome governs the brain oscillatory dynamics [17][18][19][20][21][22], which is then reflected in the phase relationship between brain regions.
Studies of networks dynamics predominantly focus on the synchronization properties, while the actual phase relationship between the oscillators is typically ignored, especially in complex networks with delays [16,23]. For the tractable case of all-to-all equally coupled phase oscillators in thermodynamic limit, the phases of each oscillator are either constantly shifted from the mean phase, or non-uniformly rotate with a speed dependent on their natural frequencies, while still preserving the overall stationary distribution [14]. For heterogeneous couplings, phases become multimodal [24] and thus imply multimodal phase shifts for stationary synchronization, and for couplings of mixed signs [25], the oscillators generally split at distance π, but for strong coupling they form a traveling wave. Glassy states with ordered, but uniformly distributed phases for each frequency also appear for distributed parameters [26], and for structured networks multiple mean fields appear with oscillating, bounded or unbounded phase differences between them [27].
A number of computational studies on brain functional connectivity use neural masses connected with connectome-defined delays and weights. This yields intermittent in-or antiphase synchronization [8] and a good agreement with experimental studies of phase relationship between local node dynamics and their degree in healthy subjects [28], and for Alzheimer's disease [29]. Time-delays have been employed as a necessary condition for modeling anti-phase spatio-temporal patterns in the brain [11,30], and for pair-wise coherence in connectome networks of phase oscillators that reproduce resting state patterns in BOLD fMRI [19], MEG [31] and EEG [21,32,33].
Oscillatory processes are particularly sensitive to delays, because shifts in phasing may render excitatory connections to inhibitory, and vice versa. The impact of time-delays on the synchronization, and indirectly to the phase-lags, has been studied for a single delay [34,35] or for homogeneously distributed delays [36,37], but so far only for all-to-all connectivities or for simple motifs. Spatially heterogeneous delays are particularly important for number of systems, foremost the brain [10,11,13]. Depending on their distribution, time-delays impose phase-shifted, in-or anti-phase clusters of oscillators, but their impact for phase-lags in largescale neural synchronization has not been properly investigated. Stronger connected network nodes have been demonstrated to lag behind the weaker for randomly distributed delays shorter than a quarter period of the oscillators [22], but this restricts a large portion of the relevant frequencies.
In the current study we identify the relationship of the brain topology and its spatio-temporal structure, with the phase lags between the brain regions at any frequency of the brain processes. Analytical insights of synchronization on networks with distributed delays and heterogeneous couplings and frequencies, are applied to in-silico large-scale brain dynamics. Phase lockings and lags are studied in consideration to the limitations of time-series analysis that depend on the regime and levels of coherence. Inhomogeneous interactions due to the connectome are shown to drive the phase relationship, whilst the regimes of synchronization are constrained by the organization of the time-delays. Besides in-phase, these include antiphase locking that for weak coherence depends on the length of the links, while for strong coupling is prevalent for the nodes of opposite hemispheres.

Model
Spatio-temporal organization of delay-coupled oscillatory networks is often studied via phase oscillators that arise for weak interactions [14,[38][39][40]. The delays are reduced to phase shifts when they are small [38,39,41], but they appear inside the state variables when they are of the order of 1/coupling-strength [38,39], yielding _ y i ðtÞ ¼ o i þ h i ðy 1 ðt À t i;1 Þ; y 2 ðt À t i;2 Þ . . . y n ðt À t i;N Þ; K i;1 ; K i;2 . . . K i;N Þ; i ¼ 1 . . . N; ð1Þ where ω i are the natural frequencies, and for each link K ij and τ ij are coupling strengths and time-delays. Phase models still exhibit rich behavior and a direct link to more complex biophysical models, while admitting analytic approaches [40][41][42][43]. A special case is the Kuramoto model, which keeps only the first sine term of the Fourier series of h i and that has been worked out to allow full analytical tractability. A large class of oscillators that are near an Andronov-Hopf bifurcation can be exactly transformed to the KM [14,44], as well as Wilson-Cowan networks [41,45]. Although the KM is not explicitly a brain model, it has been applied as one [18,19,21,[31][32][33] since it is perfectly suited to describe large-scale network synchronization. The utilized model therefore reads where network interactions are symmetric with K ij = K ji and τ ij = τ ji . The aim of the model is to show correspondence to the empirical results for phase difference. These are often captured by phase locking values [46], which are a statistical measure for similarity between phases of two signals. For 1: 1 synchronization, which is of main interest in the empirical studies, PLV is defined as where the phase difference Δθ ij (p) = θ i (p) − θ j (p) is calculated at times p = 1. . .M.

Numerical analysis and statistics
Numerical integrations utilize a Heun scheme adapted to time delays. The time-step is set as 0.01/(max([max(K), 0.05μ, D, 1])), with noise intensity D and mean of the natural frequencies μ, thus assuring that it is never larger than 0.01s and it is accordingly decreased for larger couplings, frequencies or noise. All time series are down-sampled to twice the Nyquist frequency of the fastest oscillator. Complex phase locking values, Eq (3), are calculated at sliding windows of length equal to 10 periods of the mean entrainment frequency and with 75% overlap. Qualitatively similar results are obtained for windows lengths between 5 and 10 periods, and for overlapping between 50% and 90%, although longer windows yield systematically lower level for statistical significance [47].
Signals can often be coherent just by chance and statistical testings are necessary to correctly identify the coherence due to the mutual interactions [46,47]. The level of significance for PLV is calculated as the 95th percentile of maximum values in 100 surrogate signals using two different procedures, followed by the same processing as for the original time-series. The first surrogates, which yield less strict level of significance are obtained by shuffling the time series of the phases of one of the oscillators. The second, generally stricter level is obtained by two independent uncoupled oscillators with the same frequencies, fixed or time-varying, as the original oscillators, with the same level of noise. The problem with the latter is that in empirical analysis the parameters of the uncoupled oscillators and the noise intensity can be unknown, although the noise intensity could be obtained from the variance of the signal under assumptions of stationarity.

Results
We first analyze the case with two phase oscillators with time-varying coupling strengths and natural frequencies [48]. This is valid if network interactions are either unknown, or too weak to cause synchronization, and hence assumed to be encompassed in the inherent dynamics of each oscillator, as additive noise and/or non-autonomous (NA) forcing. Then the analysis is extended to homogeneous and delay-imposed networks with bimodal δ time delays and lognormally distributed node strengths, as observed in the brain. Analytical findings are numerically validated and are used to explain phase statistics for simulated dynamics over the human connectome.

Two oscillators
The simplest case of two delay-coupled phase oscillators with constant parameters reads _ y 1;2 ðtÞ ¼ o 1;2 À K sin½y 1;2 ðtÞ À y 2;1 ðt À tÞ: Steady synchronization occurs when the oscillators start oscillating with a same adjusted frequency O, preserving a constant phase shift ϕ 1,2 = θ 1 − θ 2 , which becomes where O is described by a transcendental function and the critical coupling reads Oscillators can be locked in-or anti-phase, depending on the sign of K cos Oτ, so that for the phase shift it holds The model is made more realistic by allowing deterministic variability of the frequencies and the coupling, and additive, independent, Gaussian noise _ y 1;2 ðtÞ ¼ o 1;2 þ 1;2 sinô 1;2 t À ðK þ K sinô K tÞ sin½y 1;2 ðtÞ À y 2;1 ðt À tÞ þ Z 1;2 ðtÞ: ð8Þ Here hη i (t)i = 0 and hη 1 (t)η 2 (t 0 )i = 2Dδ(t − t 0 )δ 1,2 , with h Á i denoting time-average operator, while ω 1,2 and K are harmonically modulated. In adiabatic limit without noise [48] effective coupling K eff ðtÞ ¼ ðK þ K sinô K tÞcosOt and frequencies o eff 1;2 ðtÞ ¼ o 1;2 þ 1;2 sinô 1;2 t and Δω eff1,2 (t) = ω eff1 (t) − ω eff2 (t), can quantify the synchronization instead of fixed parameters in Eqs (5) and (6), but they give insight into the level of coherence even for stochastic dynamics, and for non-adiabatic response that occurs due to the large inherent time-scale close to incoherence. The stochastic dynamics with constant parameters is shown through the evolution of instantaneous and time-averaged phase lags, and PLVs in Fig 1. The oscillators in panel (A) are identical and the only variability is due to the noise, which causes time-varying cPLV and phases. Nevertheless, phase lags are close to zero during the periods of significant PLV, as seen by the red and magenta lines for a shorter interval in (a), and for the whole time series in (b). Thus, for a sufficient number of time-points the lags for in-phase synchronization will have a mean at 0, as seen in their histograms and estimated probability density distributions (PDF), (d, e). The mean phase difference in (B) is in the interval (−π/2, 0), as predicted by Eq (7) for in-phase locking and different natural frequencies with ω 1 < ω 2 . The results also indicate that statistical significance has no influence on the mean of the observed lags, but only impacts their variance. Hence, lower significance levels would improve the statistics of short time-series, by increasing the number of significant data points.
Non-autonomicity causes intermittent epochs of in-or anti-phase synchronization, Fig 2. These are still well captured by |2K eff | ≷ |Δω eff |, as predicted by Eq (6) for fixed parameters, with periods of insignificant coherence, (c, j), corresponding to |2K eff | ≲ |Δω eff |, (d, k). In both examples the coupling is explicitly modulated withô K , but also implicitly through the NA frequency of synchronization included in cos O(t)τ, whilst hω eff2 i is clearly larger than hω eff1 i for averaging over the times of significant coherence. The latter ensures the phase lags to be in the ranges predicted by Eqs (5) and (7) for ω 2 ! ω 1 , although they are wider distributed than for fixed parameters, due to the varying frequency mismatch. The distribution of Δθ is additionally broader due to the noise-induced variability, which gets partially averaged out for ϕ. As in Fig 1, instantaneous and averaged phases from cPLV have very similar statistics, shown through histograms for significant phase lags in regard to the both levels of significance, calculated at 50 equally spaced bins in the interval [−π, π], as it is the case in the later figures. This implies that shorter time-windows would affect observed phase shifts only indirectly, through the levels of PLV.
Results for time-varying parameters, Eq (8), in Fig 2 confirm that the theoretical insights, Eqs (5) and (7), can be also used to describe the statistics for noisy and NA parameters, which resemble the intermittent coherence observed in the real data. The distribution of phase lags depends on the time-delay compared with the frequency of synchronization, and the average ratio of the natural frequencies. The former dictates the regime of synchronization, in-or anti-phase, whereas the latter specifies in which quadrant the mean of the phases will be located.

Networks of oscillators
First we derive analytical results for two spatial configurations of time-delays in all-to-all connected oscillators with heterogeneous natural frequencies and coupling strengths. Then we Phase-lags in large scale brain synchronization generalize and numerically validate those results for a more biologically plausible scenario with stochastic inhomogeneities.
The system Eq 2 cannot be solved for general [K ij , τ ij ], such as the connectome. Still, based on certain assumptions, we characterize phase relations between different nodes, depending on their location and strength. Firstly we approximate coupling inhomogeneity by the average coupling strength of each oscillator [49], which also allows for sparse networks. The model Eq 2 is henceforth reduced to Next, global and local order parameters [13,36], are defined Here r is the global coherence or the strength of the instantaneous mean field, R i is the local coherence or the mean field strength felt by each oscillator, whilst F and C i are the phases of the global and the local mean-fields. Introducing these in Eq 9, the mean-field character of the model emerges To facilitate the analysis steady partial synchronization [24] is assumed, as opposed to the so-called standing waves [50]. We build on [13] and we derive analytical results for randomly distributed bimodal-δ delays and delay-imposed symmetrical biclusters. The PDF of the time delays with equal peaks hence reads The delays are either spatially homogeneous with the same independent probability for any link, or they are heterogeneously organized so that two identical subpopulations emerge with same internal and external time-delays, Fig 3. Besides representing distinct phenomenological structures, these topologies are motivated from the connectome. Its simplest decomposition Phase-lags in large scale brain synchronization on a left and a right hemispheres identifies the peaks in the delays distribution as intra-and inter-hemispheric links (see Fig 9(d)-9(f)), leading to the clustered organization as a first approximation. However, this division is not strict, and many links are randomly distributed, corresponding to spatial homogeneity. Due to the spatial homogeneity of the random network and of the internal links of ordered subpopulations, Fig 3, the global order parameter, Eq (10), in both cases is For the former z I, II = z can represent any proportion of nodes, while for the latter they correspond to the different delay-imposed subpopulations. Bimodal-δ spatially homogeneously distributed delays. Because of the spatial homogeneity, steady synchronization for random delays implies where the initial phase of the mean field is set to 0 without losing generality. This allows relative phases to be introduced, ϕ i = θ i − Ot, and considering that contribution of links with different delays can be separated for each oscillator [13], Eq (12) is rewritten as > 0: For steady synchronization fixed point appears in Eq (16), and considering its stability, relative phases read Thus, OΔτ and Ot are always in the same complex half-plane (left or right), as also shown by the examples in Figs 5 and 6, and the limits of synchronization are defined as Taking into account the signs of cos OΔτ and arcsin(x), Fig 4(a), relative phases read Moreover, considering the changing of arcsin(x) in those ranges, Fig 4(a), the following rules appear for the phases depending on the nodes strengths and frequencies where up/down pointing arrows indicate to the increase/decrease of the preceding parameter. Note that for any angles Oτ 1,2 , the values for OΔτ and Ot, which appear due to the summation sin(ϕ i + Oτ 1 ) + sin(ϕ i + Oτ 2 ), can be represented on the opposite side of the imaginary axis. This can be shown by adding 2π to any of the angles, so that the result of the sum is not changed, while it can be transformed to read Hence, none of the above results will change if both, sin OΔτ and cos OΔτ, get opposite signs by rotating the angles OΔτ and OΔτ by π, Fig 4(a). In addition, since for stable solutions OΔτ and OΔτ need to be on the same side of the imaginary axis, only cases when they are in the right half-plane can be considered without any loss of generality.
Numerical results in Fig 5 are for constant couplings, K i = K, and Lorentzian distributed natural frequencies with mean μ and scale γ gðoÞ ¼ g=p=½ðo À mÞ 2 þ g 2 : Regardless of angles Ot and OΔτ being in the right or left complex half-planes, panels (A) and (B), phases are as predicted by Eqs (19) and (20) with slower oscillators being closer to the mean field. Time-evolution of entrained and two closest to them unsynchronized oscillators (a, b, e, f), as well as the position of the phases in the complex plane, (c, g), are shown for the rotating reference frame O.
The same populations are simulated in Fig 6, but with equal frequencies and Gaussian noise with a same heterogeneity for the critical coupling for instantaneous interactions [15], i.e. D = γ. Node's strengths, with mean and variance given as :They are assigned to the outgoing links of each node i, whereas the symmetric connectivity is enforced by setting The stochastic frequency heterogeneity averages out over time, hence for the time-averaged phase shifts of the synchronized oscillators we approximately derive Phase-lags in large scale brain synchronization The stability condition is the same as for the deterministic case, while the criterion for synchronization is K i r cos OΔτ > |μ − O|. Predictions from Eq (20)  Time-delays imply multistable solutions for the level and frequency of synchronization, often with no analytical solutions [13,34,35]. For all-to-all equally coupled oscillators, lowdimensional evolution of the dynamics [13] reads From here the frequency of steady synchronization is given as Even though this expression is exact only for Lorentzian natural frequencies, it can be used to determine whether the frequency of synchronization is larger or smaller than the natural frequencies, i.e. O ≷ μ. This inequality depends solely on the signs of sin Ot and cos OΔτ, Eq (23), and governs the general relation of the phase shifts, Eq (21). Hence, taking also into account the stability of ϕ i , Eqs (17) and (21), the behavior of arcsin(x) , Fig 4(a), and the transformation of sine summation, Fig 4(b), we obtain the directions of change of the relative phases as nodes strength increases, Table 1. All combinations of delays are considered, because phases Oτ 1,2 might still be longer than a full cycle, despite internal delays being set shorter than the external. Phase-lags in large scale brain synchronization Results from Table 1 can be summarized by the geometrical mean of the time-delays e iOt , which can be always transformed to be in the right complex half-plane. If e iOt or e iðOtAEpÞ is in first quadrant, then O < μ and stronger nodes phase lag, whilst if it is in the fourth quadrant, then O > μ and stronger nodes lead. This is confirmed with Fig 6, where the angles are in the same quadrants as for Fig 5. Two clusters with identical bi-modally distributed delays. For fully ordered delayimposed clusters, the mean field felt by the ith oscillator is sum from the delayed mean fields of Phase-lags in large scale brain synchronization each population where the subscripts correspond to the mean field from its own (internal), and from the other (external) subpopulation. Assuming steady synchronization with frequency O for both clusters, due to the symmetry it follows that where ψ 1,2 is the phase shift of each cluster and the first is set to zero initial phase without loss of generality, ψ 1 = 0. In [13] it was shown that ψ 2 = 0 ± π depending on the sign of cos Oτ 2 , with a stability criterion Thus, from Eq (12), synchronized oscillators in the reference frame O, Eq (25), follow In phase synchronization. Stable in-phase solutions of Eq (27) are identical to the case of spatially homogeneous delays, Eq (17), with the same condition for synchronization, Eq (18). Since these are stable only for e iOt and e iOΔτ on the same side of the imaginary axis, the transformation of the sine sum and the properties of arcsin(x), lead to which is the same relation as for spatially random delays, Table 1. These results are numerically confirmed in Fig 7(A) and 7(B), where for identical coupling strengths, faster oscillators phase lead the slower ones. The limits of synchronization are also confirmed, as well as the behavior of the phases for ω i ≷ O as predicted by Eq (28). Phase-lags in large scale brain synchronization Anti phase synchronization. For Oτ 2 in the left complex half-plane, clusters are anti-phase locked with mean fields at distance π and same frequency and level of synchronization. Relative phases are the stable solutions of Eq (27), which read whilst the criterion for entrainment of single oscillators is K i r|sin OΔτ| > |μ − O|. As before, transforming e iOt and e iOΔτ to be in the right half-plane, without any loss of generality we focus Phase-lags in large scale brain synchronization only on the interval (0, π), yielding Hence, the phase shifts within the same cluster have the same dependence from the sign of ω i − O, as for in-phase synchronization, Figs 7 and 8. Next we analyze the sign of μ − O, which for deterministic bell-shaped frequency distribution or for stochastic frequencies implies the sign of ω i − O for a large majority of oscillators. For all-to-all equal couplings and Lorentzian frequencies [13], the low-dimensional dynamics Phase-lags in large scale brain synchronization reads _ z I;II ¼ ðim À gÞz I;II À K=4½ðz I;II 2 z I;II Ã tÀ t 1 À z I;II tÀ t 1 Þ þ ðz I;II 2 z II;I Ã tÀ t 2 À z II;I tÀ t 2 Þ: ð31Þ Taking also into account Eq (25), the frequency of steady synchronization becomes This allows evaluating increase/decrease of phases ϕ i as function of the nodes strength K i in Eq (27). Results for in-and anti-phase synchronization, Eqs (28) and (30), are summarized in Table 2. Note again that node-strength dependent change of the phases is identical for both models, Tables 1 and 2, when Oτ 2 is in the right complex half-plane. To compare nodes from different hemispheres during anti-phase regime, the π shift between the mean phases need to be considered in relations in Eqs (27) and (30). Henceforth, due to the periodicity of phases, the nodes whose phases lag the most in one cluster, will be closest to the leading nodes in the other. Similarly, for ω i > O stronger and slower nodes will be further away from the mean of the opposite cluster, although they will be closer to their own mean phase. This is illustrated in Fig 7(C) and 7(D) with nodes of the opposite ends in the frequency spectrum for each population being closest to each other-light red and light blue nodes are furthest apart within, but closest to each other between clusters.
Phases for distributed coupling strengths and stochastic inhomogeneities read with a same stability criterion as for deterministic case, while synchronization limits accordingly become K i r|cos OΔτ| > |μ − O| and K i r|sin OΔτ| > |μ − O|. Plots in Fig 8, show the same examples as in Fig 7, but with stochastic frequency inhomogeneity and log-normally distributed coupling strengths that dictate phase offsets. As predicted by Eqs (28) and (30)  Phase-lags in large scale brain synchronization Results for both types of networks, Figs 6 and 8, indicate that only if both delays are in the first quadrant, stronger brain regions will always phase-lag behind the weaker ones, in agreement with [22]. For this also the natural frequencies of each region are supposed to be equal on average, so that they will be all locked at a lower frequency than their natural. Hence, assuming realistic range of conduction velocities, this would hold only for lower bands of the EEG frequencies.
The above analysis and the results for bimodal δ time-delays, homogeneous or clustered, can be generalized for any multimodal δ distributed delays, or for combination of random and clustered delays. However these are not supposed to bring qualitatively new types of steady synchronization, while making the analysis more cumbersome.

Connectome based modeling (numerical results)
General analysis for networks with time-delays is practically impossible to this date. Analytic approaches exist only for certain types of complex networks [16] combined with special delay heterogeneities [13,36], but they are still limited to the thermodynamic limit or require averaging. Besides, the connectome typically consists of less than 100 nodes, rarely going above several hundreds, and the state of art large-scale brain-modeling considers personalized connectomes [51]. Therefore, numerical simulations scrutinized by the analytical insights for simpler network topologies are a reasonable direction to proceed with the analysis of the brain networks dynamics.
In-silico oscillatory neural activity is explored over connectome based architecture to better understand the phase relation between signals from distant brain areas. A human connectome, Fig 9, is randomly chosen from a list of 1200 publicly available healthy subjects part of the Human Connectome project [52]. The subject was scanned on a customized 3 T scanner at Washington University and the structural connectivity was constructed using a publicly available pipeline [53] that applies spherical deconvolution method to a probabilistic streamlines tracking algorithm [54]. The obtained connectome consists of few million tracts spatially averaged to connect 68 cortical regions defined according to Desikan-Kiliany atlas [55]. Note however, that different parcellations are also possible, for example by subdividing each of the Phase-lags in large scale brain synchronization cortical regions [56], and these can consist of several thousand nodes [57], but are not commonly used in simulations because of the computational cost.
For each link, weights are numbers of individual tracts, Fig 9(a), and lengths are their averages, Fig 9(b). Spatial distribution of the track lengths show that they are bimodaly distributed, Fig 9(d), with the modes being spatially heterogeneous, and as a first approximation corresponding to the intra-and inter-hemispheric links , Fig 9(e) and 9(f). This insight suggests that some of the aspects of the large-scale brain dynamics are expected to be explained by the results for fully ordered delays.
The propagation velocity is fixed within the realistic range [2,58] at 5m/s, and dynamics are analyzed at different frequencies and coherence levels. The latter are additionally constrained by the noise and the global coupling strength that multiplies the normalized weights of the connectome. Since the distribution of natural frequencies across brain regions is generally unknown, equal values with stochastic inhomogeneities are assumed at each node. Moreover, even band-pass filtered recordings of neural activity in most of the cases consist of several overlapping rhythms, which are time-varying and activity-dependent, henceforth equal on average for long recordings.
Time delays cause coexistence of multiple stable frequencies of synchronization, larger or smaller than the natural, and can lead to amplitude and oscillation death in more complex systems [59]. However, we observe that unlike for the networks with bi-modally distributed delays, all numerical simulations on the connectome evolve towards a state with lower frequency than the natural, as often reported for different configurations of delay-coupled phase oscillators [60][61][62][63].
Pair-wise phase lags. Even though the spatio-temporal structure of the connectome is far more complex than networks with bi-modal δ time-delays, results in Fig 10 still show in-(A, D) and anti-phase (B) clusterings between the brain hemispheres for realistic frequencies and different levels of synchronization. An intermittent state of in and anti-phase epochs is also often observed, panel (C), as seen by the mean-field parameters shown on the bottom. If the frequency is such that μhτ ext i is in the right hemisphere, then the latter regimes occurs for most of the cases with low coherence. High coherence almost exclusively leads to slowing down that pushes Ohτ ext i in the first quadrant and therefore in-phase synchronization. Nevertheless such levels of coherence are not expected to occur in a healthy brain, and these regimes are biologically implausible, Fig 10(A) and 10(D).
The good match of the simulated brain dynamics with the theoretical predictions for networks with much simpler structure, is not only limited to the regimes of synchronization. As predicted, in all examples in Fig 10 stronger nodes in each hemisphere generally lag in phase, since O < μ. This occurs for in-and anti-phase synchronization, but also during the intermittent regime. The division between the latter two is often fuzzy, since the intervals of anti-phase synchronization rarely last longer than several seconds, before being interrupted with in-phase epochs.
Anti-phase regime is assumed when the hemispheric complex order parameters are at a distance larger than π/2, and in-phase otherwise, allowing comparison with the analytical results. They capture the dynamics fairly well, even for a distance not much larger than π/2 as shown in Fig 10(C). A better approach would be intermittent intervals to be analyzed separately, since the frequency of synchronization might differ during each interval, and averaging it can lead to wrong values for the phases. Moreover, even if varying frequencies of synchronization are properly detected, averaging of the relative phases over different regimes, Fig 10(B) and 10(C), makes them to be distributed at distances smaller than π, which might be mistaken for an actual stationary clustering, rather than a mix of 0 and π clusterings. This is shown in Fig 11, where PLV and the instantaneous and time-averaged phase lags are shown for an intra and an inter hemispheric links.
The left panel of Fig 11 shows an intrahemispheric link between in-phase brain regions, and the right depicts an interhemispheric link with epochs of in-and anti-phase locking. Since Phase-lags in large scale brain synchronization K 26 > K 30 the phase difference Δϕ 30,26 2 [0, π/2) and results are very similar for the both significant levels, despite their large difference. The region 41 is stronger connected than the 14, so it is expected that during in-phase intervals Δϕ 14,41 2 (−π/2, 0], and Δθ 14,41 2 (π/2, π] for antiphase, c.f. with Fig 8(n) and 8(o). Consequently, distributed peaks appear for the histograms of phase lags in the bottom plots of Fig 11(B), but their mean is in (−π/2, −π] leading to possible wrong conclusion about the synchronization of these nodes. Whole brain phase lags statistics. Whole brain phase statistics are characterized by the mean and the standard deviation of the PLVs, and the correspondent phase-lags for each pair of brain regions. These are shown in Fig 12, where 1-standard deviation is plotted to keep the colors/coherence consistency across the images. In the upper row, the regions are arranged according to Desikan-Kiliany atlas [55], with the left hemisphere first, while in the lower row, nodes of each hemisphere are ordered increasingly according to their strength. Strengths of the tracts are reflected in PLV (first column in Fig 12), where the links with stronger direct connection show higher functional connectivity. The negative bias of the tracking techniques towards interhemispheric connections is also manifested. Consequently fewer external links have significant coherence, especially for the higher surrogates criterion, when Phase-lags in large scale brain synchronization only the strongest links survive. Phase lags within hemispheres, especially for stronger links, are around 0, revealing the in-phase synchronization. Between hemispheres, the phase lags are less informative, due to the intermittent in-and anti-phase synchronization, also seen in Fig  11(B) for the same regime. Still, many inter-hemispheric links have phases around ±π or closer to ±π/2. The latter is often hallmark of intermittent in and anti-phase regimes, as discussed for the results of Fig 11(B). The intermittency is also manifested by increased variance of the phases, visible for the inter-hemispheric links. This is much less manifested in the coherence, which stays stable during different regimes, as was also indicated by the hemispheric order parameters in Fig 10. The large variation of the overall order parameter observed there, is only due to the bursts of anti-phase ordering when it gets close to 0, whereas the coherence within each hemisphere is quite stable.
The impact of the chosen significant coherence, and the difference between instantaneous and averaged phases for the phase statistics of each link is illustrated in Fig 13 for anti-phase regime. Higher significance level causes only slightly larger variance for the phase lags, panel (A), but as seen in Fig 11, it can substantially reduce the number of accounted links, especially between hemispheres, illustrated through one such a link in panel (C). It is due to the latter mechanism that the overall distribution of the mean lags is less uniform for higher surrogates, panel (B). On contrary, time-averaging stronger decreases the variance for the links, because it diminishes the network and stochastic heterogeneities, but it does not affect the means of the phase lags for particular links, as can be also seen for the link in Fig 13(C) or in Fig 11(B).
The overall statistics for the distinctive phase regimes in the brain are illustrated in Fig 14. Phase lags and PLVs are depicted for a same subject, for various frequencies and coupling strengths, with noise proportional to the frequency to account for the frequency dependent decrease of the coherence [47]. Time-averaged coherence is shown in the first column, in addition to the mean of the PLVs of each time-window. The former speaks about the overall regime of synchronization, whilst the latter depends on the length of windows compared to the frequency and is more affected by the noise. The mean PLV alone is hence not informative, but needs to be compared with a significance level. Mean phase lags for times of significant Phase-lags in large scale brain synchronization coherence are shown in the third and fourth column for two different surrogate procedures. They produce largest difference for anti-phase regime (second row), which requires low coherence that is even smaller between hemispheres due to fewer tracts. For very low synchronization, as shown on the bottom, they are identical and therefore only one is shown, while for high overall coherence (second and third row), higher significance level discards the tails in phase lags' distribution (last column), which mainly represent links between weaker nodes, thus making the distribution sharper.
Possible paths for transition between different regimes of synchronization are also shown in Fig 14. For low frequencies, in-phase synchronization occurs (first row), which becomes intermittent/anti-phase for increased frequency (second row), at similar level of overall coherence. By increasing the coupling and henceforth the level of coherence, the brain switches to in-phase regime (third row) that can again switch to anti-phase by further increasing the frequency, but only at very low overall coherence (bottom row). Also note that the overall low coherence at the bottom row leads to spatially homogeneous values for the mean PLVs, as compared to the cases with much higher coupling and global partial coherence shown in the second and third row. Low coupling renders all links to be around a same level of coherence, without strongly coherent and incoherent links like in the middle rows, and as a result, for every pair of regions exists at least one time-window with statistically significant PLV. Henceforth the absence of links with no significant PLV.
Spatial distribution of phase-lags in Fig 14 is in agreement with the theoretical predictions. Besides being 0 centered for strong nodes within the same hemisphere regardless of the synchronization regime, lags around 0 and ±π appear between the hemispheres, resembling inand anti-phase regimes. In addition, weak regions lead the stronger, and for anti-phase hemispheres π is added. Hence, the inverted distribution of green and blue shades for the intra and inter-hemispheric links in the phase lags matrices, with darker shades corresponding to ±π/4 for internal links, and lighter for external with the values around π ± π/4. Phase-lags in large scale brain synchronization Frequency dependent spatial distribution of phase lags is illustrated in Fig 15 for intra and inter-hemispheric brain subnetworks, for two frequencies and a same global coupling. The subnetworks consist of the 10 strongest brain regions in each hemisphere based on the sum of their outgoing links. Strength of the nodes is reflected in their size, whilst links are color-coded with their phase lags taken from the upper triangle of the matrices. As predicted, strong coherence is observed during in-phase synchronization at f = 6Hz, which together with similar strengths of the nodes, renders almost zero phase lags for all the links, internal and external. During anti-phase regime observed at 20Hz, the links within the hemispheres have lags distributed around 0, but much wider than before, whilst those between them are distributed around ±π. The coherence decreases for increasing frequencies, and together with the earlier discussed Phase-lags in large scale brain synchronization non-stationarity, these cause far higher variability of phases, than during in-phase synchronization. Hence the appearance of dark shades of green and blue for in-phase, and light for antiphase synchronization.

Discussion
In this study we analyze the mechanisms by which the spatio-temporal features of the connectome impose the architecture of phase lags between distant brain regions. A general relationship is provided for how organization of time-delays drives the hemispheres to in or anti phase coherence, whereas the topology dictates the sign of the phase lags. Both aspects of connectome also determine the overall coherence, which restricts the regimes of phase organization that can be observed. The presented qualitative findings are relevant for phases in frequency decomposed neural activity.

Model set-up
Phase lags are analyzed when only pair-wise interactions are explicitly considered, and for network connectivity. The former approach is justified if interactions are too weak, and can be represented on average as a stochastic influence to the inherent dynamics at each region. This leads to mean value of the phase differences at 0 or ±π, depending on whether the delay is long enough to change the sign of the interaction.
Despite its simplicity that allows analytical tractability, the phenomenological oscillatory model resembles the non-stationary oscillations of the neural activity, which is characterized by transient synchronization. To better understand the underlying organization that regulates the large-scale brain dynamics, and henceforth the phase relationships between network nodes, we analyze synchronization for networks with bimodal delays as a first approximation of the connectome. Theoretical insights are validated numerically for more realistic frequency and couplings heterogeneities, and compared with in-silico brain dynamics, while examining the methodological limitations.
Phase lags during these epochs of coherence depend on the delays, which are constant, and on coupling strength and frequency mismatch. The latter can be different across the timeseries, but the statistics of the narrow frequency content is generally expected to be similar across the regions. Consequently, the natural frequencies are modeled as stochastic with equal means.

Frequency depression, and in-and anti-phase synchronization
Accounting for the network dynamics is a more complex approach that is more realistic, especially when the overall coherence is not insignificant. The brain network model predicts that Phase-lags in large scale brain synchronization the distribution of the phase lags will always have a peak at 0, with an additional peak at ±π appearing for anti-phase synchronization.
Of crucial importance here is whether the frequency increases or decreases during synchronization. For lowest frequencies of electro-physiological brain signals, delays cause relative phases within the first quadrant, so the frequency is depressed regardless of the topology of the delays, whilst the phase locking is in-phase. Hence stronger nodes on average phase lag behind the weaker, for any arrangements of the natural frequencies. For higher frequencies, theoretical and numerical results show that for ordered networks both directions of the frequency shift are possible. However, in-silico brain dynamics exclusively shows frequency depression, so for equal stochastic frequencies better connected brain nodes phase-lag behind the weaker, whereas for anti-phase regime the π distance should be also accounted.
The frequency depression due to delayed interactions has wider importance than the 1:1 synchronization that is discussed here. Slowing down in an anatomically constrained dynamical system with noise has been shown to induce the whole-brain FC [64], by utilizing power to phase interactions. The latter are one of the mechanisms for cross-frequency coupling [65], which besides the well-known beta-theta interactions in the hippocampus [66], are also shown to occur for cortical signals [67].
The effects of signal mixing and spread due to volume conduction cause artificial synchrony between nearby sources that are alleviated with inverse source reconstruction techniques [68]. Nevertheless, linear mixing of signals from multiple sources can still lead to wrong coherence and phase synchrony estimates and this is commonly eliminated with interaction metrics that detect exclusively lagged interactions [28,69]. This comes at the cost of an inability to detect true zero-phase lag interactions, which as we show, may be instead neurophysiologically meaningful and due to the coupling structure, as also suggested by other studies [33]. Including the actual zero-phase lagged interactions would henceforth potentially have an important impact on whole brain data analyses of M/EEG data, as it has been found that indices of Functional Connectivity sensitive to zero lag such as PLV tend to be more reliable within groups and across sessions [70,71].

Ordered versus complex networks
Although the general theoretical findings for the ordered networks still hold for the simulated dynamics over the connectome, the main contributor to their disparity is the complex spatiotemporal structure of the connectome. This is firstly reflected in the distribution of timedelays, which is neither exactly bi-modal δ, nor fully structured or random. Secondly, the weights of the links are not homogeneous for the nodes, as assumed by our approximation, but differ by several orders of magnitude. For the latter, combining other network measures, such as centrality or clustering coefficient [72], could potentially increase the predictability of the analysis. The remaining open issues of our brain network model are conceptual and are a common concern for most of the studies based on the connectomics. These are questions about the meaning of the weights and utilization of links, but also about the actual propagation velocity along tracts, which is shown to depend on large number of quantities [73].

Methodological limitations
The level of significant coherence does not impact the overall architecture of phase lags, although when it is lower it mostly increases the variance of the results by flattening the distribution. However, the stricter level of significance can fail to capture phase-locking, especially between the hemispheres where the coherence is lower due to reduced wiring.
Increased variance of the phases can also indicate an alteration between different stable states. This often causes intermittent in-and anti-phase synchronization, when several peaks appear in the distribution of pair-wise phase lags. We showed that averaging of these non-stationary dynamics leads to improper description of phase relationships and can be avoided by differentiating of the separate regimes during the analysis. However, identification of timedependent dynamics is a major challenge in analysis of biological signals [74].
Inherent variability of the frequencies or coupling strengths [48,75] is another source of non-stationarity for which we demonstrated that the observed phase-lags depend on the overall statistics of the averaged parameters. Nevertheless, dividing the time-series to different epochs for more precise identification of phase lags for different regimes is also possible in cases when the non-autonomous forcing can be recovered [74], as well as quantification of the non-autonomicity [76].
Besides the notion of synchronization, functional brain connectivity can be also described by directed information flows [77], or effective connectivity [78,79]. Bayesian frameworks [80,81], although limited to instantaneous interactions, offer another approach for studying the connectivity between neural systems, by inferring coupling functions [82] that are spatially and frequency specific [42,43].
Supporting information S1 Files. Data and scripts. The folder contains the data used in this work together with the Matlab codes (scripts and functions) necessary to perform the simulations and the analysis. Data and codes are also available at https://gitlab.thevirtualbrain.org/spase.petkoski/plos_PLSNS. (ZIP)