Human cortical dynamics reflect graded contributions of local geometry and network topography

Domain: neuroscience
Published: 2025-12-26
Editorial outcome: major_revision

Abstract
The brain is a physically embedded and heavily interconnected system that expresses neural rhythms across multiple timescales. While these dynamics result from the complex interplay of local and inter-regional factors, the relative contribution of such mechanisms across the cortex remains unclear. Our study explores geometric, microstructural, and connectome-level constraints on cortex-wide neural activity. We leverage intracranial electroencephalography recordings to derive a coordinate system of human cortical dynamics. Using multimodal neuroimaging, we could then demonstrate that these patterns are largely explainable by geometric properties indexed by inter-regional distance. However, dynamics in transmodal association regions are additionally explainable by incorporation of inter-regional microstructural similarity and connectivity information. Our findings are generally consistent when cross-referencing electroencephalography and imaging data from large-scale atlases and when using data obtained in the same individuals, suggesting subject-specificity and population-level generalizability. Together, our results suggest that the relative contribution of local and macroscale constraints on cortical dynamics varies systematically across the cortical sheet, specifically highlighting the role of transmodal networks in inter-regional cortical coordination.

Introduction
The brain generates nested and propagating rhythms from the synchronized firing of neuronal ensembles1. These dynamics signals are a core feature of brain organization, serving as established markers of brain health2–4, neural communication5–7, vigilance, and cognitive states8–15. Intracranial electroencephalography (iEEG) offers a unique window into these rhythms. Moving beyond spatially restricted recordings performed in individual patients16,17 recent data sharing initiatives have provided comprehensive maps of human neurophysiology and unlocked new avenues to study inter-regional heterogeneity in neural dynamics8,18. However, technical challenges relating to multimodal data availability, registration, and integration across and within patient cohorts have particularly hindered these investigations, which have so far been limited to single scales of inquiry. The present work fills this gap by uncovering large-scale gradients of neural dynamics and relating this complex activity landscape to underlying multiscale anatomy.

Prior work has demonstrated distinct neurophysiological signatures operating across distributed cortical territories19. These signatures generally follow overlapping spatial gradients of cortical microarchitecture involving neurotransmitter receptors and transporters, gene expression, metabolic demands, cortical thickness, and proxies of intracortical myelination20–22. Relatedly, focal lesions altering local anatomy and microarchitecture co-localize with atypical neural dynamics2,3. While this evidence supports the role of intrinsic physiological mechanisms in generating regional neural activity23,24, local structural features can only partially explain spatial variations in neural dynamics. For instance, areas made up of fundamentally different milieus can also produce rhythms with similar features: peaks in beta-range activity have been recorded in large parts of the frontal lobe, from the precentral gyrus to the frontal pole, as well as portions of the parietal lobe and the anterior insula18. As such, network-level propagation of brain activity via direct structural connections as well as polysynaptic functional interactions has also been suggested to account for this inter-regional variability25. Indeed, local dynamics systematically vary along the cortical processing hierarchy, as demonstrated by task paradigms, stimulation work, and in silico explorations26–28. Resting dynamics also propagate to nearby regions via traveling waves, scaffolded by the brain’s short- and long-range connections29, producing a complex mosaic of intrinsic and extrinsic constraints on neural dynamics.

Such multiscale inquiries can be achieved by integrating the high spatiotemporal precision of iEEG with recent methodological and analytical advances in magnetic resonance imaging (MRI) offering complementary insights into cortico-cortical network organization. At the most basic level is cortical geometry, which reflects the proximity of brain regions. The brain is a physically embedded network, with adjacent regions generally participating in similar functional processes30. Moreover, adjacent regions are generally strongly connected via intracortical axon collaterals and superficial white matter fibers31, which closely regulate local neural activity32,33. A further important component of brain network architecture is local intracortical microstructure, which has classically been studied using cytoarchitectonic mapping of post-mortem tissue34. Recently, similar approaches have been translated to in vivo MRI using myelin-sensitive contrasts, such as T1 relaxometry35–38. Both lines of investigation suggest that regions with similar microstructure generally participate in corresponding functional processes and are more likely to share a direct structural connection34,35,39,40. Beyond the study of local geometry and microstructure, techniques such as diffusion MRI tractography and functional MRI connectivity have characterized large-scale, regionally distributed structural and functional networks. While the former approximates long-range, white matter fiber tracts in vivo41–44, the latter taps into spatially distributed, and often polysynaptic intrinsic networks45–47. In line with recently developed composite wiring diagrams of the human cortex emphasizing a synergy of different MRI modalities in predicting macroscale function48–54, simultaneous consideration of these diverse modalities elucidates the contribution of short- and long-range network information to human cortical dynamics.

Our goal was (i) to chart the cortical topography of neural rhythms in the human brain interrogated with iEEG18, and (ii) identify the differential association between cortical dynamics and network organization mapped with multimodal MRI55. Applying dimensionality reduction to a cortex-wide assessment of inter-regional iEEG power spectra similarity allowed us to situate cortical areas within a compact embedding space constrained by their dynamics. These gradients, mapped here in iEEG data, captured the coordinate system of neural rhythms in the neocortex, grouping regions with similar neurophysiology independently of their anatomic allocation. We then quantified the contribution of different network level features to the topography of neural dynamics using multilinear spatial statistics at a cortex-wide scale. We complemented this approach with region-specific models capturing the relative contribution of each network metric across cortical areas, and tracked how these contributions varied along the cortical processing hierarchy56. Our assessment was based on independent datasets and systematic variations of analysis parameters to identify robust associations between human cortical dynamics and multiscale network features.

Results
Mapping cortical gradients of neural dynamics
For the main analyses, we leveraged a Discovery cohort (MNI open iEEG atlas, n=106) in which iEEG data were aggregated in stereotactic space18. We assessed correspondence to neural organization derived from two high-definition MRI datasets acquired in 50 healthy individuals at 3T and 10 healthy individuals at 7T. Our main analyses thus rest on a design comparing iEEG and imaging modalities across two different datasets. We then evaluate the generalizability of our findings using a fully independent Replication cohort (MNI clinical dataset, n=18), where both iEEG and multimodal pre-surgical MRI data were available in each participant for within-cohort comparison of modalities.

The MNI open iEEG atlas includes recordings from 106 patients with intractable epilepsy, totaling over 1700 data points taken from non-epileptic channels registered to the MNI152 reference space18. Mapping all channels to a single hemisphere as in previous work8,18 granted wide sampling of the cortical sheet (Fig. 1A). Spatial heterogeneity of spectral features was visible across channels, but could also be recovered when averaging channel power spectra within larger cortical regions (Fig. 1B)57. We captured spatial variations in macroscale neural dynamics by applying dimensionality reduction to similarity patterns in regional power spectra (Fig. 1C). After averaging vertex-wise, log-transformed power spectral densities (PSDs) within parcels, regional PSDs were cross-correlated while accounting for the cortex-wide average power spectrum using partial correlations. We converted the result to a normalized angle affinity matrix and applied diffusion map embedding, a nonlinear dimensionality reduction technique, to recover eigenvectors (gradients) of macroscale similarity in regional neural dynamics. We retained gradients explaining a cumulative sum of 50% of the variance for further analyses (G1–G14; see Fig. S1A for the topography of these gradients). The embedding space formed by the first two gradients (explaining 17.37% of the variance in PSD similarity) was used as a complementary visualization of findings.

The first gradient (G1) evolved along the anterior-posterior direction, mainly differentiating primary motor and surrounding lateral prefrontal cortices from unimodal sensory regions (Fig. 1D). G1 indeed mainly followed the structural stratification defined by Von Economo and Koskinas58–60. This suggests a close correspondence between the axis of neurophysiological differentiation captured by G1 and regional heterogeneity in laminar differentiation, notably segregating “input” sensory regions with a granular structure from “output” motor areas with an agranular structure (Fig. S1B). This gradient mirrored regional variations in spectral peaks by segregating areas with strong beta and gamma-range activity from areas showing prominent delta, theta, and alpha-range peaks18. Next, G2 differentiated unimodal sensorimotor from transmodal cortices, recapitulating the previously established principal axis of cortical connectivity and functional hierarchy (Fig. 1D)61. Indeed, G2 was specifically correlated with the unimodal-transmodal gradient of cortical organization derived from functional MRI (Spearman r=0.564; p < 0.001; Fig. S1C; correlations with all other gradients did not survive Bonferroni correction). This pattern segregated regions showing predominant lower frequency activity from other cortical areas with peaks occurring above the delta frequency range. The two-dimensional representation of neural dynamics formed by these gradients separated unimodal sensory regions, precentral motor areas, and spatially distributed transmodal systems (Fig. 1E). Cortical surface plots shown in Fig. 1B furthermore highlighted consistency of regional power spectra amplitudes with previously published normative maps of neocortical dynamics2,3,62. We observed higher delta band power in frontal, parietal, and temporal cortices relative to unimodal sensory and motor regions. High alpha- and beta-range power was observed within occipital and precentral regions, respectively, while a rostro-caudal shift was observed in gamma-range activity across the cortex. Collectively, these results map gradients of neural rhythms able to resolve fine-grained oscillatory activity of individual channels. We additionally show that these gradient recover macroscale trends in the overall similarity of spectral content across the cortex.

Control analyses were implemented to assess the robustness of described gradients and embedding space. First, results were robust to the use of different parcellation schemes (Fig. S2). Indeed, generating gradients from unparcellated PSD data (i.e., channel regions—see “Materials and methods”) or averaging PSDs within nodes defined by two different atlases (i.e., Glasser-360 and a subdivision of the aparc parcellation) yielded similar topographies for both G1 (vertexwise Spearman correlations for unparcellated: r=0.673; Glasser-360: r=0.839; aparc-200 r=0.864) and G2 (vertexwise Spearman correlations for unparcellated: r=0.524; Glasser-360: r=0.580; aparc-200 r=0.645). Furthermore, we could recover the topography of the first 10 gradients when imposing an upper bound to the number of channels included in each parcel of the Schaefer-200 parcellation (Fig. S3), suggesting that our findings could not be attributed to variable channel coverage across nodes. We also replicated our findings when preserving the bilateral distribution of channels and applying a lower resolution parcellation (Schaefer-100) to maintain broad cortical coverage (Fig. S4). The high correspondence of G1 (vertexwise Spearman correlation r_left=0.767; r_right=0.804) and G2 (vertexwise Spearman correlation r_left=0.645; r_right=0.671) topographies with gradients derived from our main analysis provides additional support for the independence of our findings to parcellation choice. Lastly, spatial correlations assessing the consistency of gradients across a range of hyperparameter combinations were generally very strong (Fig. S7).

Multiscale connectomics of neural dynamics
After establishing the topography of macroscale gradients of neural dynamics, we investigated how this embedding space co-varied with multiscale properties of brain network organization. All network measures were gathered from an independent cohort of 50 unrelated, healthy volunteers that underwent multimodal 3T MRI (23 women; age mean ± standard deviation = 30.28 ± 5.74 years)55. We first computed the Euclidean distance between all node pairs in the embedding space, which served as a proxy for inter-node similarity in neural dynamics (embedding distance, Fig. 2A). This matrix was highly robust to the number of gradients included when calculating inter-node distances (Fig. S8). Node proximity in the embedding space was cross-referenced with measures of (i) cortical geometry indexed by inter-node geodesic distance computed along the surface mesh (geodesic distance, GD), (ii) tractography-based structural connectivity (SC), (iii) intrinsic functional connectivity generated from resting-state functional MRI (rs-fMRI), and (iv) microstructural profile covariance (MPC) derived from T1 relaxometry. Matrices were computed for each participant and averaged across individuals (for details on data processing, analysis, and group-level aggregation, see “Materials and methods”). Multimodal associations were computed by correlating the upper triangular edges of each matrix. We only considered MRI data from the right hemisphere, as all iEEG data had been transposed to the right hemisphere to improve cortical coverage. For each modality, statistical significance was determined by comparing empirically observed correlations to those obtained with surrogate matrices in which each row had been shuffled using spatial autocorrelation-preserving null models63,64.

We first observed a moderate correlation (r=-0.376; p < 0.001) between GD of cortical regions and embedding distance, indicating that nodes with similar dynamics were situated closer to each other along the cortical sheet (Fig. 2B). Node-wise correlations between both modalities were particularly strong in anterior and posterior lateral regions. We also obtained a moderate correlation (r=-0.341; p < 0.001) between edge strength of the SC matrix and embedding distance, suggesting closer neural dynamic similarity between regions sharing stronger structural connections (Fig. 2C). Regional SC and embedding distances were most strongly correlated in lateral frontal and temporal regions. Polysynaptic network architectures captured by the functional connectome were weakly, yet significantly, correlated with embedding distances (r=-0.227; p < 0.001), once again indicating stronger connections between regions sharing similar neural dynamics (Fig. 2D). This relationship was particularly strong in medial and lateral prefrontal, as well as medial occipito-temporal regions. Lastly, embedding distance was significantly yet more weakly correlated with MPC (r=-0.170; p < 0.001), with significant correlations confined to heteromodal and paralimbic cortical structures such as mesial temporal and frontal regions, as well as the precuneus and lateral parietal areas (Fig. 2E). Of note, relationships between neocortical dynamics similarity and multiscale connectivity measures were consistent when binarizing network matrices across a wide range of thresholds (Fig. 2B–E). Findings were replicable when cross-referencing embedding distances with multimodal data collected in a distinct set of participants scanned at 7T MRI65 (Fig. S5), and when computing correlations between embedding distances and MRI measures from bilaterally mapped data (SC r=-0.322; FC r=-0.217; MPC r=-0.185). Together, these results demonstrate that neural dynamics similarity is most strongly associated with cortical geometry indexed by geodesic distance, with weaker yet statistically significant relationships to large-scale network metrics. Modality-specific correlations were nonetheless highly variable across regions.

These results show regional variations in the relative contributions of different network properties to cortical dynamics, while emphasizing the stronger relationship with spatial proximity between nodes. We next sought to clarify the combined contribution of network features to neural dynamic similarity by quantifying potential gains in model accuracy when considering multiscale network properties over GD alone. To do so, we constructed several multiple linear regression models, in which we varied the combination of predictors implemented in the model as well as the nature of the predicted and predictor variables. A baseline model including GD as its only predictor was used to benchmark the performance of the more complex models incorporating multimodal predictor sets. These models encompassed partial models with two (GD and another modality) or three predictors (GD and two other modalities) in addition to a full model with four predictors (all modalities; Fig. 3A).

