Long-range phase synchronization of high-gamma activity in human cortex

Inter-areal neuronal phase synchronization is thought to involve oscillations below 100 Hz while faster oscillations are presumed to be coherent only in local circuits. We show with human intracerebral recordings that 100–300 Hz high-gamma activity (HGA) can be phase synchronized between widely distributed regions. HGA synchronization was neither attributable to physiological and technical artefacts nor to epileptic pathophysiology, e.g., it was as prevalent in areas distant from the epileptic zone as within. HGA phase synchronization exhibited a reliable large-scale connectivity and cortex-wide community structures. Moreover, HGA synchronization displayed a laminar profile opposite to that of low-frequency oscillations by being stronger between the deep than the superficial cortical layers. Hence HGA synchronization constitutes a novel and temporally highly accurate form of consistent timing relationships in large-scale human cortical activity. We suggest that HGA synchronization elucidates the temporal microstructure of spiking-based neuronal communication per se in cortical circuits.


Introduction
Observations of organized neuronal population activities at frequencies above 100 Hz, such as highgamma activity (HGA, 100−300 Hz) 1 , high-frequency oscillations (HFOs, 100−200 Hz) and ripple oscillations (150−200 Hz) 2,3 are abundant in electrophysiological recordings from cortical structures of humans 1,4 and experimental animals 5,6 . Several lines of evidence link HGA with perceptual and cognitive processes 1,4,[7][8][9] , and overall indicate that high-amplitude HGA is a key signature of "active" neuronal processing. Ripple oscillations have been traditionally associated with sharp waves and off-task memory formation, but recent studies report ripples also during the performance of attention tasks 10 .
Electrophysiological HGA signals are thought to mainly arise from broad-band multi-unit spiking activity (MUA) and hence to directly reflect the local peri-electrode neuronal population activity per se 11,12 . HGA 13 and ripple oscillation signals do, however, also contain contributions from postsynaptic potentials. While the synaptic mechanisms underlying the hippocampal ripple oscillations are already well understood, it appears that genuine oscillations with presumably synapticcommunication-based mechanisms also contribute to the HGA signals 13 .
HGA has been thought to be exclusively local in terms of spike synchronization and phase coupling. For slower (< 100 Hz) neuronal oscillations, phase relationships among distributed neuronal assemblies are instrumental for coordinating neuronal communication and processing 14,15 .
Several lines of experimental and theoretical evidence have shown that these phase relations are dependent on frequency and neuroanatomical distance, and on the axonal conduction delays in particular, so that slow oscillations are generally more readily phase coupled over long distances than fast oscillations [16][17][18] . In line with this notion, measurements of long-range phase coupling in animal 15,19 and human 16

brains suggest that neuronal oscillations only in frequency bands below 100
Hz exhibit inter-areal phase synchronization 20 whereas the inter-areal cooperative mechanisms of HGA have remained poorly understood.
We hypothesized that long-range synchronization and phase coupling of HGA, even if unexpected, could nevertheless be conceivable because local synchronization and high collective firing rate can endow local pyramidal cell populations greatly enhanced efficiency in engaging their post-synaptic targets 21 , which would be experimentally observable as inter-areal HGA phase coupling. Such a finding would constitute a direct observation of the actual spiking-based long-range neuronal communication per se. Nevertheless, there is little evidence for HGA synchronization over long distances so far 22 .
In this study, we searched for long-range HGA synchronization using an extensive database of resting-state human stereo-electroencephalography (SEEG) recordings (N = 67). We used submillimeter accurate SEEG-electrode localization 23 and white-matter referencing 16 to obtain neocortical local-field potential (LFP) signals with little distortion from signal mixing with neighboring grey-matter or distant volume-conducted sources. We found that among these LFPs, long-range HGA synchronization was a surprisingly robust phenomenon and much stronger than synchronization at around 100 Hz. We rigorously excluded the possibilities of HGA synchronization being attributable to putative confounders such as the epileptiform pathophysiology or physiological and technical artefacts. The network architecture of resting-state functional connectivity achieved by HGA synchronization was split-cohort reliable and had a modular largescale architecture that were distinct from those of lower frequencies. These findings thus reveal in the human brains a novel form of spatio-temporally highly accurate neuronal coupling and thereby elucidate the cerebral organization of fast neuronal communication.

