Mesoscopic physiological interactions in the human brain reveal small-world properties.

Cognition depends on rapid and robust communication between neural circuits spanning different brain areas. We investigated the mesoscopic network of cortico-cortical interactions in the human brain in an extensive dataset consisting of 6,024 h of intracranial field potential recordings from 4,142 electrodes in 48 subjects. We evaluated communication between brain areas at the network level across different frequency bands. The interaction networks were validated against known anatomical measurements and neurophysiological interactions in humans and monkeys. The resulting human brain interactome is characterized by a broad and spatially specific, dynamic, and extensive network. The physiological interactome reveals small-world properties, which we conjecture might facilitate efficient and reliable information transmission. The interaction dynamics correlate with the brain sleep/awake state. These results constitute initial steps toward understanding how the interactome orchestrates cortical communication and provide a reference for future efforts assessing how dysfunctional interactions may lead to mental disorders.

INTRODUCTION

Cognitive computations rely on the transformation of signals and communication between brain areas. Neuronal communication between brain areas correlates with computation and behavior in non-human animals (Bastos et al., 2015; Gregoriou et al., 2009; Womelsdorf et al., 2007) and humans (Baldauf and Desimone, 2014; Michalareas et al., 2016). Furthermore, computational models of cognition require communication between different processing stages (Zhao, 2016). However, neurophysiological investigations of brain function have largely focused on individual brain areas. At least partly, the emphasis on individual brain areas has been dictated by the use of single electrodes or multielectrode arrays that target a limited region of the brain. In humans, multiple single-neuron studies (Fried et al., 2014) and field potential studies (Schwiedrzik et al., 2018) have characterized the activity of individual brain areas in isolation, while few studies have measured correlations between two areas (Fu et al., 2019; Madhavan et al., 2019).

Anatomical connectivity patterns provide essential clues about the potential flow of information between brain areas. Long-range anatomical connections between brain areas have been reported through tract-tracing studies in monkeys (Felleman and Van Essen, 1991; Gămănuţ et al., 2018; Markov et al., 2014). Studies of anatomical connections in humans are challenging (Sparks et al., 2000): only a small number of post-mortem analyses have been conducted (Burkhalter and Bernardo, 1989), revealing connectivity over limited distances. Anatomical studies constrain the potential mechanisms of communication between brain areas, but they do not reveal whether two areas interact, nor how strongly, when, and in which frequency ranges those interactions occur. Non-invasive studies capitalize on the ability to study the human brain at a macroscopic level. Diffusion tensor imaging methods produce probabilistic maps of putative white matter tracts (Assaf and Pasternak, 2008). Additionally, scalp electroencephalography (EEG) (Nolte et al., 2004) and functional magnetic resonance imaging (fMRI) (Baldauf and Desimone, 2014; Glasser et al., 2016; Yeo et al., 2011) have provided indirect measures of macroscopic correlations.

While examining correlations between brain areas has provided important insights about computation in specific tasks, the large-scale study of interactions affords the possibility of evaluating communication in the whole network (Barabási and Albert, 1999; Bassett and Bullmore, 2017; Bullmore and Sporns, 2009; Graham, 2021; Muldoon et al., 2016; Sporns et al., 2005; Watts and Strogatz, 1998). In this study, we quantify the interactions between human brain areas at the mesoscopic level. We analysed an extensive dataset comprising intracranial field potentials (IFPs) from patients diagnosed with pharmacologically intractable epilepsy; the IFPs were continuously recorded for 0.9–10.4 days. IFP recordings show high signal-to-noise ratio signals with millisecond temporal resolution and millimeter spatial resolution (Buzsáki et al., 2012; Crone et al., 2006; Dubey and Ray, 2019). IFP studies have identified key neurophysiological features that correlate with cognitive and behavioral functions (Anumanchipalli et al., 2019; Dubey and Ray, 2019). We validated the methodology by comparing both anatomical and functional interactions from data in macaque monkeys.

