Transmission time delays organize the brain network synchronization

The timing of activity across brain regions can be described by its phases for oscillatory processes, and is of crucial importance for brain functioning. The structure of the brain constrains its dynamics through the delays due to propagation and the strengths of the white matter tracts. We use self-sustained delay-coupled, non-isochronous, nonlinearly damped and chaotic oscillators to study how spatio-temporal organization of the brain governs phase lags between the coherent activity of its regions. In silico results for the brain network model demonstrate a robust switching from in- to anti-phase synchronization by increasing the frequency, with a consistent lagging of the stronger connected regions. Relative phases are well predicted by an earlier analysis of Kuramoto oscillators, confirming the spatial heterogeneity of time delays as a crucial mechanism in shaping the functional brain architecture. Increased frequency and coupling are also shown to distort the oscillators by decreasing their amplitude, and stronger regions have lower, but more synchronized activity. These results indicate specific features in the phase relationships within the brain that need to hold for a wide range of local oscillatory dynamics, given that the time delays of the connectome are proportional to the lengths of the structural pathways. This article is part of the theme issue ‘Nonlinear dynamics of delay systems’.

Transmission time delays organize the brain network synchronization Spase Petkoski and Viktor K. Jirsa Institut de Neurosciences des Systèmes (INS), Inserm, Aix Marseille Univ, Marseille, France SP, 0000-0003-4540-6293; VKJ, 0000-0002-8251-8860 The timing of activity across brain regions can be described by its phases for oscillatory processes, and is of crucial importance for brain functioning. The structure of the brain constrains its dynamics through the delays due to propagation and the strengths of the white matter tracts. We use self-sustained delay-coupled, non-isochronous, nonlinearly damped and chaotic oscillators to study how spatio-temporal organization of the brain governs phase lags between the coherent activity of its regions. In silico results for the brain network model demonstrate a robust switching from in-to anti-phase synchronization by increasing the frequency, with a consistent lagging of the stronger connected regions. Relative phases are well predicted by an earlier analysis of Kuramoto oscillators, confirming the spatial heterogeneity of time delays as a crucial mechanism in shaping the functional brain architecture. Increased frequency and coupling are also shown to distort the oscillators by decreasing their amplitude, and stronger regions have lower, but more synchronized activity. These results indicate specific features in the phase relationships within the brain that need to hold for a wide range of local oscillatory dynamics, given that the time delays of the connectome are proportional to the lengths of the structural pathways.
This article is part of the theme issue 'Nonlinear dynamics of delay systems'.

Model and methods (a) Brain network model
The model is built over connectome-based architecture that dictates the strength and the delay of the interactions between brain areas, whose inherent averaged neural activity is described with three types of self-sustained amplitude oscillators.
A healthy human connectome is chosen from the Human Connectome project [47], where a customized 3 T scanner was used for the magnetic resonance imaging (MRI). The structural connectivity is reconstructed from the diffusion tensor imaging (DTI) using a pipeline [48][49][50] that yields N = 68 cortical regions, delineated according to the Desikan-Kiliany atlas [51], figure 1a.  For each link, its weight is the numbers of individual tracts between the pair of regions, and their mean length divided by the propagation velocity that is set at 5 m s −1 , which is within the experimental range [52], gives the time delay associated with that link, figure 1b.
For all three types of oscillators, there is a variable resembling the natural frequency and it is set to ω 0 = 1 rad s −1 for all the nodes. The chaotic system has the forṁ where parameters are set to be: a = 0.2, b = 0.2 and c = 5.7, so that untrained oscillators are in chaotic regime. The nonlineary damped oscillator has the forṁ where the nonlinearity is controlled with the term m that is set to 0.75, with m = 0 corresponding to harmonic oscillator, and m 1 being the quasi-linear case. The non-isochronous oscillator has the formẊ where X i = x i + jy i , with j representing an imaginary unit, the amplitude is set to 1, and the level of non-isochronicity is set to q = 0.5, with q = 0 corresponding to the isochronous case when the phases and the amplitudes are untangled. For each link, K ij and τ ij are coupling strengths and time delays, whereas for the additive Gaussian noise η i (t) = 0 and η i (t)η j (t ) = 2Dδ(t − t )δ i,j , with · denoting time-average operator. The coupling takes the form of a linear difference, as the simplest approximation of the general coupling function, and it affects the first variables for the Rössler and VdP systems, and both variables for the LS oscillators, as is usually the case [8]. Choosing a linear additive coupling would allow obtaining the same form by rescaling the parameters of the general AF bifurcation in its realizations as LS and VdP oscillators. For the Rössler case, the difference coupling is not transformable to the additive one, and even though y and x are linearly interdependent, indicating that a rescaling should be possible in order to obtain qualitatively similar dynamics, this is not as simple as for the LS and VdP oscillators. It is worth noting that even biologically inspired neural mass models where chemical synapses lead to much more complex interactions often operate close to the AF bifurcation [21,31], allowing linearization of the coupling function. Nevertheless, an analysis of different working points of these models might also be required to fully understand the impact of the delays for the large-scale brain dynamics.
The principal component of the oscillations in each model is set to the chosen frequency for the simulated neural activity by rescaling the time. This is possible because the models are phenomenological and the parameters are chosen such that they would preserve the chaotic, nonlinear and non-isochronous regimes respectively. We fix the velocity of propagation in the connectome at 5 m s −1 , which is within the experimental range [52], and following the earlier work [19], we choose frequencies of 5 Hz and 20 Hz as values that are expected to show distinctive dynamics due to the time delays. As the natural frequencies are set at ω 0 = 1 rad s −1 , the time is accordingly rescaled with factors 5 × 2π and 20 × 2π , while the mean intra-and inter-hemispheric time delays are τ int = 6.5 ms and τ ext = 19.6 ms, respectively.