Probing human large-scale brain dynamics with SEEG
We recorded ~10 minute resting-state human cortical (local-field potential) LFP signals from 92 consecutive patients with stereo-electroencephalography (SEEG). Among them, 25 subjects were excluded from further analyses due to previous brain surgery (temporal lobotomy) or an MRIidentified cortical malformation (Suppl. Table 1, Suppl. Fig. 1). The final cohort of 67 patients yielded a total of 7068 non-epileptic grey matter contacts (113 ± 16.2 per subject, mean ± SD, range 70-152) that gave a dense sampling across all neocortical regions (Fig. 1a) and of seven canonical functional brain systems defined by fMRI intrinsic connectivity mapping 24,25 (Fig. 1b).
We assessed the phase interactions between all LFP signals recorded from non-epileptic neocortical grey-matter locations. This yielded a dense coverage of cortical interactions with 5,500 ± 1,600 (mean ± SD) contact pairs per subject (range: 2,094−9,947) and a total of 368,043 contact pairs. Out of all possible within-hemispheric connections in the 100-parcel Schaefer atlas 25 , these recordings sampled at least 80% of the left-and 90% of the right-hemispheric connections (Fig. 1c) and provided abundant sampling on the scale of functional systems (Fig. 1d). The present data thus yield comprehensive insight into large-scale dynamics and connectivity in the human brain.

Long-range phase correlations in high-gamma activity
We estimated inter-electrode-contact phase coupling with the phase-locking value (PLV) for 18 frequency bands between 3 and 320 Hz (see Methods). The inter-contact PLV estimates were averaged across subjects for three ranges of inter-contact-distance quartiles for each frequency band (Fig. 2a). The first quartile (0−2 cm) was excluded to avoid contamination from residual volume conduction. We found that the mean PLV increased from 3 to 7 Hz and then decayed rapidly from 100 Hz as expected 16 . However, from 100 Hz onward, the mean PLV increased monotonically, indicating highly significant HGA phase synchronization in all distance ranges.
To quantify the neuro-anatomical extent of HGA synchronization, we assessed the connection density, K, that was defined as the fraction of contact pairs exhibiting significant (p < 0.001 for observed > surrogate PLV) phase synchronization (Fig. 2b). Even at long ranges (> 6 cm), more than 70% of >100 Hz connections were significant and both the PLV and K findings were splitcohort reliable (Suppl. Fig. 4). HGA phase synchronization was thus a widespread phenomenon in SEEG recordings of resting-state brain activity. with these observations and readily revealed episodes of significant HGA coupling in time windows across centimeter-scale distances. Notably, synchronized HGA was observed as low-amplitude oscillation-like activity that was visible in the time series also without filtering, which provides first evidence for that HGA synchronization was not attributable to spikes or technical artefacts.
Given the novelty and unexpected nature of these observations, we performed a series of control analyses to exclude the possibilities that HGA synchronization were due to: referencing schemes and volume conduction (Suppl. Fig. 1 & 2), non-neuronal signal sources (Suppl. Fig. 3), inadequate filter attenuation or settings (Suppl. Fig. 5), amplifier noise or pathological neuronal activity, such as inter-ictal spikes, or contamination from muscular signals (Suppl. Fig. 6). The results of these analyses converged on the conclusion that the seemingly anomalous HGA synchronization can only Figure 2 High-Gamma Activity shows significant and robust long-range phase synchrony. a Mean phase synchrony strength and spatial extent (b) estimated with PLV as a function of frequency and in short (2 to 4.6 cm, pink), medium (4.6-6 cm, blue), and long (6 to 13 cm, green) distance ranges. Dashed lines represent surrogate (N=100) data level for p<0.001. Shaded areas represent confidence limits (twotail; p<0.05) for bootstrapped values (N=100). c-f Examples from two distinct subjects of long-range HGA phase synchrony with non-zero phase lag (c, d) and near-zero phase lag (e, f), and their temporal relation to fluctuations in slow rhythms; MFG (azure traces): medial frontal gyrus, IPS (red traces): intraparietal sulcus.
be explained by true long-range synchronization between local cortical neuronal assemblies (Suppl. Text). This conclusion was further consolidated by the findings, as detailed below, that HGA synchronization was predominant in specific functional systems, and had a community structure and laminar connectivity profile that were distinct from those of slower activities, which are inconceivable for technical or physiological artefacts.

