Experimental procedures
Pre-implantation neuroimaging. All participants underwent whole-brain high-resolution T1-weighted structural MRI scans before electrode implantation. In a subset of ten participants (Supplementary Table 2), RS-fMRI data were used for estimates of functional connectivity. The scanner was a 3T GE Discovery MR750W with a 32-channel head coil. The pre-electrode implantation anatomical T1 scan (3D FSPGR BRAVO sequence) was obtained with the following parameters: FOV = 25.6 cm, flip angle = 12 deg., TR = 8.50 ms, TE = 3.29 ms, inversion time = 450 ms, voxel size = 1.0 × 1.0 × 0.8 mm. For RS-fMRI, 5 blocks of 5-minute gradient-echo EPI runs (650 volumes) were collected with the following parameters: FOV = 22.0 cm, TR = 2260 ms, TE = 30 ms, flip angle = 80 deg., voxel size = 3.45 × 3.45 × 4.0 mm. In some cases, fewer RS acquisition sequences were used in the final analysis due to movement artifact or because the full scanning session was not completed. For each participant, RS-fMRI runs were acquired in the same session but non-contiguously (dispersed within an imaging session to avoid habituation). Participants were asked to keep their eyes open, and a fixation cross was presented through a projector.
iEEG recordings. iEEG recordings were obtained using either subdural and depth electrodes, or depth electrodes alone, based on clinical indications. Electrode arrays were manufactured by Ad-Tech Medical (Racine, WI). Subdural arrays, implanted in 36 participants out of 46, consisted of platinum-iridium discs (2.3 mm diameter, 5–10 mm inter-electrode distance), embedded in a silicon membrane. Stereotactically implanted depth arrays included between 4 and 12 cylindrical contacts along the electrode shaft, with 5–10 mm inter-electrode distance. A subgaleal electrode, placed over the cranial vertex near midline, was used as a reference in all participants. All electrodes were placed solely on the basis of clinical requirements, as determined by the team of epileptologists and neurosurgeons99.
No-task RS data were recorded in the dedicated, electrically shielded suite in The University of Iowa Clinical Research Unit while the participants lay in the hospital bed. RS data were collected 6.4 +/- 3.5 days (mean ± standard deviation; range 1.5–20.9) after electrode implantation surgery. In the first 15 participants (L275 through L362), data were recorded using a TDT RZ2 real-time processor (Tucker-Davis Technologies, Alachua, FL). In the remaining 31 participants (R369 through L585), data acquisition was performed using a Neuralynx Atlas System (Neuralynx Inc., Bozeman, MT). Recorded data were amplified, filtered (0.1–500 Hz bandpass, 5 dB/octave rolloff for TDT-recorded data; 0.7–800 Hz bandpass, 12 dB/octave rolloff for Neuralynx-recorded data) and digitized at a sampling rate of 2034.5 Hz (TDT) or 2000 Hz (Neuralynx). The durations of recordings were 13 +/- 11 min. In all but two participants, recording durations were between 10 and 22 min.; in one participant duration was 6 min., and in one participant the duration was 81 min.
Data analysis
Anatomical reconstruction and ROI parcellation. Localization of recording sites and their assignment to ROIs relied on post-implantation T1-weighted anatomical MRI and post-implantation computed tomography (CT). All images were initially aligned with pre-operative T1 scans using linear coregistration implemented in FSL (FLIRT)100. Electrodes were identified in the post-implantation MRI as magnetic susceptibility artifacts and in the CT as metallic hyperdensities. Electrode locations were further refined within the space of the pre-operative MRI using three-dimensional non-linear thin-plate spline warping101, which corrected for post-operative brain shift and distortion. The warping was constrained with 50–100 control points, manually selected throughout the brain, which were visually aligned to landmarks in the pre- and post-implantation MRI.
To pool data across participants, the dimensionality of connectivity matrices was reduced by assigning electrodes to one of 58 ROIs organized into 6 ROI groups (see Fig. 1; Supplementary Table 2, 3) based upon anatomical reconstructions of electrode locations in each participant. For subdural arrays, ROI assignment was informed by automated parcellation of cortical gyri102,103 as implemented in the FreeSurfer software package. For depth arrays, it was informed by MRI sections along sagittal, coronal, and axial planes. For recording sites in Heschl’s gyrus, delineation of the border between core auditory cortex adjacent non-core areas (HGPM and HGAL, respectively) was performed in each participant using physiological criteria104,105. Specifically, recording sites were assigned to HGPM if they exhibited phase-locked (frequency-following) responses to 100 Hz click trains and if the averaged evoked potentials to these stimuli featured short-latency (< 20 ms) peaks. Such response features are characteristic for HGPM and are not present within HGAL104. Additionally, correlation coefficients between average evoked potential waveforms recorded from adjacent sites were examined to identify discontinuities in response profiles along Heschl’s gyrus that could be interpreted as reflecting a transition from HGPM to HGAL. Superior temporal gyrus was subdivided into posterior and middle non-core auditory cortex ROIs (STGP and STGM), and auditory-related anterior ROI (STGA) using the transverse temporal sulcus and ascending ramus of the Sylvian fissure as macroanatomical boundaries. The insula was subdivided into posterior and anterior ROIs, with the former considered within the auditory-related ROI group5. Middle and inferior temporal gyrus were each divided into posterior, middle, and anterior ROIs by diving the gyrus into three approximately equal-length thirds. Angular gyrus was divided into posterior and anterior ROIs using the angular sulcus as a macroanatomical boundary. Anterior cingulate cortex was identified by automatic parcellation in FreeSurfer and was considered as part of the prefrontal ROI group, separately from the rest of the cingulate gyrus. Postcentral and precentral gyri were each divided into ventral and dorsal portions using the yMNI coordinate (see below) of 40 mm as a boundary. Recording sites identified as seizure foci or characterized by excessive noise, and depth electrode contacts localized to the white matter or outside brain, were excluded from analyses and are not listed in Supplementary Table 2.
Preprocessing of fMRI data. Standard preprocessing was applied to the RS-fMRI data acquired in the pre-implantation scan using FSL’s FEAT pipeline, including spatial alignment and nuisance regression. White matter, cerebrospinal fluid and global ROIs were created using deep white matter, lateral ventricles and a whole brain mask, respectively. Regression was performed using the time series of these three nuisance ROIs as well as 6 motion parameters (3 rotations and 3 translations) and their derivatives, detrended with second order polynomials. Temporal bandpass filtering was 0.008–0.08 Hz. Spatial smoothing was applied with a Gaussian kernel (6 mm full-width at half maximum). The first two images from each run were discarded. Frame censoring was applied when the Euclidean norm of derivatives of motion parameters exceeded 0.5 mm106. All runs were processed in native EPI space, then the residual data were transformed to MNI152 and concatenated.
Preprocessing of iEEG data. Analysis of iEEG data was performed using custom software written in MATLAB Version 2020a programming environment (MathWorks, Natick, MA, USA). After initial rejection of recording sites identified as seizure foci, several automated steps were taken to exclude recording channels and time intervals contaminated by noise. First, channels were excluded if average power in any frequency band (broadband, delta, theta, alpha, beta, gamma, or high gamma; see below) exceeded 3.5 standard deviations of the average power across all channels for that participant. Next, transient artifacts were detected by identifying voltage deflections exceeding 10 standard deviations on a given channel. A time window was identified extending before and after the detected artifact until the voltage returned to the zero-mean baseline plus an additional 100 ms buffer before and after. High-frequency artifacts were also removed by masking segments of data with high gamma power exceeding 5 standard deviations of the mean across all segments. Only time bins free of these artifact masks were considered in subsequent analyses. Artifact rejection was applied across all channels simultaneously so that all connectivity measures were derived from the same time windows. Occasionally, particular channels survived the initial average power criteria yet had frequent artifacts that led to loss of data across all the other channels. There is a tradeoff in rejecting artifacts (losing time across all channels) and rejecting channels (losing all data for that channel). If artifacts occur on many channels, there is little benefit to excluding any one channel. However, if frequent artifacts occur on one or simultaneously on up to a few channels, omitting these can save more data from other channels than those channels contribute at all other times. We chose to optimize the total data retained, channels × time windows, and omitted some channels when necessary. To remove shared signals unlikely to derive from brain activity, data from retained channels were high-pass filtered above 200 Hz, and a spatial filter was derived from the singular value decomposition omitting the first singular vector. This spatial filter was then applied to the broadband signal to remove this common signal.
Connectivity analysis. For RS-fMRI data, BOLD signals were averaged across voxel groupings and functional connectivity was calculated as Pearson correlation coefficients. Voxel groupings were either based on the Schaefer-Yeo 400 parcellation scheme53 in MNI-152 space, or were based on iEEG electrode location in participant space (see Fig. 1). For the latter, fMRI voxels were chosen to represent comparable regions of the brain recorded by iEEG electrodes. For each electrode, the anatomical coordinates of the recording site were mapped to the closest valid MRI voxel, E, and a sphere of 25 voxels (25 mm3) centered on E used as the corresponding recording site. This process was repeated for all N electrodes in the same ROI, and a single time series computed as the average of the fMRI BOLD signal in these N×25 voxels. These averages were used to compute an ROI-by-ROI connectivity matrix for RS-fMRI data. For comparisons between iEEG and fMRI embeddings, voxels were processed in participant space and ROI labels from the parcellation scheme illustrated in Fig. 1 and Supplementary Table 2 were applied to the fMRI data. For comparisons between fMRI embeddings derived from all cortical ROIs versus fMRI embeddings derived from just ROIs sampled in the iEEG experiments, electrode locations were transformed from participant space to MNI-152 space, then assigned to ROIs within the Schaefer-Yeo 400 scheme.
For iEEG data, two measures of functional (non-directed) connectivity were used: the orthogonalized band power envelope correlation29 and the debiased weighted phase lag index (wPLI;52), a measure of phase synchronization. Both measures avoid artifacts due to volume conduction by discounting connectivity near zero phase lag. For both measures, data were divided into 60-second segments, pairwise connectivity estimated in each segment, and then connectivity estimates averaged across all segments for that participant.
Envelope correlations were estimated for each data segment and every recording site as in29, except time-frequency decomposition was performed using the demodulated band transform107, rather than wavelets. For each frequency band (theta: 4–8 Hz, alpha: 8–13 Hz, beta: 13–30 Hz, gamma: 30–70 Hz; high gamma: 70–120 Hz), the power at each time bin was calculated as the average (across frequencies) log of the squared amplitude. For each pair of signals X and Y, one was orthogonalized to the other by taking the magnitude of the imaginary component of the product of one signal with the normalized complex conjugate of the other:
$${Y}_{orth}= \left|\text{I}\text{m}\left\{Y\times {X}^{*}/\left|X\right|\right\}\right|$$
Both signals were band-pass filtered (0.2–1 Hz), and the Pearson correlation calculated between signals. The process was repeated by orthogonalizing in the other direction and the overall envelope correlation for a pair of recording sites was the average of the two Pearson correlations.
wPLI was estimated for each data segment and every recording site pair from the sign of the imaginary part of the cross-spectrum at each frequency and averaged across frequencies within each band of interest (theta: 4–8 Hz, alpha: 8–13 Hz, beta: 13–30 Hz). The cross spectrum was calculated from the demodulated band transform as described previously108.
Prior to diffusion map embedding, connectivity matrices were thresholded by saving at least the top third (rounded up) connections for every row, as well as their corresponding columns (to preserve symmetry). We also included any connections making up the minimum spanning tree of the graph represented by the elementwise reciprocal of the connectivity matrix to ensure the graph is connected.
ROI-based connectivity analysis. Connectivity between ROIs was computed as the average envelope correlation or wPLI value between all pairs of recording sites in the two ROIs. For analyses in which connectivity was summarized across participants (Fig. 3–8), we used only a subset of ROIs such that every possible pair of included ROIs was represented in at least two participants (Supplementary Table 2). This list of ROIs was obtained by iteratively removing ROIs with the worst cross-coverage with other ROIs until every ROI remaining had sufficient coverage with all remaining ROIs.
Diffusion map embedding. See the Appendix for details about DME. In brief, the connectivity matrix K = [k(i,j)] is normalized by degree to yield a transition probability matrix P. DME consists of performing an eigendecomposition of P, with eigenvalues scaled according to the free parameter t, which sets the spatial scale of the analysis. If the recording sites are conceptualized as nodes on a graph with edges defined by K, then P can be understood as the transition probability matrix for a ‘random walk’ or a ‘diffusion’ on the graph (see Appendix;30,31). The parameter t is the number of time steps in that random walk; larger values of t correspond to exploring the structure of the data at larger spatial scales. In the analyses presented here, k(i,j) is either the orthogonalized power envelope correlations, or the debiased wPLI. We note that in recent applications of DME to fMRI functional connectivity data, K was first transformed by applying cosine similarity32. However, this additional step has not been universally implemented (e.g.,109), nor is it required within the mathematical framework of DME30,31. In the case of the dataset presented here, applying cosine similarity to the data served to smear the data in embedding space, increasing variance, and decreasing separation between ROIs, but otherwise produced qualitatively similar results.
In two sets of analyses presented here, pairs of embeddings were compared to each other: in the analysis of lateralization of speech and language networks, and in the comparison between iEEG and fMRI data. To do that, we used a change of basis operator to map embeddings into a common embedding space using the method described in Coifman et al 201431.
Dimensionality reduction via low rank approximations to P. Diffusion map embedding offers an opportunity to reduce the dimensionality of the underlying data by considering only those dimensions that contribute importantly to the structure of the data, as manifested in the structure of the transition probability matrix P. We used the eigenvalue spectrum of P to determine its ideal low rank approximation, balancing dimensionality reduction and information loss. Because P is real and symmetric, the magnitude of the eigenvalues is identical to the singular values of P. The singular values tell us about the fidelity of low rank approximations to P. Specifically, if P has a set of singular values σ1≥ σ1≥…≥ σn, then for any integer k ≥ 1,
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAALEAAAAmCAYAAACLSno+AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAAXfSURBVHhe7ZzLLyxPFMfL/QPIYGXpkRCEeEWEjYXXWjJCYknGwoYgRESEhIWFhVcsbDxjYUOwYMFCvMLeiMTW2x9wf/09+szt27dndOuenu6f+iSVqu6u7unTfeacU6dqJum3gpBIfMwvtZZIfItUYonvkUos8T1SiSW+RyqxxPdIJbbJ3d2d6OrqEqmpqSIpKYlKQ0ODGBoaEru7u3SceXp6Euvr69Rf4hyeVmIoxMnJCbWhGFNTU9T2Cjc3N6SQnZ2d4vn5WSBbGQ6HRXd3Nx2fmZkRWVlZEeVua2uj/bOzs1QzXpfT60hLbIOFhQWxsrIiioqK1D1CZGZmiqamJjE+Pi729vZIsblgu6WlRe0pcQqpxDaARU1LS1O3JIlCKrHE97iixIjxEPc5CV9TX7KzsylOxSAq3iB+NbqH8vJyGde6iG8tcV9fnwgGg9Sur6+nmPPx8VEEAgExNzcnhoeH6Vg8QTyMgRuYnJyke7i+vqbBXX9/P2UiJPHHFSWGwuEFO01JSYna+gTx6ejoKLWhyPEGnwfLrwWDvFAoRO3l5WWqJfHFkhLDssBdwo0iB4oXiG3kROG+4caxjf1IPwHUnENluB9cLgqOozhpudhCJhJkKtwmWoijLehjBbxbvCeESfprafPgicKSEiM9BNd9dnYmHh4exO3trZifnxcTExPkvnt7e8mlI2e6ublJ58AyDQwMUJvhPOni4qKoqqqi/nV1dbYnAfBAR0ZGqN3R0UG12yDfy16gtbWVajfRp/WMCvqYBUYoJydHXF1diaWlJQqXysrKItdKxBdVz7fCiYqKCkrwg/z8fKrxwiAQXCyOX15e0v5YQNGqq6up3d7eLl5eXqhtlf39fbIKsL64BuJThDBGQMn01sSo8OSDWRAD47yamhq6j52dnYhsfgZywXjBS8IgoVxcXEQ8rRfwzMAuOTlZbVmHB3Yo8A7RFBhAsbhvrGJVAXlgh3J+fk4THoBDMEyMuIGT4QQUFQYCHpZh5f34+KDajHwwCNE+k6ft7WRzPKPE/1c4BKusrFT3xBcnwwlWVG3IwGFiXl4e1Xbkwxdia2uLwkk7eEaJ39/fqbaS30WcBhCju5EX1oPPhOUHh4eHVBuBPnDDsDawWnasjptkZGRQzdYX1hbjH3gd7Uzld+XDOfCa+iwT4GuZQvlmmmZtbQ15MirKgI72KfEfbQcCAdrGfu6D/spAgI5hOxgMUp9QKBTpg+PaPrgeg+3j42NqK9/238rDozZAm6+BguNug8/U3oP2/hjIBrnD4TA9G/RRBr/q0U9wbjQ5Ew3eofbd8HtnYsmnfz5cWFYG5+hlxjbr1FdYUuKvwM1rFVQZxZJgEJCPa7e/Qiuw116uWXDP/DyUwZ6692/8LKcZ+SAb5IoGrqGXeXBwkIoZHA0npqenaRCgXJdSbT09PeLo6IhSMjygQOzkhbSMWyDMUKyZUF4ILc2Ea7aSW+WBD54fcunIz3sJu/IBhIVvb2/q1qfMyDJhJaApPnXZm+D2/G6JIQM8ELwP3CM8lZ5YcqLNLhyWKZZFSwRm5ItmibEf52vLd5DZiTijPGMaBMH7YBSuXxD/FRj4cE6+sLBQ1NbWUtsrmJEP6UqjjIhRuvM7OKLERhMICB2suhVJdJAJ2d7ejpkD/7Eo2u8IcHV8ObgWtP3o/r0ID5hRS/7FsXAiJSVFbf1JkvM+ttRmZ4okf4AFxrqUsbExcttYbCX5G8djYihrcXExrfVtbm6mfYh9lMDec/GcH4ACY0FReno6PVsza1J+Go4rsWLdqWBOXTurg1m1goKCyEyM1QU2PxUMlPiZolhZgfZTcEyJOc9nNP2L3CHyfrm5ueL+/p6sMs+9SyR2cSw7gTl1wP+toOX09JRqLFFEGgbWRGulJRI7uPL/xFjp9Pr6KkpLS8XGxgb9bAeLS37SzJ0kfjgeExtxcHBAg7rGxkb6EeXq6qpUYIljyH+Kl/geVyyxRBJPpBJLfI9UYonPEeI/dFOobPSqMfwAAAAASUVORK5CYII=)
where \(\tilde{{\mathbf{P}}_{k}}\) is the rank-k approximation to P. Thus, the magnitude of the eigenvalues corresponds to the fidelity of the lower dimensional approximation, and the difference in the magnitude of successive eigenvalues represents the improvement in that approximation as the dimensionality increases. The spectrum of P invariably has an inflection point (“elbow”), separating two sets of eigenvalues λi: those whose magnitude decreases more quickly with increasing i, and those beyond the inflection point whose magnitude decreases more slowly with increasing i. The inflection point thus delineates the number of dimensions that are most important for approximating P. The inflection point kinfl was identified algorithmically110, and the number of dimensions retained set equal to kinfl – 1.
Comparing distances in embedding space. The relative distance between points in embedding space provides insight into the underlying functional geometry. In several analyses presented here, two embeddings of identical sets of ROIs were compared as ROI distances within the two embeddings. After mapping to a common space and reducing dimensionality as described above, the two embeddings A and B were used to create the pairwise distance matrices A` and B`. The Pearson correlation coefficient r was then computed between the upper triangles (excluding the diagonal) of the corresponding elements in the distance matrices.
Signal to noise (SNR) characteristics. To measure the robustness of the embedding analysis to variability over time, an SNR was computed as follows. For each participant, a channel × channel P matrix was calculated for each 60 s segment of data. For each segment, DME analysis was applied and a channel × channel distance matrix calculated. These distance matrices were averaged across segments. The ‘signal’ of interest was defined as the variability (standard deviation) of this averaged distance matrix (ignoring the diagonals). The ‘noise’ was defined as the variability across time, estimated for each element of the distance matrix as the standard deviation across segments, then averaged across the elements of the matrix. The SNR for functional connectivity itself was computed in an analogous manner, using the original channel × channel connectivity matrix rather than the matrix of embedding distances.
Estimating precision in position and distances in embedding space. To obtain error estimates for both ROI locations in embedding space and embedding distance between ROIs, average ROI × ROI adjacency matrices were calculated. This was done by drawing each edge from an averaged bootstrap sample across participants, obtaining 10,000 such adjacency matrices, and performing diffusion map embedding for each. For locations in embedding space, these embeddings were then mapped via the change of basis procedure described above to the original group average embedding space. For each ROI, the mapped bootstrap iterations produced a cloud of locations in embedding space which were summarized by the standard deviation in each dimension. For embedding distances, no change of basis was necessary because distances were preserved across bases.
To compare the positions of STSL versus STSU relative to canonical auditory cortical ROIs (HGPM, HGAL, PT, PP, STGP, and STGM) or ROIs involved in semantic processing (STGA, MTGA, MTGP, ITGA, ITGP, TP, AGA, AGP, SMG, IFGop, IFGtr, IFGor28,38−40), we calculated the average pairwise distance from STSL or STSU to each such ROI. The difference between these averages was compared to a null distribution obtained by Monte Carlo sampling of the equivalent statistic obtained by randomly exchanging STSL/STSU labels by participant. The specific comparisons performed were chosen a priori to constrain the number of possible hypotheses to test; pairwise comparisons of all possible ROI pairs (let alone comparisons of all higher-order groupings) would not have had sufficient statistical power under appropriate corrections for multiple comparisons. Though different choices could have been made for inclusion in the “semantic processing” category, exchanging one or two of these ROIs would not strongly influence the average distance in a group of twelve ROIs.
Hierarchical clustering. Agglomerative hierarchical clustering was done using the linkage function in MATLAB, with Euclidean distance as the distance metric and Ward’s linkage (minimum variance algorithm) as the linkage method. The ordering of ROIs along the horizontal axis in the dendrogram was determined using the optimalleaforder function in MATLAB, with the optimization criterion set to ‘group’.
Comparing language dominant/non-dominant hemispheres. To test for differences in embedding location between language dominant and non-dominant hemispheres, two measures were considered: different location of individual ROIs in embedding space, and different pairwise distances between ROIs in embedding space. To calculate differences in location of individual ROIs, dominant/non-dominant average embeddings were mapped to a common space (from an embedding using the average across all participants regardless of language dominance) using the change of basis operator. The language-dominant location difference for a specific ROI was calculated as the Euclidean distance between the two locations of each ROI in this common space. For pairwise distances, the change of basis is irrelevant, so pairwise Euclidean distances were calculated in embedding space for each hemisphere and then subtracted to obtain a difference matrix. To determine whether the differences in location or pairwise distances were larger than expected by chance, random permutations of the dominant/non-dominant labels were used to generate empirical null distributions. Since this approach produces a p-value for every pair of connections, p-values were adjusted using false discovery rate (FDR) to account for multiple comparisons.
Analyses of fMRI connectivity in embedding space. Two sets of analyses were performed using fMRI data. First, iEEG and fMRI data were compared in embedding space. In this analysis, connectivity based on RS-fMRI data from voxels located at electrode recording sites was compare with the corresponding connectivity matrix derived from iEEG data. The embedding analysis was applied to the two connectivity matrices, all pairwise inter-ROI distances computed, and iEEG and fMRI data compared using the correlation of the pairwise ROI distances. The second analysis was to compare embeddings derived from all ROIs in the RS-fMRI scans to those derived from just ROIs sampled with iEEG electrodes. Here, ROI × ROI connectivity matrices were computed for all ROIs, then embeddings created from the full matrices or from matrices containing just rows and columns corresponding to the ROIs sampled with iEEG.