Geometric constraints on human brain function

Abstract

The anatomy of the brain necessarily constrains its function, but precisely how remains unclear. The classical and dominant paradigm in neuroscience is that neuronal dynamics are driven by interactions between discrete, functionally specialized cell populations connected by a complex array of axonal fibres^1,2,3. However, predictions from neural field theory, an established mathematical framework for modelling large-scale brain activity^4,5,6, suggest that the geometry of the brain may represent a more fundamental constraint on dynamics than complex interregional connectivity^7,8. Here, we confirm these theoretical predictions by analysing human magnetic resonance imaging data acquired under spontaneous and diverse task-evoked conditions. Specifically, we show that cortical and subcortical activity can be parsimoniously understood as resulting from excitations of fundamental, resonant modes of the brain’s geometry (that is, its shape) rather than from modes of complex interregional connectivity, as classically assumed. We then use these geometric modes to show that task-evoked activations across over 10,000 brain maps are not confined to focal areas, as widely believed, but instead excite brain-wide modes with wavelengths spanning over 60 mm. Finally, we confirm predictions that the close link between geometry and function is explained by a dominant role for wave-like activity, showing that wave dynamics can reproduce numerous canonical spatiotemporal properties of spontaneous and evoked recordings. Our findings challenge prevailing views and identify a previously underappreciated role of geometry in shaping function, as predicted by a unifying and physically principled model of brain-wide dynamics.

Main

The dynamics of many natural systems are fundamentally constrained by their underlying structure. For instance, the shape of a drum influences its acoustic properties, the morphology of a river bed shapes underwater currents and the geometry of a protein determines the molecules with which it can interact⁹. The nervous system is no exception, with the rich and complex spatiotemporal dynamics of anatomically distributed neuronal populations being supported by their intricate web of axonal interconnectivity^1,10. Several studies have shown correlations between various properties of brain connectivity and activity¹¹, but precisely how spatiotemporal patterns of neural dynamics are constrained by a relatively stable neuroanatomical scaffold remains unclear.

In diverse areas of physics and engineering, structural constraints on system dynamics can be understood via the system’s eigenmodes, which are fundamental spatial patterns corresponding to the natural, resonant modes of the system¹². In the linear regime, such as brain activity under normal (that is, non-seizure-like) conditions¹³, eigenmodes (hereafter also referred to as modes) offer a particularly powerful and rigorous formalism for linking brain anatomy with the physical processes that shape activity. Through this lens, spatiotemporal patterns of neuronal dynamics emerge from excitations of the brain’s structural eigenmodes, much like the harmonics of a plucked violin string arise from vibrations of its own resonant modes.

Critically, just as the resonant frequencies of a violin string are determined by its length, density and tension, the eigenmodes of the brain are determined by its structural—physical, geometric and anatomical—properties. Do any of these specific structural properties make a dominant contribution to dynamics? Here we test two influential and competing theories that make different predictions about which key elements of brain structure shape function.

One classical perspective, which represents the dominant paradigm in neuroscience, has its roots in Ramon y Cajal’s neuron doctrine¹⁴, Brodmann’s cytoarchitectonics¹⁵ and over a century of work localizing functions to specific brain regions¹⁶. According to this view, spatiotemporal patterns of neural dynamics arise from interactions between discrete, functionally specialized cell populations connected by a topologically complex array of short- and long-range axonal connections^2,3. In humans, these connections can be estimated at macroscopic scales by diffusion magnetic resonance imaging (dMRI) to yield a structural connectivity matrix, called a connectome¹⁷. This approach has been used extensively to understand brain organization and dynamics^2,17,18, and recent work has proposed that eigenmodes derived from such a discrete connectome—referred to here as connectome eigenmodes—can be used to reconstruct the spatial patterns of canonical functional networks of the human cortex mapped with functional MRI (fMRI)^19,20.

One limitation of this discrete connectome-based view is that it relies on an abstract representation of brain anatomy that does not directly account for its physical properties and spatial embedding (that is, geometry and topology). These characteristics are explicitly incorporated into a broad class of neural field theories (NFTs)^4,5,6 that describe mean-field neural dynamics on spatial scales above 0.5 mm (Supplementary Information 1). One physiologically constrained form of NFT has unified a diverse range of empirical phenomena^6,21 by treating cortical activity as a superposition of travelling waves propagating through a physically continuous sheet of neural tissue. In this theory, neural interactions between different cortical locations are approximated by a homogeneous spatial kernel that declines roughly exponentially with distance²², in line with experimental evidence that the organization of the nervous systems of numerous species is universally governed by an exponential distance rule (EDR) for connectivity^10,23,24.

Given wave-like dynamics and EDR-like connectivity, a key prediction of NFT is that the intrinsic geometry of the brain physically shapes and imposes boundary conditions on emergent dynamics^7,8,25. A remarkable corollary of this view is that, if we prioritize spatial and physical constraints on brain anatomy, we only need to consider the shape of the brain, and not its full array of topologically complex axonal interconnectivity, to understand spatially patterned activity. More formally, the theory predicts that eigenmodes derived from brain geometry—hereafter referred to as geometric eigenmodes—represent a more fundamental anatomical constraint on dynamics than the connectome^7,8,25. This view stands in stark contrast to the classical view that complex patterns of interregional anatomical connectivity shape brain activity²⁶.

Here we test these competing views of the brain with the aim of identifying the principal structural constraints on human brain dynamics. In line with theoretical predictions from NFT, we show that diverse experimental fMRI data from spontaneous and task-evoked recordings in the human neocortex can be explained more parsimoniously by eigenmodes derived from cortical geometry (geometric eigenmodes) than by those obtained from measures of brain connectivity (connectome eigenmodes). We further confirm that stimulus-evoked activity is dominated by excitations of geometric eigenmodes with long spatial wavelengths, challenging classical views that such activity is localized to focal, spatially isolated clusters. To directly link these structural constraints to the physical processes driving brain dynamics, we use a generative model to show how wave dynamics unfolding on the geometry of the cortex can explain diverse features of functional brain organization. Finally, we show that the close relationship between geometry and function captured by eigenmodes extends to non-neocortical structures, indicating that this link is a universal property of brain organization.

Geometric modes constrain cortical activity

We first examine the degree to which geometric eigenmodes can explain diverse aspects of human neocortical activity. To derive the eigenmodes, we use a mesh representation of a population-averaged template of the neocortical surface (Fig. 1a and Derivation of cortical geometric eigenmodes in Methods). We then construct the Laplace–Beltrami operator (LBO) from this surface mesh, which captures local vertex-to-vertex spatial relations and curvature, and solve the eigenvalue problem,

$${\nabla }^{2}\psi =\Delta \psi =-\,\lambda \psi ,$$

(1)

where ∇ is the gradient operator, Δ is the LBO and $\psi =\{{\psi }_{1}({\bf{r}}),{\psi }_{2}({\bf{r}}),\ldots \}$ is the family of geometric eigenmodes with the corresponding family of eigenvalues $\lambda =\{{\lambda }_{1},{\lambda }_{2},\ldots \}$. The eigenvalues are ordered sequentially according to the spatial frequency or wavelength of the spatial patterns of each mode (Fig. 1a and Extended Data Fig. 1), such that ψ₁ is the mode with the longest wavelength. The resulting eigenmodes are orthogonal, forming a complete basis set to decompose spatiotemporal dynamics unfolding on the cortex as a weighted sum of modes with varying wavelengths (Fig. 1b and Modal decomposition of brain activity in Methods). Unless otherwise specified, we use N = 200 modes throughout this study.

Fig. 1: Reconstruction of neocortical activity with geometric eigenmodes.

a, Geometric eigenmodes are derived from the cortical surface mesh by solving the eigenvalue problem, $\Delta \psi =-\,\lambda \psi$ (equation (1)). The modes ${\psi }_{1},\,{\psi }_{2},{\psi }_{3},\ldots ,\,{\psi }_{N}$ are ordered from low to high spatial frequency (long to short spatial wavelengths). Negative, zero and positive values are coloured blue, white and red, respectively. b, Modal decomposition of brain activity data. The example shows how a spatial map, $y({\bf{r}}\,,t)$, at a given time, t, can be decomposed as a sum of modes, ψ_j, weighted by a_j. c, Left, we reconstruct task-evoked data using spatial maps of activation for a diverse range of stimulus contrasts. Right, we reconstruct spontaneous activity by decomposing the spatial map at each time frame and generating a region-to-region FC matrix. d, Reconstruction accuracy of seven key HCP task-contrast maps (Supplementary Information 2.1) and resting-state FC as a function of the number of modes. Insets show cortical surface reconstructions, demonstrating the spatial scales relevant to the first 10, 100 and 200 modes corresponding to spatial wavelengths of approximately 120, 40 and 30 mm, respectively. e, Group-averaged empirical task-activation maps and reconstructions (recon.) obtained using 10, 100 and 200 modes of the seven key HCP task contrasts. Black arrowheads indicate localized activation patterns that are more accurately reconstructed when using short-wavelength modes. f, Group-averaged empirical resting-state FC matrices and reconstructions using 10, 100 and 200 modes.