Neuroanatomical localization of high-gamma synchronization
To characterize the neuro-anatomical features of HGA synchronization networks, we investigated large-scale phase coupling of HGA in two spatial resolutions. First on the level of functional systems 24 , robust HGA phase synchrony was found within and between all systems, but the strongest and most widespread phase synchrony was found between and within the default-mode (DEF) and limbic (LIM) systems ( Fig. 3a,b). This observation is in line with the hypothesis that HGA synchronization reflects healthy patterns of neuronal communication that is dominated by activity in these systems, DEF in particular 26 , during resting state. The observed systems-level connectivity pattern was split-cohort reliable (Suppl. Fig. 4), which rules out biases of individual subsets with distinct aetiologies.
To examine the architecture of HGA synchronization at a finer resolution, we used the 100-parcel Schaefer atlas 25 and pooled data across subjects to estimate the connectome of inter-regional phasesynchrony for each frequency. We found these connectomes to be split-cohort reliable (permutation test, p < 0.05, one-tailed) across the range of frequency bands investigated, including the HGA frequency band (Suppl. Fig. 4).
We then applied Louvain community detection to identify the putative communities therein and found sets of regions to represent robust modules in the HGA frequency bands (see Suppl. Text).
The numbers of regions assigned reliably to communities were much higher than expected by chance (bootstrap test, p < 0.05, one-tailed) for each HGA frequency and throughout the investigated range of the Louvain resolution parameter values γ = 1−1.5 (Suppl. Fig. 9), i.e., the parameter influencing the number of communities identified. These networks largely also demonstrated significant community structures compared to equivalent random networks (permutation test, p < 0.05, one-tailed) at a range of resolutions (Suppl. Fig. 9).
The communities in the high-gamma frequencies were similar to each other and dissimilar from communities in the lower frequency bands (Fig. 3c). For the 190 Hz connectome (e.g., in matrix form see Fig 3d), at a low spatial resolution, the network was split into four large communities that were divided into six to seven communities at higher resolutions ( Fig. 3 e,g). The majority of the communities comprised adjacent brain regions but also included spatially distant regions (Fig. 3 f,h). Demarcation of spatilly distal regions into the same modules is not only population level confirmation of long-range HGA synchrony observed at contact level, but also a strong indication of the functional relevance thereof. The whole-brain networks of HGA synchronization thus exhibited statistically significant and distinctive community structures in putatively healthy brain structures, as well as patterns of connection strengths that are stable across participants. These results thus further support the neurobiological origin of HGA synchronization.

High-gamma phase correlations have a unique laminar profile
Deep and superficial cortical laminae are known to contribute differently to inter-areal phasesynchronization at frequencies below 100 Hz 16 . We next asked whether HGA synchronization networks also would show a similar differentiation between the laminae across distances.
Leveraging our accurate localization of the electrode contacts, we divided the electrode contacts into "deep" and "superficial" by their cortical depth (see Methods), and assessed HGA phase coupling strength separately between the contacts in deep and superficial layers. This analysis reproduced the prior observation 16 of low-frequency (< 30 Hz) phase coupling being stronger among superficial sources, which is well in line with recent observations about the laminar localization of, e.g., alpha oscillation sources in human cortex 27 . In contrast, however, HGA synchronization was stronger (Fig. 4a, see inset) and more prevalent (Fig. 4b) among signals from deep cortical layers in all distance ranges (p < 0.05, Bonferroni corrected with N freq. = 18). This reverse pattern was also observed with iPLV and in bipolar recordings (Suppl. Fig. 2). Although the resolution of SEEG is insufficient for further investigating the underlying neuronal generators of these oscillations, e.g., with a current source density analysis, our results indicate that HGA Figure 4 High-Gamma phase synchrony differs between cortical layers a Mean synchrony strength (PLV between contacts) at short (2−4.6 cm, pink), medium (4.6−6 cm, blue), and long (6−13 cm, green) distance ranges shows distinct spectral profiles for contact pairs in deep cortical layers (-0.3 < GMPI ≤ 0; blue) and contact pairs in superficial cortical layers (0.5 < GMPI < 1; red). Dots represent the frequencies at which the difference between the superficial and deep contact was significantly (permutation test, N=100: p < 0.05, Bonferroni corrected along with 18 frequencies). The GMPI value represents the distance between contact and white-gray-matter border normalized by the cortical thickness. Inset: mean PLV of superficial and deep layer contact in the HGA frequency range. Shaded areas: the confidence intervals of 2.5% and 97.5% of bootstrapped population variance (100 bootstraps re-sampled with replacement); dashed lines represent confidence limit for p< 0.001 computed from 100 surrogates. b Mean fraction of significant edges (K) with p<0.001 computed from 100 surrogates.

Figure 5 Stronger
High-Gamma phase synchrony is associated with maximal amplitude correlations between cortical sites a High amplitude (normalized) samples between contact pairs is correlated with high phase consistency between these contact pairs. x and y axis represent the quintiles of normalized amplitude profiles of a pair of contact channel 1 and channel 2, respectively. Each element in the matrix is the mean of instantaneous PLV between all significant channel pairs (p <0.001 with 100 surrogates) as a function of their moment-to-moment normalized amplitudes. b Strength of moment-to-moment amplitude correlation. Distribution of instantaneous phase samples (light grey) in each amplitude bin and its distance from uniform distribution (black) of samples (i.e., null-hypothesis for the absence of moment-to-moment amplitude correlation). synchronization originates in neuronal circuitry distinct from that producing the slower LFP oscillations.