Constructing this model across all unique edges (i.e., the upper triangle of each matrix) highlighted the strong contribution of GD relative to other modalities across all models (Fig. 3B). Notably, more complex models provided minimal gains in prediction accuracy, either when implementing the model on all data points (range of adjusted r^2 = 0.147–0.171) or when bootstrapping the data (90% train and 10% test; 1000 permutations; range of median adjusted r^2 = 0.142–0.161). We assessed cross-model accuracy variations by correlating the prediction error for each edge with the number of predictors in the corresponding model. In this approach, negative correlations indicate decreasing error with more complex models (i.e., high prediction accuracy for models with more predictors), while positive correlations indicate increasing error with more complex models (i.e., including additional predictors resulted in a noisier model). Our results show overall accuracy gains with increasing model complexity in mesial and anterolateral temporal, orbitofrontal, and dorsolateral prefrontal cortices, as well as the precuneus, while opposite trends were observed primarily in occipito-parietal regions (Fig. 3B, right). These findings indicate heterogeneous contributions of network predictors across cortical regions.

We further explored region-specific associations between multiscale network features to neural dynamic similarity by fitting distinct multiple linear regression models for each node’s connections and comparing model accuracy across different predictor combinations (Fig. 3C). For instance, when comparing baseline (GD only) and full (all network features) models, the greatest gains in adjusted r^2 were seen in dorsal lateral and medial frontal cortices, orbitofrontal, mesiotemporal, lateral parietal and precuneus regions. Regional differences in accuracy gains across certain model pairs were significantly correlated with the cortical processing hierarchy represented by the principal gradient of intrinsic functional connectivity56,66. More specifically, regional gains in model accuracy were significantly correlated with this gradient in models including both GD and MPC as predictors over GD alone, and accuracy could be further increased by including either SC or FC. These findings were replicated using an alternate representation of this hierarchy mapping four functional zones co-localizing with regional changes in cortical laminar elaboration61,67 (Fig. S6). These findings suggest a complementarity of local geometry and microstructure in supporting macroscale similarity of neural dynamics.

Replicability of gradients of neural dynamics
Our work mapped the landscape of neural dynamics similarity of the cortex and decoded its multiscale network correlates using multilinear spatial statistics. We lastly leveraged an independent patient cohort to assess the replicability of macroscale gradients of neural dynamics and their association with within-sample neuroimaging metrics. The MNI clinical dataset (Replication cohort) includes 18 patients who underwent high-resolution multimodal MRI prior to clinical iEEG investigations, allowing for precise, patient-specific segmentation of cortical gray matter, co-registration of iEEG contacts, and cross-referencing to within-sample connectivity measurements. Although all data modalities were available in patients of the Replication cohort, these replication analyses study a completely distinct sample of participants relative to the Discovery cohort, for both iEEG and MRI metrics.

Following the same analytical pipeline as in our main analysis, electrophysiological data from all patients was combined to generate group-level gradients of neural dynamic similarity (Fig. 4A). The gradient space in this independent sample segregated temporo-parietal regions and frontal cortices along G1, and distinguished transmodal temporal and frontal structures from lower-order sensorimotor areas along G2 (Fig. 4B). The first 13 gradients, explaining a cumulative 50% of the variance in the input data, were all significantly correlated (p < 0.05; p < 0.003 with Bonferroni correction) with their corresponding gradient in the Discovery sample (after Procrustes alignment including translation, rotation, and scaling: mean ± SD Spearman r = 0.509 ± 0.076). After alignment, we found a significant correlation between the embedding distance matrices of Discovery and Replication cohorts (r = 0.304; p < 0.001), with variable similarity across regions (Fig. S9). We also applied a more conservative alignment strategy where Replication cohort gradients were reordered to optimally match (i.e., minimize the distance) with the Discovery cohort eigenvectors. With this approach, correlations were overall lower but also indicated generalizability of main findings in the Replication cohort (mean ± SD Spearman r = 0.196 ± 0.154; Fig. S10). Systematic within-cohort and cross-cohort correlations spanning the gradient hyperparameter space found the resulting gradients to be generally robust to parameter selection (Fig. S7). Differences in vertex-wise channel coverage between both cohorts partly contributed to node-wise displacement in the embedding space (Fig. 4B). Channel coverage in the Replication cohort was notably higher in lateral temporal regions, and lower in insular and frontal cortices. The coverage difference map across datasets was significantly correlated with average node rank changes between each gradient (r = -0.179; p = 0.038), suggesting that regional differences in channel coverage across the datasets may partially account for node shifts in the resulting embedding space.

Bivariate associations between each MRI modality and embedding distances specific to the Replication cohort again highlighted strongest correlations with GD (r = -0.323; p = 0.002), and relatively balanced correlations with SC (r = -0.233; p < 0.001), FC (r = -0.264; p < 0.001), and MPC (r = -0.256; p < 0.001) (Fig. 4C). This stronger contribution of GD was replicated in all multilinear models. The contribution of SC in these models was negligible when either FC or MPC was included as a model predictor. Of all modalities, SC was the most dissimilar across cohorts (cross-cohort correlations for each matrix: r_GD = 0.999; r_SC = 0.847; r_FC = 0.960; r_MPC = 0.947). Bivariate correlations between Replication cohort embedding distance and connectome matrices computed from our sample of healthy participants (as used in our Discovery cohort analysis) recovered the overall pattern of findings observed in the Discovery cohort, with descending correlation strength with GD (r = -0.323), SC (r = -0.268), FC (r = -0.250), and MPC (r = -0.245). These findings indicate modality-specific idiosyncrasies potentially sensitive to epilepsy-related connectome alterations when comparing findings across datasets. Leveraging within-sample data furthermore highlighted overall improved prediction of neural dynamics similarity over all edges with increasing model complexity, particularly when including FC and MPC as model predictors (Fig. 4C). When using iEEG and connectivity data extracted from distinct samples (as done in our Discovery cohort analysis), this effect was only observed in transmodal cortices and when selectively fitting models at the node-level. In additional replication analyses implemented at the individual level, we also found stronger associations with GD relative to other modalities, followed by SC and balanced correlations with FC and MPC (Fig. S11). Furthermore, reconstructing the embedding space while iteratively excluding each patient showed that embedding distances between nodes related to channel-specific PSD similarity in the held-out patient (Fig. S12). Together, these findings demonstrate the overall replicability of our main findings across cohorts as well as within-patient generalizability, while additionally identifying nuances in long-range wiring principles of neural dynamics in different analytical designs (between- versus within-subject).

Methods
The present work relies on open-sourced datasets. The studies under which these datasets were collected complied with the Research Ethics Board of McGill University and the Montreal Neurological Institute and Hospital. Written and informed consent were obtained from all participants to take part in the research MRI protocol, and they received monetary compensation for their participation in MRI studies. As iEEG data was acquired as part of each patient’s clinical workup for intractable epilepsy, patients were not compensated for this aspect of data collection.

Dataset overview
(a) Atlas of the human intracranial electroencephalogram (Discovery cohort). The MNI open iEEG atlas is a multicentre initiative providing openly available human intracranial recordings acquired during different states of vigilance8,18. This dataset includes iEEG recordings from a sample of 106 patients (52 women; mean ± standard deviation (SD) age = 33.1 ± 10.8 years) with intractable epilepsy, totaling over 1700 data points taken from non-epileptic channels mapped to a common stereotactic space. The present study focused on segments recorded during resting wakefulness with eyes closed. All data were collected during a standard clinical protocol implemented in all participating sites.
(b) MNI clinical dataset (Replication cohort). We studied a sample of 18 patients (11 women; mean ± SD age = 33.89 ± 9.02 years; mean ± SD age at seizure onset 13.28 ± 7.47) treated at the Montreal Neurological Institute and Hospital for intractable epilepsy. Patients underwent multimodal MRI at a field strength of 3 Tesla as part of a research protocol. Imaging was completed prior to clinical iEEG investigations. Patients showed heterogeneous suspected locations (frontal: n=3; temporal: n=6; parietal: n=4; occipital: n=1; insula: n=4) and laterality (right: n=11; left: n=6; bilateral: n=1) of the epileptogenic focus. An epileptologist (TA) reviewed clinical and iEEG data and labeled channels in the seizure onset zone, which were excluded from our analysis. Nine patients showed suspected focal cortical dysplasia (FCD; histopathological confirmation of FCD type II in six patients; FCD type I in one patient) and one patient presented with a ganglioglioma, while the others were deemed MRI-negative with no available surgical tissue for histopathological analysis. 13/18 patients underwent surgery to remove the suspected epileptic focus (surgical resection or radiofrequency thermocoagulation), with eight patients achieving seizure freedom after surgery (Engel Ia or Ib).

Electrophysiological data pre-processing and gradient analysis
(a) iEEG data pre-processing. All data segments from bipolar channels provided in the MNI open iEEG atlas underwent identical pre-processing18. In brief, these steps included the application of a band-pass filter at 0.5–80 Hz, downsampling to 200 Hz, application of an adaptive filter to reduce powerline interference (50 or 60 Hz depending on site) and demeaning each data segment. When no continuous and artifact-free 60-second segments were found for a given channel, multiple discontinuous segments from the channel were concatenated with a 2-second zero-padded buffer between segments. Following the procedure established in this dataset, we used Welch’s method to compute each channel’s PSD from 2-second blocks with 1-second overlap, weighted by a Hamming window. Each channel’s PSD was normalized so the total power was equal to one. As for the MNI clinical dataset, two board-certified epileptologists (BF, DM) selected data segments (60–120 second duration) free of ictal events and artifacts from periods of resting wakefulness with eyes closed. Like the MNI open iEEG atlas, selected epochs could be continuous or discontinuous, in which case data segments were concatenated with a 2-second zero-padded buffer. All recordings were acquired at a sampling rate of 2000 Hz but were downsampled to 200 Hz to match the pre-processing applied to the MNI open iEEG atlas. This cut-off allowed us to cover commonly studied frequency bands (from delta to gamma range) while eliminating high-frequency noise. However, this also filters potentially high-frequency content such as physiological ripples and fast ripple activity118.
(b) Anatomic localization and co-registration of channels. For the MNI open iEEG atlas, all channel coordinates were provided in a common stereotactic space (ICBM2009 asymmetric) alongside symmetric left and right cortical surface meshes (119k vertices per hemisphere). As this surface is symmetric, channels located in the left hemisphere were flipped to the right hemisphere to improve spatial coverage. As a result, all downstream analyses, including the computation of gradients of neural dynamics and correlations to MRI metrics, are implemented in a single hemisphere. Each channel was attributed to the closest vertex on the cortical surface determined by the shortest Euclidean distance between each vertex and the channel’s position in volumetric space. To account for the effects of volume conduction, channel PSDs were extrapolated to surrounding vertices within a 5mm radius (channel region). Distances were measured using the exact geodesic algorithm119 implemented in the Pygeodesic package (v0.1.7). Channel regions underwent surface-based registration to the fsLR surface template (~32k vertices) using nearest-neighbor interpolation.
For each patient in the MNI clinical dataset, peri-implantation scans were co-registered to pre-surgical T1w images pre-processed using micapipe (v0.2.3)120. We applied a label-based affine registration121,122 to align peri-implantation T1w (9/18 patients) and T2w (2/18 patients) to the pre-surgical native space. When only computerized tomography (CT) was available for peri-implantation imaging (6/18 patients), scans were directly aligned to pre-surgical space without prior tissue-type segmentation. We applied the resulting transform to each patient’s electrode contact coordinates, bringing them to the patient-specific native pre-surgical space. Retaining only contacts overlapping with neocortical gray matter, channel coordinates were defined as the midpoint between the two contacts making up the channel. Neocortical surface segmentations were generated from pre-surgical T1w scans using FreeSurfer 6.0123–125, and were resampled to retain the geometry of the fsLR surface template (~32k vertices). Channels were attributed to the closest vertex on the neocortical surface mesh as measured with Euclidean distance (channel vertex). As with the MNI open iEEG atlas, channel PSDs were extrapolated to surrounding vertices within a 5mm radius (channel region) as defined by the exact geodesic algorithm119 implemented in the Pygeodesic package (v0.1.7). Channels involved in seizure onset (i.e., seizure onset zone) were identified by an experienced epileptologist (TA) and excluded from further analyses.
(c) Generating gradients of neural dynamics. Channel PSDs were log-transformed and mapped to the right hemisphere labels of the Schaefer-200 parcellation57, resulting in 100 parcel-wise PSDs. This parcellation resolution was selected as it provided complete sampling of the right hemisphere. To map channel-level PSDs to this parcellation, we first identified all channel regions overlapping with a given parcel. For each frequency (0.5–80 Hz, in increments of 0.5), the amplitude of the PSD for all channels covering the parcel was averaged. More specifically, we implemented a weighted average where amplitudes were weighted according to the extent of the channel region’s coverage of the parcel (the number of vertices that channel contributed to the parcel), resulting in a parcel-level PSD. Parcel-wise PSDs were cross-correlated while controlling for the average PSD computed across all channels and underwent Fisher R-to-Z transformation. We used the BrainSpace toolbox126 to convert the Z-transformed correlation matrix to a normalized angle affinity matrix, and applied diffusion map embedding to generate eigenvectors or gradients of neural rhythm similarity across the cortex using default sparsity (keeping only the top 10% of weights) and diffusion (α=0.5) parameters. Of note, systematic correlations spanning the gradient hyperparameter space found the resulting gradients to be generally robust to parameter selection in both Discovery and Replication cohorts (Fig. S7). The same pipeline was applied in the Replication sample, which was additionally aligned to the embedding space of the Discovery dataset. We compared two complementary alignment strategies. First, using singular value decomposition, we computed a transformation matrix enabling translation, rotation, and scaling of the Replication cohort gradients to optimize their similarity to the Discovery cohort gradients. The results of this approach are presented in Fig. 4. We also applied a more conservative approach where gradients from the Replication cohort were simply reordered to minimize their distance to gradients computed in the Discovery cohort (Fig. S10). Consistency of gradient topographies across datasets was quantified using Spearman rank correlations and statistical significance was determined using spin permutation tests.