Using this decomposition we evaluate the accuracy of geometric eigenmodes in capturing both task-evoked and spontaneous brain activity (Fig. 1c) measured in 255 healthy individuals from the Human Connectome Project²⁷ (HCP; HCP data in Methods and Supplementary Information 2). For task-evoked activity, we map 47 task-based contrasts drawn from seven different tasks representing distinct evoked activation patterns. We then reconstruct each individual’s activation map using an increasing number of modes up to a maximum of 200 (Fig. 1d). For spontaneous, task-free (so-called resting-state) activity, we reconstruct the spatial map of activity at each time frame and then generate a region-to-region functional coupling (FC) matrix, describing correlations of activity among 180 discrete brain regions per hemisphere²⁸. To allow direct comparison between task-evoked and spontaneous recordings, we apply the same regional parcellation to the task-evoked data (Cortical parcellations in Methods). Finally, we quantify reconstruction accuracy by calculating the correlation between empirical and reconstructed task-evoked activation maps and spontaneous FC matrices (Fig. 1d–f).

We observe that reconstruction accuracy increases with an increasing number of modes across all task contrasts and in the resting state, with r ≥ 0.38 already achieved using just N = 10 modes (Fig. 1d). Large-scale modes are also differentially recruited across different tasks, suggesting that particular stimuli excite specific modes (Fig. 1e). Improvements in reconstruction accuracy become slow after ten modes, reaching r ≥ 0.80 at approximately N = 100 modes, with only incremental increases in reconstruction accuracy beyond this point. Beause the first 100 modes have wavelengths above around 40 mm (Supplementary Table 1), and the inclusion of shorter-wavelength modes only refines reconstruction of localized patterns (arrowheads in Fig. 1e), our findings suggest that the data are predominantly comprised of spatial patterns with long spatial wavelengths (see next section for a more detailed analysis).

These results are consistent across all 47 HCP task contrasts (Supplementary Fig. 1) and parcellations of varying resolutions (Supplementary Fig. 2), but data parcellated at higher resolution require more modes to achieve high reconstruction accuracy due to the low-pass spatial filtering effect of coarser parcellations. Our results are also not affected by the use of a population-averaged cortical surface template (rather than individual-specific surfaces) to derive the geometric eigenmodes (Supplementary Figs. 3–5 and Supplementary Information 3). Together these findings indicate that cortical geometric eigenmodes form a compact representation that captures diverse aspects of task-evoked and spontaneous cortical activity. Moreover, they show that such activity is dominated by long-wavelength, large-scale eigenmodes.

We next test the hypothesis that geometric eigenmodes provide a more parsimonious and fundamental description of dynamics than eigenmodes derived from a graph-based connectome approximation. To this end we compare the reconstruction accuracy of geometric eigenmodes against three alternative connectome-derived eigenmode basis sets (see Fig. 2a for a schematic). The first basis set is derived empirically from a connectome mapped with dMRI tractography at vertex resolution and thresholded to obtain a connection density of 0.10%, as done previously²⁹ (Derivation of connectome eigenmodes in Methods). The second basis set is derived from a connectome constructed synthetically according to a homogeneous stochastic wiring process governed by an exponential distance-dependent connection probability to mimic simple, EDR-like connectivity (Derivation of EDR eigenmodes in Methods). Because the connection densities of empirical and EDR connectomes differed, we evaluated a third basis set derived from the empirical connectome thresholded at 1.55% to match the density of the EDR connectome. The connectome, EDR and density-matched connectome eigenmodes described above are derived from the graph Laplacian (a discrete counterpart of the LBO) of their respective connectivity matrices (Fig. 2b and Extended Data Fig. 1).

Fig. 2: Geometric eigenmodes benchmarked against connectome-based eigenmodes.

a, Schematic of the anatomical properties used to derive eigenmodes for cortical geometry, the connectome and the EDR connectome. Geometric eigenmodes rely on local surface mesh information such as links (blue) between neighbouring surface mesh vertices (dots) and curvature. Connectome eigenmodes rely on local links between mesh vertices (blue) and short- and long-range connections (magenta) reconstructed empirically from dMRI. EDR eigenmodes rely on connections (red) generated from a stochastic wiring process in which the probability of connection between vertices exponentially decays as a function of their distance. b, Example connectome and EDR eigenmodes. Negative, zero and positive values are coloured blue, white and red, respectively. Despite some similarities, the spatial patterns of the modes are distinct from those derived using cortical geometry (compare with Fig. 1a). c, Reconstruction accuracy of resting-state FC matrices achieved by geometric, EDR and two variants of connectome eigenmodes: one using a connectome as defined using previous methods²⁹ and the other with the same connection density as the EDR connectome to allow fair comparison (for other densities see Supplementary Figs. 6 and 7). d, Difference in reconstruction accuracy of all 47 HCP task-contrast maps achieved by geometric eigenmodes and the other basis sets, as indicated by the text above each panel. Each row represents a different task contrast, grouped here by broad types (Supplementary Information 2.1); red indicates superior performance for geometric eigenmodes. Note that while there seems to be a performance advantage for connectome eigenmodes for reconstructions incorporating fewer than ten modes relative to geometric eigenmodes, reconstruction accuracy is generally low (average r = 0.42 across the different tasks) compared with that for 100 modes (average r = 0.71).

To summarize, geometric eigenmodes account for the intrinsic curvature of the cortical surface and local vertex-to-vertex relations in the surface mesh; connectome eigenmodes do not consider curvature but capture local spatial relations between mesh vertices, along with short- and long-range connections measured with dMRI; and EDR eigenmodes account for the effect of a homogeneous, stochastic, distance-dependent connection rule without fully capturing the cortical geometry (Fig. 2a). Contrasting these different basis sets thus allows us to disentangle the contributions to brain dynamics of cortical geometry from structural connectivity.

Direct comparison of the reconstruction accuracy of these different basis sets shows that geometric eigenmodes consistently show the highest reconstruction accuracy across both spontaneous (Fig. 2c) and task-evoked (Fig. 2d) data. EDR eigenmodes perform nearly as well as geometric eigenmodes whereas connectome eigenmodes are the least accurate. This finding holds true regardless of the parcellation used (Extended Data Figs. 2 and 3), the specific connection density used to generate the connectome eigenmodes (Supplementary Figs. 6 and 7 and Supplementary Information 4) and whether we generate the connectome using a discrete regional parcellation rather than at vertex resolution (Supplementary Fig. 8 and Supplementary Information 4). We additionally find that geometric eigenmodes show stronger out-of-sample generalization than principal components of the functional data themselves (calculated via principal component analysis (PCA); Supplementary Fig. 9, Extended Data Fig. 4 and Supplementary Information 5) and better performance than Fourier spatial basis sets (Extended Data Fig. 5, Supplementary Information 6 and Comparisons with statistical basis sets in Methods).

Taken together, these results demonstrate the parsimony, robustness and generality of geometric eigenmodes as a basis set for brain function. They also support the prediction of NFT that brain activity is best represented in terms of eigenmodes derived directly from the shape of the cortex, thus emphasizing a fundamental role of geometry in constraining dynamics.

Long wavelengths dominate cortical activity

Reconstructions of both spontaneous and task-evoked data with geometric eigenmodes show that the spatial organization of brain activity is dominated by patterns with spatial wavelengths of about 40 mm or longer (Fig. 1d–f). This result counters classical approaches to analysis of neuroimaging data, in which stimulus-evoked activations are mapped by thresholding statistical maps to identify focal, isolated areas of heightened activity. This classical approach rests on the assumption that focal loci represent discrete brain regions putatively engaged by the stimulus and that subthreshold activity in other regions is of negligible interest. The surprisingly long-wavelength content of task-activation data (Fig. 1d–e) suggests that classical procedures focus only on the tips of the iceberg and obscure the underlying spatially extended and structured patterns of activity evoked by the task (see Extended Data Fig. 6 for an explanation of the reasons involved). These observations accord with the theoretical predictions of NFT and previous analyses of task-evoked electroencephalography (EEG) signals^30,31.

Here we leverage the modal decomposition described in Fig. 1b to characterize the complete spatial pattern—the entire iceberg—of task-evoked activation. To this end we analyse the spatial power spectrum obtained using a geometric mode decomposition of group-averaged unthresholded activation maps from the 47 task contrasts in HCP^27,32 (Modal power spectra of task-evoked activation maps in Methods). As an independent replication, we also analyse 10,000 unthresholded activation maps from 1,178 independent experiments available in the NeuroVault repository³³, thus providing a comprehensive picture of the diversity of stimulus-evoked activation patterns mapped in the human brain.

Despite the wide range of stimuli, paradigms and data-processing approaches used to acquire these activation maps, we observe that a large fraction of power in the maps is concentrated in the first 50 modes, corresponding to spatial wavelengths greater than around 60 mm (Fig. 3a; similar results are found separately for each of the key HCP task-contrast maps; Extended Data Fig. 7). Using surrogate data, we confirm that these findings cannot be explained by the spatial smoothing induced by typical fMRI processing pipelines, which can filter out short-wavelength spatial patterns of activity (Extended Data Fig. 8 and Supplementary Information 7). We further observe that incremental, sequential removal of long-wavelength modes has a much greater impact on reconstruction accuracy than removal of short-wavelength modes (Fig. 3b and Contributions of long- and short-wavelength modes in Methods). For instance, across the seven key HCP task contrasts, removal of the top 25% long-wavelength modes (modes 1–50) yields a drop in reconstruction accuracy of around 40–60% whereas removal of the top 25% short-wavelength modes (modes 151–200) yields a drop of only around 2–4% (Fig. 3b, insets). These results indicate that, on temporal and spatial scales accessible with fMRI, evoked cortical activity comprises large-scale, nearly brain-wide spatial patterns, challenging classical views that such activity should be described in terms of discrete, isolated and anatomically localized activation clusters.