Long-range HGA synchronization is associated with local synchronization indexed by HGA amplitude
To further delve into the physiological plausibility of HGA synchronization, we asked whether it was related to the moment-to-moment variability in the amplitudes of local HGA. HGA amplitude is likely to tightly reflect the number of neurons in the local assembly and the degree of their local synchronization, both of which are central to the ability of a local assembly to engage its postsynaptic targets effectively. We thus hypothesize that the moments of strongest inter-areal HGA phase synchrony would coincide with the moments when both contacts recorded high HGA amplitudes. We selected electrode contacts exhibiting significant HGA synchronization and for each contact pair, distributed the data sample-by-sample into a 2D matrix according to quintiles of the normalized local amplitudes (see Methods).
First, inspecting the variability in the numbers of samples among the cells of these 2D matrices we found that there was a slight positive correlation between the amplitudes so that the coincidence of the largestquintile amplitudes was ~1 % more prevalent than the coincidence of the smallest and largest amplitudes. Second, after equalizing the sample counts in amplitude quintile pairs, we estimated the PLV from samples in each quintile pair, and averaged the PLVs across electrode pairs and subjects for each frequency separately. We found that while HGA phase coupling was significant across a range of local amplitudes, it was the strongest in those moments when the HGA amplitudes were the largest in both contacts and essentially insignificant when either location exhibited the lowest HGA amplitudes (Fig. 5a, b). Quantitatively, PLV was very dependent on the local oscillation amplitudes and exhibited a ~400 % difference between smallest and largest values in data averaged across the HGA band and a 200-800% difference in Figure 6 High-Gamma amplitude is modulated by the phase of slower oscillations a Left: Among putatively healthy regions (nEZ), beta to low gamma band amplitude is coupled with the phase of low frequencies (theta to alpha). nPAC: population mean of individual phase-amplitude-coupling (PAC) PLV estimates normalized by within-subject surrogate mean. Right: fraction of significant edges (K) for phase-amplitude coupling on the population level. b Within the epileptogenic zone (EZ), prominent PAC is observed between delta, theta and alpha and faster brain rhythms. c Connections between EZ and nEZ recording sites also show PAC coupling between theta to alpha and faster brain rhythms.
individual HGA frequencies with the strongest effects observed at highest frequencies (Suppl. Fig.   8). These data thus strongly support the notion of local HGA synchronization being instrumental for long-range HGA phase coupling. observed within the EZ (Fig. 6b), where instead of predominant narrow-band theta/alpha coupling, the amplitudes of oscillations from the lowest frequencies to HGA were locked to the phase of delta-band (here 3 Hz) oscillations. This putatively pathological pattern was also reflected in PAC evaluated between the healthy and epileptic areas, which exhibited both the delta and theta/alpha nesting of faster activities (Fig. 6c).
Is high-gamma synchrony a feature of epileptic brain tissue?
A crucial question for the functional implications of HGA synchronization is whether it is mainly a pathological property of the epileptogenic network or a feature of healthy brain activity and neuronal communication therein. The earlier finding of systematic HGA community structures across subjects with diverse aetiologies already provided evidence in support of the latter alternative. We further addressed this question first by asking whether the strength of HGA synchronization was significantly stronger in the EZ or between the EZ and nEZ. Controlling for the neuroanatomical distance of the electrode-electrode interactions, we found the epileptogenic zone to exhibit significantly stronger phase synchronization than healthy areas at low frequencies (3−10 Hz) but no significant differences were observed in the HGA frequency band (Fig. 7, see also Suppl. Fig. 7e).
We then tested whether the strength of HGA synchronization overall was correlated with the frequency of inter-ictal spikes; a proxy for the magnitude of epileptic pathology. No such correlation was observed (Suppl. Fig. 6). These data thus converged onto the conclusion that HGA synchronization is a property of healthy brain dynamics that is preserved in these patients.

Figure 7 Phase synchronization is stronger within the epileptogenic zone (EZ) in low frequencies Top:
At low frequencies (3-10Hz) and short-to-medium distances, mean phase synchrony among EZ regions (pink) is significantly stronger than for nEZ-nEZ (blue) or EZ-nEZ (green) connections. At higher frequencies (2 nd to bottom row), this difference between EZ-EZ and nEZ-nEZ and EZ-nEZ recording site becomes insignificant across all distance ranges. Shaded areas represent 5 th and 95 th percentile of PLV values (bootstrapping, N=500). Distance ranges: short (2 cm ≤ x < 4.6 cm); medium (4.6 cm ≤ x < 6 cm), long (6 ≤ x < 13 cm) In either case, little is known about the large-scale phase synchrony of HGAs or HFOs, and they have been considered exclusively as local phenomena. HGA amplitude fluctuations may be coupled inter-areally in human cortical networks during reading 33 , and similarly, ripple oscillations bursts co-occur between the rat hippocampus and association cortices during memory formation 34 .