(b) Analysis of the phase dynamics (i) Phases from time-series
The complexity of the model makes the analytical derivations of phases at each region a cumbersome, and probably impossible task. Instead, in order to quantify the phase dynamics, angle variables [8] are defined for each node as These represent protophases that are then used to obtain the phases [53], which by definition need to grow linearly, although this step is often skipped [8]. Alternatively, in the time-series analysis the protophase is often estimated from an oscillatory variable obtained by Hilbert or Wavelet Transform [8,54]. The phases are calculated from protophases as [53] Note that for all the analysis, a steady synchronization arises at frequency Ω =Ψ =θ i | sync [18,19], where the mean field phase Ψ [55] is obtained from the complex order parameter of each hemisphere where r defines the level of synchrony. Thus, the analysis of the phases is completed in the rotating frame Ω where the relative phases are rewritten as θ i → θ i − Ωt.
(ii) In-and anti-phase coherent network dynamics: theoretical background oscillators [19]. Assuming that the synchronization of each node depends on its strength K i , which is defined as the sum of weights of all links connecting the node i, it reads where τ = (τ ext − τ int )/2 andτ = (τ ext + τ int )/2 for the means of the weighted inter-(external) and intra-hemispheric (internal) delays. Anti-phase synchronization appears for Ωτ ext in the left complex half-plane [19], and the condition for synchronization of each region reads

(iii) Statistical analysis of the phase locking
Phase relationships between brain regions are quantified using phase locking values (PLV) [7], which are a statistical measure for similarity between phases of two signals, frequently employed in the analysis of empirical data. Although synchronization can occur between different frequencies [54,56], most commonly studied is the one-to-one entrainment, for which the complex phase locking values (cPLV) are defined as Here, we calculate cPLV at sliding windows of length equal to 5 periods of the calculated mean frequency of the entrainment, Ω, [55] and with 25% overlap. To identify only the phase coherence due to mutual interactions [7,19,57], we calculate a level of statistical significance for PLV as the 95th percentile of maximum values in 100 surrogate signals that are obtained by shuffling one of the phases, and which follow the same processing as the original signals.