Fig. 3: Task-evoked activity excites long-wavelength modes.

a, Normalized mean power spectra of the 47 HCP task-contrast maps (left) and 10,000 contrast maps from the NeuroVault database (right). Insets show cortical surface reconstructions demonstrating spatial scales relevant to the first 50, 100 and 200 modes corresponding to spatial wavelengths of approximately 60, 40 and 30 mm, respectively. Contrast-specific spectra for the seven key HCP task contrasts are presented in Extended Data Fig. 7. b, Reconstruction accuracy of the seven key HCP task-contrast maps as a function of the percentage of modes (among 200) removed in the reconstruction process. Solid and dashed lines correspond to removal of the top long- and short-wavelength modes, respectively. Insets show group-averaged empirical activation maps (data) and their reconstructions after removal of 25% of modes. Negative, zero and positive values are coloured blue, white and red, respectively.

Wave dynamics bridge geometry and function

Geometric eigenmodes of the cortex are obtained by solving the eigenvalue problem of the LBO, which is also known as the Helmholtz equation (equation (1)). In physically continuous systems, the solutions of the Helmholtz equation correspond to the spatial projections of the solutions of a more general wave equation, such that the resulting eigenmodes inherently represent the vibrational patterns, or standing waves, of the system’s dynamics³⁴. This equivalence implies that the superior efficacy of geometric eigenmodes in the reconstruction of diverse patterns of brain activity results from a fundamental role of wave dynamics in shaping these patterns, as predicted by NFT. This prediction has been confirmed through models of EEG recordings^21,35, but waves across the whole brain have only recently been observed in fMRI signals^36,37 and thus far lack a theoretical explanation. Here we use NFT and geometric eigenmodes to show that wave dynamics can provide a unifying account of diverse empirical and physiological phenomena observed at scales accessible with fMRI.

We model neural activity using an isotropic damped NFT wave equation without regeneration⁶ (Fig. 4a and NFT wave model in Methods). Under this model, activity propagates between points on the neocortex through their white-matter connectivity with a strength that decays approximately exponentially with distance (Supplementary Fig. 10 and Supplementary Information 1 and 8). To simulate resting-state neural activity we use a white noise input to mimic unstructured stochastic fluctuations²¹ (Modelling resting-state dynamics in Methods). We compare the performance of this simple wave model with a biophysically based neural mass model (balanced excitation–inhibition (BEI) model) that has been used extensively to understand resting-state fMRI signals³⁸ (Fig. 4a and Neural mass model in Methods). The neural mass model is closely aligned with the classical, connectome-centric view of brain function, representing dynamics as the result of interactions between neuronal populations in discrete anatomical regions, coupled according to an empirically measured connectome.

Fig. 4: Wave dynamics shape patterns of spontaneous and stimulus-evoked activity.

a, For the wave model, the activity $\varphi ({\bf{r}},t)$ at location r and time t is governed by a wave equation with damping rate γ_s, spatial length scale r_s and input Q. For the neural mass model the activity ${S}_{i}\left(t\right)$ of region i is described by the function f, which depends on the activity of other regions S, local population parameters θ_i and connectome C scaled by the global coupling parameter G. The model dynamics are used to calculate a simulated FC matrix (Methods). M_fixed and M_free correspond to the number of fixed and free parameters of each model, respectively. b, Comparison of data and model simulations based on various metrics which, from left to right, are: FC matrix (for visual purposes), edge FC, node FC and FCD. For edge FC and node FC, red lines represent linear fits with Pearson correlation coefficient r; for FCD, the probability density function (pdf) of the similarity of interregional synchrony in data and model dynamics is compared using the Kolmogorov–Smirnov (KS) statistic. c, Wave propagation of activity after 1 ms stimulation of V1 for t = 1–2 ms. Arrows indicate the direction of propagation (Supplementary Video 1). d, Activity profile of 17 regions in the visual cortical hierarchy. Insets show spatial locations of regions on the cortical surface coloured according to their activation profiles. e, Relationship of ranked activity profile time to peak and ranked T1w:T2w value of the regions in d. Red line represents a linear fit of the ranked variables, with Spearman correlation coefficient r and one-sided P_spin, estimated from 10,000 permutations.

We first compare the efficacy of the two models in capturing distinct and commonly studied properties of spontaneous, task-free FC: namely, static pairwise FC (edge FC), static node-level average FC (node FC) and time-resolved dynamic properties of FC (FCD) (Modelling resting-state dynamics in Methods). Across all FC-based benchmark measures, the wave model shows comparable or superior performance in reconstruction of empirical data relative to the neural mass model (Fig. 4b). The wave model also captures time-lagged properties^36,37,39 of empirical resting-state activity more accurately than the mass model (Extended Data Fig. 9 and Measuring time-lagged properties of resting-state dynamics in Methods). This strong performance of the wave model is remarkable given its relative simplicity: the wave model only requires the geometry of the cortex (that is, the surface mesh) as input and includes one fixed parameter and one free parameter (r_s) for fitting to data (Extended Data Fig. 10) whereas the neural mass model requires a dMRI-derived interregional anatomical connectivity matrix and comprises 15 fixed parameters and four free parameters (Supplementary Information 9). These considerations indicate that wave dynamics offer a more accurate and parsimonious mechanistic account of macroscale, spontaneous cortical dynamics captured by fMRI.

We next consider stimulus-evoked cortical activity in the wave model. We analyse cortical responses to sensory stimulation of primary visual cortex (V1), because it elicits a well-defined hierarchy of regional cortical responses^40,41 (Modelling stimulus-evoked dynamics in Methods). A 1 ms pulse input to V1 yields a propagating wave of activity that rapidly splits along the dorsal and ventral visual processing streams (Fig. 4c (arrows) and Supplementary Video 1), consistent with the mainstream understanding of hierarchical visual processing⁴². Remarkably, this result indicates that geometric constraints on travelling waves of evoked activity are sufficient for the segregation of the dorsal and ventral processing streams, which have traditionally been thought to be driven primarily by complex patterns of layer-specific connectivity^40,42,43. Furthermore, the temporal profile of evoked responses across the visual system follows a well-defined timescale hierarchy, with higher-order association areas showing peak responses that are delayed and prolonged compared with lower-order visual areas (Fig. 4d). These findings thus indicate that this hierarchical ordering, previously identified in experimental and modelling studies^41,44,45, emerges naturally from waves of excitation propagating through the cortical medium. Critically, this hierarchical temporal ordering of areal responses strongly correlates with an independent anatomical measure of the cortical processing hierarchy based on non-invasive estimates of myeloarchitecture (T1-weighted (T1w) and T2-weighted (T2w) ratio)^46,47. This correlation is particularly strong within the visual processing hierarchy (r = −0.72, one-sided spin-test P value (P_spin) = 0.003; Fig. 4e) but is also present when considering all cortical areas (r = −0.44, P_spin = 0.037; Supplementary Fig. 11). Together, our modelling results show how simple wave dynamics unfolding on the geometry of the cortex provide a unifying generative mechanism for capturing complex properties of spatiotemporal brain activity.

Geometry constrains subcortical activity

Our analyses thus far have focused on the strong coupling of geometry and dynamics in the neocortex. We next investigate this coupling in non-neocortical areas, focusing on the thalamus, striatum and hippocampus, because these structures have geometries easily captured using MRI data and their functional organization has been extensively studied⁴⁸.

We first generalize our eigenmode analysis to three-dimensional (3D) volumes (Estimating the geometric eigenmodes of non-neocortical structures in Methods), yielding geometric eigenmodes that extend spatially through the three spatial dimensions of each structure. Next, to fully capture the macroscale functional organization of these non-neocortical regions, we apply a widely used manifold learning procedure to voxel-wise FC data to obtain the key functional gradients in each structure⁴⁹ (Mapping the functional organization of non-neocortical structures in Methods). These functional gradients describe the principal axes of spatial organization dictated by similarities in FC, thus representing the dominant modes of variation in functional organization, ordered according to the percentage of variance in FC similarity that they explain.

The spatial profiles of the first three functional gradients of the thalamus, striatum and hippocampus (accounting for 24, 50 and 47% of the variance in FC similarity, respectively) show a near-perfect match to the first three geometric eigenmodes (Fig. 5a–c; spatial correlations r ≥ 0.93). This tight correspondence generalizes out to the first 20 modes and first 20 gradients of each structure (with the first 20 gradients respectively accounting for 49, 70 and 68% of total variance in FC similarity), with all absolute spatial correlations |r| > 0.5, except for the 20th gradient and 20th mode in the striatum and hippocampus (Fig. 5d–f). This strong relationship is striking given that the functional gradients are generated via a complex processing pipeline applied to fMRI-derived FC measures whereas the eigenmodes are derived simply from each structure’s geometry, independent of the functional data. These findings suggest that the functional organization of non-neocortical structures derives directly from their geometric eigenmodes.