Discussion
However, amplitude correlation does not have the temporal precision that phase-coupling has for carrying out important computational functions during neuronal processing. Nonetheless, inter-areal phase-coupling has been thought to exist only in slow brain rhythms but not in HGA frequency bands. To the best of our knowledge, there is only one reported HGA inter-areal phase coupling 22 .
In this study, Yamamoto et al (2014) reported phase synchronization of ripple oscillations within the entorhinal-hippocampal circuit in behaving rats.
We report here that in the human brains, neuronal activity in the 100−300 Hz frequency range may be phase synchronized across several centimeters during awake resting-state brain activity. Longrange synchronization was highly reliable and not explained by physiological or technical artefacts in the recorded data. HGA synchronization was also primarily a physiological phenomenon and not

Micro-and macroscale cortical architecture of HGA synchronization
HGA synchronization exhibited systematic organization at several spatial and temporal scales. At the level of cortical systems determined by fMRI intrinsic connectivity 24 , HGA synchronization was most pronounced within and between the limbic-and default-mode systems. These areas exhibit high level of activity during resting-state 26 and thus this finding is in line with the notion that HGA synchronization reflects neuronal processing or communication per se. Interestingly, in a subset of subjects who had recordings both during resting-state and during sleep, we found differences in HGA synchronization patterns between these two conditions (Suppl. Fig. 7). This suggests that HGA synchronization is state-dependent similar to fMRI resting-state connectivity in humans 35 and non-human primates 36 .
In addition to this insight into HGA synchronization among fMRI-derived brain systems, we examined the intrinsic community structure in the whole-brain connectome of HGA synchronization. We found evidence for significant and split-cohort reliable community structures in HGA synchronization. Together with the split-cohort reliability of the connectome itself, these findings show that HGA synchronization has a group-level stable cortical topology that is independent of individual subjects and thus also a structure that is not dictated by the individual pathogenesis or electrode placement. Moreover, this ties HGA synchronization with phase and amplitude interaction networks in slower frequencies, which as known to have salient community structures [37][38][39][40] . The HGA communities, however, were dissimilar with those at slower frequencies, which may be attributable to distinct cortical generator mechanisms of the slow and HGA signals.
This conclusion was supported by the finding that at the scale of cortical laminae, HGA synchronization was significantly stronger among electrode contacts in deep than superficial layers whilst the opposite laminar organization was found for theta, alpha, and beta oscillations in SEEG here and in prior work 16 . HGA synchronization thus exhibits reliable neuroanatomical organization at scales from cortical laminae to brain systems with phenomenological differences with the oscillations in slower frequencies. HGA synchronization is thus not a "by-product" of neuronal interactions coupling the slower oscillations, but rather a hitherto undiscovered component in the organization of large-scale brain dynamics.

HGA synchronization is not attributable to artefactual sources
The connectivity and community structures of HGA synchronization as well as its laminar organization strongly suggest that it may not be attributable to physiological or technical artefacts such as signals from muscles or extra-cranial sources. Nonetheless, we corroborated this notion with a number of controls. First, in addition to white-matter referencing, HGA synchronization was observable with bipolar referencing, indicating that its current sources are millimeter-scale local in cortical tissue and very unlikely to originate from extracranial or muscular sources during inter-ictal events 41 . Second, HGA synchronization was also observable with linear-mixing insensitive interactions metric and hence not attributable to volume conduction 42

HGA activity is not explained by pathological epileptic activity
The role of epileptic pathology is overall a concern for research carried out with epileptic patients.
Although HGA within the epileptogenic zone show peculiar spectral and amplitude contents 43,44 , and are temporally predictive of upcoming seizures, most recent work has questioned its pathological-only origin 45,46 . We addressed question of physiological vs. pathophysiological genesis of HGA synchronization through several lines of analyses. First, only the putatively healthy brain areas and time-windows with no epileptic spikes were included in the primary analyses of HGA synchronization. Second, if HGA synchronization was attributable to epileptic pathophysiology, one would expect highly individual large-scale patterns, but instead of such, we observed that both the connectivity and community structures of HGA were split-cohort reliable at the group level and thus not driven by individual pathology. Third, in direct comparisons of HGA synchronization among healthy areas, between healthy areas and the epileptic zone, and within the epileptic zone, we did not find significant differences between these conditions. We found the epileptic zone to be distinct from healthy brain areas both in delta-frequency phase synchronization and in the phaseamplitude coupling of delta and theta oscillations with HGA. Hence, if HGA synchronization were a property of the epileptic zone, it would likely have been observed in this analysis. Finally, we observed no phase coupling of HGA with epileptic spikes. Taken together, HGA synchronization appears to be a property of healthy brain dynamics rather than being attributable to epileptic activity or pathophysiology.