MRI data acquisition
(a) Participants. Data were collected in a sample of 50 healthy volunteers (23 women; age mean ± SD = 30.28 ± 5.74 years) denying any history of neurological and psychiatric illness. The Ethics Committee of the Montreal Neurological Institute and Hospital approved the study. Raw and pre-processed connectomes data are openly available as part of the MICA-MICs dataset55 (osf.io/j532r/, note that updated outputs generated with micapipe v0.2.3120 were used in the present study). Patients included in the MNI clinical dataset (Replication cohort) underwent the same multimodal imaging protocol as healthy volunteers.
(b) MRI acquisition. Scans were acquired at the McConnell Brain Imaging Centre of the Montreal Neurological Institute and Hospital on a 3T Siemens Magnetom Prisma-Fit equipped with a 64-channel head coil. We acquired two T1-weighted (T1w) scans with identical parameters using a 3D magnetization-prepared rapid gradient-echo sequence (MPRAGE; 0.8mm isovoxels, TR = 2300ms, TE = 3.14ms, TI = 900ms, flip angle = 9°, FOV = 256×256mm2). T1w scans were inspected to ensure minimal head motion before undergoing further processing. Quantitative T1 (qT1) relaxometry data were acquired using a 3D-MP2RAGE sequence (0.8mm isotropic voxels, 240 sagittal slices, TR = 5000ms, TE = 2.9ms, TI 1 = 940ms, T1 2 = 2830ms, flip angle 1 = 4°, flip angle 2 = 5°, iPAT=3, bandwidth = 270Hz/px, echo spacing = 7.2ms, partial Fourier = 6/8). We combined two inversion images for qT1 mapping to minimize sensitivity to B1 inhomogeneities and optimize reliability127,128. A 2D spin-echo echo-planar imaging sequence with multi-band acceleration was used to obtain DWI data, consisting of three shells with b-values 300, 700, and 2000 s/mm2 and 10, 40, and 90 diffusion weighting directions per shell, respectively (1.6mm isotropic voxels, TR = 3500ms, TE = 64.40ms, flip angle = 90°, refocusing flip angle = 180°, FOV = 224×224mm2, slice thickness = 1.6mm, multi-band factor = 3, echo spacing = 0.76ms). b0 images acquired in reverse phase encoding direction were used for distortion correction. One 7min resting-state function MRI (rs-fMRI) scan was acquired using multiband accelerated 2D-BOLD echo-planar imaging (3mm isotropic voxels, TR = 600ms, TE = 30ms, flip angle = 52°, FOV = 240×240mm2, slice thickness = 3mm, mb factor = 6, echo spacing = 0.54ms). Participants were instructed to keep their eyes open, look at a fixation cross, and not fall asleep. Distortion correction of rs-fMRI data was performed using two spin-echo images with reverse phase encoding (3mm isotropic voxels, TR = 4029ms, TE = 48ms, flip angle = 90°, FOV = 240×240mm2, slice thickness = 3mm, echo spacing = 0.54ms, phase encoding = AP/PA, bandwidth = 2084Hz/Px).
(c) Replication 7T MRI dataset. Data were collected in a sample of 10 healthy volunteers (6 women; 26.60 ± 4.60 years) denying any history of neurological and psychiatric illness. Raw and pre-processed connectome data are openly available as part of the MICA-PNI dataset. Only the first available study session for each participant was analyzed as part of this study. Detailed acquisition and pre-processing parameters are available in the MICA-PNI descriptor article65.

MRI data processing
Imaging data were processed with micapipe (v0.2.3; micapipe.readthedocs.io/), an openly available pipeline for multimodal MRI processing, data fusion, and multiscale connectome generation120. This pipeline is notably adapted to sample features in modality-specific spaces, therefore limiting data blurring that can occur with co-registration to standard template space.
(a) T1w pre-processing and geodesic distance measurements. Structural processing was carried out using several software packages, including tools from AFNI, FSL, and ANTs129. Each T1w scan was deobliqued and reoriented to standard neuroscience orientation (LPI: left to right, posterior to anterior, and inferior to superior). Both scans were then linearly co-registered and averaged, automatically corrected for intensity nonuniformity130, and intensity normalized. Cortical surface segmentations were generated from native T1w scans using FreeSurfer 6.0123–125. We computed individual GD matrices along each participant’s native cortical midsurface using workbench tools131,132, specifically the -surface-geodesic-distance function. First, a centroid vertex was defined for each cortical parcel by identifying the vertex with the shortest summed Euclidean distance from all other vertices within its assigned parcel. The GD between centroid vertices and all other vertices on the native midsurface mesh was computed using Dijkstra’s algorithm. Notably, this implementation computes distances not only across vertices sharing a direct connection, but also across pairs of triangles which share an edge to mitigate the impact of mesh configuration on calculated distances. Vertex-wise GD values were then averaged within parcels. Of note, longer distances measurements may be obtained here due to our use of high-resolution meshes, compared to techniques using surfaces resampled according to a lower-resolution template’s geometry.
(b) DWI pre-processing and tractography-derived SC. DWI data were pre-processed using MRtrix133,134. DWI data were denoised135,136, underwent b0 intensity normalization130, and were corrected for susceptibility distortion, head motion, and eddy currents using a reverse phase encoding from two b=0 s/mm2 volumes. Required anatomical features for tractography processing (e.g., tissue type segmentations, parcellations) were non-linearly co-registered to native DWI space using the deformable SyN approach implemented in ANTs137. Diffusion processing was performed in native DWI space. Structural connectomes were generated with MRtrix from pre-processed DWI data133,134. We performed anatomically-constrained tractography using tissue types (cortical and subcortical gray matter, white matter, cerebrospinal fluid) segmented from each participant’s pre-processed T1w images registered to native DWI space138. We estimated multi-shell and multi-tissue response functions139 and performed constrained spherical deconvolution and intensity normalization140. We generated a tractogram with 40M streamlines (maximum tract length = 250; fractional anisotropy cutoff = 0.06). We applied spherical deconvolution informed filtering of tractograms (SIFT2) to reconstruct whole brain streamlines weighted by cross-sectional multipliers141. The reconstructed cross-section streamlines were mapped to the Schaefer-200 parcellation, which was also warped to DWI space using the previously generated transforms137. The connection weights between nodes were defined as the weighted streamline count. Group-level SC was computed using distance-dependent thresholding to preserve the distribution of within- and between-hemisphere connection lengths found in individual subjects142. However, SC weights of the resulting matrix remain highly skewed, with few very strong connections and many weak or zero connections. As this property of the matrix violates assumptions of normality, which is problematic for the statistical tests applied in the present work, connection weights in the resulting group-level were log-transformed to reduce connectivity strength variance.
(c) Resting-state functional connectivity (FC). After discarding the first five volumes of each subject’s rsfMRI run to ensure magnetic field saturation, images underwent reorientation, as well as motion and distortion correction. Motion correction was performed by registering all timepoint volumes to the mean volume, while distortion correction leveraged main phase and reverse phase field maps. Nuisance variable signal was removed using an ICA-FIX143 classifier and by regressing timepoints with high motion detected from motion outlier outputs provided by FSL. Volumetric time series were averaged, and the mean volume was registered to each subject’s T1w image using a label-based affine registration121,122. We applied the resulting transform to previously computed surface files to sample volumetric time series using trilinear interpolation. Time series were averaged within nodes defined by the Schaefer-200 parcellation and were cross-correlated to generate individual-specific functional connectomes. Correlation values (Pearson’s r) subsequently underwent Fisher-R-to-Z transformations. Group-level functional connectivity (FC) was computed by averaging matrices across participants.
(d) qT1 pre-processing and MPC. A series of 14 equivolumetric surfaces were constructed for each participant between pial and white matter boundaries. We aligned qT1 and native FreeSurfer space T1w volumes using a label-based affine registration121,122. We applied the resulting transform to previously computed equivolumetric surface files to sample image intensities in native qT1 space. This resulted in vertex-wise microstructural intensity profiles, which were used to calculate individual microstructural profile similarity matrices35,37. Intensity profiles were averaged within parcels, cross-correlated across using partial correlations controlling for the average cortex-wide intensity profile, and log-transformed37,48. Resulting matrices thus represented participant-specific similarity matrices in myelin proxies across the cortex. Matrices were subsequently averaged for group-level analyses.