Fig. 5: Geometry shapes non-neocortical function.

a–c, First three geometric eigenmodes and FC-based functional gradients in the thalamus (a), striatum (b) and hippocampus (c). The modes and gradients are shown in 3D coordinate space with negative, zero and positive values coloured blue, white and red, respectively. D, dorsal; P, posterior; R, rightward. Scatter plots show the relationship between modes and gradients, with red lines representing linear fits with Pearson correlation coefficient r. d–f, Bottom, absolute correlation (|r|) of the first 20 geometric eigenmodes and functional gradients in the thalamus (d), striatum (e) and hippocampus (f). Top, highest |r| value obtained for each functional gradient (grey bars), taking into account order flips in geometric eigenmodes, along with the percentage of variance explained by each functional gradient (blue lines). Max., maximum.

Discussion

The dynamics of many physical systems are constrained by their geometry and can be understood as excitations of a relatively small number of structural modes¹². Here we show that structural eigenmodes derived solely from the brain’s geometry provide a more compact, accurate and parsimonious representation of its macroscale activity than alternative connectome-based models. This mode-based view of the brain further shows that both spontaneous and evoked brain activity captured by fMRI are dominated by large-scale eigenmodes with relatively long wavelengths, whose dynamics are derived from a biophysically motivated wave equation. These findings challenge the classical neuroscientific paradigm in which topologically complex patterns of interregional connectivity between discrete, specialized neuronal populations are viewed as a critical anatomical foundation for dynamics. Instead, our results indicate that a physically grounded approach that treats the brain as a continuous, spatially embedded system offers a unifying framework for understanding structural constraints on diverse aspects of macroscopic neuronal function.

The extensive comparisons of geometric eigenmodes with other anatomical (connectome and EDR eigenmodes) and statistical (PCA and Fourier) basis sets show that the superior performance of geometric eigenmodes in capturing macroscale neocortical activity is not trivially driven by generic mathematical properties of basis set expansions. Rather, this result indicates that geometry represents a fundamental anatomical constraint on dynamics. Additionally, the strong performance of the EDR eigenmodes derived from a synthetic network suggests that a homogeneous, distance-dependent connectivity with near-exponential form represents another important anatomical constraint on activity. EDR-like connectivity is mathematically embedded within the Helmholtz equation in equation (1)^6,7 (Supplementary Information 8), so the role of such connectivity is implicitly captured by the geometric eigenmodes.

The comparatively poor performance of connectome eigenmodes indicates that topologically complex connections that exist beyond a simple EDR afford minimal further benefit in obtaining eigenmodes that can accurately explain spatiotemporal patterns of cortical activity as measured with fMRI. Our findings thus counter traditional views that emphasize intricate patterns of anatomical connections as the primary driver of coordinated dynamics^26,50,51,52. Indeed, recent estimates indicate that long-range cortical connections are rare⁵³—they may therefore represent a relatively minor perturbation of the dominant effect imposed by EDR-like connectivity. Nonetheless the topological centrality, metabolic cost and tight genetic control of such connections^54,55 suggest that they provide important functional and evolutionary advantages beyond wave-like dynamics⁵⁶ (Supplementary Information 11). The limited resolution and sensitivity to preprocessing pipelines of dMRI and fMRI data^57,58 complicate attempts to fully uncover the functional role of these connections. High-quality animal tract-tracing and electrophysiological data may be helpful in this regard.

The close coupling between geometry and dynamics is apparent in neocortical and non-neocortical structures alike, suggesting that the functional organization of regions outside the neocortex is also dominated by distance-dependent anatomical connectivity and wave dynamics, as found in recent experiments^36,59,60. These observations indicate that geometric eigenmodes offer a simpler, more parsimonious and mechanistically informative account of putative gradients of functional organization in non-neocortical structures than the complex manifold learning procedures currently used in the literature⁴⁹. This is because such procedures are phenomenological, providing statistical descriptions of dominant sources of variances in the data, whereas the study of structural eigenmodes derives from a generative process.

Geometric mode decomposition offers unique insights into the spatial properties of brain activation maps. Classical brain-mapping analyses focus on responses in isolated clusters of spatial locations that exceed a statistical threshold⁶¹. By contrast, our approach aligns with rigorously established results from physics and engineering in which perturbations of spatially continuous systems elicit system-wide responses, just as the musical notes of a violin string result from vibrations across its entire length rather than from the behaviour of a restricted segment⁶². Accordingly, the use of geometric eigenmodes indicates that, across more than 10,000 diverse maps from task-based fMRI studies, task engagement is associated predominantly with the excitation of modes with wavelengths of roughly 60 mm and longer. This result coincides with similar observations of long-wavelength excitations in empirical EEG and evoked response potential data^30,31 and suggests that classical analyses reliant on thresholding of pointwise statistical maps obscure the spatially extended and complex patterns of activity actually evoked by a task.

Our modelling results offer insight into the physical processes that underlie the close link observed between geometry and function. In particular, the relative simplicity and superior performance of the wave model in capturing diverse aspects of spontaneous fMRI dynamics indicate that the model provides a more parsimonious account than a complex neural mass model, which treats the brain as a graph of discrete anatomical regions (nodes) coupled via the connectome (edges). This finding is consistent with experimental observations of wave dynamics in both human and animal fMRI data^63,64. Future work could explore whether the introduction of spatial heterogeneities⁶⁵ or complex structured input⁶⁶ into the wave model further improves its accuracy in explaining diverse empirical phenomena.

Application of the wave model to mimic visual stimulation shows that waves propagating from the stimulation site segregate along the classical dorsal and ventral visual pathways, and that regional responses to the perturbation conform to a well-described hierarchy of timescales that range from rapidly responding unimodal areas to slower transmodal regions^41,44,45. These canonical properties of hierarchical visual processing have been extensively studied for decades and are conventionally thought to be driven by complex patterns of layer-specific interregional connectivity^40,42,43, but our analysis shows that waves travelling through the cortical geometry are sufficient for the emergence of segregated, hierarchical processing streams. Thus, while our findings cannot rule out a role for complex interregional connectivity they do indicate that such connectivity is not necessary for the emergence of these macroscale dynamics.

The superior performance of geometric eigenmodes offers an immediate practical benefit because the modes can be estimated using only a mesh representation of the structure of interest, which can easily be derived using well-established, automated processing pipelines for T1w anatomical images⁶⁷. By contrast, connectome eigenmodes require a graph-based model of macroscopic interregional connectivity generated via complex data-processing pipelines applied to both T1w and dMRI images^58,68; the definition of graph nodes, which is a topic of contention⁶⁹; and the application of a thresholding procedure to remove putatively spurious connections, which our own analysis shows can affect the findings (Supplementary Figs. 6 and 7). The fact that such choices are not required to obtain the geometric eigenmodes means that they can be applied robustly and flexibly across different experimental contexts in both humans and other species^70,71, opening new avenues of research. For example, one can investigate how geometric eigenmodes vary through neurodevelopment or are disrupted in clinical disorders. Indeed, the close link we identify between geometry and function implies that interspecies differences in spatiotemporal dynamics may largely be driven by differences in brain shape. The characterization of how variations in brain geometry, both within and between species, shape brain function will be essential for understanding physical and anatomical constraints on neuronal activity.

Methods

Derivation of cortical geometric eigenmodes

If brain structure can be approximated as being constant in time, the resulting spatial and temporal dynamics can be treated separately via eigenmode decomposition^1,7, similar to the treatment of other physical systems¹². In particular, the spatial aspect satisfies the Laplacian eigenvalue problem, which is also known as the Helmholtz equation, defined in equation (1).

For the cerebral cortex, which we consider as a two-dimensional (2D) model embedded within 3D Euclidean space, the LBO in equation (1) captures intrinsic geometry, which includes the curvature of the cortical surface⁷² and is defined generally as^73,74,

$$\varDelta {\rm{:=}}\frac{1}{W}\sum _{i,j}\frac{\partial }{\partial {x}_{i}}\left({g}^{ij}W\frac{\partial }{\partial {x}_{j}}\right),$$

(2)

where x_i,x_j are the local coordinates, g^ij is the inverse of the inner product metric tensor ${g}_{ij}:\,=\langle \frac{\partial }{\partial {x}_{i}},\frac{\partial }{\partial {x}_{j}}\rangle$, W := √det(G), det denotes the determinant and $G:\,=(\,{g}_{ij})$.