Putative generative mechanisms HGA synchronization
What mechanisms could underlie the signal generation of long-range phase coupled HGA? Whereas synaptic mechanisms are known to contribute to the generation of synchronized ~200 Hz oscillations during hippocampal sharp-wave events 47 , it has remained disputed whether broad-band HGA in the 100-300 Hz range is associated with synaptic mechanisms generating neuronal population oscillations with rhythmicity in specific time scales 19 . A recent study 13 argues that putatively-spiking-related broadband and genuinely oscillatory components with presumably synaptic-communication-based mechanisms are dissociable in the HGA signals, even though several studies suggest HGA to arise from neuronal spiking activity 48,49 per se. We speculate that the HGA signals observed in this study may reflect both local cortical population spiking activity and the consequent downstream post-synaptic potentials resulting from these volleys of spikes.
Long-range MUA synchronization has not been previously reported, but it is important to note that unlike the micro-electrodes used in animal research, the larger surface area of the electrodes heavily predisposes SEEG specifically to detecting population spiking activity 50 that is already locally synchronized in a sizeable assembly and thus able to achieve post-synaptic impact in distant targets 21 . Hence, by construction, SEEG may be effectively filtering out asynchronous multi-unit activity that is unlikely to achieve well-timed downstream effects.

Data acquisition
We recorded SEEG data from 67 subjects affected by drug resistant focal epilepsy and undergoing pre-surgical clinical assessment for the ablation of the epileptic focus. We acquired monopolar (with shared reference in the white matter far from the putative epileptic zone) local field potentials (LFPs) from brain tissue with platinum-iridium, multi-lead electrodes. Each penetrating shaft has 8 to 15 contacts, and the contacts were 2 mm long, 0.8 mm thick and had an inter-contact border-to-

Signal pre-processing
We excluded electrode contacts (1.3±1.2, range 0−10, contacts) that demonstrate non-physiological activity from analyses. We employed a novel referencing scheme for SEEG data where electrodes in grey-matter were referenced by the contacts located in the closest white-matter (CW) 2 . This referencing scheme is proven optimal for preserving phase relationship between SEEG contact data 2 . Since one same white-matter contact can be used for referencing multiple cortical contacts, we rejected derivations with shared reference. The final size of channels analysed is on average 110±25 for each subject and 7491 in total.
Prior to the main analysis, SEEG time series were filtered with 18 Finite Impulse Response (FIR) band pass filters with central frequency (F c ) ranging from 2.50 to 320Hz. We used a relative bandwidth approach for filter banks such that pass band (W p ) and band stop (W s ) were defined as 0.5×F c and 2×F c , respectively. We excluded all 50 Hz line-noise harmonics using a band-stop equiripple FIR filter with 1 % of maximal band-pass ripples and 3 up to 8Hz width for the stop band parameters.
Epileptic events such as interictal spikes are characterized by high-amplitude fast temporal dynamics as well as by widespread spatial diffusion. Due to possible filtering artefacts around epileptic spikes and the resultant increase in synchrony, we discard periods of 500ms containing Interictal Epileptic Events (IIE). We defined such periods as the temporal windows where at least 10% of cortical contacts demonstrated abnormal concurrent sharp peaks in more than half of the 18 frequency bands. Such episodes were excluded from within-and cross-frequency synchrony analysis. To identify such periods, we divided the signal in multiple 500 ms non-overlapping windows and detected IIE events in amplitude envelopes as the time point exceeding 5 times the standard deviation greater than the channel mean amplitude.

Defining the epileptic zones based on seizure activities in SEEG signals
The epileptogenic and seizure propagation zone were identified by clinical experts by visual analysis of the SEEG traces 1,3 . Epileptogenic areas are the hypothetical brain areas that are necessary and sufficient for the origin and early organization of the epileptic activities 4 , from where contacts recording often show low voltage fast discharge or spike and wave events at seizure onset.
Seizure propagation area are recruited during the seizure evolution, but they do not generate seizures 5,6 , from where contact recording show delayed, rhythmic modifications after seizure initiation in the epileptogenic areas. In this study, we combined epileptogenic and propagation areas as the epileptogenic zone (EZ) to distinguish from the rest of brain areas that are referred to as putative healthy zones (nEZ).