Multiscale connectomics of neural dynamics
(a) Bivariate associations between neural dynamics and multiscale network metrics. The similarity of neural dynamics was quantified by computing the Euclidean distance between all node pairs in the gradient space. For this purpose, we retained gradients explaining a cumulative sum of 50% of the variance (G1-G14; note that results are displayed in a 2D representation composed of G1 and G2). Node proximity in the embedding space was cross-referenced with measures of (i) cortical geometry indexed by GD, (ii) tractography-based SC, (iii) resting-state FC and (iv) microstructural similarity (MPC) by correlating the upper triangular edges of each matrix. For each modality, statistical significance was determined by comparing empirically observed correlations to those obtained with surrogate matrices in which each row had been shuffled using spatial autocorrelation-preserving null models63,64. Consistency of cross-modality associations were also assessed by comparing embedding distance distributions across thresholded network matrices preserving the 5–95% strongest edges (in increments of 5%). We also assessed node-level bivariate associations by repeating correlations independently for each matrix row, and cross-referencing resulting correlation coefficients with those obtained spatial null models63,64.
(b) Multiple linear regression models. Multiple linear regression models were implemented with distinct sets of modality predictors (i.e., a baseline model including GD only; partial models including GD and one other modality, or two other modalities from: SC, FC, or MPC; a full model including all four modalities). Models were additionally implemented using either all edges of the upper triangular part of each matrix or on a region-specific basis. Goodness-of-fit of each model was assessed using adjusted R-squared. In models implemented over all edges, we compared the importance of each feature from its t-statistic, i.e., each feature’s coefficient divided by its standard error. We assessed the generalizability of model parameters by bootstrapping, specifically subsampling 90% of edges to estimate model parameters and testing the goodness-of-fit over the remaining 10% of edges. Region-specific models were also implemented to quantify node-level gains in adjusted R-squared, which were spatially correlated with the principal gradient of resting-state functional connectivity derived from an independent sample as provided in the BrainStat toolbox56,66.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper builds on a large body of prior studies that aim to decipher how the brain's structure shapes its dynamic activity. The authors present compelling evidence that geometric and microstructural properties jointly influence cortex-wide neural rhythms based on intracranial EEG (iEEG) recordings. Their results suggest that transmodal cortices exhibit uniquely complex dependencies on cortical geometry and microstructural similarity. At the core of the analysis is the use of diffusion embedding applied to a similarity matrix derived from power spectral measures of iEEG. The study of iEEG data in this context is challenging and novel, and the authors demonstrate meticulous data curation and processing steps to ensure the accuracy of their results. Sufficient detail is provided in the methods section for the study to be replicated in future work. The authors claim that the obtained embedding represents an \"intrinsic\" coordinate system of human cortical dynamics. The study is comprehensive and includes various complementary analyses, e.g., relating the embedding distance to other connectomes—analyses that help to support the claims. The obtained embedding is, to a rather reasonable degree, replicated in an independent cohort, which makes it a bit difficult to fully accept the term \"intrinsic\" from this reader's perspective. Overall, the work presents novel findings, and I enjoyed reading it. I have two semi-major comments, followed by some minor notes.\n[1] (Robustness of derived gradients) Given the high spatial resolution of the mapping between iEEG electrode positions and the cortical surface, the choice of a single atlas to parcellate the cortical surface in this context is not well justified. To demonstrate the robustness of your results, please show how your main finding (the color-coded 2D gradient embedding comparing G1 to G2, along with its projection on the parcellated cortical surface, i.e., Figure 1E) would appear for a lower-resolution atlas (Schaefer-100) and a higher-resolution atlas (Schaefer-300/400 and/or Glasser360). Please feel free to argue against this analysis if you deem it unnecessary and believe it is sufficient to address this point only in the discussion section of the paper.\n[2] (Replication results) The authors argue that the gradients are highly replicable. However, when comparing Fig. 2A to Fig. 4B, the G1-G2 space embedding and its surface projection appear to deviate substantially from each other—e.g., see the lateral frontal and occipital lobes. Furthermore, the authors show that after Procrustes alignment, the G1-G13 gradients from the discovery and replication cohorts are significantly correlated, although the actual correlation values are moderate (average correlation: ~0.5); please also see the first minor comment below regarding Procrustes alignment. These observations make it difficult for the reader to accept the claim that the obtained embedding represents an \"intrinsic\" coordinate system of human cortical dynamics, as stated in the abstract. For a more holistic comparison, it would be helpful to compare the embedding distances (ED) for the discovery and replication cohorts; that is, reproduce the middle and right-hand plots from Figure 2A (discovery cohort) for the replication cohort. Instead of the right-hand plot, you may simply report the correlation between the ED matrices and test whether the correlation is significant, using a null model similar to that used in Figure 2B–E. Regarding the results shown in Figure 3 and Figure 4C: are the same connectivity matrices used for the analyses in both the discovery and replication cohorts? Or are the four connectivity matrices separately computed using data from the discovery and replication cohorts? If they are separately computed, please clarify and discuss how your results might be affected by mere differences between the two versions of the connectivity matrices—independent of the differences in the derived gradients for the discovery and replication cohorts.\nMisc comments:\nRe: Procrustes alignment; Line 541: \"The same pipeline was applied in the replication sample, which was additionally aligned to the embedding space of the discovery dataset using Procrustes alignment.\" What exactly was aligned using the Procrustes transform (PT)? The resulting gradients? Please clarify what you mean by \"align,\" since with PT, given two sets of vectors, there are two common approaches. One approach is to only reorder the indices so that vectors in the two sets optimally match. The second approach is to actually transform (i.e., modify) the vectors in one set—via scaling, translation, and rotation—to optimally match the vectors in the other set. If you have used the second strategy, then the resulting transformed gradients are more difficult to interpret. However, if you have used the first approach, that is much more conservative, and the correlation matrix shown in Figure 4B is quite neat. In either case, please explicitly state which of the two approaches you have used.\nAbstract: \"our results support gradual shifts in the balance of local and macroscale constraints on cortical dynamics\" Please revise this sentence as it is not clear what is meant by \"gradual shifts\", e.g., gradual along what axis? time? space? both?\nAlthough the discovery cohort iEEG atals is from prior work, please provide some minimal demographic details about the 106 patients that were used to build the atlas. If sex and age distribution is noteably different from that of your replication cohort, this can potentially help to further argue why the gradients obatined from the two settings do not highly correlate.\npage 15, line 509: \"As this surface is symmetric, channels located in the left hemisphere were flipped to the right hemisphere to improve spatial coverage.\" Maybe helpful to explicutely state that the gradients are are derived in one hemisphere, as in prior related work.\nline 535: \"This parcellation resolution (Schaefer200) was selected as it provided complete sampling of the right hemisphere.\" It can be insightful if robustness of results was also verfified with Scahefer100 and Schaefer300, see comment [1] above.\nline 534: \"Parcel-wise PSD amplitudes at a given frequency were derived from a weighted-average considering each channel’s parcel coverage.\" This is a bit confusing. How does this operation relate to the extrapolation that you describe earlier above?, on line 512: \"To account for effects of volume conduction, channel PSDs were extrapolated to surrounding vertices within a 5mm radius (channel region).\"\nline 538: \"applied diffusion map embedding to generate eigenvectors or gradients of neural rhythm similarity across the cortex using default sparsity (keeping only the top 10% of weights) and diffusion (alpha=0.5) parameters\" To what degree are the gradients that you derive stable relative to changes in these two parameters?\nFigure 2: It is not very clear what the histograms shown as insets in B-E (related to null models tests for statiscal significance) are showing; please label teh axes and specify what the dashed red is.\nFigure 2: In the right-hand plots, it is not clear what \"mean\" refers to on the y-axis. It is also not clear what is meant by \"connected\" and \"not connected\".\nFigure 3A-B, Figure 4C: its it not straightforward to understand the color-encoding you are using to link the barplots (T-statistic) to the colormap used for the four connectovity matrices. Please enhance the legend.\nHarry H. Behjat\n(Remarks on code availability)\nI only briefly checked the code and I can confirm that it includes necessary scripts that provide details for replicating the main findings in combination with the description of methods in the paper.", "reviewer_2": "The authors use a dimensionality reduction technique to evaluate the contribution of brain geometry and connectivity to regional dynamics. While the paper is well written and the analysis is performed carefully, I have many concerns about the approach to integrate data across subjects that leave me less enthusiastic of the interpretation of the findings.\nI fully understand the difficulties of working with iEEG data that is clinically acquired. In each patient, electrodes will be placed for clinical reasons, and it is therefore difficult if not impossible to achieve similar coverage of brain activity recordings across patients – each patient will have recordings from different brain areas. In order to achieve full coverage of brain areas, here the authors employ a technique of mapping the data from all patients onto a single common space. While they are not the first to use this approach, I am not convinced that it is appropriate for the types of questions that the authors want to answer about the data. By looking at group level data, the authors are fundamentally choosing to ignore individual differences across the population that I believe would be essential to preserve in order to truly ask about questions of how brain geometry, connectivity, and dynamics are related. We know that both structural and functional brain connectivity can be used to perform fingerprinting of individuals. I have the following very specific concerns about how individual variability is ignored in this data.\n1. In the Discovery Cohort, not only are all subjects mapped onto a single common space, but also, the left hemisphere is flipped to be on the right hemisphere and only one side of the brain is considered. This adds further distortion to the data as it assumes perfect symmetry between both left and right dynamics and connectivity, on top of ignoring individual differences between subjects. Although there is some mention of using patient specific data in the Replication Cohort, from what I can tell, for the iEEG data, ultimately this same group level mapping was performed.\n2. Structural connectivity analysis was also performed at the group level, and additionally, the connection weights in the group level matrix were log-transformed to reduce variance in connectivity strength. However, presumably, the differences in connectivity strength are important! Again, these are related to individual differences.\n3. Functional connectivity for rsfMRI data was also group averaged, ignoring patient-specific differences.\n4. At the microscale, it seems that patient-specific matrices were obtained. I don’t understand how this analysis can be considered comparable to the others that were performed at the group level.\n5. The data used in this study all comes from patients with epilepsy. As each patient has a different focus of their seizures and potential lesions in different locations, even when the channels are selected from “healthy” regions, it is to be expected that the underlying brain connectivity and resulting dynamics are affected. This should lead to increased individual variability within each cohort, making it even less reasonable to perform a group level analysis.\nI assume that the reason that the analysis was performed on two cohorts was to alleviate some of these concerns about group level analysis. However, I still believe that while the “average” analysis might be similar in 2 cohorts, the actual findings could simply not hold in either. In any sort of population that has sub-groups within the population or even just simply lots of variability across the population, one can detect relationships between the average geometry and structure with average dynamics that simply do not exist for any of the individuals in the population.\nWhile I understand the desire to fully cover the brain with the iEEG data, I think that to do this study properly, one would have to perform the gradient analysis separately for each individual, and accept that between individuals, the spatial coverage would be different. I realize that this makes it hard to collapse the analysis across individuals and perform regional analysis, but at least one could look at the correlations between gradients and structure at the individual level and later try to find where there was spatial overlap between patients to do a limited spatial analysis.\nIn summary, because this work is based on group averaged data, making it hard to know where true correlations might exist at the individual level, I do not feel that it has the ability to be of high impact to the field.\n(Remarks on code availability)\nBecause I have issues with the fundamental assumptions about the data set, I did not review any code.", "reviewer_3": "In this work, Royer et al. investigate how regional dynamics, as measured by iEEG, are constrained by geometric, microstructural, and connectomic features. They find that geodesic distances strongly constrain brain dynamics, but this relationship is weakest in transmodal cortices, consistent with other findings showing that neuronal organisation in these regions clearly departs from that in the rest of the brain. The authors further strengthen the reliability of their findings by using a replication dataset, where they observe similar results. This is a very well-written paper, with excellent figures and good science! (cid:1) I have no major concerns, and my comments are largely restricted to suggestions on how the authors could improve the manuscript and enhance clarity where needed.\n1. As the authors will be oh so very well aware, there are many papers now correlating gradients with all the other gradients and all the other brain related things, so it is easy to get lost in what has been done before and where something clearly novel is actually being done. To the best of my knowledge, this is the first study using iEEG, correct? If so, I think this is worth highlighting more explicitly as a key novelty. I suggest making it very clear in the Introduction how this technique provides insights that previous methods could not. At present, the importance of using iEEG is somewhat underemphasised.\n2. All iEEG data used in this study were obtained from patients. Given that these individuals have epilepsy, and would be expected to exhibit altered patterns of brain activity, how applicable are the findings to non-clinical populations? Or is it assumed that, outside of seizure episodes, patients’ brain activity patterns are comparable to those of non-epileptic individuals?\n3. In the main analysis using the Discovery Cohort, iEEG data from patients were combined with SC/FC/MPC measures calculated in healthy individuals. However, for the replication cohort, were both the iEEG and MRI-derived measures calculated in the same patient cohort (i.e., was SC/FC/MPC calculated for the MNI clinical cohort and used in the analysis)? I believe this is what was done, but if so, could it be made more explicit? More broadly, the main concern here is that data from healthy and clinical populations are being combined and used to make general inferences about brain organisation. This issue is briefly acknowledged on lines 445–446 (“as well as epilepsy-related changes in SC that could only be detected in a within-patient design”), but I think a clearer discussion or acknowledgement of this limitation should be made.\n4. All analyses are conducted in a single hemisphere (due to understandable technical constraints). Especially when we consider that the pattern of inter-hemispheric connections is very different to intra-hemispheric connections and virtually all such inter-hemispheric connections are long-range, which as the authors note seem to have some special role in shaping communication across brain networks (thus inter-hemispheric dynamics might be very different!). At a minimum, would be worth noting as a limitation.\n5. I found it interesting that the first two gradients explain relatively little variance, whereas most other decompositions of brain features based on regional similarity or covariance explain much more. Do the authors have any intuition as to why this might be?\n6. I assume that the 50% cutoff for selecting gradients was chosen because it's a nice big round number. Do the results change meaningfully if only the first 3–5 gradients are used? I realise re-running all analyses would be burdensome (though it could serve as a nice reproducibility check!). However, if you can show that the correlation between the full embedding matrix (computed with G1–G14) and a reduced matrix (e.g., computed with just G1–G2, G1–G3, etc.) remains high, that would suggest the results are robust to the number of gradients included.\n7. In Figure 2A, the matrix is labelled as indicating “Euclidean distance,” but in the following figure it is described as “embedding distance.” While the embedding distance is indeed a Euclidean distance, for consistency I suggest updating the colourbar label to “Embedding distance.” Also (as I am sure you are aware) “Euclidean distance” is typically understood to refer to physical distances between parcels, so to some people (i.e., me) labelling a matrix of the distances in high dimensional gradient space as “Euclidean distance” may be confusing.\n8. Also for Figure 2A: in the brain projection, could the seed parcel be made more visually prominent or central? The selected hippocampal parcel is somewhat obscured, and I missed it at first glance.\n9. Is the geodesic distance measured in millimetres? The colour scale suggests a maximum distance of 550 mm, which seems quite high. In my own data, the maximum geodesic distance across individuals tends to be around 250–300 mm. Could you please double-check?\n10. For Figure 3B, it wasn’t immediately clear that the colours of the bars correspond to the colours of the matrices for each feature. I think there’s enough space in the legend to explicitly state what each colour represents.\n11. Also for Figure 3B: are the ridge plots showing the distribution of bootstrap results? This isn’t stated clearly. Is the vertical line the mean or median of the bootstrap distribution, and is the dot the non-bootstrapped result?\n12. For Figure 3C: in the far-right plot showing correlations with FC G1, is the variable shown the correlation between nodal FC G1 values and the full model–baseline difference from the previous plot? If so, could you ensure the name of this variable is consistent across figures to help with interpretation?\n13. Figure 4A: what does the rainbow line overlaid on the brain slices represent?\n14. The paragraphs in the Discussion section are quite long. Breaking them up would improve readability.\nStuart Oldham, MCRI\n(Remarks on code availability)\nI quickly checked the code, seems to all run and reproduce the results in the paper"}, "author_response": "RESPONSE TO REVIEWERS (NCOMMS-25-05787-T)\nWe would like to thank the Editors for giving us the opportunity to submit a revised manuscript, and the\nReviewers for their constructive evaluations. We are grateful to the Reviewers for their thoughtful\ncomments and found the suggestions very helpful to improve our paper. Assessments were overall positive\nand mainly asked for minor additional analyses to ensure robustness of main results, as well as clarifications\nof terminology, methods, and expanded discussions on implications of presented findings. Dr. Behjat’s\n(Reviewer 1) concerns focused on ensuring the replicability and robustness of results across processing\nparameters and cohorts, which we address in supplementary analyses and figures. Reviewer 2 questioned\nour focus on group-level inference to uncover general principles of brain organization. While we\nacknowledge and discuss the limitations of our approach, we address the Reviewer’s concerns by including\nadditional analyses implemented at the within-patient level in our Replication cohort, where iEEG and\nmultimodal MRI data were available in the same group of participants. Dr. Oldham (Reviewer 3) also raised\nconcerns on robustness of findings to parameter variations, clarifying contents of the figures, and\ncontextualizing findings with the existing literature. We addressed all suggestions point-by-point and\nhighlighted corresponding manuscript changes in gray.\nREVIEWER #1\nSigned Review (Harry H. Behjat)\n1.1 This paper builds on a large body of prior studies that aim to decipher how the brain's structure\nshapes its dynamic activity. The authors present compelling evidence that geometric and microstructural\nproperties jointly influence cortex-wide neural rhythms based on intracranial EEG (iEEG) recordings.\nTheir results suggest that transmodal cortices exhibit uniquely complex dependencies on cortical\ngeometry and microstructural similarity. At the core of the analysis is the use of diffusion embedding\napplied to a similarity matrix derived from power spectral measures of iEEG. The study of iEEG data in\nthis context is challenging and novel, and the authors demonstrate meticulous data curation and\nprocessing steps to ensure the accuracy of their results. Sufficient detail is provided in the methods\nsection for the study to be replicated in future work. The authors claim that the obtained embedding\nrepresents an \"intrinsic\" coordinate system of human cortical dynamics. The study is comprehensive\nand includes various complementary analyses, e.g., relating the embedding distance to other\nconnectomes—analyses that help to support the claims. The obtained embedding is, to a rather\nreasonable degree, replicated in an independent cohort, which makes it a bit difficult to fully accept the\nterm \"intrinsic\" from this reader's perspective. Overall, the work presents novel findings, and I enjoyed\nreading it. I have two semi-major comments, followed by some minor notes.\nWe thank Dr. Behjat for the thoughtful and positive assessment of our manuscript.\n1.2 (Robustness of derived gradients) Given the high spatial resolution of the mapping between iEEG\nelectrode positions and the cortical surface, the choice of a single atlas to parcellate the cortical surface\nin this context is not well justified. To demonstrate the robustness of your results, please show how your\nmain finding (the color-coded 2D gradient embedding comparing G1 to G2, along with its projection on\nthe parcellated cortical surface, i.e., Figure 1E) would appear for a lower-resolution atlas (Schaefer-100)\nand a higher-resolution atlas (Schaefer-300/400 and/or Glasser360). Please feel free to argue against\nthis analysis if you deem it unnecessary and believe it is sufficient to address this point only in the\ndiscussion section of the paper.\nand G2 (individual gradients and 2D embedding space) when using unparcellated data, or parcellating the\ndata with the Glasser-360 atlas as well as a 200-node subdivision of the aparc parcellation. These\nparcellations were selected to assess the replicability of findings across atlases defined by functional\n(Schaefer-200), multimodal (Glasser-360), and anatomical (aparc) criteria. Moreover, we show in FIGURE\nS4 that the topography of G1 and G2 is preserved when applying the Schafer-100 parcellation to bilaterally\ndistributed channels. We expanded sections referring to these analyses in the revised manuscript, and added\nquantitative comparisons of G1/G2 obtained from each supplementary analysis (p.6):\n“A series of control analyses were implemented to assess the robustness of the described gradients and embedding\nspace. First, results were robust to the use of different parcellations schemes (FIGURE S2). Indeed, generating\ngradients from unparcellated PSD data (i.e., channel regions – see MATERIALS AND METHODS) or averaging PSDs\nwithin nodes defined by two different atlases (i.e., Glasser-360 and a subdivision of the aparc parcellation) yielded\nsimilar topographies for both G1 (vertexwise Spearman correlations for unparcellated: r=0.673; Glasser-360:\nr=0.839; aparc-200 r=0.864) and G2 (vertexwise Spearman correlations for unparcellated: r=0.524 ; Glasser-360:\nr=0.580; aparc-200 r=0.645). Furthermore, we could recover the topography of the first 10 gradients when imposing\nan upper bound to the number of channels included in each parcel of the Schaefer-200 parcellation (FIGURE S3),\nsuggesting that our findings could not be attributed to variable channel coverage across nodes. We also replicated\nour findings when preserving the bilateral distribution of channels and applying a lower resolution parcellation\n(Schaefer-100) to maintain broad cortical coverage (FIGURE S4). The high correspondence of G1 (vertexwise\nSpearman correlation r =0.767; r =0.804) and G2 (vertexwise Spearman correlation r =0.645; r =0.671)\nleft right left right\ntopographies with gradients derived from our main analysis provides additional support for the independence of our\nfindings to parcellation choice.”\n1.3 (Replication results)\n1.3.1 The authors argue that the gradients are highly replicable. However, when comparing Fig. 2A to\nFig. 4B, the G1-G2 space embedding and its surface projection appear to deviate substantially from each\nother—e.g., see the lateral frontal and occipital lobes. Furthermore, the authors show that after\nProcrustes alignment, the G1-G13 gradients from the discovery and replication cohorts are significantly\ncorrelated, although the actual correlation values are moderate (average correlation: ~0.5); please\nsee the first minor comment below regarding Procrustes alignment. These observations make it difficult\nfor the reader to accept the claim that the obtained embedding represents an \"intrinsic\" coordinate\nsystem of human cortical dynamics, as stated in the abstract.\nAs suggested, we toned down our descriptions of the magnitude of replicability across cohorts. Our initial\nsubmission discussed potential reasons for these differences, specifically looking at iEEG channel coverage\nspecific to each cohort. We reported that the coverage difference map across datasets was significantly\ncorrelated with average node rank changes between each gradient (r=-0.179; p =0.038), suggesting that\nspin\nregional differences in channel coverage across the datasets may partially account for node shifts in the\nresulting embedding space.\nTo avoid confusion on the meaning of this term, we removed references to “intrinsic” activity in the paper\nwhen referring to the iEEG gradients or embedding space.\n2\n1.3.2 For a more holistic comparison, it would be helpful to compare the embedding distances (ED) for\nthe discovery and replication cohorts; that is, reproduce the middle and right-hand plots from Figure 2A\n(discovery cohort) for the replication cohort.\nWe added a new supplementary analysis (FIGURE S9) to quantify the similarity of embedding distances\nacross both cohorts. We reproduce the visualization in FIGURE 2A for the replication cohort, selecting two\nseeds for visualization (one high/one low replicability) and show the histogram of node-specific correlations.\nThis figures also shows matrix-wide correlations of Discovery/Replication cohort embedding distance and\nsignificance testing with the same null model procedure as used in Figure 2A. We refer to this new figure\nin the revised RESULTS (p.11) and append it below:\n“After alignment, we found a significant correlation between the embedding distance matrices of Discovery and\nReplication cohorts (r = 0.304; p < 0.001), with variable similarity across regions (FIGURE S9).”\nFigure S9. Replication of embedding distances across Discovery and Replication cohorts. (A) Embedding\ndistance matrices for each cohort. (B) Cross-correlating corresponding rows of each matrix showed\nvariable correspondence of inter-node distances, with the distribution of correlations coefficients\n(Pearson’s r) shown in a histogram. Correlation coefficients for each region are also displayed on the\ncortical surface, alongside example seeds with poor correspondence (10th percentile; temporal pole) and\ngood correspondence (90th percentile; central sulcus) across cohorts. (C) Correlating the upper triangular\npart of each embedding distance matrix showed a moderate (r = 0.304; p<0.001) correlation across cohorts.\nStatistical significance was determined by comparing empirically observed correlations to those obtained\nwith surrogate matrices in which each row had been shuffled using spatial autocorrelation-preserving null\nmodels (as in FIGURE 2 of the main manuscript).\n3\n1.3.3 Instead of the right-hand plot, you may simply report the correlation between the ED matrices and\ntest whether the correlation is significant, using a null model similar to that used in Figure 2B–E.\nResults are now reported in the caption of supplementary FIGURE S9 (see response to point 1.3.2)\n1.3.4. Regarding the results shown in Figure 3 and Figure 4C: are the same connectivity matrices used\nfor the analyses in both the discovery and replication cohorts? Or are the four connectivity matrices\nseparately computed using data from the discovery and replication cohorts? If they are separately\ncomputed, please clarify and discuss how your results might be affected by mere differences between the\ntwo versions of the connectivity matrices—independent of the differences in the derived gradients for the\ndiscovery and replication cohorts.\nDifferent connectivity matrices were used in analyses involving the Discovery vs. Replication cohorts. We\nused connectivity matrices from a cohort of healthy controls (n=50) in all Discovery cohort analyses\n(FIGURE 2 and FIGURE 3), while Replication cohort analyses used iEEG and MRI data from the same set\nof patients (n=18; FIGURE 4). All 3T MRI data was acquired on the same scanner. We added a clarifying\nsentence to the RESULTS (p.11):\n“Although all data modalities were available in patients of the Replication cohort, these replication analyses study a\ncompletely distinct sample of participants relative to the Discovery cohort, for both iEEG and MRI metrics.”\nAs suggested, differences in the architecture of the connectivity matrices between the two cohorts could\ninfluence the iEEG-MRI associations. We found lower replicability of SC findings, which we attribute to\nstronger differences across cohorts for this modality. We quantify and discuss connectivity matrix\ndifferences across cohorts in the RESULTS (p.11). While this analysis was included in our original\nsubmission, we have rephrased some sections to clarify:\n“The contribution of SC in these models was negligible when either FC or MPC were included as model predictors.\nOf all modalities, SC was the most dissimilar across cohorts (cross-cohort correlations for each matrix: r =0.999;\nGD\nr =0.847; r =0.960; r =0.947). Bivariate correlations between Replication cohort embedding distance and\nSC FC MPC\nconnectome matrices computed from our sample of healthy participants (as used in the Discovery cohort) recovered\nthe overall pattern of findings observed in the Discovery cohort, with descending correlation strength with GD (r=\nr=-0.323), SC (r=-0.268), FC (r=-0.250), and MPC (r=-0.245). These findings indicate modality-specific\nidiosyncrasies potentially sensitive to epilepsy-related connectome alterations when comparing findings across\ndatasets.”\nThese findings are also revisited in the DISCUSSION (p.15), alongside a broader discussion on limitations\nof studies relying on cross-dataset contextualization approaches:\n“Some differences were nevertheless observed across cohorts. First, bivariate correlations in the Replication cohort\nshowed a stronger relationship between neural dynamic similarity and microstructural similarity, and a lower\ncorrelation with structural connectivity relative to the Discovery cohort. In our multilinear approach, the contribution\nof structural connectivity was generally minimal in the Replication cohort. This may be due to a multitude of factors,\nincluding more precise co-registration and cross-modal alignment in the Replication cohort, differences in channel\ncoverage that may reflect different clinical features of both Discovery and Replication cohorts (e.g., the higher\ncoverage in temporal, relative to frontal, insular, and parietal channels, in the Replication cohort; slight\noverrepresentation of women in the Replication cohort), as well as epilepsy-related changes in SC that could only be\ndetected in a within-patient design.”\n4\n1.4 Re: Procrustes alignment; Line 541: \"The same pipeline was applied in the replication sample, which\nwas additionally aligned to the embedding space of the discovery dataset using Procrustes alignment.\"\nWhat exactly was aligned using the Procrustes transform (PT)? The resulting gradients? Please clarify\nwhat you mean by \"align,\" since with PT, given two sets of vectors, there are two common approaches.\nOne approach is to only reorder the indices so that vectors in the two sets optimally match. The second\napproach is to actually transform (i.e., modify) the vectors in one set—via scaling, translation, and\nrotation—to optimally match the vectors in the other set. If you have used the second strategy, then the\nresulting transformed gradients are more difficult to interpret. However, if you have used the first\napproach, that is much more conservative, and the correlation matrix shown in Figure 4B is quite neat.\nIn either case, please explicitly state which of the two approaches you have used.\nWe clarified that the method used here corresponds to the second option described by the Reviewer. For\ncompleteness, we now report the results of the more conservative approach that was suggested.\nEach approach is now described in the MATERIAL AND METHODS (p.17):\n“We compared two complementary alignment strategies. First, using singular value decomposition, we computed a\ntransformation matrix enabling translation, rotation, and scaling of the Replication cohort gradients to optimize their\nsimilarity to the Discovery cohort gradients. The results of this approach are presented in FIGURE 4. We also applied\na more conservative approach where gradients from the Replication cohort were simply reordered to minimize their\ndistance to gradients computed in the Discovery cohort (FIGURE S10).”\nThe RESULTS have also been updated accordingly (p.11):\n“We also applied a more conservative alignment strategy where Replication cohort gradients were reordered to\noptimally match (i.e., minimize the distance) with the Discovery cohort eigenvectors. With this approach, correlations\nwere overall lower but also indicated generalizability of main findings in the Replication cohort (mean±SD Spearman\nr=0.196±0.154; FIGURE S10).”\nFigure S10. Gradient alignment across Discovery and Replication cohorts. (A) The absolute value of\nSpearman r of gradients across cohort (without alignment). We show the absolute value of correlation\ncoefficients as signs of eigenvectors can be flipped across datasets. (B) Reordering gradients in the\nSpearman r values across the correlation matrix diagonal.\n1.5 Abstract: \"our results support gradual shifts in the balance of local and macroscale constraints on\ncortical dynamics\" Please revise this sentence as it is not clear what is meant by \"gradual shifts\", e.g.,\ngradual along what axis? time? space? both?\nWe modified this sentence in the revised manuscript to better reflect the link we aimed to describe between\nspatial gradients of neural dynamics and the association with multiscale network metrics (p.2):\n“Together, our results suggest that the relative contribution of local and macroscale constraints on cortical dynamics\nvaries systematically across the cortical sheet, specifically highlighting the role of transmodal networks in inter-\nregional cortical coordination.”\n1.6 Although the discovery cohort iEEG atlas is from prior work, please provide some minimal\ndemographic details about the 106 patients that were used to build the atlas. If sex and age distribution\nis notably different from that of your replication cohort, this can potentially help to further argue why\nthe gradients obtained from the two settings do not highly correlate.\nAs suggested, we now report the demographic information of patients included in the Discovery cohort\n(MNI open iEEG atlas) in the MATERIALS AND METHODS (p.16):\n“This dataset includes iEEG recordings from a sample of 106 patients (52 women; mean  standard deviation (SD)\nage = 33.110.8 years) with intractable epilepsy, totalling over 1,700 data points taken from non-epileptic channels\nmapped to a common stereotactic space.”\nThe average age of patients composing each cohort was comparable (Discovery: 33.110.8 years vs.\nReplication: 33.99.0 years). However, the Replication cohort slightly overrepresented women (61.1%\nwomen), while the Discovery cohort included a more balanced sex ratio (49.1% women). We now mention\nthis slight cohort difference as a potential factor that may have influenced replicability (p.15):\n“More broadly, differences in socio-demographic profiles (e.g., slight overrepresentation of women in the Replication\ncohort vs. a more balanced sex ratio in the Discovery cohort) may have influenced resulting iEEG gradients as well\nas MRI-derived measures.”\n1.7 page 15, line 509: \"As this surface is symmetric, channels located in the left hemisphere were flipped\nto the right hemisphere to improve spatial coverage.\" Maybe helpful to explicutely state that the gradients\nare are derived in one hemisphere, as in prior related work.\nWe added this information to both MATERIAL AND METHODS (p.16) and RESULTS (p.6):\n“As a result, all downstream analyses, including the computation of gradients of neural dynamics and correlations to\nMRI metrics, are implemented in a single hemisphere.”\n“We only considered MRI data from the right hemisphere, as all iEEG data had been transposed to the right\nhemisphere to improve cortical coverage.”\n1.8 line 535: \"This parcellation resolution (Schaefer200) was selected as it provided complete sampling\nof the right hemisphere.\" It can be insightful if robustness of results was also verfified with Scahefer100\nand Schaefer300, see comment [1] above.\nresults of supplementary analyses computing gradients of neural dynamics similarity using several different\nparcellations (Schaefer-100, Glasser-360, 200-node subdivision of aparc), as well as in unparcellated data.\n1.9 line 534: \"Parcel-wise PSD amplitudes at a given frequency were derived from a weighted-average considering each channel’s parcel coverage.\" This is a bit confusing. How does this operation relate to the extrapolation that you describe earlier above?, on line 512: \"To account for effects of volume conduction, channel PSDs were extrapolated to surrounding vertices within a 5mm radius (channel region).\"\nWe agree that additional detail on the different steps of this analysis will further clarify our approach. In\nessence, after attributing each channel to a specific vertex, we define a region of the cortical surface\nreflecting the channel’s putative sensitivity profile to co-localizing neural potentials (corresponding to this\n5mm-radius “channel region”). This also considers the spatial uncertainty involved in the successive\nregistration of each channel to stereotaxic space. For the present paper, we implement additional\nmanipulations to improve cortical coverage (i.e., mapping to a single hemisphere and parcellation of PSDs)\nand facilitate comparisons to MRI-derived metrics (p.17):\n“ Channel PSDs were log-transformed and mapped to the right hemisphere labels of the Schaefer-200 parcellation\n(56), resulting in 100 parcel-wise PSDs. This parcellation resolution was selected as it provided complete sampling\nof the right hemisphere. To map channel-level PSDs to this parcellation, we first identified all channel regions\noverlapping with a given parcel. For each frequency (0.5-80Hz, in increments of 0.5), the amplitude of the PSD for\nall channels covering the parcel were averaged. More specifically, we implemented a weighted average where\namplitudes were weighted according to the extent of the channel region’s coverage of the parcel (the number of\nvertices that channel contributed to the parcel), resulting in a parcel-level PSD.”\n1.10 line 538: \"applied diffusion map embedding to generate eigenvectors or gradients of neural rhythm\nsimilarity across the cortex using default sparsity (keeping only the top 10% of weights) and diffusion\n(alpha=0.5) parameters\" To what degree are the gradients that you derive stable relative to changes in\nthese two parameters?\nWe thoroughly explored the stability of gradients across a range of sparsity and diffusion parameters in both\nthe Discovery and Replication cohorts. Spatial correlations for each gradient across parameter sets were\noverall strong, as shown in supplementary FIGURE S7. We refer to this analysis in the RESULTS (p.6):\n“Lastly, spatial correlations assessing the consistency of gradients across a range of hyperparameter combinations\nwere generally very strong (FIGURE S7).”\nWe also added a reference to this figure to the MATERIALS AND METHODS (p.17):\n“Of note, systematic correlations spanning the gradient hyperparameter space found the resulting gradients to be\ngenerally robust to parameter selection in both Discovery and Replication cohorts (FIGURE S7).”\n1.11 Figure 2: It is not very clear what the histograms shown as insets in B-E (related to null models\ntests for statiscal significance) are showing; please label teh axes and specify what the dashed red is.\nWe clarified the nature of the data plotted in the inset histograms in the figure legend (p.9). We also added\na label to the axes of the histograms in FIGURE 2.\n7\n“Statistical significance was determined by comparing empirically observed correlations to those obtained with\nsurrogate matrices in which each row had been shuffled using spatial autocorrelation-preserving null models (62,63).\nThe distribution of correlations obtained from these surrogate matrices are displayed in the inset histogram for panels\nB-E, in which the dashed red line represents the true correlation obtained in the data.”\n1.12 Figure 2: In the right-hand plots, it is not clear what \"mean\" refers to on the y-axis. It is also not\nclear what is meant by \"connected\" and \"not connected\".\nWe added this information to the legend of FIGURE 2 (p.9):\n“For each threshold i, edges were labelled as “connected” if their weight was within the i% highest values for the\nanalysed modality. Corresponding embedding distances at each edge labelled “connected” were then averaged, and\ncompared to the average embedding distances of edges labelled “not connected”. Error bars denote the standard\nerror of the mean.”\n1.13 Figure 3A-B, Figure 4C: its it not straightforward to understand the color-encoding you are using\nto link the barplots (T-statistic) to the colormap used for the four connectovity matrices. Please enhance\nthe legend.\nA sentence explaining the colour correspondence between the cited panels in FIGURE 3 (p.10) and FIGURE\n4 (p.12) was added.\n“The colour of the bars corresponds to the different MRI modalities shown in panel A (blue: GD; orange: SC; pink:\nFC; green: MPC).”\n“The colour of the dots (correlations to ED) and bars (T-statistic) correspond to the different MRI modalities shown\nin the left-hand portion of the panel (blue: GD; orange: SC; pink: FC; green: MPC).”\n1.14 (Remarks on code availability): I only briefly checked the code and I can confirm that it includes\nnecessary scripts that provide details for replicating the main findings in combination with the\ndescription of methods in the paper.\nWe thank the Dr. Behjat for verifying the replication scripts provided in the repository.\nREVIEWER #2:\n2.1 The authors use a dimensionality reduction technique to evaluate the contribution of brain geometry\nand connectivity to regional dynamics. While the paper is well written and the analysis is performed\ncarefully, I have many concerns about the approach to integrate data across subjects that leave me less\nenthusiastic of the interpretation of the findings.\nWe thank the Reviewer for their positive feedback on the writing and analytical precision of the paper. We\naddress their concerns on the group-level nature of our approach in the points below.\n2.2 I fully understand the difficulties of working with iEEG data that is clinically acquired. In each\npatient, electrodes will be placed for clinical reasons, and it is therefore difficult if not impossible to\nachieve similar coverage of brain activity recordings across patients – each patient will have recordings from different brain areas. In order to achieve full coverage of brain areas, here the authors employ a\ntechnique of mapping the data from all patients onto a single common space. While they are not the first\nto use this approach, I am not convinced that it is appropriate for the types of questions that the authors\nwant to answer about the data. By looking at group level data, the authors are fundamentally choosing\nto ignore individual differences across the population that I believe would be essential to preserve in\norder to truly ask about questions of how brain geometry, connectivity, and dynamics are related. We\nknow that both structural and functional brain connectivity can be used to perform fingerprinting of\nindividuals. I have the following very specific concerns about how individual variability is ignored in\nthis data.\nAssessing interindividual consistency of iEEG-MRI associations was beyond the scope of the present paper,\nwhich rather aimed to uncover fundamental principles of brain organization expressed across individuals.\nWe emphasize the leaps enabled by group-level iEEG studies in functional human brain mapping, making\nit a well-suited technique for our objective (discussed in response to point 2.6). Importantly, we developed\ntwo supplementary analyses that directly assess the question of individual-level replicability of our findings.\nThese analyses and corresponding supplementary figures (FIGURE S11 and FIGURE S12) are described in\ndetail in our responses to point 2.7 and 2.8. In summary, these supplementary analyses show that (1) the\nproximity of nodes in the embedding space captures the similarity of channel PSDs in held-out patients,\nand (2) the relative differences in modality-specific iEEG-MRI correlations generalize to the single-patient\nlevel (i.e., strongest associations with GD, followed by SC, and similar associations to FC and MPC).\nFor completeness, we more clearly acknowledge the limitations of our approach in the revised DISCUSSION.\nIntracranial EEG indeed presents challenges for interindividual comparisons compared to techniques with\nconsistent brain coverage. However, it is a powerful technique both in its unique spatiotemporal precision\nand in its ability to directly record neural activity. We thus balance discussions of technical limitations by\nhighlighting how such approaches can still bring valuable insights into our fundamental understanding of\nbrain organization (p.15):\n“Nonetheless, investigations of human brain organization using iEEG are limited by the technique’s unique challenges,\nincluding its restricted accessibility, exclusive use in clinical settings, as well as sparse and patient-specific sampling\n(111). It also remains an open question to what extent focal epileptic activity in one brain region might affect activity\nin the rest the brain and skew the generalizability of iEEG findings to healthy brain function. Recent data sharing\ninitiatives, such as the resource used here (18) and others (8, 112-115), have been instrumental in addressing these\nlimitations. Although iEEG data is collected in a clinical context, pooling carefully curated single-patient data\nsignificantly boosts spatial coverage, thus enabling novel explorations of brain cartography and facilitating the\nidentification of overarching principles of neural organization. These innovations add on to the recognized advantages\nof iEEG regarding its high spatiotemporal precision to study human brain function, as well as its specificity to neural\nsignals over methods offering full brain coverage such as fMRI (113). Studies focusing on group-level inference, such\nas the present work, thus offer new perspectives for human brain mapping to which other modalities may be blind.\nGroup-level paradigms are a valuable first step towards more targeted individual-level inquiries able to clarify\npopulational variations. For example, by testing the consistency of these principles across individuals, extreme\ndeviations from the population norm can be identified and related to clinical presentations (3).”\n2.3 In the Discovery Cohort, not only are all subjects mapped onto a single common space, but also, the\nleft hemisphere is flipped to be on the right hemisphere and only one side of the brain is considered. This\nadds further distortion to the data as it assumes perfect symmetry between both left and right dynamics\nand connectivity, on top of ignoring individual differences between subjects. Although there is some\nmention of using patient specific data in the Replication Cohort, from what I can tell, for the iEEG data,\nultimately this same group level mapping was performed.\nWhile we expand our discussion of group-level analyses in response to other points, we indeed mapped all\nchannels to a single hemisphere to further improve spatial coverage of the cortex. Nonetheless, control\nanalyses in FIGURE S4 demonstrate that highly similar gradients are obtained when preserving the bilateral\ndistribution of channels (i.e., when no flipping was applied). These findings are also reported in the main\nmanuscript (p.6), suggesting that gradient findings were not significantly affected by this decision to map\nchannels to a single hemisphere:\n“We also replicated our findings when preserving the bilateral distribution of channels and applying a lower\nresolution parcellation (Schaefer-100) to maintain broad cortical coverage (FIGURE S4). The high correspondence of\nG1 (vertexwise Spearman correlation r =0.767; r =0.804) and G2 (vertexwise Spearman correlation r =0.645;\nleft right left\nr =0.671) topographies with gradients derived from our main analysis provides additional support for the\nright\nindependence of our findings to parcellation choice.”\nWe also edited the DISCUSSION to encourage explorations of asymmetries in iEEG signals across\nhemispheres in future work (p.15):\n“Future studies based on even larger datasets or providing a high proportion of bilaterally implanted cases may\nfurthermore help to elucidate principles of asymmetry in neural dynamic across the two hemispheres.”\n2.4 Structural connectivity analysis was also performed at the group level, and additionally, the\nconnection weights in the group level matrix were log-transformed to reduce variance in connectivity\nstrength. However, presumably, the differences in connectivity strength are important! Again, these are\nrelated to individual differences. Functional connectivity for rsfMRI data was also group averaged,\nignoring patient-specific differences.\nRegarding structural connectivity (SC), we followed established procedures for group-level analyses of this\nmodality as described in previous work (Betzel et al., 2019a). Specifically, distance-dependent consensus\nthresholding leverages edge length data to preserve topological parameters of subject-level networks by\nreliably representing both short- and long-range connections at the group-level. Resulting SC weights are\nhowever highly skewed, with few very strong connections and many weak or zero connections (as is also\nthe case at the individual-level). This property of SC matrix weights violates assumptions of normality but\ncan notably be mitigated by log-transforming edge values (Hansen et al., 2023; Sebenius et al., 2023; van\nden Heuvel et al., 2015). Given that statistical methods used here assume normality of the underlying data\n(i.e., Pearson correlation, multilinear regression), this log-transformation is an adequate transformation to\noptimize analytical validity (Nelson et al., 2023). We have added this explanation to the revised manuscript\n(p.19):\n“Group-level structural connectivity (SC) was computed using distance-dependent thresholding to preserve the\ndistribution of within- and between-hemisphere connection lengths found in individual subjects (140). However, SC\nweights of the resulting matrix remain highly skewed, with few very strong connections and many weak or zero\nconnections. As this property of the matrix violates assumptions of normality, which is problematic for the statistical\ntests applied in the present work, connection weights in the resulting group-level were log-transformed to reduce\nconnectivity strength variance.”\nWe agree with the Reviewer that individual differences in SC, FC, and other modalities leveraged in our\nwork are important topics of investigation. However, we respectfully disagree that individual-focused\napproaches should necessarily be conducted at the expense of discovering general principles of brain\norganization expressed at the group-level. While we argue in favour of the importance and relevance of\ngroup-level inference throughout this response (notably in points 2.6, 2.7, and 2.8), several steps in our\nprocessing pipeline were taken to sample individual-specific data as precisely as possible using state-of\nthe-art techniques. For example, surface-based features such as the timeseries used to generate the FC\nmatrices and the intensity profiles used to generate MPC matrices are sampled in the participant’s native\nspace for each modality, limiting data interpolation with co-registration to a standard space. Geodesic\ndistance measurements also relied on high-resolution, native space surfaces. These procedures are\nimplemented in addition to the SC consensus thresholding approach described above. This information is\nnow more clearly highlighted in the revised manuscript (p.18):\n“Imaging data were processed with micapipe (v0.2.3; micapipe.readthedocs.io/), an openly available pipeline for\nmultimodal MRI processing, data fusion, and multiscale connectome generation (119). This pipeline is notably\nadapted to sample features in modality-specific spaces, therefore limiting data blurring that can occur with co-\nregistration to a standard template space.”\n2.5 At the microscale, it seems that patient-specific matrices were obtained. I don’t understand how this\nanalysis can be considered comparable to the others that were performed at the group level.\nWe believe the Reviewer is referring to the microstructural profile covariance (MPC) matrices, which were\ncomputed at the patient-level and subsequently parcellated and averaged across individuals for all analyses\ninvolving the iEEG data. This is the case for all MRI-derived modalities, which we clarify in the revised\nmanuscript (p.6):\n“Matrices were computed for each participant and averaged across individuals (for details on data processing,\nanalysis, and group-level aggregation, see MATERIALS AND METHODS).”\n2.6 The data used in this study all comes from patients with epilepsy. As each patient has a different focus\nof their seizures and potential lesions in different locations, even when the channels are selected from\n“healthy” regions, it is to be expected that the underlying brain connectivity and resulting dynamics are\naffected. This should lead to increased individual variability within each cohort, making it even less\nreasonable to perform a group level analysis.\nThe Reviewer is correct that patients referred for iEEG often show a complex clinical presentation (e.g.,\ninconclusive or contradictory localization findings using non-invasive methods). This may increase some\ndegrees of variation across participants due to different seizure loci locations and clinical etiologies.\nHowever, the stringent data selection procedure of the MNI Open iEEG Atlas and Replication cohort\nremoving channels involved in epileptogenic processes (i.e., excluding channels located in abnormal tissue;\nselection of segments days after implantation; selection of segments several hours away from a seizure; no\nsegments selected after stimulation; absence of interictal discharges or slow-wave anomalies; consistent\nrecording conditions) helps mitigate these concerns. Furthermore, the replication analyses presented in our\npaper indicate good generalizability of our findings across two distinct cohorts, despite variability in clinical\npresentations in these samples.\nWe respectfully disagree with the Reviewer that this limitation of iEEG studies pertaining to inter-patient\nvariability make group-level analysis inappropriate for the objective of our study, that is the discovery of\ngeneral principles of brain organization. Historically, iEEG recordings performed in patients with epilepsy\nhave led to significant contributions to our understanding and mapping of brain function. For example,\nintracranial stimulation studies have produced representative maps of the cortical organization of function,\nnamely the motor and sensory homunculi (Penfield & Boldrey, 1937), and there have been further exciting\ndevelopments in our understanding of brain network organization based on recent applications of\nstimulation techniques (Momi et al., 2025). Aggregation of functional connectivity maps in predefined\nfrequency bands across patients implanted with ECoG can replicate network topographies found using rs-\nfMRI in healthy cohorts (Betzel et al., 2019b). Similarly, iEEG recordings allow the measurement of neural\ntimescales from power spectrum density characteristics (Gao et al., 2020; Honey et al., 2012), with longer\ntimescales in association cortex relative to sensory and motor regions replicating fMRI findings in healthy\nindividuals (Golesorkhi et al., 2021). Moreover, recordings in putatively normal tissue are increasingly used\nto build normative maps of cortical dynamics and can accurately identify abnormal brain dynamics in\nunseen data on a region-specific basis (Bernabei et al., 2022; Taylor et al., 2022). This suggest that datasets\nsuch as the MNI Open iEEG Atlas provide a valid normative repertoire of inter-individual variability in\ncortical dynamics amenable to investigations of group-level characteristics of these maps. Patients with\ndifferent epileptogenic zone locations can also undergo the same task to clarify the topography and\ndynamics associated with different experimental conditions. For example, signals originating from deep\nstructures such as the insula (Das & Menon, 2024; Huang et al., 2025), medial temporal lobe (Solomon et\nal., 2019), and amygdala (Wang et al., 2017) can be precisely captured with invasive techniques such as\niEEG. For significant group-level effects to be detected, individual-level effects must present a reasonable\nlevel of consistency, despite the clinical heterogeneity of the study sample. Similarly, many features of\ngeneral brain organization can be replicated in patients: in epilepsy, we can find similar sets of canonical\nnetwork communities and macroscale gradient patterns (Cabalo et al., 2024; Maneshi et al., 2014; Royer et\nal., 2023; Xie et al., 2024), and the spread of abnormal dynamics can be accurately modelled using\nconnectome information derived from healthy samples (Xia et al., 2025) . While quantitative differences\nhave been observed in patients vs. controls on metrics derived from these modalities (e.g., connectivity\nstrength, balance of segregation/integration of networks (Coito et al., 2016; Larivière et al., 2020a)), the\npreservation of these global principles supports the utility of patient studies to improve our understanding\nof human brain dynamics despite the presence of patient-specific disease-related alterations.\nNonetheless, even with careful selection of data segments, it remains an open question to what extent focal\nepileptic activity in one brain region might affect activity in the rest the brain. We now acknowledge this in\nthe revised DISCUSSION (p.15), in the segment cited in point 2.2.\n2.7 I assume that the reason that the analysis was performed on two cohorts was to alleviate some of\nthese concerns about group level analysis. However, I still believe that while the “average” analysis\nmight be similar in 2 cohorts, the actual findings could simply not hold in either. In any sort of population\nthat has sub-groups within the population or even just simply lots of variability across the population,\none can detect relationships between the average geometry and structure with average dynamics that\nsimply do not exist for any of the individuals in the population.\nThe Reviewer is correct that the analysis performed in the Replication cohort (Figure 4) served to assess\nthe robustness of findings reported in the main cohort to a distinct sample of participants. We now report\nthis clearly in the RESULTS (p.11):\n“We lastly leveraged an independent patient cohort to assess the replicability of macroscale gradients of neural\ndynamics and their association with within-sample neuroimaging metrics.”\nWhile group-level findings can be difficult to generalize to individual-level inference, the objective of the\npresent work was rather to identify general principles of brain organization from neurophysiological data\nand their network-level correlates. For example, the sensory-association axis is known to be repeatedly\nexpressed across structural, functional, and molecular modalities (Sydnor et al., 2021), yet its exact\ntopography varies throughout development and aging (Bethlehem et al., 2020; Larivière et al., 2020b), as\nwell as across disease states (Hong et al., 2019). We do not deny that there is likely inter-individual\nheterogeneity in the expression of these principles and their relationships across modalities, and that aspects\nof these dimensions and their coupling likely relate to symptoms and behaviour (Baum et al., 2020; Paquola\net al., 2019; Royer et al., 2023). These questions, however, are beyond the scope of the present study.\nNonetheless, our work represents a first step in identifying these overarching dimensions of brain\norganization and sets the foundation for future work implementing these techniques at an individual level.\nTo address the Reviewer’s concerns regarding individual-level replicability, we designed two additional\nanalyses (see also point 2.8). We first sought to assess if the group-level embedding space could capture\nindividual-level similarity in PSD architectures. To do so, we regenerated the group-level gradient space\nwhile iteratively excluding each subject, and calculated embedding distances in each group-level space. For\nthe left-out subject at each iteration, we evaluated if the proximity of sampled nodes in the embedding space\nwas related to the similarity of corresponding channel PSDs covered by the gradient parcels. PSD similarity\nwas assessed using a partial correlation of channel PSDs, while controlling for the within-subject average.\nWe specifically compared partial correlation values at short embedding distances (embedding distance <20th\npercentile) versus long embedding distances (embedding distance >80th percentile) for each patient. In all\npatients, channel PSDs were more highly correlated across nodes situated closer in the embedding space\ncomputed from the remaining 17 patients, although the strength of this effect varied across individuals\n(FIGURE S12).\nWe report this analysis in the revised RESULTS (p.11-12) and include this additional figure below:\n“Furthermore, reconstructing the embedding space while iteratively excluding each patient found that embedding\ndistances between nodes related to channel-specific PSD similarity in the held-out patient (FIGURE S12).”\ngradient space was regenerated while iteratively excluding each subject, and we calculated embedding\ndistances in each group-level space. For the left-out subject at each iteration, we evaluated if the proximity\nof sampled nodes in the embedding space was related to the similarity of corresponding channel PSDs\ncovered by the gradient parcels. PSD similarity was assessed using a partial correlation of channel PSDs,\nwhile controlling for the within-subject average. We specifically compared partial correlation values at\nshort embedding distances (embedding distance <20th percentile) versus long embedding distances\n(embedding distance >80th percentile) for each patient. In all patients (each represented by their own\nboxplot), PSDs were more highly correlated across nodes situated closer in the embedding space computed\nfrom the remaining 17 patients, although the strength of this effect varied across individuals (mean T-\nstatistic=16.004; S.D.=13.152). All p-values are FDR-corrected.\n2.8 While I understand the desire to fully cover the brain with the iEEG data, I think that to do this study\nproperly, one would have to perform the gradient analysis separately for each individual, and accept that\nbetween individuals, the spatial coverage would be different. I realize that this makes it hard to collapse\nthe analysis across individuals and perform regional analysis, but at least one could look at the\ncorrelations between gradients and structure at the individual level and later try to find where there was\nspatial overlap between patients to do a limited spatial analysis.\nThe Reviewer’s point is well taken that individual-level investigations would strengthen presented group-\nlevel findings. For this reason, we provide a new supplementary analysis and figure in which gradients of\nneural dynamics were generated for each patient independently, over covered regions. To improve coverage,\nwe also mapped all contacts to a single hemisphere as in our group analysis (although only 4/18 patients\nhad bilateral implantations). We computed embedding distance matrices for each patient and calculated\nedgewise correlations with MRI measures derived from each patient’s own data. The distribution of\ncorrelation coefficients is presented in FIGURE S11:\nFigure S11. Individual-level correspondence of iEEG-MRI associations in Replication cohort. We\ncomputed patient-specific gradients from covered regions and computed the embedding distance between\navailable nodes for gradients explaining up to 50% in the input data. Bilateral implantations were mapped\nto the right hemisphere (n=4), left-sided implantations were flipped to the right hemisphere (n=5), and\nright hemisphere implantation were analysed without flipping (n=9). For each patient, embedding\ndistances were correlated with MRI features of corresponding nodes, and resulting correlation coefficients\n(Pearson r) are represented as black markers. Each black dot represents a single patient, with the group\ndistribution shown in ridge plots for each modality (black line: median). Replicating group-level findings\nin the Discovery cohort, similarity of neural dynamics in the cortex were most strongly correlated with\ngeodesic distance (GD; median r=-0.325), followed by structural connectivity (SC; median r = -0.249),\nand lower associations with functional connectivity (median r = -0.151) and microstructural profile\ncovariance (median r = -0.156).\nThese findings are also reported in the revised RESULTS (p.11) and DISCUSSION (p.14-15):\n“Additional replication analyses implemented at the individual level also found stronger associations with GD relative\nto other modalities, followed by SC and balanced correlations with FC and MPC (FIGURE S11).”\n“Indeed, translating our approach at an individual-level revealed considerable inter-individual variability in the\nmagnitude of multimodal association. However, given the technical limitations of iEEG, it is difficult to disentangle\nthe relative contributions of individual differences in brain organization versus spatial sampling effects in explaining\nthis variability.”\n15\n2.9 In summary, because this work is based on group averaged data, making it hard to know where true\ncorrelations might exist at the individual level, I do not feel that it has the ability to be of high impact to\nthe field.\nWe are confident that our answers and new individual-level analyses in the Replication cohort (point 2.8)\nhave helped address the Reviewer’s concerns. We agree with the Reviewer that there is certainly high value\nin identifying sources of individual uniqueness. However, and as outlined in our response, numerous\nimportant findings in neuroscience focus on group-average data to identify general principles of human\nbrain organization, an endeavour to which to which we contribute through the present work.\n2.10 (Remarks on code availability): Because I have issues with the fundamental assumptions about the\ndata set, I did not review any code.\nThe code remains openly available for independent review.\nREVIEWER #3:\nSigned Review (Stuart Oldham)\n3.1 In this work, Royer et al. investigate how regional dynamics, as measured by iEEG, are constrained\nby geometric, microstructural, and connectomic features. They find that geodesic distances strongly\nconstrain brain dynamics, but this relationship is weakest in transmodal cortices, consistent with other\nfindings showing that neuronal organisation in these regions clearly departs from that in the rest of the\nbrain. The authors further strengthen the reliability of their findings by using a replication dataset,\nwhere they observe similar results. This is a very well-written paper, with excellent figures and good\nscience! 😊 I have no major concerns, and my comments are largely restricted to suggestions on how the\nauthors could improve the manuscript and enhance clarity where needed.\nWe thank Dr. Oldham for the positive assessment and thoughtful suggestions to improve our paper.\n3.2 As the authors will be oh so very well aware, there are many papers now correlating gradients with\nall the other gradients and all the other brain related things, so it is easy to get lost in what has been done\nbefore and where something clearly novel is actually being done. To the best of my knowledge, this is the\nfirst study using iEEG, correct? If so, I think this is worth highlighting more explicitly as a key novelty.\nI suggest making it very clear in the Introduction how this technique provides insights that previous\nmethods could not. At present, the importance of using iEEG is somewhat underemphasised.\nTo our knowledge, our paper is indeed the first to derive gradients from iEEG data. In accordance with the\nReviewer’s suggestion, we now further emphasize the novelty of our approach in INTRODUCTION (p.3-4):\n“The present work established a novel approach to fill this gap by uncovering large-scale gradients of neural dynamics\nand relating this complex activity landscape to underlying multiscale anatomy.”\n“These gradients, mapped for the first time in iEEG data in the present work, captured the coordinate system of neural\nrhythms in the neocortex, grouping regions with similar neurophysiology independently of their anatomical location.”\n3.3 All iEEG data used in this study were obtained from patients. Given that these individuals have\nepilepsy, and would be expected to exhibit altered patterns of brain activity, how applicable are the\nfindings to non-clinical populations? Or is it assumed that, outside of seizure episodes, patients’ brain\nactivity patterns are comparable to those of non-epileptic individuals?\nWhile iEEG presents unique advantages in terms of spatiotemporal precision, data collection with this\ntechnique is indeed limited to clinical populations, and particularly people with epilepsy. Without direct\naccess to equivalent data in healthy populations, it is difficult to ascertain how generalizable the present\nfindings are to non-clinical cohorts. However, different sets of findings allow us to speculate that results\nmay at least partially replicate in healthy samples. First, using resources such as the MNI Open iEEG atlas\n(as in the present work) as a normative dataset, previous studies could accurately identify abnormal brain\ndynamics in unseen data (Bernabei et al., 2022; Taylor et al., 2022). This suggests relatively homogeneous\nregional properties at the group-level, perhaps reflecting typical neurophysiological processes allowing the\ndetection of regional divergence from a normative distribution of neural dynamics. Second, previous work\napplying dimensionality reduction to a large number of timeseries features derived from\nmagnetoencephalography data collected in healthy participants showed a principal component notably\nseparating unimodal and transmodal cortical systems (Shafiei et al., 2023), similar to the present study.\nWe revised the DISCUSSION to include this information (p.15):\n“Nonetheless, investigations of human brain organization using iEEG are limited by the technique’s unique challenges,\nincluding its restricted accessibility, exclusive use in clinical settings, as well as sparse and patient-specific sampling\n(113). It also remains an open question to what extent focal epileptic activity in one brain region might affect activity\nin the rest the brain and skew the generalizability of iEEG findings to healthy brain function. Recent data sharing\ninitiatives, such as the resource used here (18) and others (8, 114-117), have been instrumental in addressing these\nlimitations. Although iEEG data is collected in a clinical context, pooling carefully curated single-patient data\nsignificantly boosts spatial coverage, thus enabling novel explorations of brain cartography and facilitating the\nidentification of overarching principles of neural organization. These innovations add on to the recognized advantages\nof iEEG regarding its high spatiotemporal precision to study human brain function, as well as its specificity to neural\nsignals over methods offering full brain coverage such as fMRI (113). Studies focusing on group-level inference, such\nas the present work, thus offer new perspectives for human brain mapping to which other modalities may be blind.\nGroup-level paradigms are a valuable first step towards more targeted individual-level inquiries able to clarify\npopulational variations. For example, by testing the consistency of these principles across individuals, extreme\ndeviations from the population norm can be identified and related to clinical presentations (3). Future studies based\non even larger datasets or providing a high proportion of bilaterally implanted cases may furthermore help to\nelucidate principles of asymmetry in neural dynamic across the two hemispheres.”\n3.4 In the main analysis using the Discovery Cohort, iEEG data from patients were combined with\nSC/FC/MPC measures calculated in healthy individuals. However, for the replication cohort, were both\nthe iEEG and MRI-derived measures calculated in the same patient cohort (i.e., was SC/FC/MPC\ncalculated for the MNI clinical cohort and used in the analysis)? I believe this is what was done, but if\nso, could it be made more explicit? More broadly, the main concern here is that data from healthy and\nclinical populations are being combined and used to make general inferences about brain organisation.\nThis issue is briefly acknowledged on lines 445–446 (“as well as epilepsy-related changes in SC that\ncould only be detected in a within-patient design”), but I think a clearer discussion or acknowledgement\nof this limitation should be made.\nThe main analyses performed in the Discovery cohort indeed focus on cross-dataset comparisons, while the\nReplication cohort uses within-sample data. As requested, we clarified this in the revised RESULTS (p.4):\n“For the main analyses, we leveraged a Discovery cohort (MNI open iEEG atlas, n=106) in which iEEG data were\naggregated in stereotaxic space (18). We assessed correspondence to neural organization derived from two high-\ndefinition MRI datasets acquired in 50 healthy individuals at 3T and 10 healthy individuals at 7T. Our main analyses\nthus rest on a design comparing iEEG and imaging modalities across different datasets. We then evaluate the\ngeneralizability of our findings using a fully independent Replication cohort (MNI clinical dataset, n=18), where both\niEEG and multimodal pre-surgical MRI data were available in each participant for within-cohort comparison of\nmodalities.”\nWe also expanded our DISCUSSION to acknowledge potential limitations of this approach (p.14):\n“Cross-referencing iEEG data collected in patients with MRI measures acquired in healthy participants allowed us\nto maximize spatial coverage and sample sizes in our Discovery cohort analyses. This approach is in line with the\nwidespread use of data contextualization techniques in neuroscience, which measures and integrates distinct maps of\nhuman brain organization (107). However, this remains a limitation of the present work, as it is unclear to what extent\nfindings from these distinct cohorts generalize to overall principles of brain organization.”\n3.5 All analyses are conducted in a single hemisphere (due to understandable technical constraints).\nEspecially when we consider that the pattern of inter-hemispheric connections is very different to intra-\nhemispheric connections and virtually all such inter-hemispheric connections are long-range, which as\nthe authors note seem to have some special role in shaping communication across brain networks (thus\ninter-hemispheric dynamics might be very different!). At a minimum, would be worth noting as a\nlimitation.\nWe explore this suggestion by using data mapped to the Schaefer-100 parcellation. We use this parcellation\nin our supplementary analysis assessing the replicability of our findings in a lower resolution parcellation\npreserving the bilateral distribution of channels in the MNI Open iEEG atlas (FIGURE S4). When computing\nthe embedding distances on this bilaterally mapped data and correlating edge values with corresponding\nSC, FC, and MPC entries, we find a similar sequence of correlation strengths across modalities: SC r=-\n0.322; FC r=-0.217; MPC r=-0.185 (as in the main analyses where SC |r| > FC |r| > MPC |r|). As GD is only\ncomputed within hemisphere, we could not include it in this analysis.\nWe report these findings in the revised RESULTS (p.7):\n“Findings were replicable […] when computing correlations between embedding distances and MRI measures from\nbilaterally mapped data (SC r = -0.322; FC r = -0.217; MPC r = -0.185).”\n3.6. I found it interesting that the first two gradients explain relatively little variance, whereas most other\ndecompositions of brain features based on regional similarity or covariance explain much more. Do the\nauthors have any intuition as to why this might be?\nThe variance explained by the different gradients is indeed quite balanced here. This observation notably\nguided our choice of using all gradients explaining up to 50% of the variance in the input data, rather than\nfocusing on the first one to three gradients as is often done in the literature.\nDetermining the unique contributions of factors with the potential to shift each gradient’s explained\nvariance (e.g., modality and its specific topographic properties, pre-processing, hyperparameters used to\ngenerate the gradients) would require careful analysis and comparisons across a multitude of conditions.\nHowever, we note that in the existing literature, the scaled eigenvalue as a measure of variance/information\nexplained of each gradient is usually computed on a smaller number of components. Previous work has\noften restricted the number of computed gradients to a much lower number relative to the dimensionality\nof their input data (parcellated with a higher resolution parcellation or using vertexwise data). In the present\ncase, we computed a relatively high number of gradients (99; restricted by the dimensionality of our\nparcellated input data). For example, scaling the eigenvalues of the iEEG gradients according to the first\n10 components increases the scaled eigenvalues of G1-G3 to 21.3%, 16.9%, and 14.2%. As a result, we\nbelieve that the interaction between the dimensionality of the input data and the number of gradients\nretained in the analysis affected the reported variance explained.\n3.7 I assume that the 50% cutoff for selecting gradients was chosen because it's a nice big round number.\nDo the results change meaningfully if only the first 3–5 gradients are used? I realise re-running all\nanalyses would be burdensome (though it could serve as a nice reproducibility check!). However, if you\ncan show that the correlation between the full embedding matrix (computed with G1–G14) and a reduced\nmatrix (e.g., computed with just G1–G2, G1–G3, etc.) remains high, that would suggest the results are\nrobust to the number of gradients included.\nWe found that the embedding distance matrix was highly correlated across many different combinations of\ngradients from 2 to 14. The correlation matrix is shown below and has been added to our Supplementary\nmaterials. We refer to this analysis in our revised RESULTS (p.6):\n“This matrix was highly robust to the number of gradients included when calculating inter-nodes distances (FIGURE\nS8)”\nFigure S8. Robustness to variations in number of gradients included in embedding distance matrix. We\nsystematically varied the number of gradients included when computing inter-node embedding distances\nand cross-correlated each matrix (excluding the diagonal). The resulting correlation matrix (min r = 0.902;\nmean r = 0.965) shows high reproducibility of the embedding distance matrix when successively including\nmore gradients to compute embedding distances.\n3.8 In Figure 2A, the matrix is labelled as indicating “Euclidean distance,” but in the following figure\nit is described as “embedding distance.” While the embedding distance is indeed a Euclidean distance,\nfor consistency I suggest updating the colourbar label to “Embedding distance.” Also (as I am sure you\nare aware) “Euclidean distance” is typically understood to refer to physical distances between parcels,\nso to some people (i.e., me) labelling a matrix of the distances in high dimensional gradient space as\n“Euclidean distance” may be confusing.\nWe thank Dr. Oldham for this helpful suggestion and updated the labelling in FIGURE 2A accordingly.\n3.9 Also for Figure 2A: in the brain projection, could the seed parcel be made more visually prominent\nor central? The selected hippocampal parcel is somewhat obscured, and I missed it at first glance.\nWe updated the figure following the Reviewer’s suggestion.\n3.10 Is the geodesic distance measured in millimetres? The colour scale suggests a maximum distance\nof 550 mm, which seems quite high. In my own data, the maximum geodesic distance across individuals\ntends to be around 250–300 mm. Could you please double-check?\nWe use the Workbench Command function -surface-geodesic-distance to compute individual\ngeodesic distance matrices, which we believe to output distance in millimeters (although this does not seem\nto be specified in the documentation). We speculate that the reason output distances are higher here might\nbe related to our use of high-resolution, native surfaces to compute the distances (coastline paradox).\nIn detail, we used micapipe (https://github.com/MICA-MNI/micapipe/) to align selected parcellations to the\nnative surface space, select a centroid vertex for each parcel, compute the geodesic distance of this centroid\nvertex to all other surface vertices, and average distances within each parcel. We added this clarification to\nthe METHODS (p.18-19):\n“Of note, longer distances measurements may be obtained here due to our use of high-resolution meshes, compared\nto techniques using surfaces resampled according to a lower-resolution template’s geometry.”\n3.11 For Figure 3B, it wasn’t immediately clear that the colours of the bars correspond to the colours of\nthe matrices for each feature. I think there’s enough space in the legend to explicitly state what each\ncolour represents.\nThis information was added to the revised legend (p.10):\n“The colour of the bars corresponds to the different MRI modalities shown in panel A (blue: GD; orange: SC; pink:\nFC; green: MPC).”\n3.12 Also for Figure 3B: are the ridge plots showing the distribution of bootstrap results? This isn’t\nstated clearly. Is the vertical line the mean or median of the bootstrap distribution, and is the dot the non-\nbootstrapped result?\nThis is indeed correct, and the information was added to the revised legend (p.10):\n“Each ridge plot shows the distribution of bootstrap results, with their median value indicated by the black line. The\nblack dot marker indicates the true adjusted r2 for each model.”\n3.13 For Figure 3C: in the far-right plot showing correlations with FC G1, is the variable shown the\ncorrelation between nodal FC G1 values and the full model–baseline difference from the previous plot?\nIf so, could you ensure the name of this variable is consistent across figures to help with interpretation?\nThis specific correlation (full model – baseline model difference) corresponds to the rightmost marker in\nthis plot (dark red). The other markers show the correlations values for each combination of partial model\n– baseline model difference. We have changed the label of the y-axis to better reflect this.\n3.14 Figure 4A: what does the rainbow line overlaid on the brain slices represent?\nWe added this information to the revised caption of FIGURE 4 (p.12):\n“The coloured markers overlaid on the T1-weighted images indicate the locations of contacts of a single electrode\nimplanted in this patient (deeper contacts in blue, superficial contacts in orange).”\n3.15 The paragraphs in the Discussion section are quite long. Breaking them up would improve\nreadability.\nWe followed the Reviewer’s suggestion and have broken up the paragraphs of the DISCUSSION.\n3.16 (Remarks on code availability): I quickly checked the code, seems to all run and reproduce the\nresults in the paper.\nWe thank the Reviewer for verifying the replication scripts provided in the repository.\nBaum, G. L., Cui, Z., Roalf, D. R., Ciric, R., Betzel, R. F., Larsen, B., et al. (2020). Development of\nstructure–function coupling in human brain networks during youth. Proceedings of the National\nAcademy of Sciences, 117(1), 771-778.\nBernabei, J. M., Sinha, N., Arnold, T. C., Conrad, E., Ong, I., Pattnaik, A. R., et al. (2022). Normative\nintracranial EEG maps epileptogenic tissues in focal epilepsy. Brain, 145(6), 1949-1961.\nBethlehem, R. A. I., Paquola, C., Seidlitz, J., Ronan, L., Bernhardt, B., Consortium, C.-C. A. N., et al.\n(2020). Dispersion of functional gradients across the adult lifespan. Neuroimage, 222, 117299.\nBetzel, R. F., Griffa, A., Hagmann, P., & Mišić, B. (2019a). Distance-dependent consensus thresholds for\ngenerating group-representative structural brain networks. Network neuroscience, 3(2), 475-496.\nBetzel, R. F., Medaglia, J. D., Kahn, A. E., Soffer, J., Schonhaut, D. R., & Bassett, D. S. (2019b). Structural,\ngeometric and genetic factors predict interregional brain connectivity patterns probed by\nelectrocorticography. Nature Biomedical Engineering, 3(11), 902-916.\nCabalo, D. G., Leppert, I. R., Thevakumaran, R., Dekraker, J., Hwang, Y., Royer, J., et al. (2024).\nMultimodal Precision Neuroimaging of the Individual Human Brain at Ultra-high fields. bioRxiv,\n2024.2006. 2017.596303.\nCoito, A., Genetti, M., Pittau, F., Iannotti, G. R., Thomschewski, A., Höller, Y., et al. (2016). Altered\ndirected functional connectivity in temporal lobe epilepsy in the absence of interictal spikes: a high\ndensity EEG study. Epilepsia, 57(3), 402-411.\nDas, A., & Menon, V. (2024). Electrophysiological dynamics of salience, default mode, and frontoparietal\nnetworks during episodic memory formation and recall: A multi-experiment iEEG replication. In:\neLife Sciences Publications, Ltd.\nGao, R., van den Brink, R. L., Pfeffer, T., & Voytek, B. (2020). Neuronal timescales are functionally\ndynamic and shaped by cortical microarchitecture. eLife, 9, e61277.\nGolesorkhi, M., Gomez-Pilar, J., Tumati, S., Fraser, M., & Northoff, G. (2021). Temporal hierarchy of\nintrinsic neural timescales converges with spatial core-periphery organization. Communications\nbiology, 4(1), 277.\nHansen, J. Y., Shafiei, G., Voigt, K., Liang, E. X., Cox, S. M., Leyton, M., et al. (2023). Integrating\nmultimodal and multiscale connectivity blueprints of the human cerebral cortex in health and\ndisease. PLoS biology, 21(9), e3002314.\nHoney, C. J., Thesen, T., Donner, T. H., Silbert, L. J., Carlson, C. E., Devinsky, O., et al. (2012). Slow\ncortical dynamics and the accumulation of information over long timescales. Neuron, 76(2), 423-\n434.\nHong, S.-J., De Wael, R. V., Bethlehem, R. A., Lariviere, S., Paquola, C., Valk, S. L., et al. (2019). Atypical\nfunctional connectome hierarchy in autism. Nature communications, 10(1), 1-13.\nHuang, W., Lyu, D., Stieger, J. R., Gotlib, I. H., Buch, V., Wagner, A. D., et al. (2025). Direct interactions\nbetween the human insula and hippocampus during memory encoding. Nature Neuroscience, 28(8),\n1763-1771.\nLarivière, S., Rodríguez-Cruces, R., Royer, J., Caligiuri, M. E., Gambardella, A., Concha, L., et al. (2020a).\nNetwork-based atrophy modeling in the common epilepsies: A worldwide ENIGMA study. Science\nAdvances, 6(47), eabc6457.\nLarivière, S., Vos de Wael, R., Hong, S.-J., Paquola, C., Tavakol, S., Lowe, A. J., et al. (2020b). Multiscale\nstructure–function gradients in the neonatal connectome. Cerebral cortex, 30(1), 47-58.\nManeshi, M., Vahdat, S., Fahoum, F., Grova, C., & Gotman, J. (2014). Specific Resting-State Brain\nNetworks in Mesial Temporal Lobe Epilepsy [Original Research]. Frontiers in Neurology, Volume\n5 - 2014.\nMomi, D., Wang, Z., Parmigiani, S., Mikulan, E., Bastiaens, S. P., Oveisi, M. P., et al. (2025). Stimulation\nmapping and whole-brain modeling reveal gradients of excitability and recurrence in cortical\nnetworks. Nature Communications, 16(1), 3222.\nbrain connectome weighted by the myelin content and total intra-axonal cross-sectional area of\nwhite matter tracts. Network Neuroscience, 7(4), 1363-1388.\nPaquola, C., Vos De Wael, R., Wagstyl, K., Bethlehem, R. A., Hong, S.-J., Seidlitz, J., et al. (2019).\nMicrostructural and functional gradients are increasingly dissociated in transmodal cortices. PLoS\nbiology, 17(5), e3000284.\nPenfield, W., & Boldrey, E. (1937). Somatic motor and sensory representation in the cerebral cortex of man\nas studied by electrical stimulation. Brain, 60(4), 389-443.\nRoyer, J., Larivière, S., Rodriguez-Cruces, R., Cabalo, D. G., Tavakol, S., Auer, H., et al. (2023). Cortical\nmicrostructural gradients capture memory network reorganization in temporal lobe epilepsy. Brain,\n146(9), 3923-3937.\nSebenius, I., Seidlitz, J., Warrier, V., Bethlehem, R. A. I., Alexander-Bloch, A., Mallard, T. T., et al. (2023).\nRobust estimation of cortical similarity networks from brain MRI. Nature Neuroscience, 26(8),\n1461-1471.\nShafiei, G., Fulcher, B. D., Voytek, B., Satterthwaite, T. D., Baillet, S., & Misic, B. (2023).\nNeurophysiological signatures of cortical micro-architecture. Nature Communications, 14(1), 6000.\nSolomon, E. A., Stein, J. M., Das, S., Gorniak, R., Sperling, M. R., Worrell, G., et al. (2019). Dynamic\nTheta Networks in the Human Medial Temporal Lobe Support Episodic Memory. Current Biology,\n29(7), 1100-1111.e1104.\nSydnor, V. J., Larsen, B., Bassett, D. S., Alexander-Bloch, A., Fair, D. A., Liston, C., et al. (2021).\nNeurodevelopment of the association cortices: Patterns, mechanisms, and implications for\npsychopathology. Neuron, 109(18), 2820-2846.\nTaylor, P. N., Papasavvas, C. A., Owen, T. W., Schroeder, G. M., Hutchings, F. E., Chowdhury, F. A., et al.\n(2022). Normative brain mapping of interictal intracranial EEG to localize epileptogenic tissue.\nBrain, 145(3), 939-949.\nvan den Heuvel, M. P., de Reus, M. A., Feldman Barrett, L., Scholtens, L. H., Coopmans, F. M. T., Schmidt,\nR., et al. (2015). Comparison of diffusion tractography and tract-tracing measures of connectivity\nstrength in rhesus macaque connectome. Human Brain Mapping, 36(8), 3064-3075.\nWang, S., Yu, R., Tyszka, J. M., Zhen, S., Kovach, C., Sun, S., et al. (2017). The human amygdala\nparametrically encodes the intensity of specific facial emotions and their categorical ambiguity.\nNature Communications, 8(1), 14821.\nXia, J., Yang, S., Li, J., Meng, Y., Niu, J., Chen, H., et al. (2025). Normative structural connectome\nconstrains spreading transient brain activity in generalized epilepsy. BMC Medicine, 23(1), 258.\nXie, K., Royer, J., Larivière, S., Rodriguez-Cruces, R., Frässle, S., Cabalo, D. G., et al. (2024). Atypical\nconnectome topography and signal flow in temporal lobe epilepsy. Progress in Neurobiology, 236,\n102604.\n23"}, {"round": 2, "reviewer_comments": {"reviewer_1": "the authors have done an excellent revision and have fully addressed the two main issues in my initial review, complemented with satisfactory new results that strengthen their conclusions. specifically, they have provided (i) clearer justification of the methodological rationale and analysis steps, and (ii) additional evidence supporting the robustness of their findings across parcellation schemes. both concerns have been convincingly resolved, in part through a toning down of some claims. overall, this is a clear, careful, and thoughtful revision. (note: i came across few very minor typos in the quotes within response letter; please double-check in the revised manuscript: Figure S9, caption: “distribution of correlations coefficients” -> “distribution of correlation coefficients”; Section 1.3.4 (quoted passage): “(r=r=-0.323)” -> “(r=-0.323)”; Section 1.9 (quoted passage): “were averaged” -> “was averaged”.)\n(Remarks on code availability)\ni checked the code during the first review round, and no new commits have been added to the repository since then.", "reviewer_2": "While the authors and I somewhat disagree on the importance of group level analysis in this context, I appreciate that they have performed additional analysis and toned down their messaging by more directly addressing the limitations of the approach. This better allows readers to make their own judgement about how or if the results might generalize at the individual level.\n(Remarks on code availability)", "reviewer_3": "The authors have addressed all of my comments, and I have no further issues to raise. I am happy to recommend the paper for publication :)\nStuart Oldham\n(Remarks on code availability)\nI was able to run the code"}, "author_response": "RESPONSE TO REVIEWERS (NCOMMS-25-05787A)\nWe would like to thank the Editors for accepting our paper for publication in Nature Communications under\nan open access license. We are grateful to the Reviewers for their thoughtful comments and found the\nsuggestions very helpful to improve our paper. As Reviewers 2 and 3 had no further comments following\nthe previous round of revisions, we only address the last few suggestions of Reviewer 1 below.\nREVIEWER #1\n1.1 the authors have done an excellent revision and have fully addressed the two main issues in my initial\nreview, complemented with satisfactory new results that strengthen their conclusions. specifically, they\nhave provided (i) clearer justification of the methodological rationale and analysis steps, and (ii)\nadditional evidence supporting the robustness of their findings across parcellation schemes. both\nconcerns have been convincingly resolved, in part through a toning down of some claims. overall, this is\na clear, careful, and thoughtful revision. (note: i came across few very minor typos in the quotes within\nresponse letter; please double-check in the revised manuscript: Figure S9, caption: “distribution of\ncorrelations coefficients” -> “distribution of correlation coefficients”; Section 1.3.4 (quoted passage):\n“(r=r=-0.323)” -> “(r=-0.323)”; Section 1.9 (quoted passage): “were averaged” -> “was averaged”.)\nWe thank Dr. Behjat for their careful revision of the manuscript. We made all suggested corrections in the\nrevised paper and supplementary materials.\n1"}]

Decision letter