Previous studies of mesoscopic networks examined resting state interactions over short periods spanning several seconds to minutes (Betzel et al., 2019; Fox et al., 2020; Keller et al., 2014; Solomon et al., 2017; Trebaul et al., 2018), except for one study that followed six subjects for almost a whole day each (Kramer et al., 2011). By examining interactions between 148,404 electrode pairs in 48 subjects during week-long continuous recordings, we provide an extensive neurophysiological characterization of long-range interactions between brain areas in the human brain. The resulting network shows that interactions in the human brain follow a log-normal distribution, are correlated with anatomical connections, and demonstrate small-world properties, whereby information could potentially travel from any one area to another rapidly and flexibly. Additionally, these brain interactions correlate with the cognitive state of the subject (sleep versus awake). This work is accompanied by a supplemental website where it is possible to scrutinize interactions in every electrode pair and every pair of brain areas, both in individual subjects and a common brain (http://braininteractome.com).

RESULTS

Definition of pairwise interactions

We examined IFP recordings from 48 subjects with pharmacologically resilient epilepsy (3–47 years old, 42% female, 90% right-handed, Table S1). We analyzed continuous IFPs recorded for 0.9–10.4 days (5.2 ± 1.8, here and throughout, mean ± SD) from 4,142 electrodes (86 ± 28 electrodes per subject, Table S1). The electrode positions in each subject are illustrated in Figure 1, and the anatomical areas covered are listed in Table S3.

Figure 1.

Location of electrodes in each of the 48 subjects

Locations with respect to the reconstructed pial surfaces of 4,142 electrodes (black dots) in 48 subjects. Arrow markers join the two physical electrodes that were used to compute each bipolar referenced signal. An average of 86 ± 28 electrodes was implanted per subject (minimum, 26; maximum, 176). The number of electrodes per subject is shown in Table S1, and a summary of electrode locations is shown in Table S3.

We implemented rigorous and conservative pre-processing steps to ensure that the recordings were free from potential artifacts (STAR methods): (1) bipolar referencing to remove volume conduction effects (Shirhatti et al., 2016) (we show the two electrodes that compose the signal throughout our figures); (2) custom algorithm to remove artifacts (Figures S1A–S1C; Table S2); (3) focus on seizure-free epochs; and (4) restriction to non-adjacent electrode pairs (Figures S1D and S1E) (Nunez et al., 1997). A detailed description of each of these steps and their importance to ensure the high fidelity of the recordings is presented in STAR methods.

We studied 148,404 electrode pairs during a total of 6,024 h. To explain the methodology step by step, we start by showing an illustrative example 10-s segment of IFP signals from two electrodes, one in the right superior temporal gyrus and one in the right pars opercularis (Figure 2A). We defined the relationship between the IFP signals from these two electrodes by calculating the degree of coherence, which ranges from 0 to 1 (STAR methods). Throughout this article, we refer to the coherence as a metric of interaction between the two signals (Bastos et al., 2015; Madhavan et al., 2019; Zhou et al., 2016) (see Discussion). We focus on interactions in the broadband frequency range (0.5–125 Hz, STAR methods; see supplemental website for results in other frequency bands). The coherence between the two signals in Figure 2A was 0.36. We computed coherence in non-overlapping windows over the whole duration of continuous recordings of up to 10.4 days (Figures 2B and 2C).

Figure 2.

Example recordings and definition of interactions

(A) Example pair of bipolar referenced intracranial field potentials (IFPs) during a 10-s time window from subject 3. One bipolar electrode pair was located in the right superior temporal gyrus (top), and the other bipolar electrode pair was located in the right pars opercularis (bottom). The coherence between these two signals was 0.36.

(B and C) Coherence was calculated in consecutive 10-s windows (represented by dots). A zoomed-in segment of 60-min duration (gray rectangle in C around the 10-s segment in A) is shown in (B). (C) Coherence was calculated for the entire ~5.5 days of recordings. (C) Segment of 24 h around the interval shown in (B) (note the two cuts in the x axis). The dark arrows in (B) and (C) mark the segment in (A). The horizontal dashed lines indicate the coherence under the null hypothesis (STAR methods).

Validation of pairwise interactions

We assessed the statistical significance of coherence measurements by comparing them with a null model whereby we randomly shifted one of the signals in time (Figure S2; STAR methods). Inserting random temporal shifts in one of the two electrodes ensured that the statistical properties of the individual signals remained intact while removing any temporal structure leading to interactions. The coherence in Figure 2A decreased from 0.36 to 0.05 following the shifting procedure (Figures S2A–S2D). This approach defined a significance threshold for the coherence for each electrode pair (horizontal dashed lines in Figures S2E and S2F).

We validated the use of coherence as a measure of neural interactions using anatomical and neurophysiological data from monkeys. We focused on monkey measurements for validation to leverage tract tracing studies of long-range anatomical connectivity (Markov et al., 2014) (Figure 3A) and known inter-areal neurophysiological interactions (Bastos et al., 2015; Zhou et al., 2016). We analyzed resting state neurophysiological data during the awake state from 128 intracranial electrodes implanted throughout the macaque cortical surface (Yanagawa et al., 2013), computing pairwise interactions for every pair of electrodes using the same analyses as for the human dataset (Figure 3B; STAR methods). Electrode positions were localized to 39 areas in the Markov-Kennedy parcellation (Markov et al., 2014).

Figure 3.

Validation using intracranial recordings in macaque monkeys

(A and B) Comparison between neurophysiological recordings and neuronal tracing studies in monkeys. (A) Known anatomical cortico-cortico connections between 39 macaque monkey brain areas based on the M132 atlas (Markov et al., 2014). Entry (i, j) in this matrix shows the log of the symmetrized normalized fraction of neurons (log10 FLNe) projecting from area i to area j (see Markov et al. [2014] for details). The color scale is shown on the right. (B) Summary of physiological interactions in the macaque monkey brain based on intracranial recordings from Yanagawa et al. (2013) for the same 39 brain areas shown in (A). The methodology to compute coherence is the same as the one described in the text for human data (STAR methods). Electrode positions were mapped onto the 39 corresponding brain areas used in the Markov-Kennedy connectivity atlas. Entry (i, j) in this matrix shows the degree of coherence between electrodes placed in area i and area j; the color scale is shown on the right. White entries denote location pairs where we did not have coverage in either anatomy (A) or neurophysiology (B). Gray entries denote location pairs where there were no significant anatomical connections (A) or physiological interactions (B). The Spearman correlation coefficient between the upper triangular matrix in (A), anatomical connections, and in (B), physiological interactions, was 0.28 (p = 1 × 10−8), controlling for interelectrode distances. The dendrogram in (B) displays the degree of similarity between areas according to hierarchical clustering on the coherence values (see Figure 5 and v). A shorter horizontal fork length indicates higher similarity.

(C–E) Mapping between human brain locations and macaque brain locations. (C) Lateral view (top) and medial view (bottom) of macaque cortical areas (reproduced from Figure 1 in Markov et al., 2014). Each color denotes a different brain area; some of the area names are labeled. (D) The pial surface from human subject 3 was mapped to the 91-area Markov-Kennedy parcellation for interspecies comparisons. (E) The surface in (D) cut and flattened. Colors in (D) and (E) correspond to the brain areas in (C). Inter-species registration was generated using the caret package (Van Essen et al., 2001), and flat maps were generated using the FreeSurfer package (Dale et al., 1999) (STAR methods).

(F–I) Known physiological interactions in the macaque monkey were found in human and monkey neural recordings. (F) Example of macaque intracranial electrode locations from two studies that reported physiological interactions between a pair of brain areas (Bastos et al., 2015; Zhou et al., 2016). The locations shown are the macaque electrocorticogram (ECoG) electrodes nearest to the originally reported locations in the corresponding studies. (G) Example of macaque intracranial electrode locations from three pairs of areas that are not anatomically connected according to the tracing study of Markov et al. (2014). Lines indicate the locations of bipolar references. (H) For each of the regions in the two studies in (F), we considered all possible electrodes in the macaque ECoG data (Yanagawa et al., 2013). We report here the proportion of electrode pairs that showed significant interactions (dark gray bars; see number of significant electrode pairs and total number of electrode pairs above each bar). Additionally, using the maps described in (C)–(E), we considered all of the electrodes in our human intracranial field potential recordings. We report here the proportion of electrode pairs in the human data that showed significant interactions (light gray bars). (I) Number of electrode pairs showing significant interactions for the pairs of areas with no known anatomical connection as shown in (G). A larger fraction of interacting electrode pairs was found in electrodes with previously reported interactions (H) versus in electrode pairs with reported absences of anatomical tracts (I), both in macaques and humans. Pairs of areas with a reported absence of anatomical tracts showed lower average coherence compared to pairs of areas that were known to interact.

As a first validation test, we compared the pairwise interactions with anatomical connectivity (Markov et al., 2014) (Figure 3A). There was a weak but significant correlation between the anatomical connectivity matrix (Figure 3A) and the physiological interaction matrix (Figure 3B, Spearman correlation coefficient = 0.28; p = 10−8, two-sided t test) after regressing out inter-electrode distances. Despite the differences between structural connectivity and physiological interactions, as well as the limited behavioral repertoire of the physiological recordings, the network described using the coherence metric approximated known long-range anatomical connections derived from tracing studies.

The second validation test examined transient interactions between brain areas while monkeys performed cognitive tasks, as positive controls to evaluate whether we could find evidence of those interactions using our metrics. We focused on two studies that documented interactions between visual areas V1 and V4 and between V1 and TEO (Bastos et al., 2015), and also between V4 and TEad/TEpd (Zhou et al., 2016). In the resting state data, we found 44, 54, and 15 electrode pairs that mapped onto V1-V4, V1-TEO, and V4-TEad/TEpd, respectively. Of these, 100%, 94%, and 100%, respectively, showed interactions according to our metrics (Figures 3F and 3H). These high percentages were remarkable given the methodological differences across labs and studies, and especially given that we compared transient interactions over hundreds of milliseconds during cognitive tasks (Bastos et al., 2015; Zhou et al., 2016) against interactions evaluated over hours of resting-state recordings.

In the third validation test, we evaluated these same interactions between monkey visual cortical areas in the continuous human neurophysiological data. We used putative homologies to map the locations of electrodes on the monkey cortical surface onto human cortical surface coordinates (Van Essen, 2004) (STAR methods; Figures 3C–3E). We found 248, 32, and 67 electrode pairs that mapped onto putative homologs of the V1-V4, V1-TEO, and V4-TEad/TEpd regions, respectively. Of these electrode pairs, 43%, 19%, and 4% demonstrated interactions in the human brain (Figures 3F and 3H). These percentages are lower than the ones reported in monkeys in the previous paragraph, especially the last value. Despite the extrapolation across temporal scales, across species, and across experimental conditions, the coherence metric was able to detect interactions that were consistent with those previously reported in the literature.

As a final validation test, we evaluated whether the coherence metric was valid when reporting the absence of interactions. Negative controls are more challenging to find due to potentially infrequent or weak interactions that may be difficult to detect (Markov et al., 2014), and due to reporting biases favoring positive interactions. As a coarse approximation to a negative control, we considered pairs of areas that do not show anatomical connectivity (Markov et al., 2014), that is, V1-area 2 (not to be confused with area V2), V1-area 5 (not V5), and V4-F5. We found 72, 27, and 20 pairs of electrodes, respectively, in the monkey resting-state neural data that mapped onto those locations. Of those electrode pairs, 4%, 4%, and 0% showed significant interactions (Figures 3G and 3I); these percentages are well below the 100%, 94%, and 100% reported in Figures 3F and 3H. We also found 27, 34, and 6 electrode pairs in the human data that mapped onto these locations. None of those electrode pairs showed significant interactions (Figures 3G and 3I), in stark contrast with the 43%, 19%, and 4% values reported in Figures 3F and 3H.

In sum, coherence constitutes a reasonable metric to assess interactions between brain areas: the values were unlikely to arise by chance (Figure S2), the interactions correlate with anatomical measures (Figures 3A and 3B), and interactions measured from both human and monkey recordings are consistent with positive (Figures 3F and 3H) and negative (Figures 3G and 3I) neurophysiological and neuroanatomical examples from previous studies.

Mesoscopic interactions change over time and are broadly distributed

We return to the 10-s IFP segment for the example electrode pair in Figure 2A. Other 10-s segments (Figure S3) for the same electrode pair show coherence values ranging from no statistical significance (Figure S3A) to typical coherence values (Figure S3D), to strong interactions (Figure S3C). Coherence measurements for a 1-h segment are shown in Figure 2B, and coherence measurements for a 24-h segment within the whole period of 5.5 days are shown in Figure 2C. We considered the fraction of all 10-s windows that showed statistically significant coherence values as a measure of temporal consistency (STAR methods). Throughout the 1-h segment shown in Figure 2B, this electrode pair showed significant coherence 39% of the time. During 5.5 days, the temporal consistency was 41%, with a coherence averaged over significant segments of 0.22 ± 0.05. Here, we focus on this time-averaged coherence value as a summary metric and we return to the temporal variations in coherence in Figure 7. We define a pair of electrodes to show significant interactions when the temporal consistency was above 5%.

Figure 7.

Brain interactions change with sleep and wake states

Cross-validated accuracy of a nearest neighbor classifier trained to decode whether the subject was awake or asleep using coherence values from pairwise interactions in each of 14 subjects (STAR methods). Error bars denote standard deviation. Dotted lines indicate the permutation significance threshold. Accuracies were significantly above chance (chance = 0.5) in all 14 subjects. **p < 0.001, permutation test, corrected for multiple comparisons.

After examining the temporal consistency of physiological interactions, we evaluated the degree of spatial specificity. We asked whether the interaction illustrated in Figure 2A was spatially restricted by considering IFPs from other locations. Figures S3E and S3F show three other electrodes that were simultaneously recorded with the electrode pair in Figure 2. None of the other six possible coherence measurements between the two 10-s segments shown in Figure 2 and the three simultaneous 10-s segments in Figure S3F reflected significant interactions. The maximum coherence for these six non-significant pairs was 0.13, compared to 0.36 for the interaction between the example electrode pair in Figure 2. Of the 91 bipolar electrodes for this subject, there were 44 electrodes where we could evaluate potential interactions with the two electrodes in Figure 2 after applying the exclusion criteria for both electrodes. Of these 44 electrodes, 8 (18%) showed significant coherence with one but not the other electrode, 22 (50%) showed significant coherence with both, and 14 (32%) did not show significant coherence with either (including the 3 electrodes in Figure S3F). In sum, both electrodes showed spatially restricted interactions with only a specific subset of the other simultaneously recorded electrodes. This subset included broadly distributed interactions.

We extended the study of spatial specificity from these example 10-s segments to time-averaged coherences calculated over the entire duration of recordings. In the same electrode pair as in Figures 2A and S4A, we illustrate in Figure S4B all interactions for the electrode in the superior temporal gyrus (Figure 2A, top). This electrode interacted with 20 of the 75 electrodes remaining after applying the distance exclusion criterion (27%). Similarly, Figure S4C shows the 23 significant interactions out of 60 electrodes (38%) for the electrode in the pars opercularis (Figure 2A, bottom). In sum, both electrodes showed spatially selective interactions, typically targeting a broad range of other simultaneously recorded electrodes.

A mesoscopic network of interactions between human brain areas

Next, we considered all possible electrode pairs. Figure S4D (right) depicts all significant interactions between the 91 electrodes for this example subject. Entry i, j in this symmetric coherence matrix denotes the interaction between electrode i and electrode j. For example, row 35, column 84 shows the example electrode pair in Figure 2A, the whole row 35 shows the interactions in Figure S4B for the superior temporal gyrus electrode, and the whole column 84 shows the interactions in Figure S4C for the pars opercularis electrode. White entries denote neighboring electrode pairs excluded from analyses (Figures S1D and S1E). Gray entries denote pairs that did not reach statistical significance (Figure S2). Of the 4,095 possible electrode pairs (91 choose 2), 3,193 pairs (78%) remained after applying the distance exclusion criterion (Figures S1D and S1E), and 557 of these pairs (17%) showed significant interactions.

We mapped each of the 91 electrodes (Figure S4D, left) to its corresponding brain area in the 36-area Desikan-Killiany parcellation (Desikan et al., 2006). The symmetric coherence matrix in Figure S4E shows the average coherence for all pairs that map onto the corresponding brain areas (18 areas with coverage). For example, row 5, column 13 shows the average of all pairs between the pars opercularis and the superior temporal gyrus, including the one in Figure 2A. The hierarchical clustering of the coherence matrix (Figure S4E, right) reflects the similarity between brain areas based on coherence values, where shorter horizontal distances mean higher similarity (STAR methods).

We extended the analyses to the full cohort of 48 subjects. Using the same procedure illustrated for subject 3 in Figures 2 and S4, we present 48 maps showing all of the interactions in each individual subject in the supplemental website (STAR methods). Of the 148,404 electrode pairs across all subjects, we analyzed 124,865 pairs (84%) after applying the distance exclusion criterion. We found a total of 11,692 significant pairwise interactions (9.4%). Across all subjects, the average fraction of significant pairwise interactions was 10% ± 15%.

The variation in the fraction of interactions across subjects arises from the heterogeneity in electrode locations dictated by clinical needs (Figure 1; STAR methods). To compare results between subjects, we considered two complementary approaches to map electrode locations across brains: (1) mapping all electrode positions onto 36 Desikan-Killiany areas (Table S3) (Desikan et al., 2006) as illustrated for subject 3 in Figure S4E, and (2) mapping all of the electrodes onto a custom parcellation based on spatial clustering of electrode locations (Figures S5A and S5B; STAR methods).

To illustrate the methodology, we follow the two areas in Figure 2, the superior temporal gyrus and the pars opercularis. Thirty out of 48 subjects had electrode coverage in both areas. We found at least one pair of interacting electrodes in these areas in 10 out of these 30 subjects (33%). All interacting pairs between these two areas are shown in Figure 4. Of the total of 448 electrode pairs spanning the two areas in these 10 subjects, 49 (11%) showed physiological interactions. The average coherence between the superior temporal gyrus and pars opercularis was 0.30 ± 0.13. Note that these parcellations reflect large brain areas, and we expect heterogeneity in the electrode positions within each area across subjects. Despite the variance in electrode positions across subjects, we consistently found significant interactions between these two regions.

Figure 4.

The interaction between superior temporal gyrus and pars opercularis was consistent across subjects

Ten subjects (subject numbers shown next to each brain) where there was at least one significant pair of electrodes showing significant interactions between electrodes mapping onto the superior temporal gyrus (cyan) and pars opercularis (beige). Bipolar electrode pairs are indicated by the droplet-shaped markers. The black markers denote the bipolar electrodes participating in significant interactions. The color of the line joining any two dots indicates the degree of coherence (see color scale on bottom). The arrow on the color bar shows the average degree of coherence across all of these electrode pairs. The white arrows denote the electrodes shown in Figure 2A for subject 3.

We applied the same analyses to the rest of the electrode pairs. The analyses were restricted to 31 of the 36 possible areas due to lack of coverage in 5 areas (Table S3; STAR methods), resulting in 465 area pairs (31 choose 2). There was adequate electrode coverage in 378 area pairs, of which 193 (51%) showed significant interactions after applying thresholds for consistency, as illustrated for the example interaction between superior temporal gyrus and pars opercularis (STAR methods). These interactions are summarized in the symmetric matrix in Figure 5. For example, row 3, column 10 in this matrix reflects the interactions between electrodes in the pars opercularis and superior temporal gyrus, after averaging over all interacting electrode pairs in these two areas as shown in Figure 4.

Figure 5.

Mesoscopic interactome of the human brain

This symmetric matrix shows all significant subject-averaged coherence values between areas in the Desikan-Killiany parcellation. Pairwise interactions between anatomical areas defined in the Desikan-Killiany parcellation (Desikan et al., 2006) are shown. There was electrode coverage in 31 of the 36 defined areas (STAR methods). Of the 31-choose-2 = 961 pairs of areas, 180 (19%) did not have electrode coverage in our data (white squares). Of the 781 covered pairs of areas, 345 (44%) interactions were statistically significant, and the color shows the average coherence (see color scale on the right). Interactions shown are significant in at least two subjects, and across at least 10 bipolar electrode pairs. Gray squares denote interactions that were not significant. The interaction between the superior temporal gyrus and pars opercularis (Figures 2 and 3A) showed an average coherence of 0.30 ± 0.02 and can be seen at row 3, column 10. The hierarchical clustering dendrogram on the right shows the relationship between Desikan-Killiany areas based on their similarities in coherence with other areas (STAR methods).

Of those area pairs that showed significant interactions, the weakest interaction was between the inferior temporal gyrus and superior parietal gyrus (#7 and #24 in Figure 5 and 0.145 ± 0.003 coherence), and the strongest interaction was between the parahippocampal gyrus and caudal middle frontal gyrus (#1 and #13, 0.563 ± 0.072 coherence).

The distribution of interaction strengths followed a log-normal distribution (Figure S6A, p = 0.59, Kolmogorov-Smirnov test) rather than a normal distribution (p = 0.04, Kolmogorov-Smirnov test). The supplemental website (Video S1) renders all of the interactions in a 3D brain that can be rotated to scrutinize the location and coherence of each electrode pair. The paracentral gyrus had the least number of interactions (#18 in Figure 5): only electrode pairs where both coordinates mapped within the paracentral gyrus showed significant interactions. The fusiform gyrus showed the most interactions (#5 in Figure 5): 23 of 29 (79%) areas with electrode coverage interacted with the fusiform gyrus.

The areas in the interaction matrix shown in Figure 5 reflected a common parcellation of the human brain based on anatomical landmarks (Desikan et al., 2006). We also implemented a custom parcellation dictated by electrode coverage. We mapped electrodes from each of the 48 subjects onto a common average surface (Figure S5) and used K-means clustering to form 150 non-overlapping areas (Figure S5B; STAR methods). The pair of electrodes in Figure 2 mapped onto areas 111 and 148 in this custom parcellation. Area 111 was one of 10 areas in the custom parcellation that overlapped with the superior temporal gyrus. Area 148 was one of 3 areas that overlapped with the pars opercularis. We described the interaction in Figure 4 in terms of the custom parcellation: within the 10 × 3 = 30 possible area pairs between these two sets of custom parcellation areas, 25 had electrode coverage. Of these 25 custom area pairs, 18 (72%) showed significant interactions, with an average coherence value of 0.30 ± 0.10 that was comparable to the one reported for Desikan-Killiany areas in Figure 4 (0.30 ± 0.13).

We described the interactions from the whole set of 11,175 custom area pairs (150 choose 2). Of the 4,644 custom area pairs with electrode coverage, 2,387 (51%) showed significant interactions (Figure S5C). Interaction strengths ranged from 0.13 (between areas 121 and 143, both overlapping with the precentral regions in the Desikan-Killiany map) to 0.79 (between areas 137 and 140, overlapping with the precentral and rostral middle frontal regions in the Desikan-Killiany map). As described for the anatomical map, the distribution of interaction strengths followed a log-normal distribution (Figure S6C). The area that showed the most interactions with other areas was area 19, which overlaps with the middle temporal gyrus in the Desikan-Killiany map. Area 19 interacted with 82 of the 111 regions with coverage.

Our dataset contained subjects from 3 to 47 years old (Table S1). We evaluated whether there was a difference in the brain interactome between young and older subjects by splitting the data using the median age of 17 into two subgroups of ages 11.2 ± 3.9 and 26.9 ± 10.2. We re-calculated physiological interaction networks in each subgroup (Figure S7). Of the original 193 significant brain area pairs, there were 140 (73%) brain area pairs that were found to be significant in both age subgroup datasets. We found a significant difference in coherence values between the two groups (p = 2 × 10−17, two-sided rank-sum test). Similarly, for the gamma-band coherence, of the original 184 significant brain area pairs, 125 (68%) were found to be significant in both age subgroups. There was also a difference in gamma band coherence values between the two groups (p = 2 × 10−19, two-sided rank-sum test).

The interactions reflect small-world network properties

The interaction matrix (Figure 5) can be thought of as a network graph composed of nodes (brain areas), connected by weighted edges (coherence values). As noted above, the number of edges per node (i.e., the number of other areas that any brain region interacted with) was not uniformly distributed. Therefore, we sought to describe the network structure of the interaction matrix by asking whether it exhibited small-world properties (Watts and Strogatz, 1998). In a small-world network, connections between nodes are neither distributed uniformly nor randomly; instead, the connectivity includes “neighborhoods” with high local connectivity and sparser interconnectivity between neighborhoods, enabling traversing from one node to any other node in a small number of steps. Using the reciprocal of the coherence as a proxy for the length between two nodes, we measured the mean shortest path length (L) and the clustering coefficient (C) of the interaction matrices (Bassett and Bullmore, 2017) (STAR methods). Small-world networks are characterized by short path lengths and high clustering coefficients. We calculated these metrics in both the human (Figure 6) and monkey (Figure S8) brains.

Figure 6.

The mesoscopic interactome showed small-world network properties

(A) Graph of neural interactions based on 20 nodes corresponding to areas in the Desikan-Killiany parcellation after removing areas with missing values (see Table S3 for abbreviations). The color indicates the coherence value (see scale on the right).

(B) Random graph of neural interactions obtained by re-assigning coherence weights from the network in (A).

(C) Lattice graph of neural interactions. The networks in (B) and (C) have the same number of nodes and edges as the actual network in (A). The small-worldcoefficient (Humphries and Gurney, 2008) was σ = 1.2 (95% CI: 1.1, 1.3), satisfying the range σ > 1. The small-world measurement (Telesford et al., 2011) was ω = −0.41 (95% CI: −0.42, −0.39), satisfying the range −0.5 < ω < 0.5 (STAR methods).

(D) Schematic rendering of each of the Desikan areas, adapted from the Desikan-Killiany parcellation (Desikan et al., 2006).

We compared the L and C values against two null models: random graphs (rand, Figure 6B), and lattice graphs (latt, Figure 6C). Both null models have the same number of nodes and edges as the original network but they show different topology. We computed the small-world coefficient, σ, defined as the ratio between the clustering coefficient and the mean shortest path length, both normalized by the corresponding values for the random graph (σ = C/C /L/L) (Humphries and Gurney, 2008). For the human brain, we obtained σ = 1.20 (95% confidence interval (CI): 1.14, 1.26); for the macaque brain, we obtained σ = 1.15 (95% CI: 1.09, 1.20). Both species revealed small-world networks, which are characterized by σ > 1. Another commonly used small-world metric, u, is defined as L/L – C/C (Telesford et al., 2011), and it takes values between −0.5 and +0.5 for small-world networks. For the human brain, we obtained ω = −0.406 (95% CI: −0.417, −0.394); for the macaque brain, we obtained ω = −0.188 (95% CI: −0.195, −0.182). A summary of all computed small-world measures is shown in Table S4. Both the human and monkey interaction networks deviated from random and lattice graphs and were consistent with small-world network properties found in other neural systems (Bassett and Bullmore, 2017). Furthermore, we evaluated small-world properties in 13 individual subjects who had sufficient coverage: all 13 subjects showed small-world properties based on σ, and 11 of 13 subjects showed small-world properties based on ω (Table S5). In sum, efficient communication via small-world networks was evident at the individual subject level as well as in the whole interactome, in humans and monkeys, and in broadband and gamma-band signals.

The interactions between brain areas correlate with cognitive state

The interactions described in Figures 3, 4, 5, and 6 show a temporal average that removes the rich dynamic fluctuations in coherence. However, the coherence values between any two given electrodes changed over time (Figures 2B, 2C, and S3). Such dynamic fluctuations in brain interactions might reflect ongoing changes in cognitive state and function. As a proof-of-principle evaluation of whether the brain interactions reported here correlate with ongoing computations, we analyzed sleep, an important and salient aspect of cognitive state and function (Figure 7).

We considered a random subset of 14 subjects, and annotated whether the subjects were awake or asleep. These annotations were performed continuously over the entire dataset for those subjects based on the video, and independently validated by manual annotations from two sleep experts using EEG data on three subjects (STAR methods). We used a nearest neighbor machine learning classifier with 20-fold cross-validation to decode whether subjects were awake or asleep using the coherence values between all pairs of electrodes as features. On average, the cross-validated accuracy was 0.66 ± 0.11 (chance = 0.5). This analysis was separately performed for each of the 14 subjects. We found that we were able to distinguish sleep versus awake states above chance levels in all subjects (Figure 7). There was variation in the decoding performance across subjects, ranging from 0.53 ± 0.003 (subject 39) to 0.85 ± 0.034 (subject 11). The variation among subjects is largely governed by the different electrode locations (Figure 1).

DISCUSSION

Cognition emerges through rapid communication between brain areas. In this study, we provide an extensive, population-based, and dynamic description of the mesoscopic interactome in the human cortex over a continuous period of up to 10.4 days. Interactions were defined by coherence measurements between IFPs recorded from pairs of bipolar-referenced electrodes (Figure 2) in 48 subjects (Figure 1; Table S1). We achieved ample coverage with 4,142 electrodes and 148,404 electrode pairs mapping onto 31 of 36 areas in the Desikan-Killiany parcellation (Figure 1; Table S3) (Desikan et al., 2006). Interactions were spatially localized (Figures 3B, 4,5, S3E, and S3F), fluctuated over time (Figures 2A and S3), and were consistent across subjects (Figure 4). We found 11,682 pairwise interactions between electrodes, summarized by 193 interactions between 31 Desikan-Killiany areas (Figure 5), and 2,387 interactions between 150 areas in a custom parcellation (Figure S5). The physiological interactome is accompanied by a comprehensive supplemental website showing each interaction, in individual and average brains, in different frequency bands, using multiple visualization tools (http://braininteractome.com). Coverage of the human cortex through the 4,142 electrodes was extensive (Figure 1; Table S3), but not exhaustive. Because mapping electrodes onto an average brain necessarily reduces the precision of reported locations, the supplemental website depicts all interactions in each subject’s own brain without any distortions. The networks in humans and monkeys showed small-world properties (Figures 6 and S8) and correlated with brain state (Figure 7).

We analyzed continuous recordings during the course of 5.2 days per subject, on average. We observed a significant difference in the interactome when considering only 10-min recordings (STAR methods), demonstrating that the inclusion of coherence measurements over long timescales was necessary to arrive at the coherence network results.

Stringent precautions were taken to avoid contamination from electrical noise, movement artifacts, seizures, and volume conduction, leading to conservative estimates of the number and magnitude of interactions (STAR methods). A caveat of these restrictions, especially the bipolar referencing and 17-mm minimum spacing (Figures S1D and S1E), is the exclusion of short-distance interactions.

We used the coherence between IFPs to measure interactions (Nunez et al., 1997). This metric has been used extensively to assess neural communication (Baldauf and Desimone, 2014; Bastos et al., 2015; Bullmore and Sporns, 2009; Fries, 2005; Gregoriou et al., 2009; Madhavan et al., 2019; Zhou et al., 2016). We validated the coherence metric against anatomical connections and physiological interactions (Figure 3). Additionally, the interaction matrix and small-world network properties were robust to the choice of the interaction metric (Table S6) and were also detected in individual subjects (Table S5).

Importantly, note that all of the work herein reflects recordings in patients with epilepsy and all electrode positions were dictated by clinical needs. Seizure-related activity was not included in the analyses (STAR methods) and the overall methodology was independently compared against monkey neuroanatomy and neurophysiology (Figure S3). Despite these precautions and validation, it is conceivable that some (or none, or all) of the conclusions could be different in non-epilepsy subjects. With current techniques, it is not possible to perform this type of long-term continuous monitoring for a week of invasive neurophysiological data at high spatiotemporal resolution from healthy subjects.

Physiological interactions did not reflect global activity changes but instead were locally specific (Figures S3E and S3F). Coherence measurements revealed spatially restricted subnetworks of interactions (Figures 3B, 4, 5, and S4D). We refer to the spatiotemporal scale studied here as mesoscopic, distinguishing these measurements from finer resolution in anatomical connections between individual neurons (Kasthuri et al., 2015) and coarser-resolution interactions evaluated by EEG (Nunez et al., 1997), magnetoencephalography (MEG) (Michalareas et al., 2016), or fMRI (Baldauf and Desimone, 2014; Glasser et al., 2016; Yeo et al., 2011).

The broadband coherence included signals between 0.5 and 125 Hz. Multiple studies described interactions in narrower frequency bands, particularly in the gamma range (Bastos et al., 2015; Gregoriou et al., 2009; Madhavan et al., 2019). We repeated all analyses for the following four bands: 3–8 Hz (theta), 8–12 Hz (alpha), 12–30 Hz (beta), and 30–100 Hz (gamma) (supplemental website). Filtering the data in these different frequency bands does not assume oscillatory activity to be the origin of our measures (Miller et al., 2009), but rather it describes potentially different bands of communication between areas. The general principles of spatial specificity, temporal fluctuations, and small-world properties were consistent across the different frequency bands. However, the fine details of the interaction networks depended on the frequency band. These differences open up the intriguing possibility of multiplexing different routes of communication within the network (Buzsaki, 2006).

The millisecond temporal resolution of IFP recordings was critical to uncovering rapid interactions. Indeed, most previous studies have focused on short-lived interactions during behavioral tasks (Bastos et al., 2015; Gregoriou et al., 2009; Madhavan et al., 2019; Zhou et al., 2016). We observed a wealth of interesting temporal features in our coherence metric including a baseline coherence level that fluctuated slowly over the course of many hours, gradual shifts in the baseline over the course of days, accompanied by brief modulations in coherence (Figure 2). The brief modulations are correlated with whether the subject was in a sleep versus awake state (Figure 7). The coherence results in Figures 3, 5, and S5 focused on averaged interactions over long periods, which encompassed a diverse spectrum of temporal features (Figure 2) (see also discussion in Liégeois et al. [2017] and Zalesky et al. [2014]).

Subjects were bedridden; awake behaviors included rest (Glasser et al., 2016; Yeo et al., 2011), social interactions, head and body movements, sitting, and watching TV. Most of the “housekeeping” and pervasive cognitive functions, including visual/auditory processing, language, and internal thought processes, are reflected in the interactome. Infrequent brain interactions contribute to a lesser extent to the current results. As an initial proof of principle to evaluate the relationship between the dynamic changes in brain interactions and cognitive function, we compared the sleep versus awake states (Figure 7). Using only the brain interactions assessed by coherence measurements, it was possible to decode the cognitive state of the subject (awake versus sleep). The initial results introduced in Figure 7 suggest that the brain interactions that we quantify in this study can discriminate the aggregate functions involved during different cognitive states.

Coherence values do not necessarily reflect direct interactions. Anatomical connectivity correlated well with interactions as measured by coherence (Figures 3A and 3B). However, significant coherence values do not imply direct anatomical connectivity (Petersen and Sporns, 2015). Conversely, anatomical connectivity does not imply physiological interactions either. The coherence values are not directional, and no claim is made here about the causality of the interactions (Reid et al., 2019). One exciting possibility to assess causality is to quantify the consequences of current injection via electrical stimulation followed by examination of the evoked cortico-cortical potentials (Fox et al., 2020; Trebaul et al., 2018). One of these direct stimulation studies revealed a network that partially exhibited small-world characteristics (Keller et al., 2014), consistent with the results reported herein.

The interaction strengths follow a log-normal distribution (Figure S6), in line with log-normal distributions in anatomical tract-tracing network weights (Buzsáki and Mizuseki, 2014; Ercsey-Ravasz et al., 2013; Gámánuţ et al., 2018; Markov et al., 2014). Such distributions are a key property of hierarchical networks where few weaker long-distance projections are amplified by local circuits that make up most of the connections (Markov et al., 2013b). The current results extend this notion to mesoscopic physiological interactions.

In small-world networks, it is possible to rapidly convey information from any point in the network to any other point in a small number of steps, and edges are configured economically to facilitate both global and local integration (Bullmore and Sporns, 2012; Telesford et al., 2011). The coherence weights used for our network analyses are distinct from measures used in other modalities and in anatomical networks (Muldoon et al., 2016). Despite these differences, the combination of local and long-range connectivity gave rise to a network topology revealing small-world characteristics (Table S4; Figures 6 and S8), consistent with properties of anatomical connectomes (Bassett and Bullmore, 2017). Small-world networks efficiently organize information flow by facilitating both local and global connectivity (Watts and Strogatz, 1998). We speculate that the coherence network characterized herein may facilitate information transmission across the cortex (Telesford et al., 2011). The existence of weaker, long-range interactions has been shown to be critical for anatomical network specificity (Markov et al., 2013a). The coherence metric captured these critical “shortcut” interactions.

The results presented in the present study provide initial steps toward the mesoscopic characterization of the flow of information in the human brain. The physiological interactome should not be thought of as static, let alone immutable. Interactions between brain areas can be modified over developmental timescales, and much more rapidly during cognitive tasks. The network structure described herein provides a reference map for future studies of brain physiology and pathology.

STAR★METHODS

Detailed methods are provided in the online version of this paper and include the following:

RESOURCE AVAILABILITY

Lead contact

Further information and requests should be directed to and will be fulfilled by the lead contact, Gabriel Kreiman (Gabriel.Kreiman@childrens.harvard.edu).

Materials availability

This study did not use or generate any reagents.

All source code is publicly available on Github (https://github.com/kreimanlab/MesoscopicInteractions).

All the coherence data have been deposited at http://braininteractome.com.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Human subjects

We studied 48 subjects with pharmacologically intractable epilepsy (20 female, 43 right-handed, 24 right hemisphere coverage, aged 3–47 years old, 19 ± 11, here and throughout we report mean ± standard deviation, Table S1). All procedures were performed with the subjects’ informed consent. Recruitment for the study was conducted according to protocols approved by the Institutional Review Boards at Boston Children’s Hospital, Harvard Medical School (35 of 48 subjects), and Brigham and Women’s Hospital, Harvard Medical School (13 of 48 subjects).

METHOD DETAILS

Neurophysiological recordings

We studied 48 subjects with pharmacologically intractable epilepsy (20 female, 43 right-handed, 24 right hemisphere coverage, aged 3–47 years old, 19 ± 11, here and throughout we report mean ± standard deviation, Table S1). Because the majority of subjects were right-handed, we did not have enough data to directly compare brain interactions in right-handed versus left-handed subjects. The recordings were conducted as described previously (Liu et al., 2009). Briefly, subjects with pharmacologically resilient epilepsy were recruited based on their need for long-term invasive monitoring by surgically implanted electrodes. All procedures were performed with the subjects’ informed consent. Recruitment for the study was conducted according to protocols approved by the Institutional Review Boards at Boston Children’s Hospital, Harvard Medical School (35 of 48 subjects), and Brigham and Women’s Hospital, Harvard Medical School (13 of 48 subjects).

Subjects were implanted with platinum electrodes (Ad-Tech, Racine, WI, USA; 2.3 mm diameter exposed area, with inter-electrode distances of either 5 mm (64 of 4,142 electrodes) or 10 mm (4,078 of 4,142 electrodes), impedance < 1kOhm). The location and number of electrodes varied across patients as determined by clinical needs. There were 26 to 176 electrodes implanted per subject (86 ± 28, Table S1). Electrodes consisted of either grids or strips and were implanted in the subdural cavity. We refer to the electrical signals from these electrodes throughout the text as the intracranial field potential (IFP). IFP signals were continuously recorded for 0.9 to 10.4 days (5.2 ± 1.8 days, Table S1), with a sampling rate of 250 Hz, 256 Hz, 500 Hz, 512 Hz, 1000 Hz, 1024 Hz, 2000 Hz, or 2048 Hz using Natus NeuroWorks software (Oakville, Ontario, Canada).

Electrode localization

We collected T1-weighted magnetic resonance imaging (MRI) images of the head before electrode implantation for each subject. The pial surface of each subject was then reconstructed from image sets using the freesurfer package (Dale et al., 1999). We followed the recon-all pipeline using the fully-automated option (“-autorecon-all” flag). Post-operation computed tomography (CT) images of the implanted electrodes were localized to the MRI reconstruction using the iELVis package (Groppe et al., 2017). CT gantry tilt was corrected using the dcm2niix package from https://www.nitrc.org. Electrode grid and strip orientation was identified in the CT scan based on pre-surgical sketches and platinum marker guides (Ad-Tech, Racine, WI, USA). Figure 1 shows the resulting electrode locations for all subjects. Electrode placement and coverage varied across individuals, and this variation can be appreciated from the reconstructed coordinates in Figure 1. All subjects except subject 26 had unilateral electrodes. All analyses were limited to within-hemisphere interactions only. The accompanying supplemental website provides three-dimensional renderings of the electrode locations for each subject (http://braininteractome.com).

In some cases (Figure S5; supplemental website), we rendered electrodes from different subjects onto an average brain. Registration of electrodes from individual subjects onto the average brain was done using iELVis, followed by cross-hemispheric registration using freesurfer to register electrodes from the fsaverage right and left hemisphere surfaces onto the symmetric fsaverage_sym right hemisphere surface. Electrodes from all 48 subjects on a common pial surface (Figure S5A) revealed extensive coverage over a large expanse of the cerebral cortex surface. Electrode positions were measured on a three-dimensional coordinate system and were used to calculate inter-electrode distances for neighborhood identification (Figure S1E).

We summarized the electrode locations by mapping them onto one of 36 brain areas based on the parcellation of the Desikan-Killiany atlas (Desikan et al., 2006) using the recon-all function in freesurfer. Table S3 provides summary information about the number of electrodes in each of these 36 areas. We also present results under a custom parcellation based on electrode coverage (Figure S5, described next).

Custom brain parcellation

The Desikan-Killiany parcellation is based on anatomical landmarks denoted by brain sulci and gyri. When describing electrode locations, we also used an alternative custom parcellation map based on the density of electrode locations. We constructed a 150-area parcellation of locations that best matched the electrode coverage (Figure S5C; supplemental website). This parcellation was constructed by first considering the 3,615 bipolar electrodes across all subjects and their coordinates on the fsaverage_sym pial surface (Figure S5A). Distances between these coordinates were calculated and used to define regions based on electrode coverage (Figure S5B). We defined 150 areas, a number approximating that of the multimodal parcellation of brain regions in (Glasser et al., 2016). Following the reconstruction procedure computed by the freesurfer package (Dale et al., 1999), 3,615 bipolar electrodes were given coordinates in the corresponding subject spaces. These coordinates were then transformed into the fsaverage_sym space using the iELVis package (Groppe et al., 2017). The 3,615 by 3,615 distance matrix was calculated by taking the arclength of all pairwise combinations of the spherical electrode coordinates u and v, using freesurfer spherical reconstructions (Dale et al., 1999). The freesurfer default radius of 100 mm for spherical surfaces was used as the value for r to calculate the distance: Using this distance matrix, each of the 3,615 bipolar electrodes was assigned to one of 150 clusters to form sets of coordinates delineated into custom brain areas. The k-means clustering algorithm was computed using the function kmeans in MATLAB. Multiple runs of the algorithm were used to overcome the noise introduced by random initializations. The final version was taken as the clustering with the minimum total within-cluster variance after running the algorithm for 900 random initializations (Figure S5).

In addition to clustering bipolar electrodes into 150 regions, we further defined region membership for the other vertices making up the fsaverage_sym surface. Each vertex of the fsaverage_sym surface was given membership in one of the 150 regions by considering the 3-nearest-neighboring bipolar electrodes. The majority cluster out of 3 was assigned to that vertex. The distance was computed using the same arclength definition as before. The rendering of the custom parcellation as shown in Figure S5B revealed 150 regular and non-overlapping regions. Regions contained 5 to 78 bipolar electrodes, with an average of 24 ± 14. Regions measured between 8.8 mm to 28 mm at their widest dimension, with a mean distance of 7.7 ± 1.9 mm across.

Signal preprocessing

Several precautions were taken to remove potential artifacts from raw IFP recordings. Here we describe each of these preprocessing steps.

We aimed to measure local signals to investigate mesoscopic interactions between brain regions. Volume conduction effects could potentially lead to inflated measures of interactions (Nunez et al., 1997). We erred on the side of caution and removed signals that may arise from volume conduction. We maximized the locality of the signal by re-referencing the IFP following a bipolar montage, which measured voltages with respect to adjacent electrodes. The bipolar montage has been shown to reduce the cortical spread of underlying neural signals (Dubey and Ray, 2016; Kajikawa and Schroeder, 2011). Additionally, bipolar referencing helps remove sources of global artifacts by subtracting signals coming from the ground reference. The bipolar IFP response measured local neural activity, with a reported spatial resolution of up to 3 mm (Bansal et al., 2014; Dubey and Ray, 2019; Vidal et al., 2010). Coherence measurements from bipolar-referenced IFP signals were less prone to false-positive reports of physiological interactions (Kajikawa and Schroeder, 2011). We refer to the resulting voltage signal, i.e., the difference between voltage in two physical electrodes, as the signal coming from a “bipolar electrode.” Throughout the main text and figures, we used the term “electrodes” to refer to “bipolar electrodes” for simplicity.

Bipolar IFP signals were defined by subtracting a pair of raw IFPs from adjacent electrodes within a strip or grid. The definition of bipolar electrodes from physical electrodes followed the clinical convention. Let V(t) represent the voltage measured at time t at the physical electrode i (i = 1, …, n). Voltage at bipolar electrode j was defined as V(t) = V(t) – Vi+1(t) for j = 1, …, n−1. Grid electrodes were first split into individual strips before following the same procedure as defined for strips. In total, the number of bipolar electrodes per subject equaled the number of physical electrodes minus the number of strips, after converting grids into strips (Table S1).

We assumed that neighboring electrodes before implantation remained neighbors post-implantation. We checked this assumption by measuring the distances between physical electrodes using their electrode localization coordinates. For the 46 out of 48 subjects with a pre-implantation interelectrode distance of 10 mm, the average distance measured between physical electrodes defined as a bipolar electrode pair was 10 ± 2 mm (n = 3,440, Figure S1E). The variance in the inter-electrode distance was mainly attributed to the flexibility of the electrode housing material and precision of electrode localization. The interactions between pairs of bipolar electrodes reported throughout this study encompass 4 physical electrodes (two bipolar electrodes); we refer to this arrangement as a bipolar electrode pair. The average number of physical electrodes per subject was 86 (range: 26 – 176). After applying the bipolar reference, the average number of bipolar electrodes per subject was 75 (range: 20 – 154).

When mapping electrodes to a parcellation of the brain (e.g., Figure S4), bipolar electrodes were localized based on the positions of the two physical electrodes. If both physical electrodes mapped onto the same area, the bipolar electrode was assigned to that area (2,101 out of 3,615 bipolar electrodes). Otherwise, the bipolar electrode was mapped based on the primary physical electrode.

The recording system removed DC biases. Interference from power line noise was removed by using a fifth-order Butterworth notch filter at 60 Hz, 120 Hz, and 180 Hz with a bandwidth of 3 Hz (Liu et al., 2009). In subjects where IFPs were recorded with a sampling rate of 250 Hz, the power line interference of the third harmonic at 180 Hz was apparent as the aliased frequency of 70 Hz. For subjects with sampling rates of 256 Hz, the apparent frequency of the third power line oscillation was 76 Hz. Notch filters with centers of 180 Hz were adjusted to their corresponding apparent center frequencies according to the sampling rates of the signals being filtered. Interference at 20 Hz from medical devices was detected in some of the IFP recordings (Dzwonczyk et al., 2009), and was also removed by applying the same notch filter with a center frequency of 20 Hz and a bandwidth of 4 Hz. To ensure that this filtering did not impact the coherence, frequency samples between (17, 23) Hz were excluded. The analyses in Figure 3 refer to data recorded from macaque monkeys (Yanagawa et al., 2013) (see “Comparisons with macaque monkey neurophysiological recordings”). In those recordings, the power line noise frequency was 50 Hz, and the notch filter centers were adjusted accordingly.

We maintained a standardized sampling rate throughout analyses. We used the lowest sampling rate available for analysis over the whole dataset: 250 Hz. All recordings were down-sampled to 250 Hz after the data filtering steps. In 10 of 48 subjects, sampling rates were integer multiples of 256 Hz. In these cases, the target for down-sampling was 256 Hz. The difference between subjects with 250 or 256 Hz sampling rates was accounted for throughout our analyses and did not impact any of the results. We found these sampling rates to be adequate in achieving a high temporal resolution (4 ms).

To simulate the effect of coarser temporal resolution, we repeated the analysis of the interaction network for subject 3 after applying a Gaussian filter (scipy.ndimage.gaussian_filter1d, sigma = 0.1 s (Jones et al., 2001)) to bipolar IFP recordings and found that only 58 (2%) out of 3,193 bipolar electrode pairs showed significant interactions. This number was much smaller than the original 557 (17%) pairs found to be significant, demonstrating that 0.1 s is too coarse to adequately assess interactions and emphasizing the necessity of recordings with high temporal resolution.

We removed artifacts from our dataset

The IFP recordings may contain artifacts of non-physiological origin, which may affect the analyses, especially when considering such long recording sessions as the ones used in this study. Some of those artifacts may persist even after bipolar referencing. Bipolar IFPs were separated into one-second segments to determine the presence of potential artifacts. We used a series of conservative criteria to define 6 types of potential artifacts (three of which are illustrated in Figure S1):

Large amplitude. The strongest evoked IFP signals that other groups and we have characterized span several hundred μV (Liu et al., 2009; Tang et al., 2014). A one-second segment was considered to contain a potential artifact if the voltage range exceeded 2,000 μV: max(()) − min(()) > 2000 μ.

Large slope. A one-second segment was considered to contain a potential artifact if the magnitude of the slope of voltage versus time exceeded 100 μV/ms, |() / | > 100(μ / = ), at any time during the segment. This slope was much larger than the range of neurophysiological slopes observed in previous IFP studies (Liu et al., 2009; Tang et al., 2014).

Small amplitude. A one-second segment with a voltage range of less than 10 μV was also marked as a potential artifact: max(()) − min(()) < 10 μ

Examples of these types of artifacts are shown in Figure S1. As a conservative measure, any segment with a potential artifact was removed from further analyses. These three criteria led to the removal of 0.72%, 0.44%, and 0.1% of the data, respectively.

Seizure events and interictal discharges. Clinical experts defined seizure onset and offset events as well as interictal discharges. Recordings from all channels starting 15 minutes before seizure onset and ending 15 minutes after seizure offset were removed from analyses. An average of 0.68 ± 1.81% of the data was marked as seizures.

Electrical stimulation events. Subjects received external cortical stimulation as part of the procedure to functionally map areas relevant for surgery. All the cortical stimulation events were removed from the analyses. An average of 0.41 ± 0.17% of the data was marked as stimulation events.

Electrodes with many artifacts. Individual bipolar electrodes that had > 50% of the data segments marked as potential artifacts were not considered for further analyses. This criterion removed a total of 0.34% of the 3,615 bipolar electrodes (Table S2).

In total, 2.1% of the data were marked as potential artifacts (out of 1.6 × 109 one-second segments). Individual fractions of artifacts marked for each criterion are listed in Table S2. Coherence values (see “Definition of interactions”) were computed over 10 s windows. Thus, there were 20 one-second segments for each coherence window, considering the two bipolar electrodes. We only computed coherence measurements for a pair of bipolar electrodes when none of these 20 one-second windows contained potential artifacts.

We used bipolar referencing to improve spatial resolution and minimize artifacts (Dubey and Ray, 2016; Kajikawa and Schroeder, 2011). To further reduce potential false positives, we removed interactions between neighboring bipolar electrodes from analyses.

First, we considered same-grid neighbor electrodes. Two electrodes were neighboring if they were adjacent on strips and grids, or diagonally adjacent on grids (Figure S1D). A bipolar electrode pair involved a total of 4 physical electrodes, i.e., two pairs of physical electrode pairs. We represented this relationship schematically in Figure S1D by drawing the relationships between the midpoint of physical electrode pairs. A bipolar electrode pair was determined to be neighboring if any of the 4 possible physical electrode pairs were neighboring. After applying this definition to the 148,404 bipolar electrode pairs across 48 subjects, 16,451 (11%) same-grid pairs were removed from analyses.

Next, we considered different-grid neighbor electrodes. For bipolar electrodes residing on different strips or grids, neighbors were defined as follows. We considered the 16,451 same-grid neighbor electrodes. A histogram of the electrode distances for this set is shown in Figure S1E. This distance was used to determine an optimal threshold for neighborhood classification. A threshold of 17 mm was sufficient to achieve a true-positive rate of 92.1% and a false positive rate of 1.8%. After applying this threshold to bipolar electrode pairs from different grids, 7,088 (4.8%) additional bipolar electrode pairs were removed from analyses (shaded gray area under the solid line in Figure S1E).

The combined fraction of bipolar electrode pairs removed due to both same- and different-grid electrodes was 15.9%, leaving 124,865 pairs for subsequent analyses. This distance constraint impacted our study of physiological interactions by removing short-distance interactions while ensuring high sensitivity in reported coherence measurements (see Discussion).

Comparisons with monkey neurophysiology

This study focused on human neurophysiological data. As an independent validation, we also analyzed similar IFP recordings obtained in a macaque monkey from the study of Yanagawa and colleagues (Yanagawa et al., 2013). These recordings consisted of resting intracranial field potentials from 128 electrodes similar to the ones used in humans, covering 42 of the 91 areas in the Markov-Kennedy atlas (Markov et al., 2014; Yanagawa et al., 2013).

The same data preprocessing steps described above for the human data were used on the monkey IFP recordings. To account for the 5 mm interelectrode distance of the monkey intracranial electrodes, the electrode neighborhood exclusion threshold was decreased from 17 mm to 8 mm. The periods including and between the administration of anesthesia (time = 40 minutes), and the administration of the anesthesia antagonist (time = 85 minutes) were removed from analyses. Method parameters were taken as originally documented (Yanagawa et al., 2013). Signal preprocessing excluded 12% of the recordings, while the anesthesia condition removed an additional 39% of the recordings. In total, we analyzed 450 minutes across 112 bipolar electrodes and 6,216 bipolar electrode pairs. Electrodes were localized by referencing locations provided by the study of Yanagawa and colleagues (Yanagawa et al., 2013) to the Markov-Kennedy areas (Markov et al., 2014).

We compared human and macaque brain locations using the following steps. The M132 parcellation of the macaque monkey cortex (Markov et al., 2014) was applied to human subjects using the surface-based interspecies registration method (Van Essen and Dierker, 2007) from the software package caret (Van Essen et al., 2001). The interspecies algorithm transformed each of the 91 areas in the M132 macaque parcellation into areas represented on the coordinate space of the human Population Average Landmark and Surface-based (PALS) atlas (Van Essen, 2005). These areas were then mapped from the human PALS atlas to the human fsaverage atlas (Dale et al., 1999) by caret_command from the caret package. For each human subject, the putative macaque area homologs were mapped onto the individual coordinate space by mris_apply_reg from the freesurfer package (Dale et al., 1999). The calculations for registering each subject to the average surface was conducted with default parameters set by freesurfer and using automatic reconstructions of surfaces as described in the electrode localization methods (see “Electrode localization”). Figures 3C and 3E shows a comparison between original M132 areas on the surface of the monkey brain (Figure 3C) and the putative interspecies homologs mapped on the surface of the brain of human subject 3 (Figures 3D and 3E).

We compared physiological interactions calculated from monkey IFPs to tractography connectivity information from the study of Markov et al., 2014) (Figures 3A and 3B). We used the extrinsic fraction of labeled neurons (FLNe) metric to quantify anatomical connections. We collapsed the directed matrix by combining FLNe values in the forward and reverse directions, averaging the directed values to create the symmetric connectivity matrix. A missing value connection was defined as the absence of connectivity between pairs of areas where only one direction was tested. These missing values were omitted from the symmetric matrix averaging. Cortical distances were predictive of connectivity strength (Markov et al., 2014). Thus, we regressed out the Euclidean interelectrode distances when computing the partial correlation between the anatomical connectivity network (Figure 3A) and the coherence network (Figure 3B). The correlation coefficient obtained without taking into account the interelectrode distances was 0.31; the correlation coefficient after incorporating the interelectrode distances was 0.28 (Figure 3). Both of these correlation coefficients were significantly different from zero (p = 10−10 and p = 10−8, respectively).

We compared the results with known physiological interactions. Electrode locations from two studies of known physiological interactions (Bastos et al., 2015; Zhou et al., 2016) (Figures 3F and 3H) were taken from their respective sources directly. Electrodes were localized to areas defined by the Markov-Kennedy parcellation (Markov et al., 2014).

Definition of physiological interactions

We considered pairwise interactions between two bipolar electrodes in a segment of duration T = 10 s. We used segments of T = 10 s duration; however, the millisecond temporal resolution of IFP signals was essential for the analyses. Similar results were observed when using T = 5 s or T = 15 s to calculate coherence.

Let V(t) and V(t) represent the voltage in bipolar electrodes a and b at time t (0 ≤ t ≤ T). The autocorrelation of each voltage signal and the associated auto-spectral density are defined as: and (similarly for bipolar electrode b). We omit spelling out the duration T in the variables for simplicity. The cross-correlation between the voltage signals and the associated cross-spectral density are defined as: and . From these variables, we define the coherence between the two signals as (Nunez et al., 1997): The real-valued coherence ranges from 0 (no relationship between the two signals) to 1 (for two identical signals).

Several studies have documented physiological interactions in narrow frequency bands (Bastos et al., 2015), such as the gamma band (Gregoriou et al., 2009; Madhavan et al., 2019). Coherence was separately integrated over five conventional frequency bands: broadband (0.5 ≤ f < 125Hz), theta (3 ≤ f < 8Hz), alpha (8 ≤ f < 12Hz), beta (12 ≤ f < 17Hz) and (23 ≤ f < 27Hz), and gamma (30 ≤ f < 100Hz). To ensure that power line filtering did not have any effect on coherence calculations, a bandwidth of ± 4 Hz (56 – 64 Hz), which was larger than the notch filter used in preprocessing of ± 3 Hz, was used to exclude confounding frequency samples. This exclusion removed artifacts arising from imperfect notch filtering. When the frequency interval included 60 Hz or harmonics (50 Hz and harmonics for the monkey data), a window of ± 4 Hz was removed from the analyses. Throughout the main text, we focus on the broadband coherence analyses. Results for the other frequency bands are shown in the supplemental website.

The coherence was computed using Welch’s method of overlapping segments. The input parameters were chosen to minimize the mean of the null distribution of the magnitude of coherence computed over permuted time-shifted IFP signals. The input parameters were 0.2 s segment width and 80% segment overlap. These parameters resulted in a spectral resolution of 5 Hz and were used for computing the broadband and gamma coherence. For the theta, alpha, and beta coherence, a segment width of 2 s and 50% overlap was used to achieve 0.5 Hz spectral resolution. Given the large amounts of data processed in this study, computations were accelerated by writing a graphics processing unit (GPU) coherence function based on Welch’s method (Welch, 1967) using the package reikna (Opanchuk, 2014). The 95-percentile of the difference in the magnitude of coherence between the GPU-accelerated implementation and scipy.signal.coherence (Jones et al., 2001) was 0.00073 (n = 15,336 randomly sampled time segments). There was a 16-fold decrease in processing time compared to scipy.signal.coherence, bringing down the analysis time from months to weeks per run.

The coherence value does not necessarily reflect direct interactions between the two areas. Anatomical connectivity correlates with physiological interactions (Figures 3A and 3B). However, a significant coherence value does not imply direct anatomical connectivity between the two areas (Petersen and Sporns, 2015). Conversely, direct anatomical connectivity does not imply a physiological interaction either.

The coherence calculations reported here are not directional, and we do not make any statements about the causality of the pairwise interactions. The incidence of reciprocity in axonal tracts in animal models is variable, with near-complete reciprocity in mice (Gam anut ‚ et al., 2018), and lower reciprocity in macaque monkeys (Markov et al., 2014). Functional differences due to anatomical directionality have been reported in the macaque monkey visual cortex (Bastos et al., 2015; Nassi et al., 2013). However, the construction of directional anatomical connectivity maps and the extent to which these asymmetries may impact function in the human cerebral cortex is the subject of ongoing investigation (Teillac et al., 2019, Organization for Human Brain Mapping, conference; Van Essen et al., 2013).

Hierarchical clustering

We defined a distance metric for hierarchical clustering. For each area i, the n-dimensional vector (i) = {1(i), 2(i), …, (i)} contained the coherence values between that area and each of the other areas of the parcellation of interest. The distance d(j,k) between two areas j and k was calculated by taking the Euclidean distance between vectors c(j) and c(k). Missing coherence values with insufficient coverage were omitted from the distance calculation.

Hierarchical clustering was computed according to Ward’s method (Ward, 1963). The resulting hierarchical clustering was visualized as a dendrogram. The position on the horizontal axis corresponded to the relative value of the within-group variance metric, where shorter distances indicated better group fit. Hierarchical clustering was performed using the following functions in MATLAB: linkage, optimalleaforder, dendrogram.

Hierarchical clustering provided an initial step for examining interaction networks (Figures 3, 5, S4E, and S5). High-level clusters mapped approximately onto the temporal-parietal, lateral prefrontal, medial prefrontal, and occipital lobe areas. Brain areas with similar functional localization also had similar physiological interactivity patterns. The clustering algorithm was not given anatomical information a priori, yet we found correspondence between the clustered order and anatomic proximity. For example, the Precentral gyrus and Postcentral gyrus were clustered together (see areas 14, 15 in Figure 5). These neighboring areas are related both functionally and anatomically (Van Essen et al., 2013).

Network graph analyses

We considered the graph formed by brain areas as nodes and interactions as edges. We performed several calculations described in this section to characterize the structure of this interaction network graph, both in monkeys (Figure S8A) and in humans (Figure 6A).

In monkeys, the Markov-Kennedy areas (Markov et al., 2014) served as the nodes. Of the 91 areas in the atlas, 39 areas were covered by intracranial electrodes (Figure 3B). In humans, the Desikan-Killiany areas served as the nodes (Figure 5). Of the 36 areas in the atlas, 31 areas were covered by intracranial electrodes.

We converted networks into edge-complete networks by removing nodes with incomplete coverage. The node with the highest number of missing values was identified and removed from the adjacency matrix. This operation removed all edges associated with the identified node. The procedure was repeated iteratively until no missing values were present. An 18 by 18 matrix remained in the macaque physiological interaction network (Figure S8), and a 20 by 20 matrix remained in the human network (Figure 6).

To investigate whether the network graphs showed small-world properties, we calculated the characteristic path length (L), and the clustering coefficient (C) (Humphries and Gurney, 2008; Muldoon et al., 2016; Telesford et al., 2011; Watts and Strogatz, 1998). Two small-world indices were derived from the values of L and C for both the original networks and generated random and lattice networks described below. Physiological interaction network edges were constructed from the magnitude of coherence, ranging from 0 to 1. We used the extended definitions of L and C suitable for weighted networks as defined by (Muldoon et al., 2016), as implemented using the accompanying MATLAB code, with modifications. The weighted clustering coefficient C was computed using the definition given by (Onnela et al., 2005). To compute the characteristic path length L, the coherence values c between nodes i and j were converted to path lengths w = 1 / c, following the definition from (Newman, 2001). The resulting distance matrix was the input to compute path lengths using the MATLAB function graphallshortestpaths. The magnitude of the path length between any two nodes was calculated as the sum of the individual distances of the edges making up the shortest path between them. The path length was calculated for all pairs of nodes, excluding self-paths, and the average of the resulting values represented the weighted characteristic path length. If a node was disconnected from the network (no edge connected it to any other node), it was excluded from the average calculation.

We defined small world network indices as follows. We normalized the L and C metrics using the following definitions. The small-world coefficient (Humphries and Gurney, 2008) was defined as σ = γ/λ, where γ = C/C and λ = L/L. σ > 1 indicate networks with small-world properties. We also considered the following alternative small-world measurement (Telesford et al., 2011), defined as ω = L/L – C/C. Values −1≤ω≤1 indicate networks with small-world properties.

We computed small-world coefficients in individual subjects when individual networks had adequate coverage and the network was reachable (13 of 48 subjects). Coherence values between Desikan-Killiany areas, as shown for subject 3 in Figure S4E, were used to compute small-world network coefficients. Following the same missing-value removal procedure as applied to the 31-node subject-averaged coherence value matrix (as shown in Figure 5), individual subject networks were converted to edge-complete networks. There were 5 of 48 subjects where this conversion resulted in an empty matrix and these subjects were omitted from further analysis. After omitting subjects without sufficient electrode coverage, an additional 30 subjects did not have sufficient interactions for calculating both σ and ω. The path length, clustering coefficient, s and u was calculated for the remaining 13 individual subject networks, where the number of nodes ranged from 4 to 12 (Table S5).

Sleep annotations and decoding

We manually annotated parts of the continuous video recordings of the patients in order to assess whether the brain interactions correlated with cognitive state and behavior (Figure 7). We selected a random subsample of 14 out of the 48 subjects to continuously annotate whether subjects were awake or asleep. These annotations were based on the video recordings. We also selected a subset of 4 out of these 14 subjects who had scalp EEG recordings to validate the annotations. Two certified sleep clinicians performed sleep and wake annotations based on the scalp EEG recordings in Subjects 5, 8, 14, and 21 for a total of 6 nights of sleep, 47.3 hours. We annotated periods of sleep and wake based on audiovisual recordings in the same subset of data and found good agreement with clinical experts with an overall accuracy of 0.90, true positive rate of 0.96, true negative rate of 0.64, and an inter-rater reliability of 0.66 (Cohen’s kappa). We then annotated sleep based on audiovisual recordings for a total of 1362.4 hours in the 14 subjects (97.3 ± 36.0 hours per subject).

QUANTIFICATION AND STATISTICAL ANALYSIS

Statistical analyses and null hypothesis

To assess whether random voltage fluctuations could explain the coherence values between two signals, we defined a null hypothesis whereby one of the signals was randomly shifted in time, thus preserving all its structure except for the synchronous interactions. A permutation test for each bipolar electrode pair was performed by randomly shifting IFP signals in time: V(t) V(t + τ). The time shift τ was chosen uniformly without replacement between 2 minutes and the upper bound set by the length of the recording (Table S1). The time window used for time-shifted coherence was the same window used in the other analysis, that is, 10 s. This permutation procedure was run 10,000 times per bipolar electrode pair. The resulting null distribution was fitted to a t-distribution using the MATLAB fitdist function (The MathWorks, Inc., Natick, MA). A coherence threshold was determined for each bipolar electrode pair by taking the inverse cumulative distribution function of the fitted t-distribution, using the MATLAB icdf function. The p value used for icdf was corrected for multiple hypotheses testing by controlling the family-wise error rate (FWER = 0.05). Coherence values above this threshold were marked as significant interactions (Figures 2B and 2C). Further constraints were imposed by the temporal consistency and between-subject consistency, described next. In the total of 124,865 bipolar electrode pairs, the resulting coherence thresholds had a median of 0.22 ± 0.18 (broadband), 0.42 ± 0.19 (theta band), 0.25 ± 0.30 (alpha band), 0.28 ± 0.26 (beta band), and 0.15 ± 0.16 (gamma band).

Temporal consistency

Coherence values were computed for every segment of size T over the entire dataset (Table S1). For each pair of bipolar electrodes, given the number of segments where coherence was above the significance threshold, n, and the total number of segments, N, we defined the temporal consistency as n/N. The coherence between two electrodes was defined to be consistent over time if the temporal consistency was above 5%. This threshold was well above (i.e., more conservative), than the one obtained by considering a Bonferroni-corrected binomial distribution based on each time segment. In a simulated dataset of coherence values generated by randomly sampling from the fitted null distribution, only 1 bipolar electrode pair was temporally consistent out of 124,865 bipolar electrode pairs across the 48 subjects. Thus, the method was confirmed to be robust at controlling the Type-I error rate.

To summarize the interactions occurring through time for each bipolar electrode pair, we defined the time-averaged coherence value to be the average coherence of those n segments where coherence was above the significance threshold (Figures 2B and 2C). The exclusion of coherence values below the threshold ensured a high-quality metric at the cost of potentially missing weak interactions with low coherence magnitudes.

Between-subject consistency

To summarize physiological interactions across subjects, it was necessary to reconcile electrode maps from one subject to another. We considered two types of anatomical maps for these analyses: (i) the locations of all electrodes were registered onto a standardized cortical parcellation (Desikan et al., 2006) (see “Electrode localization”; e.g., for subject 3, the 91 bipolar electrodes mapped onto 18 areas, Figure S4E); (ii) electrodes were mapped onto our custom 150-area parcellation (see “Custom brain parcellation”).

Once each bipolar electrode was assigned to a brain area, we combined results from multiple bipolar electrode pairs mapping onto the same pair of areas. In the example shown for subject 3 in Figure 4, there were 8 bipolar electrodes in the Superior Temporal Gyrus and 3 bipolar electrodes in the Pars Opercularis, yielding a total of 24 bipolar electrode pairs between these two areas. For a given pair of areas i and j, the symmetric matrix M(i,j) indicates the total number of bipolar electrode pairs. The size of M was 36 by 36 for the Desikan-Killiany parcellation (Figure 5), and 150 by 150 for the custom parcellation (Figure S5). Only 31 of the 36 areas had electrode coverage (Figure 5; Table S3). We followed the same procedure for the monkey data, resulting in a 39 by 39 matrix (Figure 3B).

Bipolar electrode pairs mapping onto the same area were also considered. In some cases, significant interactions occurred within brain areas, as shown on the diagonal of Figure 5. These interactions were possible when an area pair involved an electrode pair where both electrodes localized to the same area. The distance criterion was applied to all analyzed electrodes (Figure S1D–E); thus, these intra-areal interactions mostly occurred in areas with a large surface. The interactions reported here represented interactions between sub-areas and do not necessarily represent a type of internal correlation. The interaction networks calculated using the custom parcellation (Figure S5C) did not contain these intra-areal interactions. The larger number of areas resulted in smaller surfaces that did not accommodate electrode pairs that could pass the distance criterion (Figure S5B).

We combined the area-wise interaction matrices across subjects. Given the number of bipolar electrode pairs with significant interactions m(i,j) for a given pair of areas i and j, and the total number of bipolar electrode pairs M(i,j), the between-subject consistency was defined as m(i,j)/M(i,j). The significance of the between-subject consistency was determined in a similar way as for the temporal consistency. A threshold of between-subject consistency ≥ 10% was used to define significant interactions for a pair of brain areas. An analysis of simulated coherence values using the null hypothesis was conducted to validate the methodology. The method reported no pairs of brain areas showing significant interactions in this null dataset.

The placement of intracranial recording electrodes varied from subject to subject (Figure 1; Table S3). To minimize potential sampling bias when combining multiple subjects, we restricted the analyses to M(i,j) ≥ 10. We further restricted analyses to cases where bipolar electrode pairs came from at least 2 subjects.

The coherence values from significant electrode pairs mapping onto each brain area pair were averaged to give the area-pair coherence (color-mapped values in Figure 5). We performed this averaging for each of the 193 significant between-area interaction sets, each of which contained 401 ± 518 electrode pairs, and 14 ± 9.2 unique subjects. The standard deviations, numbers of significant and total pairs, and the numbers of subjects spanned by these pairs is shown in the supplemental website. For the gamma frequency band, in each of the 184 significant interactions between areas, there were 377 ± 515 electrodes pairs and 13 ± 10 unique subjects. For each of the 2,387 significant interactions between areas in the custom parcellation (Figure S5C), there were 23 ± 18 electrode pairs and 7.4 ± 3.9 unique subjects. For the gamma frequency band, in each of the 2,110 significant interactions between areas in the custom parcellation, there were 23 ± 14 electrode pairs and 7.4 ± 4.0 unique subjects. In the simulated dataset where coherence values were generated by randomly sampling from the fitted null distribution, no significant interactions were found in the broadband coherence or other frequency bands.

We assessed whether the coherence values from the interaction networks shown in Figures 5 and S17 for the Desikan-Killiany areas (Desikan et al., 2006), and the networks shown in Figure S5C for our custom parcellation, were fit by a log-normal distribution. These distributions are shown in Figures S6A and S6B for the Desikan-Killiany areas and in Figures S6C and S6D for the custom parcellation. Areas with insufficient electrode coverage or nonsignificant interactions were not considered for distribution fits. A one-sample Kolmogorov-Smirnov test was performed for the empirical observations using the MATLAB function kstest after normalizing coherence values according to a standard normal distribution. The resulting p value was interpreted as significantly deviating from a standard normal distribution if p < 0.05. In cases where we performed the test against a log-normal distribution, we performed the kstest algorithm after taking the logarithm of coherence values. We found no significant interactions (as shown in Figure 5) when we sampled coherences from the null distribution, thus there was no opportunity for log-normal assessment for the distribution of coherence values under the null distribution.

Robustness of the subject-averaged coherence

We assessed the robustness of our subject-averaged coherence metric by splitting our dataset in half based on the electrode pairs, calculating the physiological interaction network individually for the two halves, and comparing the results. Of the 193 significant brain area pairs reported from the original 124,865 bipolar pairs (see Results), there were 186 (96%) brain area pairs that were found to be significant in both half datasets. The difference in coherence values between the two halves was 0.005 ± 0.048 and was not significantly different (p = 0.78, two-sided rank-sum test, n = 186). The gamma-band frequency coherence was also found to be robust. Of the 184 significant brain area pairs, there were 174 (95%) brain area pairs that were found to be significant in both half datasets. The difference in gamma coherence values between the two halves was 0.001 ± 0.055 and was not significantly different (p = 0.75, two-sided rank-sum test, n = 186).

As an additional control, instead of splitting the data by electrode pairs, we further assessed the robustness of our subject-average coherence metric by splitting the dataset in half by time windows. For each of 48 subjects, coherence time windows were segmented into 4-hour epochs. These epochs were split evenly to form two derived datasets, each receiving half the amount of data from the full dataset. Of the 193 significant brain area pairs reported for the full dataset, there were 186 (96%) brain area pairs that were found to be significant in both half datasets. There was no significant difference between the coherence values in the two halves (p = 0.88, two-sided rank-sum test, n = 186, absolute difference of 5.4% ± 6.1%). The gamma frequency coherence was also found to be robust. Of the 184 significant brain area pairs, there were 179 (97%) brain area pairs that were found to be significant in both half datasets. There was no significant difference between the gamma band coherence values in the two halves (p = 0.79, two-sided rank-sum test, n = 179, absolute difference of 3.3% ± 3.7%).

The current study examines a large dataset encompassing 6,024 hours of neurophysiological recordings. To assess whether short recordings could be sufficient to elucidate brain interactions, we considered 130 random segments of 10-minute duration (encompassing 60 segments of 10 s duration) for each of the 48 subjects. Of the 193 significant brain area pairs reported for the full dataset, there were 142 (74%) brain area pairs that were found to be significant in both the short segment and in the complete data. The coherence values in the complete data and short samples were significantly different (p = 5·10−4, two-sided rank-sum test, n = 142, absolute difference of 17.2% ± 14.2%). The gamma frequency band coherence showed similar results: of the 184 significant brain area pairs, there were 120 (65%) brain area pairs that were found to be significant in both the short segments and in the complete data. The gamma band coherence values in the complete data and short segments were significantly different (p = 8·10−6, two-sided rank-sum test, n = 120, absolute difference of 19.5% ± 13.3%).

Both the number of electrodes (Table S1) and the regions covered (Table S3) varied widely across subjects. To assess the effect of individual subjects on the results, we considered the 193 between-area interactions shown in Figure 5 and randomly removed n subjects from the analysis. After removing all bipolar electrode pairs of the selected subjects, each of the 193 area pair coherence averages were re-calculated using data from the remaining subjects. We then examined the interactome matrix (Figure 5), the log-normal distribution, and the small-world network properties.

After 50,000 random removals and re-calculations, there was a high correlation between the original 193 significant coherences and the mean of the random re-calculations (r = 1.00, p = 0.00, n = 193) for n = 1. This correlation remained high (r = 0.99, p = 0.00, n = 193) for n = 24 when half of the subjects were randomly removed. For the broadband signals, the maximum absolute difference between the original 193 significant coherences and the mean of the random re-calculations was 0.003 for n = 1 (this maximum occurred for the interaction between inferior Parietal Gyrus and Pars Orbitalis and represented 0.9% of the original coherence value of 0.36 between these two areas), and 0.06 for n = 24 (this maximum occurred for the interaction between the Cuneus and Parahippocampal Gyrus and represented 16.5% of the original value of 0.37 between these two areas). In the gamma frequency band, the maximum absolute difference between the 184 significant gamma coherence interactions and the mean of the random re-calculations was 0.002 for n = 1 (this maximum occurred for the interaction between the Pars Orbitalis and Frontal Pole and represented 0.9% of the original coherence value of 0.23 between these areas), and 0.04 for n = 24 (this maximum occurred for the interaction between the Pars Orbitalis and Supramarginal gyrus and represented 9.7% of the original coherence value of 0.37 between these areas).

The distributions of coherence values remained log-normal after re-calculation. The deviation from a normal distribution at n = 1 was significant (p = 0.033, Kolmogorov-Smirnov test), while a deviation from a lognormal distribution was not (p = 0.59, K-S test). For n = 24, the deviation from a normal distribution was no longer significant (p = 0.064, K-S test), but the deviation from a lognormal distribution was not significant (p = 0.56, K-S test). Lognormality was less robust in the gamma coherence values, where for n = 1, the K-S test p value was 0.052 for deviation from a normal distribution, and 0.12 for deviation from a lognormal distribution. For n = 24, the K-S test p value was 0.091 for deviation from a normal distribution, and 0.23 for deviation from a lognormal distribution. Thus, taking these analyses together, our metric of inter-areal coherence was robust to variations across individual subjects.

Robustness to the choice of coherence metric

As an alternative to the coherence metric, we computed the Pearson correlation coefficient for each pair of 10 s IFP signals and repeated the analyses. The interactome characterized by the Pearson correlation coefficients revealed significant interactions between 182 pairs of Desikan areas. Of these area pairs, 138 (76%) overlapped with the 193 area pairs that showed significant interactions using the coherence metric (Figure 5). The correlation between the interaction matrices calculated using the Pearson metric and the coherence metric was significant (r = 0.27, p = 0.001, n = 138). The correlation between the Pearson metric and the gamma band coherence metric was also significant (r = 0.25, p = 0.003, n = 138). The standard deviation of the Pearson correlation values was comparable to that of the coherence metric. We also found small-world network properties in the Pearson correlation metric (Table S6). The distribution of 182 Pearson correlation coefficients did not deviate from a log-normal distribution (p = 0.71, n = 182, Kolmogorov-Smirnov test), or from a normal distribution (p = 0.36, n = 182, Kolmogorov-Smirnov test). The numbers of significant bipolar pairs and numbers of significant subjects using the Pearson correlation metric was similar to those using the coherence metric. Therefore, the main findings were robust to the choice of the interaction metric.

Small-world graph analysis statistics

We compared the characteristic path length L and the clustering coefficient C of physiological interaction networks (Figure 6A) with two classes of null-hypotheses networks: (i) random network, assuming that the network arose from a random edge assignment (Figure 6B), and (ii) lattice network, assuming that the network arose from a regular lattice edge assignment (Figure 6C)(Watts and Strogatz, 1998). Random and regular lattice networks were generated as described in (Muldoon et al., 2016). Randomization was performed by preserving weights from the edges of the original network. Additional networks were randomly initialized with the same number of edges and non-edges, where the order was determined by random permutation (randperm MATLAB function). To compute network indices, 10,000 random networks were generated to estimate confidence intervals. The same procedure was repeated for the macaque monkey interaction network (Figure S8A), with corresponding constructions of the random network (Figure S8B) and the lattice network (Figure S8C).

Lattice networks were generated while preserving coherence weights from the original networks. Network edge weights were sorted in descending order to prioritize assigning higher weights closer to the lattice center. The innermost radius was filled first, containing n edges across n nodes. The initial network was the ring network where all nodes were connected to every other node via one circular path. The order of weights for filling the ring was determined by random permutation using the n highest weights. Following this construction, the next innermost radius was filled containing up to n edges with the n next highest edge weights. These edges connected every node to its second-neighboring node. The procedure was repeated until all edge weights had been depleted. The final assignment at level k required the random permutation of mod(m,k) weights, where m is the total number of edges in the network. Randomization was performed k times, with either n or mod(m,k) number of elements per batch. The lowest edge weight from each batch was greater than or equal to the highest edge weight from the next batch.

For each network graph, we created the random and lattice null network models and calculated the respective L, L, C, and C values. We sampled the random and lattice graphs with 10,000 initializations (Tables S4 and S6).

Robustness of small-world properties

Following our analysis of randomly re-calculated coherence values after removing subjects, we found small-world properties in the re-calculated networks. For random removal of n = 1 subjects, there were 20 nodes remaining (of the original 20 nodes, as shown in Figure 6) that formed an edge-complete network suitable for network analysis. The re-calculated small-world coefficient σ (Humphries and Gurney, 2008) and small-world measurement ω (Telesford et al., 2011) showed small-worldedness: σ = 1.20 (95% CI: 1.14 to 1.26), ω = −0.40 (95% CI: −0.42 to −0.39) for the broadband coherence network, and σ = 1.28 (95% CI: 1.19 to 1.36), ω = 0.44 (95% CI: −0.45 to −0.43) for the gamma coherence network. For random removal of n = 24 subjects, there were again 20 nodes remaining that formed an edge-complete network. The resulting network measures also showed small-worldedness: σ = 1.21 (95% CI: 1.14 to 1.27), ω = −0.41 (95% CI: −0.42 to −0.40) for the broadband coherence network, and σ = 1.28 (95% CI: 1.19 to 1.37), ω = −0.44 (95% CI: −0.45 to −0.42) for the gamma coherence network.

Classification of sleep versus wake states

We used a machine learning approach to assess whether it was possible to decode the annotation labels from the brain interaction data (Kriegeskorte and Kreiman, 2011). First, we used a nearest neighbor classifier to decode whether the subjectwas asleep or awake using only the coherence measurements. Each 10 s segment was assigned a binary label: sleep or awake. Transitions from sleep to wake and viceversa created ambiguous windows andwere removed. Since there weremoreperiodsof wakethansleepinevery subject,the wake labels were randomly subsampled to maintain an equal balance of wake and sleep labels. We used 20-fold cross-validation to separate the data into a training set and a test set, using 1000 iterations. We considered all pairwise coherence measurements for each subject. Using only the training data, we selected the top 10 most significant coherence features, based on the p value of a two-sided ranksum test between awake and sleep states. The results shown in Figure 7 represent cross-validated test accuracy values ranging from 0 to 1, where chance = 0.5. Statistical significance was determined by comparing the accuracies against those obtained after randomly shuffling the sleep/wake labels. The statistical analyses were corrected for multiple comparisons across subjects by Bonferroni’s method.

ADDITIONAL RESOURCES

Supplemental website

Additional materials and interactive rendering of the results can be obtained from the following supplemental website: http://braininteractome.com (Video S1)

The web materials include:

Maps of all physiological interactions in the average brain

Maps of all physiological interactions in each of 48 subjects

Maps of all physiological interactions in additional frequency bands

Interactive maps that allow the user to rotate the brains, highlight specific areas and electrodes, and visualize individual pairs of interacting areas and electrodes.

Note about the displays in the supplemental website: we use line objects to show interactions between two positions in the supplemental website. Bipolar electrodes positions were first defined on a three-dimensional coordinate system using the iELVis package (Groppe et al., 2017). These positions were represented by spheres plotted on reconstructed pial surfaces. Next, a series of line segments were constructed to connect pairs of spheres. The lines connecting different electrodes do not represent an anatomical basis for interactions, nor are they meant to represent axonal tracts. These representations were used only for visualization purposes. To construct these line segment paths, electrode positions were first localized on spherical surfaces using freesurfer tools (Dale et al., 1999). The path between these two electrode vertices was then determined by tracing the great-circle that passed through these two points, taking the shorter of the two possible paths. The path was mapped back onto the pial surface of the corresponding subject following a nearest-neighbor assignment of points to the great-circle path. Meshes were generated from each path using the tubeplot package (Wesenberg, 2016). Following these preparatory steps, graphical elements were combined into one common view in the supplemental website. The website was made using blender (The Blender Foundation, Amsterdam, Netherlands) and blend4web (Triumph LLC, Moscow, Russia).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Raw IFP data recorded from a monkey during rest	Neurotycho (Yanagawa et al., 2013)	http://www.neurotycho.org/anesthesia-and-sleep-task
Functional Interactome Results	This paper	http://braininteractome.com/
Software and algorithms
MATLAB R2020b	The MathWorks, Inc., Natick, MA	https://www.mathworks.com/; RRID: SCR_001622
Natus NeuroWorks EEG 8	Natus Medical Inc., Oakville, Ontario, Canada	https://neuro.natus.com/products-services/natus-neuroworks-eeg-software
FreeSurfer 6	Dale et al., 1999	https://surfer.nmr.mgh.harvard.edu/; RRID:SCR_001847
Intracranial Electrode Visualization (iELVis) Toolbox	Groppe et al., 2017	https://github.com/iELVis/iELVis; RRID:SCR_016109
Computerized Anatomical Reconstruction and Editing Toolkit (Caret) 5.616	Van Essen Lab (Van Essen et al., 2001)	http://brainvis.wustl.edu/wiki/index.php/Caret:Download; RRID:SCR_006260
Coherence algorithm	This paper	https://github.com/kreimanlab/MesoscopicInteractions