Functional connectivity estimates
We where T is the sample number of the entire signal (i.e., ~10 minutes), and * is complex conjugate.
We computed cPLV for the entire recording excluding 500 ms time windows showing epileptic or artefactual spikes (see below). The PLV is the absolute value of complex cPLV ( = | |), and it is a scalar measure bounded between 0 and 1 indicating absence of phase and full phase synchronization, respectively.
Additionally, we used imaginary part of cPLV ( = ( )), metric insensitive to zero-lag interactions caused by volume conduction [8][9][10] , for verification. For both PLV and iPLV connectivity, the fraction of significant edges (K) is the number of significant edges divided by the total possible edge number.

Statistical hypothesis tests
We estimated the null-hypothesis distributions of interaction metrics with surrogates that preserve the temporal autocorrelation structure of the original signals while abolishing correlations between two contacts. For each contact pair, we divided each contact's narrow band time series of into two blocks with a random time point k so that 1 ( ) = (1 … ) and 2 ( ) = ( … ), and constructed the surrogate as ( ) = [ 2 , 1 ]. We computed surrogate PLV across all channel pairs and assembled the surrogate interaction matrix, and its means and standard deviations were later used in hypothesis testing.

Correlation estimates post-processing
To demonstrate the changes of interaction strength as a function of spatial distance between recording sites, we divided the inter-contact distances into three ranges (short (SH) 2 cm ≤ x < 4.6 cm; medium (MD) 4.6 cm ≤ x < 6 cm, and long-range (LG) 6 ≤ x < 13 cm) with same number of edges falling into each distance bin. Therefore, we averaged across subjects all the edges falling within each distance bin and for each frequency band separately (N=48702/range). The confidence intervals for PLV and iPLV, were expressed relatively to the surrogate means (SM) for PLV (3.42*SM corresponding to p < 0.001, Rayleigh distribution), and the surrogate standard deviations (SD) for iPLV (3.58*SD corresponding to p < 0.001, normal distribution).
To compare signals from superficial and deep layers in the grey matter ( Fig. 4 and Suppl. Fig. 3), we divided the contacts into "shallow" and "deep" based their Grey Matter Proximity Index (GMPI) 2 that is defined as the relative distance between the contact location and the nearest whitegrey surface, normalized by the grey matter thickness at that location: where P(x, y, z), W(x, y, z) and C(x, y, z) are the vertices on the pial, white-matter surface and contact coordinates in 3D individually reconstructed brain from MRI scan, respectively. Values 0< GMPI < 1 indicate that the contact midpoint is located in grey matter whereas a negative GMPI indicates that the contact midpoint is in the white matter.
We used -0.3<GMPI<0 and 0.5<GMPI<1.2 to classify deep and superficial layer contacts respectively. Next, PLV and iPLV estimates were averaged across subjects between deep-deep (DD) and superficial-superficial (SS) contact-pairs. We tested for between groups difference with a paired permutation test (100 random samples created by shuffling SS and DD labels within subjects; threshold for significance corrected for multiple comparisons with Bonferroni p < 0.05/N, with N=18).

Anatomical co-localization of SEEG implant and Functional System characterization
To show that high-gamma phase coupling is associated with a common systems-level mechanism, we assessed the fraction of significant edge (K) between functional systems 11 . We extracted cortical parcels from pre-surgical T1 MRI 3D-FFE (used for surgical planning) using Freesurfer 12  discarded from this analysis. We thus assigned each cortical contact to a cortical parcel that belongs to a functional system. We then computed the fraction of significant edges (K) between 7 functional systems.
Note that we initially conducted this analysis in the Destrieux 148-parcel atlas 14 , but we reanalysed the whole dataset with the Schaefer atlas 13 after its release for both verifying the observations in the Destrieux atlas and for achieving optimal parcel-to-system morphing quality.

Cross-frequency coupling of slow rhythm phase and fast rhythm amplitude
Two signals of distinct rhythms are cross-frequency phase-amplitude coupled (PAC) if the phase of a slow neuronal oscillation modulates the amplitude fluctuations of the faster neuronal oscillations.
The rationale is that if the power fluctuations of fast rhythms are modulated by the phase of the slow oscillations, the fluctuations of these two time-series should be synchronized. In this study, we estimated PAC with the phase locking value (PLV) as: where ( ) is the phase time series of the power envelope of fast rhythm while ℎ ( ) is the narrow band phase time series of the slow rhythm. When there is a consistent relationship between these two time-series, the vector length of the mean phase differences (in the polar coordinate across all n samples) should be greater than zero, and a maximum value of 1 indicates perfect coupling. The significance of PAC PLV value was determined in the same manner in individual subjects as we conducted for 1:1 phase synchrony described earlier.