Results
Even though all three oscillatory systems making up the BNM have more complex dynamics than the harmonic phase oscillators [19], they also synchronize such that brain hemispheres are in-and anti-phase locked, depending on the frequency, as in figure 2. Here, scatter plots of phases and nodes strengths are depicted, together with theirs probability density distributions (PDF) and the time-evolution of the phases in the rotating frame Ω of the mean-field. Left panels of figure 2 are for a frequency f = 5 Hz, so that even the long inter-hemispheric time delays do not cause phase shifts Ωτ larger than π/2, while in the right panels f = 20 Hz and the phase shift is in the left complex half-plane. As with the KM, within the hemispheres the stronger nodes generally lag in phase behind the weaker, while the π shift needs to be accounted for when comparing contralateral regions during anti-phase regime. Moreover, the analytical result for the relative phases, equation (2.7), holds quite well, with the only significant deviation appearing for the chaotic oscillators.
For the isochronous case of LS oscillators, their phase dynamics is fully captured by the KM [34]. However, even for the non-isochronous case chosen here, as with the nonlinear VdP, the phase shifts during different regimes of synchronization are similar to the case of KM [19]: in general, weaker nodes drift forward, and the stronger lag behind, while being locked to the mean field for most of the time. Contrary to the KM, the anti-phase synchronization is not necessarily intermittent for BNM using amplitude oscillators, and therefore the arrangement of the phases within the hemispheres is not affected by the type of locking (in-or anti-phase). The only anomalous behaviour is for small delays and chaotic oscillators, when some of the stronger connected nodes lag behind the mean field [55], which is more aligned to the weaker synchronized nodes. This slowing down is not reflected in the peak frequencies of the power spectrum, which are identical for synchronized oscillators of different strength, figure 3, and they   coincide with the frequency of the complex order parameter. The power spectrum also shows that the entrainment for chaotic low-frequency activity is at a higher frequency than the natural, i.e. the frequency of an uncoupled oscillator. Nevertheless, this does not seem to reverse the dependence of the relative phases of the nodes strength, unlike for phase oscillators [19], suggesting that an artefact in the recovery of phases might be the cause of this anomaly. Phase portraits in figure 3 show that network interactions inevitably cause distortions of limitcycles, beyond the stochastic nature of the added noise. This gets more pronounced when delays are comparable with the inherent time-scale of the oscillators, as is demonstrated in figure 4 for     all three BNMs at frequency f = 20 Hz. The strongest nodes are mostly affected, and in general the amplitude decreases with the strength of the nodes, so the oscillations become less nonlinear and more stochastically driven. In the case of chaotic oscillators, they become more regular and harmonic, the same as for the nonlinearly damped model.
Time-series of the neural activity from nodes with a different strength, figure 4, depict the inand anti-phase synchronization within and between the hemispheres, respectively. Three of the strongest, the weakest and medium strength nodes within the same hemisphere are shown for each of the models, as well as the activity of strongest contralateral nodes. In all the cases, the weakest nodes are unsynchronized with the rest, and the strongest and the medium nodes are rather coherent between each other, while being anti-correlated with the opposite hemisphere. This is visible even though the oscillations are quite distorted, as is also reflected in their power spectra. These show that the strongest nodes, despite losing their harmonics, still have the most complex dynamics during anti-phase arrangement, with the largest peak at the synchronization frequency that is always lower than the natural. On the contrary, the weakest, unsynchronized nodes are speeding up compared to the isolated nodes for the LS and the Rössler models and are in general faster than the rest of the network, as is also visible in the time-series. It is also worth noting that the frequency depression is highest in the LS system, while chaotic oscillators are the least influenced by this.

(iv) Pair-wise phase lags statistics
The most typical pair-wise correlation and phase coherence are illustrated in figure 5, and besides the in-and anti-phase synchronization, they include the varying case, when both are intermittently appearing. The left panel shows an intra-hemispheric link between in-phase synchronized brain regions, and the middle and the right panels depict inter-hemispheric links.  In the first case, both nodes are weakly connected to the rest of the network, and hence epochs of in-and anti-phase locking exist, leading to multi-modal distribution of the phase lags. Unlike phase oscillators, when even the strongest nodes slip from anti-to in-phase arrangement, this only occurs between weak nodes for amplitude oscillators. As predicted by equation (2.7) [19], since K 24 < K 44 and the link is inter-hemispheric, the phase difference φ 24,44 ∈ (π/2, π ). On the other hand K 22 < K 34 and both nodes are in the left hemisphere, so φ 22,34 ∈ (−π/2, 0). Finally, K 5 < K 65 , so the inter-hemispheric phase difference should be φ 5,65 ∈ (π/2, π ) during the antiphase epochs, as it is during the beginning of the highlighted interval, before slipping in (−π/2, 0) for the rest of the interval when both regions are in-phase. (v) Whole-brain phase-lags statistics The whole brain statistics for the distinctive phase regimes observed during oscillatory brain dynamics are shown in figure 6. It depicts the mean PLV for each pair of brain regions, and the correspondent mean and the standard deviation of the phase-lags calculated from cPLV. To keep the colours/coherence consistency across the images 1−standard deviation is shown. All links in the plots have statistically significant coherence detected during at least five time windows, because the chosen method yields a low level of significance for PLV [19] (see also figure 5). The latter two metrics generally mirror the strength of the structural links and can be used to describe the functional connectivity, but are still differently affected by the underlying dynamics. For example, underrepresentation of the inter-hemispheric connections by the tracking techniques is clearly reflected in the PLV values for in-phase regimes, although lower levels of PLV do not necessarily invoke high variability of the phases. This is especially the case for the anti-phase dynamics of the chaotic BNM, but also the absence of intermittent in-and anti-phase dynamics for all the models, contrary to the Kuramoto oscillators [19], induces small spread of the phases due to their uni-modality, which is mainly affected by the strengths of the nodes. Thus, the high variance in figure 5 for the two weakest VdP oscillators is due to their very low overall phase coherence, despite their quite high average PLV. The consistency of the regimes of synchronization, despite the weak level of statistical significance that renders all the links to be significantly coherent, still induces clearly pronounced peaks at 0 and ±π radians. The inter-hemispheric (external) links during the anti-phase regime are around ±π rad, contrary to the intra-hemispheric (internal) links that have always 0 centred phase lags, the same as the external links at low frequencies. Hence, the spatial distribution of phase-lags in figure 6 is in agreement with the theoretical predictions for Kuramoto oscillators, [19]. Moreover, strong regions lag behind the weaker, hence the green and blue shades for links in the phase lags matrices during in phase regimes. These get inverted for the anti-phase interhemispheric links, with darker shades corresponding to ±π/4 for internal links, and lighter for external with the values around π ± π/4.