We employed the LaPy python library^72,75 installed in the MASSIVE high-performance computing facility⁷⁶ to derive the geometric eigenmodes of the human cortex. Specifically, we used a triangular surface mesh representation of the midthickness human cortical surface, comprising 32,492 vertices in each hemisphere, obtained from a downsampled, left–right symmetric version of the FreeSurfer’s fsaverage population-averaged template⁷⁷ (https://github.com/ThomasYeoLab/CBIG/tree/master/data/templates/surface/fs_LR_32k). This template is independent of the data sample used in all our analyses, thus obviating any concerns about circularity.

Note that the continuous LBO operates on the underlying Riemannian manifold of the surface and not directly on mesh vertices. LaPy uses the cubic finite element method on the surface mesh to achieve numerically tractable solutions of equation (1) on an interpolated smooth manifold. This distinguishes it from the discrete graph Laplacian⁷⁸, which does not encode spatial relations between points. All our analyses were focused on unihemispheric eigenmodes, but our approach can easily be extended to the whole brain because bihemispheric eigenmodes can be represented as symmetric or antisymmetric combinations of the eigenmodes derived from each hemisphere⁷; symmetric combinations correspond to mirror symmetry across the sagittal midplane and asymmetric combinations correspond to cases in which the hemispheres have the same spatial structure but with flipped signs.

The eigenvalue solutions of equation (1) are ordered sequentially according to the spatial frequency or wavelength of the spatial patterns of each eigenmode—that is, $0\le {\lambda }_{1}\le {\lambda }_{2}\le \ldots$. Note that the first eigenvalue, λ₁, is approximately equal to zero (wavelength much greater than brain size) and the corresponding eigenmode, ψ₁, is a constant function with no nodal lines (zero sets of the function). Throughout our study, we used the first 200 modes (including the constant mode, ψ₁) in our analyses given the diminishing improvements in reconstruction accuracy observed when using an increasing number of modes (Fig. 1d).

Each eigenmode comprises spatial patterns with a specific spatial wavelength. Following ref. ⁷, we approximate the eigenmode wavelengths using an idealized spherical case because it is topologically comparable to the human cortex. By solving equation (1) on a sphere, degenerate solutions exist such that certain eigenmodes have the same eigenvalue and spatial wavelength—this is analogous to spherical harmonics in quantum physics. In fact, because the eigenmodes will approach the spherical harmonics in the limit of vanishing cortical folding⁷, the former can be grouped together into an eigengroup with spatial wavelength⁷¹,

$${\rm{wavelength}}=\frac{2\pi {R}_{{\rm{s}}}}{\sqrt{l(l+1)}},$$

(3)

where R_s is the radius of the sphere (for the fsaverage population-averaged template used in this study, R_s ≈ 67.0 mm) and l is the eigengroup number (the angular momentum quantum number in atomic physics). The wavelengths of the first 15 eigengroups and the eigenmodes included in the eigengroup are shown in Supplementary Table 1.

Modal decomposition of brain activity

We used the geometric eigenmodes to decompose spatiotemporal fMRI data, measured at spatial location r and time t, for each individual, as a weighted sum of modes,

$$y({\bf{r}},t)=\mathop{\sum }\limits_{j=1}^{N}{a}_{j}(t){\psi }_{j}({\bf{r}}),$$

(4)

where a_j is the amplitude of mode j in explaining the data, ψ_j is the jth mode and N is the number of modes used; we used N = 200 for our analyses. For spatiotemporal data—that is, recordings of spontaneous dynamics from task-free fMRI—each time frame of the data was substituted into equation (4), resulting in a time-dependent amplitude ${a}_{j}\left(t\right)$ for each mode ψ_j. For purely spatial data—that is, task-evoked activation maps—amplitudes are independent of time such that ${a}_{j}\left(t\right)\to {a}_{j}$. In both cases, the amplitudes can be obtained by integrating over the cortical surface,

$${a}_{j}(t)=\int y({\bf{r}},t){\psi }_{j}({\bf{r}})d{\bf{r}},$$

(5)

which can be derived from equation (4) using the orthogonal property of the eigenmodes^62,79. If there are insufficient measurements to evaluate the integral, the amplitudes can also be estimated via a statistical general linear model.

After obtaining the amplitudes, equation (4) was used to calculate the reconstructed data. We quantified the accuracy of this reconstruction by calculating the correlation between empirical and reconstructed data. For the spatiotemporal task-free data, we first parcellated the empirical and reconstructed data by taking the average of the data within discrete parcels/regions as per standard practice in the field⁸⁰, and then constructed a matrix of interregional FC by calculating the Pearson correlation coefficient of pairs of parcel time series. For task-evoked data, we applied the same parcellation on the activation maps to allow direct comparison. Finally, reconstruction accuracy for task-free data was calculated by taking the correlation of the upper triangular elements of the empirical and reconstructed FCs. For task-evoked data, reconstruction accuracy was calculated from the spatial correlation of parcellated empirical and reconstructed maps. We then took the average reconstruction accuracy across all participants.

HCP data

We used preprocessed fMRI data from HCP²⁷. We did not perform any additional preprocessing steps, such as global signal removal, because the first eigenmode (considered as the global, constant mode) already explicitly captures global deviations in the data, allowing the other modes to capture functionally relevant non-global activity. We analysed data from 255 unrelated healthy individuals (aged 22–35 years, 132 females and 123 males), which is the largest HCP sample excluding twins or siblings, and with all participants having completed task-evoked and task-free resting-state data. All participants were volunteers and provided informed consent. The open-access data were acquired by the WU–Minn HCP consortium with local overseeing ethics committee approval and were shared with the authors according to HCP’s data use terms. All our procedures were carried out in accordance with protocols set by these data use terms. For a detailed account of the image acquisition protocol, preprocessing pipelines and ethics oversight, see refs. ^27,32.

Within the HCP dataset we analysed task-evoked fMRI measured in seven task domains (Supplementary Table 2) and task-free resting-state fMRI, which were already preprocessed by HCP. Both datasets were mapped onto the fsLR-32k CIFTI space with 32,492 vertices in each hemisphere. Further details on data and preprocessing can be found in Supplementary Information 2.1 and 2.2. We also analysed individual connectomes derived from dMRI data, as provided in ref. ⁸¹. The connectomes represent a high-resolution weighted matrix of size 32,492 × 32,492 for each hemisphere. Further details on data and connectome construction can be found in Supplementary Information 2.3.

Cortical parcellations

The main results of the study analysed task-evoked activation maps and task-free FC data parcellated into discrete regions. We presented results using the HCP–MMP1 parcellation with 180 regions per hemisphere (we term this the Glasser360 parcellation), which reflects sharp areal boundaries based on the combination of cortical architecture, function, connectivity and topography²⁸. The parcellation was on the fsLR-32k space as provided by HCP. To test the robustness of our results (Supplementary Fig. 2 and Extended Data Figs. 2 and 3) we also performed our analysis using parcellations provided by Schaefer et al.⁸² on the fsLR-32k space at varying resolution (100, 200, 400, 600, 800 and 1,000 parcels across both hemispheres; we refer to these as, respectively, Schaefer100, Schaefer200, Schaefer400, Schaefer600, Schaefer800 and Schaefer1000 parcellations).

Derivation of connectome eigenmodes

Connectome eigenmodes were derived according to previous methods²⁹ to enable comparison with previous findings. Note that, in previous studies^19,20,29,83, connectome eigenmodes have been referred to as connectome harmonics but we use the term eigenmodes here because the term harmonics implies integer frequency ratios, which is not necessarily guaranteed for brain-derived modes.

We obtained high-resolution maps of connectivity measured with dMRI tractography as described in ref. ⁸¹, in which the connectivity of each of the 32,492 vertices in the cortical surface mesh was estimated by tracing streamlines from each point until they terminated at some other point. Connection weights between vertices (considered as nodes) were estimated as the number of interconnecting streamlines without the need for normalization⁸⁴. Tractography was performed on individuals from HCP (see Supplementary Information 2.3 for further details on the data and tractography method). From the tractography data we combined the individual weighted connectivity matrices of size 32,492 × 32,492 to generate a group-averaged connectome, W_connectome, with weights representing the average number of streamlines. We then generated a binary adjacency matrix, A_local, that captures the discrete representation of local spatial relations between points in the cortical surface mesh model constructed by connecting two vertices that are direct neighbours in the mesh. These links are intended to capture local, very-short-range connectivity that cannot be resolved by traditional dMRI tractography²⁹.

Following ref. ²⁹, the group-averaged weighted connectome, W_connectome, was thresholded to remove the smallest weights such that the number of connections was fourfold greater than A_local. The resulting thresholded matrix was binarized to obtain the group adjacency matrix, A_connectome. Finally we generated a merged adjacency matrix ${A}_{{\rm{C}}}={A}_{{\rm{local}}}| | {A}_{{\rm{connectome}}}$ (with matrix of size 32,492 × 32,492; || is the logical OR operator), which captures both local vertex-to-vertex connectivity and complex short- and long-range connections as measured empirically. Note that the connection density of the adjacency matrix, A_C, that resulted from the thresholding process described above was 0.10%.

The connectome eigenmodes were obtained by solving the eigenvalue problem,

$${L}^{{\prime} }\psi =-\lambda \psi ,$$

(6)

where L′ is the normalized graph Laplacian, a discrete counterpart of the LBO. The normalized graph Laplacian is related to the unnormalized graph Laplacian, L, as ${L}^{{\prime} }={D}^{-1/2}L{D}^{-1/2}$, with L defined according to previous work³⁴,

$$L=\frac{1}{2}\left[\left(D-{A}_{{\rm{C}}}\right)+{(D-{A}_{{\rm{C}}})}^{T}\right],$$

(7)

where D is the diagonal degree matrix and superscript T denotes matrix transpose. As with geometric eigenmodes, the eigenvalue solutions of equation (6) form the sequence $0\le {\lambda }_{1}\le {\lambda }_{2}\le \ldots$. Note also that the use of a high-resolution, vertex-level connectome results in connectome eigenmodes spanning a space with dimensions (number of modes) equal to the number of vertices, allowing fair comparison with geometric eigenmodes.

Derivation of EDR eigenmodes

Exponential distance rule eigenmodes were derived also by solving equation (6). However, in this case the non-normalized graph Laplacian, L, was defined using a synthetically constructed 32,492 × 32,492 adjacency matrix, A_E, that follows a stochastic EDR. To construct A_E we used the group-averaged, unthresholded weighted connectome, W_connectome, from the previous section. We then fitted the variation of weights as a function of Euclidean distance, d, between vertices in the cortical surface by an exponential function of the form ${e}^{-\alpha d}$, where α is a scale exponent parameter. The fitting was performed in MATLAB 2019b using a nonlinear least-squares method, resulting in an optimal empirical parameter value of ${\alpha }_{{\rm{empirical}}}=0.12$, consistent with previous estimates based on the connection probability versus distance function⁸⁵. We then generated a random, binary adjacency matrix following a stochastic EDR wiring process in which two vertices were connected with probability ${e}^{-\alpha d}$ with $\alpha ={\alpha }_{{\rm{empirical}}}=0.12$.

The adjacency matrix, A_E, generated by instances of the EDR model had a connection density of 1.55%, which was substantially higher than the 0.10% density of A_C. We therefore constructed another version of connectome eigenmodes that thresholded the group-averaged empirical connectome to achieve a final A_C with a density that matched A_E to allow fair comparison (we termed this density-matched connectome eigenmodes). Supplementary Figs. 6 and 7 present a more thorough evaluation of the performance of connectome eigenmodes as a function of network connection density.

Comparisons with statistical basis sets

The geometric, connectome and EDR eigenmode basis sets were all derived from a generative model that accounts for how brain function emerges from anatomy. This approach contrasts with the statistical basis sets commonly used in the literature, which can efficiently summarize the data but offer no insights into the underlying generative process. We evaluated the performance of geometric eigenmodes with respect to two statistical basis sets, one derived from PCA of the functional data themselves and the other based on a Fourier spatial basis set. Further details of these analyses are provided in Supplementary Information 5 and 6, respectively.

Modal power spectra of task-evoked activation maps

To investigate the spectral content of task-evoked activation maps we calculated their modal power spectra using the modal decomposition in equation (4) and taking the absolute square of the amplitudes, a. This is analogous to calculation of the temporal power spectral density from Fourier analysis. We then normalized the power in mode j with respect to the total power in all modes, such that

$${P}_{j}=\frac{{|{a}_{j}|}^{2}}{{\sum }_{j=1}^{N}{|{a}_{j}|}^{2}}.$$

(8)

We calculated the modal power spectra of two sets of task-evoked data. The first set comprises unthresholded activation maps from the HCP task-evoked data (Supplementary Information 2.1). Before perfoming spectral analysis we took the group average of the activation maps across 255 individuals and then analysed the power spectrum of the group-averaged activation map of each task contrast. The results in Fig. 3a and Extended Data Fig. 8a show the mean power spectrum across the 47 HCP task-contrast maps. Contrast-specific power spectra for the seven key HCP task contrasts are shown in Extended Data Fig. 7.

The second set comprises 10,000 activation maps from 1,178 independent experiments from the NeuroVault repository (https://neurovault.org/)³³. We used the python module Nilearn⁸⁶ to retrieve activation maps from NeuroVault that were unthresholded and with a modality tag of fMRI-BOLD. We projected the activation maps from volume onto the fsLR-32k CIFTI space to match the HCP data using Nilearn (via the function nilearn.surface.vol_to_surf). We then analysed the power spectrum of each activation map. The results in Fig. 3a and Extended Data Fig. 8a show the mean power spectrum across the 10,000 NeuroVault maps.

We then compared the power spectra of the empirical maps with those of surrogate random maps with varying levels of smoothing to further investigate the relevance of long-wavelength modes. In particular, we generated 10,000 random maps in volume space, which we smoothed at kernel sizes with full-width at half-maximum ranging from 0 to 50 mm, and projected onto the fsLR-32k CIFTI space. We then calculated the mean square logarithmic error (MSLE) of the power spectra between the empirical and surrogate maps. Further details of this analysis and the MSLE measurements are provided in Supplementary Information 7.

Contributions of long- and short-wavelength modes

To understand the contributions of long- and short-wavelength geometric eigenmodes in the reconstruction of task-evoked activation maps, we sequentially removed modes before performing the reconstruction process. Specifically, we started by reconstructing the seven key HCP contrast maps using 200 geometric eigenmodes and calculated the reconstruction accuracy (that is, correlation between empirical and reconstructed maps), serving as our baseline. We then performed an incremental, sequential removal of modes starting from long-wavelength modes (that is, removal of mode 1, modes 1–2, modes 1–3, modes 1–4, …, modes 1–200) and calculated reconstruction accuracy at every increment. We repeated the same procedure but starting from short-wavelength modes (that is, removal of mode 200, modes 199–200, modes 198–200, modes 197–200, …, modes 1–200).

NFT wave model

As stated above, an implication of the superior performance of geometric eigenmodes is that neuronal activity is dominated by wave dynamics as predicted by NFT. To investigate whether wave dynamics can explain complex spatiotemporal patterns of neuronal activity, we implemented a simple NFT wave model in which dynamics are described by an isotropic damped wave equation without regeneration^5,6,87,

$$\left[\frac{1}{{\gamma }_{s}^{2}}\frac{{{\rm{\partial }}}^{2}}{{\rm{\partial }}{t}^{2}}+\frac{2}{{\gamma }_{s}}\frac{{\rm{\partial }}}{{\rm{\partial }}t}+1-{r}_{s}^{2}{{\rm{\nabla }}}^{2}\right]\varphi ({\bf{r}},t)=Q({\bf{r}},t),$$

(9)

where $\varphi ({\bf{r}},t)$ is the neural activity at location r and time t, Q is an external input, γ_s is the damping rate and r_s is the spatial length scale of the wave propagation (conceptually related to the α-length scale of the stochastic EDR wiring process in the derivation of EDR eigenmodes). This form tells us that an impulse input will produce an activity that dissipates at a rate of γ_s and propagates at a velocity of γ_sr_s. Here we treated γ_s as a fixed parameter with the value of 116 s⁻¹ taken from electrophysiological estimates²¹ and r_s as a free parameter. We applied the model on the cortical midthickness surface mesh with 32,492 vertices per hemisphere to solve the activity at each vertex. Note that the propagation of activity between points is governed by their white-matter connectivity, with strength that decays approximately exponentially with distance (Supplementary Fig. 10). This distance dependence is more apparent when equation (9) is converted into its equivalent integral form (see Supplementary Information 8 for details).

Neural mass model

We compared the performance of our simple NFT wave model with a biophysical large-scale neural mass model in which mesoscopic dynamics emerge from the interactions of neural populations (that is, neural masses) coupled via an empirical anatomical connectivity¹⁸. In a typical neural mass model, each brain region i has its own mean-field population dynamics and its temporal activity, S_i, is defined by the general equation,

$$\frac{{\rm{d}}{S}_{i}}{{\rm{d}}t}=f({\bf{S}},{\theta }_{i},C,G),$$

(10)

where f is a function describing the evolution of the region’s activity. The function depends on the activity of other regions S, local population parameters θ_i, anatomical connectivity between regions C (a parcellated version of the vertex-resolution connectome, W_connectome, used to derive the connectome eigenmodes) and global coupling parameter G that scales the connectivity between regions.

There are several whole-brain neural mass models available in the literature that we can use⁸⁸, from simple phase oscillator models (for example, the Kuramoto model⁸⁹) to more complex biophysical population models (for example, the Wilson–Cowan model⁹⁰). All these models follow the form of equation (10), especially their reliance on an anatomical interregional connectivity matrix, C. Here we focus on one widely used neural mass model, the BEI model^38,65,91,92. The BEI model uses a mean-field approach to approximate local population dynamics in each brain region, which are coupled via an anatomical connectivity matrix derived from dMRI. Each brain region i comprises interacting populations of excitatory (E) and inhibitory (I) neurons governed by the following nonlinear stochastic differential equations,

$$\frac{{\rm{d}}{S}_{i}^{(E)}}{{\rm{d}}t}=\,-\,\frac{{S}_{i}^{(E)}}{{\tau }_{E}}+\left(1-{S}_{i}^{(E)}\right)\gamma {r}_{i}^{(E)}+\sigma {\nu }_{i}(t),$$

(11)

$$\frac{{\rm{d}}{S}_{i}^{(I)}}{{\rm{d}}t}=\,-\,\frac{{S}_{i}^{(I)}}{{\tau }_{I}}+{r}_{i}^{(I)}+\sigma {\nu }_{i}(t),$$

(12)

$${r}_{i}^{(E)}={H}^{(E)}\left({I}_{i}^{(E)}\right)=\frac{{a}_{E}{I}_{i}^{(E)}-{b}_{E}}{1-\exp \left[-{d}_{E}\left({a}_{E}{I}_{i}^{(E)}-{b}_{E}\right)\right]},$$

(13)

$${r}_{i}^{(I)}={H}^{(I)}\left({I}_{i}^{(I)}\right)=\frac{{a}_{I}{I}_{i}^{(I)}-{b}_{I}}{1-\exp \left[-{d}_{I}\left({a}_{I}{I}_{i}^{(I)}-{b}_{I}\right)\right]},$$

(14)

$${I}_{i}^{(E)}\left(t\right)={I}^{{\rm{ext}}}+{W}_{E}{I}_{0}+{w}_{EE}{S}_{i}^{(E)}\left(t\right)+GJ\sum _{j}{C}_{ij}{S}_{j}^{(E)}\left(t\right)-{w}_{IE}{S}_{i}^{(I)}\left(t\right),$$

(15)

$${I}_{i}^{(I)}\left(t\right)={W}_{I}{I}_{0}+{w}_{EI}{S}_{i}^{(E)}\left(t\right)-{S}_{i}^{(I)}\left(t\right),$$

(16)

where ${S}_{i}^{(E,I)}$, ${r}_{i}^{(E,I)}$ and ${I}_{i}^{\left(E,I\right)}$ represent the synaptic gating variable, firing rate and total input current, respectively, for E and I populations. The parameters τ_E,I are the time constants, γ is a kinetic rate constant and v_i(t) is a time-varying random Gaussian input with standard deviation σ. The functions H^(E,I) are sigmoidal neuronal response functions transforming total input currents ${I}_{i}^{(E\,,I)}$ into firing rates ${r}_{i}^{(E,I)}$, which are parametrized by the gain factors a_E,I, threshold currents b_E,I and curvature parameters d_E,I. In equations (15) and (16), ${I}_{i}^{{\rm{ext}}}$ is the external input, I₀ is the local input current scaled by W_E and W_I for the excitatory and inhibitory populations, respectively, w_EE is the excitatory–excitatory strength, w_EI is the excitatory–inhibitory strength, w_IE is the inhibitory–excitatory strength, G is the global coupling parameter, J is the effective N-methyl-d-aspartate (NMDA) conductance and C_ij is the structural connectivity strength between regions i and j estimated from dMRI.

Overall, the BEI model has 15 fixed parameters and four free parameters, as detailed in Supplementary Table 3. The values of the fixed parameters were taken from ref. ⁹¹. To allow direct comparison we also used the discretized connectome data provided by ref. ⁹¹ to define the interregional structural connectivity matrix, C. These connectome data were derived from minimally preprocessed dMRI data of 334 unrelated HCP subjects and constructed via the probabilistic tractography tool of FMRIB Software Library (FSL)⁹³. We refer readers to ref. ⁹¹ for further details of data processing and connectome construction. Our results did not change when the parcellated version of our connectome data was used. We note that numerical solutions of NFT equations (such as in equation (9)) also spatially discretize the cortex, but into a very fine array of points, before integrating whatever temporal dynamics have been chosen at each point with connectivity that is correct in the continuum limit. As such, neural mass models are approximations to the more general NFT approach^94,95. We also reiterate that local dynamics do not affect the spatial eigenfunctions of the NFT wave model.

Haemodynamic model

To simulate fMRI data we transformed the neural activity generated by the NFT wave and BEI neural mass models to a blood oxygen-level dependent (BOLD) signal using the well-established Balloon–Windkessel haemodynamic model⁹⁶. Note that, although this model is a simple approximation to more detailed models of the physiological haemodynamic processes underlying the BOLD signal^97,98,99, we use this approximation here to allow direct comparison with the vast majority of modelling studies in the literature^{65,91,92,100,101}. The BOLD–fMRI signal in each vertex or brain region i is governed by the following differential equations,

$$\frac{{\rm{d}}{z}_{i}}{{\rm{d}}t}={N}_{i}(t)-\kappa {z}_{i}(t)-\gamma [\,{f}_{i}(t)-1],$$

(17)

$$\frac{{\rm{d}}{f}_{i}}{{\rm{d}}t}={z}_{i}(t)\,,$$

(18)

$$\frac{{\rm{d}}{v}_{i}}{{\rm{d}}t}=\frac{1}{\tau }[\,{f}_{i}(t)-{v}_{i}^{1/\alpha }(t)],$$

(19)

$$\frac{{\rm{d}}{q}_{i}}{{\rm{d}}t}=\frac{1}{\tau }\left\{\frac{{f}_{i}(t)}{\rho }[1-{(1-\rho )}^{1/{f}_{i}(t)}]-{v}_{i}^{1/\alpha -1}(t)\right\},$$

(20)

$$\frac{{\rm{d}}{y}_{i}}{{\rm{d}}t}={V}_{0}\left\{{k}_{1}\left[1-{q}_{i}(t)\right]+{k}_{2}\left[1-\frac{{q}_{i}(t)}{{v}_{i}(t)}\right]+{k}_{3}\left[1-{v}_{i}(t)\right]\right\},$$

(21)

where z_i, f_i, v_i, q_i and y_i are the vasodilatory signal, blood inflow, blood volume, deoxyhaemoglobin content and BOLD signal variables, respectively. The variable N_i represents the neural activity generated by the wave and BEI neural mass models. For the wave model, N(t) is ϕ(t), and for the BEI neural mass model, ${N}_{i}\left(t\right)$ is ${S}_{i}^{\left(E\right)}(t)$. The model parameters and their values taken from previous work^96,102 were as follows: κ = 0.65 s⁻¹ is the signal decay rate, γ = 0.41 s⁻¹ is the flow-dependent elimination rate, τ = 0.98 s is the haemodynamic transit time, α = 0.32 is the Grubb’s exponent, ρ = 0.34 is the resting oxygen extraction fraction, V₀ = 0.02 is the resting blood volume fraction and k₁ = 3.72, k₂ = 0.53 and k₃ = 0.53 are 3T fMRI parameters.

Modelling resting-state dynamics

We used the NFT wave and BEI models, together with the haemodynamic model, to estimate various spontaneous FC properties. For the wave model we solved equation (9) with a white noise input to mimic the absence of any structured stimulus, following previous studies^21,103. We combined the resulting solution with the haemodynamic model in equations (17–21) to simulate the BOLD–fMRI signal. The simulated BOLD signal was downsampled to a sampling interval of 0.72 s with 1,200 time frames to match the resting-state HCP data described in Supplementary Information 2.2. Finally we parcellated the simulated BOLD signal using the HCP–MMP1 parcellation and calculated the correlation matrix representing the model FC.

For the BEI model we can use the method described above to calculate the model FC starting from the numerical solutions of equation (11) for each brain region. However, it has been shown that one can linearize the equations of the BEI and haemodynamic models to obtain an analytic approximation of the model FC⁹¹. Due to the large number of model parameters in the BEI model we used this analytic approximation in this study because it allows more comprehensive and computationally efficient model fitting. See ref. ⁹¹ for further details.

We fitted the model FCs to empirical FCs (also parcellated using the HCP–MMP1 parcellation) of the same HCP participants used in all our analyses by fitting each model’s free parameters. We divided the participant sample into training and test sets, each with 125 individuals, to enable out-of-sample evaluation of model performance and to avoid overfitting. In particular, we fitted the model parameters on the training set data and used the fitted model parameters to predict the test set data. Model fitting and performance evaluation were based on three widely used FC-related metrics^65,92: edge FC, node FC and FCD.

The edge FC metric was calculated by taking the Pearson correlation of the upper triangular elements (that is, the strength of FC edges) of the z-transformed model and empirical FCs. We only took the upper triangular elements because the FC values are symmetric with respect to the diagonal. Higher correlations represent a better fit between model and data.

The node FC metric was calculated by taking the Pearson correlation of the average FC strength of each brain region in the model and empirical FCs. The average FC strength of a brain region i was defined as $\frac{1}{n}{\sum }_{j=1}^{n}{{\rm{F}}{\rm{C}}}_{ij}$, where n is the total number of brain regions. Once again, higher correlations indicate better model fits.

The FCD metric captures the spatiotemporal statistics of resting-state activity and was calculated as follows⁹². The time series at each region i was filtered between 0.04 and 0.07 Hz using a second-order Butterworth filter; this band was based on ref. ¹⁰⁴ and its inclusion was motivated by its functional relevance to the brain^105,106. The time series was then Hilbert (H) transformed to calculate the quantity ${y}_{i}\left(t\right)={x}_{i}\left(t\right)+j{H}_{i}(t)$, where j is the imaginary number. The instantaneous complex argument ${\theta }_{i}\left(t\right)={\tan }^{-1}\left[{H}_{i}(t)/{x}_{i}(t)\right]$ was then computed. The level of synchrony between regions i and j at time t, $\Delta (i,j,t)$, was calculated as

$$\Delta \left(i,j,t\right)=\cos \left[{\theta }_{i}\left(t\right)-{\theta }_{j}(t)\right].$$

(22)

Note that we do not term ${\theta }_{i}\left(t\right)$ a phase, as is sometimes done in the literature, because our signals are broadband whereas this interpretation is applicable only to signals that are nearly monochromatic. We computed this quantity for comparison with previous work^65,92, then we calculated the similarity of global synchrony, φ_uv, between two time instances, τ_u and τ_v, with φ_uv defined as

$${\phi }_{uv}=\frac{1}{{d}_{u}{d}_{v}}\sum _{i > j}\Delta (i,j,{\tau }_{u})\Delta (i,j,{\tau }_{v}),$$

(23)

where

$${d}_{x}=\sqrt{\sum _{i > j}{[\Delta (i,j,{\tau }_{x})]}^{2}}.$$

(24)

Here, φ_uv is the FCD, which is a symmetric time × time matrix interpreted as the similarity of global synchrony between times τ_u and τ_v. We then compared the distributions of the upper triangular elements of the model and empirical FCD estimates, concatenated across all individuals or model realizations as per ref. ⁹², to better capture the general fluctuations of the data. The FCD distributions of the model and the data were compared using the KS statistic, with lower KS statistic indicating better model fit.

The wave model has one free parameter—the spatial length scale of wave propagation, r_s—that was optimized to fit the empirical fMRI data. Specifically, we constructed a vector of 20 values of r_s evenly spaced between 10 and 100 mm. For each r_s value the FC metrics described above were calculated. The optimization landscapes are shown in Extended Data Fig. 10. We took the value of r_s that minimized the FCD KS statistic because this metric has been found to be the most stringent benchmark for model–data comparisons among the three metrics discussed above⁹². This procedure resulted in an optimized value of r_s = 28.9 mm, which is smaller than those obtained previously in analyses of EEG data²¹, possibly due to either the haemodynamic processes limiting neural activity propagation to neighbouring regions or our focus here on cortico-cortical dynamics^97,98.

The BEI neural mass model has 15 fixed parameters and four free parameters—that is, w_EE, w_EI, w_IE and G—that were optimized to fit the data, following ref. ⁹¹. This procedure resulted in optimized parameters w_EE = 9.80, w_EI = 1.48, w_IE = 7.13 and G = 6.87. Further details of the optimization process can be found in Supplementary Information 9.

Measurement of time-lagged properties of resting-state dynamics

In addition to the FC-based metrics used for model fitting and evaluation of model performance, we investigated whether the wave and neural mass models could also capture the temporal properties of propagated activity. In particular we analysed the lag structure (or lag threads) of resting-state BOLD–fMRI time courses, as proposed by refs. ^39,107. We briefly discuss the algorithm for calculation of the lag structure of both empirical and simulated fMRI data in Supplementary Information 10 and refer readers to previous articles^39,107 for further details. Note that there are several other ways of characterizing the spatiotemporal properties of resting-state activity^37,108,109, but the method currently chosen is sufficient for our purposes.

Modelling stimulus-evoked dynamics

We also evaluated the degree to which the wave model could capture classical properties of evoked neural responses. Specifically we used the optimized wave model (with r_s = 28.9 mm) to investigate neural dynamics in response to a stimulus applied to the left V1. In particular, the stimulus Q(r,t) was a 1 ms pulse (t = 1–2 ms) with magnitude 20 s⁻¹ (the results are robust to changes in amplitude) restricted to vertices in the cortical surface that fall within the V1 region as defined by the HCP–MMP1 parcellation. We performed the simulation over a 100 ms time period with 0.1 ms resolution. Finally we parcellated the activity $\varphi ({\bf{r}},t)$ into 180 regions using the HCP–MMP1 parcellation.

We compared the activity profile (amplitude versus time) of each brain region, focusing on the time for the activity to reach peak amplitude (that is, time to peak). Specifically we investigated whether the temporal precedence of activity follows the human visual cortical hierarchy from visual to frontal cortices, which includes the following 17 brain regions: V1, V4, 7m, 7Am, TE1p, 7AL, 24dd, 2, 24dv, 8BM, 10r, 10v, 8BL, 10pp, 10d, 9-46v and 9-46d. These regions closely resemble areas previously identified in the macaque neocortex visual hierarchy using tract-tracing data and nonlinear network modelling⁴¹.

We also compared the regional estimates of peak response times to T1w:T2w values, which is a non-invasive measure sensitive to intracortical myelin content⁴⁶ and a good proxy for cortical hierarchy rank¹¹⁰. The myelin map in fsLR-32k space was obtained from the HCP dataset^27,28 and then parcellated using the HCP–MMP1 parcellation. We quantified the relationship between the regional values of time to peak and myelin content via Spearman rank correlation. The statistical significance of the correlation was assessed by comparison with a null distribution of 10,000 correlation values obtained using a spatially constrained spin-test approach^111,112. This approach calculates the correlation between one map and random spatially rotated versions of the other map, thereby preserving the spatial relationship of the parcels in the map. The resulting P value, P_spin, is the fraction of null correlation values greater than the empirical correlation value. Finally we repeated the above statistical test including all brain regions (that is, not restricted to visual hierarchy brain regions) to ensure that our findings were not driven by our particular selection of regions of interest (ROIs). This is a more conservative test because not all brain regions show strong evoked responses to visual stimulation.

Estimation of the geometric eigenmodes of non-neocortical structures

We extended our eigenmode analysis to regions outside the neocortex, focusing in particular on the subcortex (thalamus and striatum) and archicortex (hippocampus). Unlike the cortical ribbon, which can be modelled as a 2D sheet, these structures are solid 3D objects. We therefore calculated the geometric eigenmodes using a tetrahedral mesh rather than the surface-based triangular mesh to account for the full 3D geometry of the non-neocortical structures⁷², as outlined below.

We first used the probabilistic Harvard–Oxford subcortical atlas (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases) to generate a volumetric binary mask in each hemisphere for the thalamus, striatum and hippocampus. Voxels with a probability of 25% or more of belonging to these structures were included in the masks. We used FreeSurfer’s mri_mc function, which implements a marching-cubes algorithm, to construct a 2D surface by tessellating the volumetric masks, followed by the Gmsh software (https://gmsh.info/) to convert the 2D surface into a 3D tetrahedral mesh. We then used the LaPy python library to solve equation (1) on the tetrahedral mesh of each non-neocortical structure to obtain the eigenmodes and their corresponding eigenvalues. Finally we projected the eigenmodes in tetrahedral space back into the natural volumetric space via interpolation. Hence, the resulting eigenmodes spatially vary through the 3D voxels comprising each non-neocortical structure’s volume. For this part of the study we discarded the first constant mode and used the next 20 modes in our analyses.

Mapping the functional organization of non-neocortical structures

The signal-to-noise ratio of fMRI is generally weaker in non-neocortical structures than in the neocortex¹¹³. Moreover, fine-grained task activations are often hard to resolve in these smaller structures due to the limited resolution of fMRI and the spatial smoothing induced by common fMRI processing methods. Thus, to efficiently map the functional organization of the thalamus, striatum and hippocampus we used connectopic mapping⁴⁹ of the resting-state fMRI signals in each structure to obtain their dominant functional modes (often called gradients, although the patterns capture spatial variations in pointwise similarities of FC profiles). This technique, and related procedures, have been extensively used in past work to study the functional organization of these structures^114,115,116.

We applied connectopic mapping to the volumetric voxel-wise resting-state fMRI data of the HCP individuals following the procedure described in ref. ⁴⁹. Specifically, for each individual and non-neocortical ROI (that is, thalamus, striatum and hippocampus) we constructed an ROI time-series data matrix, A, of size T × N, where T is the number of time frames and N is the number of voxels in the ROI. Similarly, for each individual we constructed the grey matter time-series data matrix, B, of size T × M, where M is the number of grey matter voxels outside the ROI. Because M is generally large (10³ or above), we reduced the dimensionality of B using singular value decomposition to construct the matrix $\widetilde{B}$ of size T × (T – 1). The connectivity fingerprint of every voxel within the ROI was then calculated using the Pearson correlation of A and $\widetilde{B}$ to obtain the matrix, C, of size N × (T – 1). Similarity in the connectivity fingerprints between each pair of ROI voxels was calculated using the η² coefficient¹¹⁷, resulting in the matrix, S, of size N × N. The η² coefficient represents the fraction of variance in one connectivity profile accounted for by the variance in another. We then took the average S across all individuals. A nonlinear manifold learning procedure using the Laplacian Eigenmaps algorithm¹¹⁸ was applied to the S matrix of each ROI to calculate its eigenvectors and eigenvalues. The eigenvectors represent the functional patterns, termed functional gradients, of the ROI, ordered according to the variance in FC similarity they explain. We only analysed the first 20 non-constant gradients to enable direct comparisons with the geometric eigenmodes. Typical applications of FC-based gradient analysis rarely consider more than the first five gradients^49,119.

We analysed the correspondence between the geometric eigenmodes and functional gradients in each non-neocortical structure by taking their absolute spatial correlations, given that the signs of the modes and gradients are arbitrary. Moreover, the ordering of the geometric eigenmodes within the same eigengroup can change (Supplementary Table 1; for example, the first eigengroup can have order flips among modes 2–4) and hence one-to-one correspondence between the indices of geometric eigenmodes and functional gradients is not guaranteed. We therefore examined all possible pairwise correlations between the modes and gradients and evaluated correspondence with respect to the mode that maximally correlated with each gradient. The maximal order differences observed between geometric modes and functional gradients were one, eight and five out of 20 for the thalamus, striatum and hippocampus, respectively.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.