Detecting the modular structures
If observed individual level long-range phase synchrony between SEEG contacts are truthful, we next ask whether on the population level, the observed synchrony networks are functionally meaningful, i.e., the networks across frequencies, especially the HGA band, contains well defined modules of brain regions. To answer this question, we first unambiguously assigned each SEEG contact located in the cortex to one of the 100 Schaefer parcels 13 (subcortical contacts were not analyzed). For each frequency, we pooled PLV values of contact-pair originating from homologous parcel-pair across population, and then collapsed those PLV values in parcel-pair using the median of these contact-pair PLV values (using median due to the non-normality of PLV distribution). By doing so, we collapsed sparsely sampled patient contact-pair networks into well-sampled population-level phase synchrony networks between brain regions (see Fig. 1c). Due to low interhemispheric contact coverage, we limit our investigation to intra-hemispheric modular structures.
Moreover, given the sparse and non-homogenous sampling of cortical space across patients, 10% to 20% of all possible parcel-pairs have not been sampled, i.e., missing values (Fig. 1c).
Next, to identify modular structures in these synchrony networks, we employed a consensus clustering approach while missing values were accounted for 20 . Thereby, the observed modules of brain regions are less likely confounded by distributed local fragmentation of the network topology due to missing values. Briefly, the assumption is that the PLV values of the missing part of the networks follow the same probability distribution function of observed networks. For each observed phase synchrony network at any frequency, we generated 5,000 variants of the original network, wherein each missing edge was filled in with the PLV value of a randomly and independently selected existing edge. We next applied Louvain algorithm [21][22][23] to identify modules in each of the 5,000 variants. The range of the resolution parameter γ, which influences the number of communities identified, was set from 1 to 1.5 with 0.05 interval. The resolution parameter (γ) weights the importance of the null model (i.e., absence of modular structure) against which the original network is compared, when identifying the network partition or modules which maximise the modularity value: where ( ) is the modularity value at resolution parameter , is the total strength of PLV values, Next to assess the consensus in community assignments 20 across the 5,000 variants at any given resolution (γ), we derived community co-assignment (binary) matrix C for each of the variants, wherein M ij equals one when brain region i and j are demarcated to the same module otherwise zero. We then collapsed these 5,000 community co-assignment matrices into an average matrix and detected the modules therein with Louvain method with the same γ value used earlier to create the partitions. This consensus clustering method assumes that the set of 5,000 variants of the functional network are noisy versions of the true underlying network, and by averaging their partition assignments, the graph noise at individual level can be mitigated.
We also evaluated similarity between two community assignments solutions m and n across frequencies and Louvain resolutions. Here m and n are two vectors with each element associated with a cortical parcel, and the value of the element is the community ID assigned to that parcel.
With any given m and n, we could assemble the community co-assignment ( ) and ( ) as mentioned earlier. Thus, the partition similarity between m and n can be computed as 24 : . Since the dot product 〈 , 〉 satisfies the Cauchy-Schwartz inequality such that 〈 , 〉 ≤ �〈 , 〉〈 , 〉 , the value of partition similarity is 0 for uncorrelated partitions and 1 for identical partitions. We did not compare the similarity between left and right hemisphere of the brain because in the Schaefer atlas, the brain region demarcation are asymmetric across hemispheres.
We also assessed the range of resolution parameter values at which modules could be reliably identified. This was done both at the level of the entire network as well as at the level of individual regions (see Suppl. text).

PLV in amplitude bins
To assess whether larger values of phase synchrony were correlated with higher amplitude values, we estimated instantaneous amplitude and phase profiles of the filtered time-series by means of Hilbert transform. We then normalised each amplitude time-course by its median and divided the normalized amplitude samples in quintiles. For each frequency, we built an amplitude-amplitude correlation matrix containing the instantaneous phase difference between channel pairs of each time-samples at that amplitude bins. We discarded amplitude samples larger than twice the amplitude median to remove effects of subthreshold spikes. We quantified the number of timesamples falling in each bin as a simple amplitude correlation measure. We hypothesized that that, if a given contact pair is amplitude correlated, the time-samples would not be randomly distributed over amplitude bins. Indeed, it will result in highly skewed distribution of time-samples towards larger amplitude bins. On the other hand, in the absence of real amplitude correlation that distribution would be undistinguishable from a uniform distribution. Hence, to test for a moment-tomoment amplitude correlation, we quantified the distance of the time-sample distribution from a uniform distribution for each amplitude-amplitude bin under the above null-hypotheses of no correlation. To quantify whether phase consistency was correlated with moment-to-moment amplitude modulation, we quantified PLV in each amplitude-amplitude bin. The PLV is a measure sensitive to the sample number used 25 , hence we quantified the minimum number of samples falling in each bin and then quantified PLV in amplitude bin with matched time-samples. Specifically, for each contact pair, we computed instantaneous phase difference across the entire time course. By grouping instantaneous amplitude samples falling in same bin, we averaged phase differences in each amplitude bin and later averaged this quantity across channel pairs and subjects.