Conclusion
In this paper, we computationally analyse the architecture of phases and amplitudes of large-scale neural oscillations, as they are shaped by the connectome. We identify exact features in the activity of the brain regions, given that the inherent neural dynamics is oscillatory and that the time delays due to the propagation in the connectome are proportional to the lengths of the structural pathways. In this way, we offer a framework for experimentally verifying the BNMs that utilize self-sustained oscillatory systems, which becomes more important in light of recent efforts towards building a probabilistic atlas of human cortical connections [52], specifically their time delays, which will eventually allow improving the utilized BNMs. In particular, the results can be used to test whether the inter-hemispheric time delays are significantly longer than the intrahemispheric ones, which in the case of detectable coherence and large enough frequency should inevitably lead to anti-phase arrangement between some of the stronger-connected contralateral regions. In addition, we demonstrate that the delays necessarily reduce the activity at stronger regions for increasing frequency in linearly coupled oscillatory brain nodes.
We have extended previous studies of phases of brain regions that were derived from Kuramoto oscillators [19,28], and we have shown that the lagging of the stronger regions is universal for oscillatory BNMs, as has been experimentally observed [28,33]. The other consistent characteristic across the models is the decrease of the amplitude with an increase in the frequency and the coupling strength, especially pronounced for stronger nodes. This is in contrast to the findings about LS oscillators with time delays reduced to phase shifts [28], suggesting that besides the anti-phase locking, reduction of delays to phases also disregards other important aspects in the emerging dynamics. Furthermore, the amplitude reduction, as observed here, for the systems near the AF bifurcation is associated with a shift of the working point towards the bifurcation. This is the same mechanism responsible for the amplitude death [35,58], which for identical oscillators can be due to delayed interactions [59]. The amplitude death is amplified for distributed delays [60], such as those due to the connectome, but it is also increasing with the coupling strength, which for heterogeneous networks is stronger for stronger connected nodes.
With this study, we also confirm that phase reductions of time-delayed couplings, specifically the Kuramoto model, are sufficient for capturing the synchronization related aspects of the dynamics in networks of oscillators, especially those near AF bifurcation. Nevertheless, besides allowing for the study of the impact of the heterogeneous delays on the amplitude of the brain activity, increasing the dimensionality of the models as expected brings new peculiarities in the oscillatory dynamics. For instance, the phase-amplitude entanglement stabilizes the anti-phase ordering, which is otherwise sporadic for the delays of the connectome [19], while the frequency spectra become more complex due to nonlinearities. In addition, although delays can cause period doubling even for the mean-field of the KM [61], amplitude oscillators allow for coexistence of different types of local oscillations, such as chaotic, quasi-periodic and harmonic, within the same network. However, the current model assumes identical parameters for each region, and a more realistic approach would require this to be specified based on the data [23] or on some clinical hypothesis [25,[62][63][64]. In this way, the specific characteristics of each region that translate in various spectral, isochronous or chaotic properties would also impact the overall dynamics, thus increasing the model's authenticity.