Functional harmonics reveal multi-dimensional basis functions underlying cortical organization.

The human brain consists of specialized areas that flexibly interact to form a multitude of functional networks. Complementary to this notion of modular organization, brain function has been shown to vary along a smooth continuum across the whole cortex. We demonstrate a mathematical framework that accounts for both of these perspectives: harmonic modes. We calculate the harmonic modes of the brain's functional connectivity graph, called "functional harmonics," revealing a multi-dimensional, frequency-ordered set of basis functions. Functional harmonics link characteristics of cortical organization across several spatial scales, capturing aspects of intra-areal organizational features (retinotopy, somatotopy), delineating brain areas, and explaining macroscopic functional networks as well as global cortical gradients. Furthermore, we show how the activity patterns elicited by seven different tasks are reconstructed from a very small subset of functional harmonics. Our results suggest that the principle of harmonicity, ubiquitous in nature, also underlies functional cortical organization in the human brain.

Introduction

The topographic organization of the brain into functionally specialized areas is one of its fundamental properties (Felleman and Van Essen, 1991), suggested to have been present in evolution as early as the last common ancestor of vertebrates (Krubitzer and Kaas, 2005; Eickhoff et al., 2018). The individuality of each brain area is determined by its functional specification, its microstructure (cyto- and myeloarchitecture) (Eickhoff et al., 2018), and its inter- and intra-areal connectivity (Glasser et al., 2016). Significant effort in neuroscience has been directed toward subdividing the brain into adjoining parcels based on functional characteristics and inter-areal connectivity (Eickhoff et al., 2018; Glasser et al., 2016). However, in parallel to this modular view of brain organization, where separate, adjoining brain areas with uniform functionality and homogeneous structural connectivity form distinct functional units, there is also extensive evidence that gradually varying boundaries between brain areas exist, suggesting a degree of transition (Bailey and Von Bonin, 1951) and context-dependence (Salehi et al., 2020) instead of sharply separated brain areas. Furthermore, topographic mappings including retinotopy (Sereno et al., 2001), somatotopy (Penfield and Rasmussen, 1950), and tonotopy (Perrone-Capano et al., 2017), show that representations of our visual field, body, and auditory frequency spectrum are spatially continuous across the areas of the primary visual, somatomotor, and auditory cortices, respectively, exemplifying an organizational principle that is complementary to the finding of relatively uniform functionality within some brain areas.

In order to allow for the immense complexity of human brain function (Tononi et al., 1994), a multitude of functionally distinct brain areas coordinate through synchronous fluctuations in their activity (Varela et al., 2001). Coherent oscillations among distinct brain areas have been shown to be another evolutionarily conserved aspect of brain activity (Vincent et al., 2007). However, more complex or elaborate mental processes are hypothesized to result from the convergence of information from sensory modalities onto association cortices (Mesulam, 1998). This convergence is assumed to increase with spatial distance on the cortex from the highly functionally specialized primary cortices (Buckner and Krienen, 2013). As a consequence, gradiental organization might be a general organizational principle throughout the cortex and not only in primary sensory areas. In line with this hypothesis, a principal connectivity gradient of cortical organization in the human connectome has been identified, where the functional networks of the human brain are located according to a functional spectrum from perception and action to more abstract cognitive functions (Margulies et al., 2016; Huntenburg et al., 2018). This suggests that in addition to coherent oscillations between brain regions (“modular view”), macroscopic networks can also be described in terms of “gradiental” organization. It is important to note that both organizational principles are complementary and not mutually exclusive (i.e., brain regions may differ in the degree to which their boundaries are sharp) (Tian and Zalesky, 2018; Bajada et al., 2020). However, the principles underlying the functional organization of macroscopic networks from specialized brain areas remain largely unknown.

Here, we propose that the human brain’s functional organization is governed by the natural principle of harmonic modes. In particular, we demonstrate how harmonic modes of the resting state functional connectivity of the human cortex can explain both gradiental organization and the presence of distinct functional areas (parcels). The principle of harmonic modes underlies a multitude of physical and biological phenomena including harmonic waves encountered in acoustics (Chladni, 1802), optics (Bedzyk et al., 1988), electron orbits (Schrödinger, 1926; Moon et al., 2008), electro-magnetism (Roos, 2012; Britton et al., 2012), and morphogenesis (Murray, 1988; Xu et al., 1983). The principle of harmonicity is also respected in the human brain across multiple scales, ranging from the ocular dominance patterns of the early visual areas (Swindale, 1980), visual hallucinations (Ermentrout and Cowan, 1979; Bressloff et al., 2002; Billock and Tsou, 2007; Rule et al., 2011), to the organization of cortical and thalamic tissues (Wilson and Cowan, 1973; Atasoy et al., 2016). On the macroscopic scale, harmonic modes of the circle (Nunez and Srinivasan, 2006) and of the sphere have been proposed to underlie cortical communication observed in magneto-/electroencephalography (M/EEG) (Nunez and Srinivasan, 2006; Tokariev et al., 2019) and in functional magnetic resonance imaging (fMRI) (Robinson et al., 2016).

In the graph domain, it was shown that harmonic modes of the structural connectome can explain functional connectivity, in particular, resting state networks (Atasoy et al., 2016). More generally, harmonic modes of the structural connectome are useful for our understanding of how functional activity is variably shaped by underlying white matter connectivity (Preti and Van De Ville, 2019; Glomb et al., 2020). Moreover, harmonic modes of the structural connectivity have been found to predict disease progression in dementia (Raj et al., 2012). Harmonic modes of the “functional” connectome have been studied for topographic organization by Haak et al. (2013, 2018). By using predefined regions of interest, they revealed up to two principal fine-grained gradients along the visual and somatomotor hierarchy. In the seminal work by Margulies et al. (2016), the authors provided a detailed analysis and interpretation of the first two gradients of cortical functional connectivity. In a recent publication, Tian et al. (2020) applied the procedure described in Haak et al. (2018) in order to obtain a multi-scale parcellation of the subcortex by considering the first three gradients of its functional connectivity. The link between structural and functional connectivity is made by dynamical models, which suggest that space and time are linked via the oscillatory frequencies of certain brain networks (Atasoy et al., 2018). Further evidence for a link between spatial patterns and oscillations comes from applications of harmonic modes of the structural connectivity to faster timescales (i.e., M/EEG) (Glomb et al., 2020; Tokariev et al., 2019; Raj et al., 2020).

In this work, we uncover the spatial shapes of the first 11 harmonic modes that result from synchronous hemodynamic fluctuations in large-scale brain activity as measured in fMRI, by solving the time-independent (standing) wave equation (Levy, 2006; Atasoy et al., 2016) on the functional connectivity (FC) of the human brain. These harmonic modes, called “functional harmonics,” decompose the FC of the human brain into a hierarchical set of (graph) frequency-specific smooth activity patterns, or gradients. We demonstrate how functional harmonics capture spatial organization of the cortex on several spatial scales, spanning from within-area topographic mappings to combinations of macroscopic networks. Thereby, the functional harmonics unveil both, the principal connectivity gradient (Margulies et al., 2016), as well as cortical parcellations (Glasser et al., 2016), while also exploring, for the first time, higher-order gradients revealed by the harmonic decomposition of the dense FC of the human cortex.

Results

Estimation of functional harmonics

We computed functional harmonics of the dense functional connectivity matrix (dense FC) of the Human Connectome Project (HCP) (Glasser et al., 2013; Van Essen et al., 2012, 2013; Moeller et al., 2010; Feinberg et al., 2010; Setsompop et al., 2012; Xu et al., 2012; Jenkinson et al., 2002). Mathematically, the patterns of harmonic modes of a dynamical system are estimated by the eigendecomposition of the Laplace operator. In the present case, we are interested in the Laplace operator of the dense FC. This dense FC is estimated from the pairwise temporal correlations between all pairs of vertices on the cortical surface (V = 59.412 vertices in total) (Figures 1A–1C). Correlations can be interpreted as measuring the functional similarity of two cortical regions. Thus, the dense FC can be understood as encoding functional integration and segregation in the cortex.

Figure 1

Workflow for the estimation of functional harmonics

(A) Brain activity measured with functional magnetic resonance imaging (fMRI) in resting state for 812 subjects provided by the Human Connectome Project (HCP; 900 subjects data release).

(B) Illustration of brain activity time series of three representative vertices on the cortex (x1, x2,…, x).

(C) The dense functional connectivity (FC) matrix computed from the temporal correlations between the time courses of each pair of vertices as shown in (B) averaged across 812 subjects.

(D) Representation of the dense FC as a graph, where the edges indicate strong correlations between the corresponding vertices. The anatomical locations of the vertices are color-coded (Glasser et al., 2016).

(E) Functional harmonics are estimated by the eigenvectors of the graph Laplacian computed on the graph representation of the FC. The first five functional harmonics ordered from the lowest to higher spatial frequencies are illustrated on the FC graph representation (top), in the eigenvector format as 59,412 × one-dimensional vectors (middle), and on the cortical surface (bottom). Note that here we show the patterns on the left hemisphere for illustrative purposes, yet the entire cortex was used in the analysis. Likewise, the graph representations in (D) and (E) are shown for a parcellated version of the FC matrix using the HCP parcellation (Glasser et al., 2016), i.e., each node represents an HCP parcel, but the computation of the functional harmonics were performed on the dense FC using 59,412 × 59,412 without any parcellation.

We utilized the discrete counterpart of the harmonic modes that are defined on the graph built from the dense FC (i.e., the eigenvectors of the graph Laplacian) (Figures 1D and 1E). This is a suitable approach because brain connectivity is commonly conceptualized as a graph in which each node is a location in the brain and the edges are given by the connectivity matrix; this is illustrated in Figure 1D (but note that although in the figure, for illustrative purposes, each node corresponds to a brain region, in the following analyses, the nodes of the graph are formed by the full set of vertices on the cortical surface). In our case, in order to build the graph, we used a binary adjacency matrix derived from the dense FC matrix: connections were set between each vertex i and the k = 300 vertices j ≠ i with the largest correlations to vertex i, such that the entries a, j ≠ i were set to 1 where connections between vertices exist and to 0 otherwise. The graph Laplacian is defined as , where is a diagonal matrix that holds the degree (number of connections) of each vertex i in its diagonal entries d.

By definition, the first functional harmonic with eigenvalue 0 is constant over the whole cortex and is discarded. Figure 2 shows the first 11 non-constant functional harmonics (referred to as ψ1, ψ2, ⋯, ψ11), ordered starting from the lowest eigenvalue >0, and illustrating that each harmonic is a smoothly varying pattern on the cortex between a positive and a negative polarity (i.e., a gradient). To give an intuitive interpretation of functional harmonics, we emphasize that the actual magnitudes of the gradients plotted on the cortex are not per se meaningful, but the difference between the values assigned to two vertices reflects how different they are in terms of their “functional connectivity profile” (i.e., their pattern of connectivity to the rest of the cortex).

Figure 2

Functional harmonics capture existing characterizations of functional anatomy

The first 11 non-constant functional harmonics plotted on the cortical surface. It is clearly visible that the first two functional harmonics (A and B) constitute global gradients over the entire cortex, whereas subsequent maps (C–K) include increasingly more local details. In each functional harmonic, known functional regions (e.g., C), processing streams (e.g., E) or networks (e.g., B) are discernible, and we have annotated the most conspicuous ones. In order to illustrate that similarly colored patches of cortex correspond to known functional regions, borders of HCP parcels have been added (white lines). V1–V4, visual areas 1 to 4; MT, middle temporal visual area; 24 dd, an area that contains a higher order representation of the hand; fusiform face complex, an area that responds specifically to images of human faces. The functional harmonics were derived from the HCP’s dense functional connectivity matrix, which is an average of over 812 subjects.

With increasing eigenvalue, the functional harmonics become increasingly more complex and segregate the cortex into an increasing number of nodal areas (Levy, 2006) (contiguous areas of the cortex with similar colors in the surface plots in Figure 2). Mathematically, as the eigenvalue increases, (spatial) graph frequency also increases (Shuman et al., 2013). This intrinsic link between the graph frequency and cortical scale implies that functional harmonics yield not only a multi-dimensional, but also a multiscale description of the cortex (Tian et al. 2020). Note that the ordering by graph frequency is a property that emerges from the Laplace eigenfunctions and therefore is not present in function bases proposed by other methods, such as principle component analysis (PCA) or independent component analysis (ICA), which, as a result, do not implicitly possess this multiscale property.

Functional harmonics capture sub-areal topographic organization

We first tested whether functional harmonics capture cortical organization on a sub-areal scale (i.e., within a cortical parcel). Specifically, we investigated whether functional harmonics capture somatotopy (Penfield and Rasmussen, 1950) and retinotopy (Sereno et al., 2001), two major topographic mappings found in the brain. Topographic mappings represent sensory input on the cortical surface such that the relative positions of the receptors that receive these inputs (e.g., the photo receptors in the retina) are preserved. In a broader sense, it refers to any mapping in which neurons that are close together on the cortical surface in one functional area project to neurons that are close together in another functional area. This organization leads to a functional gradient within a specialized brain area (Haak et al., 2013).

In the HCP parcellation (Glasser et al., 2016), the sensorimotor cortex is, in each hemisphere, divided into five areas (1, 2, 3a, 3b, and 4) that are defined cytoarchitechtonically and using microstructural measures such as cortical thickness and myelination (see “Supplementary neuroanatomical results” in Glasser et al. [2016]). However, functional connectivity is not homogeneous within these regions: five somatotopic sub-areas are defined by the HCP (Glasser et al., 2016) and form a topographic map of the surface of the body on the cortex, which is oriented orthogonally to the anatomically defined areal boundaries. These sub-areas correspond to representations of the face, upper limbs (mapped by moving the hands), eyes, lower limbs (mapped by moving the feet), and trunk. We observed somatotopic mappings within functional harmonics 3 (ψ3), 7 (ψ7), and 11 (ψ11) (Figures 2C, 2G, and 2K). Figure 3A illustrates the two-dimensional subspace formed by functional harmonics 3 (ψ3) and 11 (ψ11), which accounts for the mapping of the human body onto the somatotopic regions of the cortex. This example also illustrates how specialized functional regions can be explained by the interaction of functional harmonics across multiple dimensions within the functional harmonics framework instead of a single gradient (see Figures S1A–S1C for further examples).

Figure 3

Functional harmonics capture somatotopy and retinotopy

(A) Functional harmonics 3 (ψ3) and 11 (ψ11) in their own space. Vertices are color-coded according to their anatomical locations (see Figure 1D), and the location of 4 somatotopic areas in this space is annotated.

(B and C) Retinotopies of functional harmonics 4 (B, ψ4) and 8 (C, ψ8). Each panel shows, on the right, the colors of the respective functional harmonic in early visual areas V1–V4 on a polar plot of eccentricity (distance in degree from the fovea) and angle on the visual field (see legend at the bottom of the figure). On the left, the respective functional harmonic is shown on a flat map of early visual cortex (left hemisphere). V1, V2, V3, and V4: visual areas 1, 2, 3, and 4. The shown figures are derived from the functional harmonics obtained from the HCP’s dense functional connectivity matrix, which is an average over 812 subjects.

We derived a measure that we term “somatotopic silhouette value,” which quantifies the degree to which each somatotopic sub-area is delineated within these three functional harmonics (see STAR Methods for details). We compared our results to those obtained from spherical rotations (Alexander-Bloch et al., 2018) of functional harmonics 3 (ψ3), 7 (ψ7), and 11 (ψ11). We found that for functional harmonic 3 (ψ3), the face and foot areas were significantly separated from the rest of the cortex as well as other somatotopic areas; for functional harmonic 7 (ψ7), we found the face and hand areas, and for functional harmonic 11 (ψ11), we found the hand areas to be significantly separated (see Figure S2A; pcorr < 0.05 after Bonferroni correction, Monte Carlo tests with 300 rotated versions of the functional harmonics). This finding indicates that functional harmonics capture somatotopic organization in the cortex.

We next investigated the presence of retinotopic mapping of early visual regions (V1–V4), where cortical representations of the visual field reflect the positions of the receptors in the retina of the eye such that each vertex within the patterns of functional harmonics is assigned an eccentricity (distance from the fovea) and a polar angle (position in the visual field, i.e., top, bottom, left, right), according to the HCP retinotopy dataset (Benson et al., 2018). It has previously been demonstrated that topographic organization can be evaluated using functional connectivity gradients, either within the visual system (Glasser et al., 2016; Haak et al., 2018) or on a whole-brain level (Yeo et al., 2011). Examples of polar plots of the retinotopic gradients are shown in Figures 3B and 3C (all polar plots are shown in Figure S2B). To investigate the degree of agreement between functional harmonics 1–11 (Figure 2) and the retinotopic mappings, we measured the correlation between eccentricity as well as polar angle maps and functional harmonic patterns in V1–V4. We found significant correlations (pcorr < 0.05 after Bonferroni correction) between the retinotopic eccentricity map and all functional harmonics 1–11 except functional harmonic 9 and also between the retinotopic angular map and functional harmonics 1–4, 7–9, and 11.

These results, obtained with whole-brain connectivity, demonstrate that retinotopic organization of the early visual areas is implicitly present in the resting state brain activity (Yeo et al., 2011) and is revealed by the functional harmonic basis. We draw the reader’s attention to the fact that retinotopic organization has been shown in rich detail elsewhere, using either retinotopic mapping techniques based on specific visual stimuli (Benson et al., 2018; Sereno et al., 2001) or functional connectivity within the visual system (Glasser et al., 2016; Haak et al., 2018). Thus, although functional harmonics capture some prominent aspects of retinotopic organization, there are also details that require more dedicated approaches.

It is also important to note that the borders of visual areas V1–V4, which were identified by the existence of orderly topographic maps with clear visual field reversals (Sereno et al., 1995), correspond to a sign reversal between the positive and negative polarity of the harmonic pattern in various functional harmonics (e.g., ψ5–ψ10). Yet in functional harmonics, the delineation of the visual areas occurs in a hierarchical fashion (i.e., although the sign changes in ψ1 and ψ2 follow the borders of the visual cortex as a whole, higher order functional harmonics respect various borders occurring between V1–V4 within the visual cortex).

Functional harmonics reveal specialized brain areas

By visual inspection of Figure 2, one notices that patches of cortex that are colored in the same shade seem to correspond to a certain degree to known specialized brain regions (parcels) delineated by the HCP parcellation (Glasser et al., 2016) or groups of such regions (all parcel borders are shown in Data S1, page 1). This suggests that the gradients described by functional harmonics are flat within many specialized cortical areas, indicating a near-homogeneous functional connectivity profile. In order to quantify the homogeneity within cortical areas of the functional harmonics shown in Figure 2, we compared the within-area-variability to the average between-area-variability of each parcel. The resulting measure is similar to the silhouette value used in quantifying the quality of cluster solutions (de Amorim and Hennig, 2015), but ranges between 0 and 1, where a value close to 1 indicates that the functional harmonic is homogeneous within that region. We refer to this measure as “modified silhouette value” (see STAR Methods for details).

We applied this analysis to four alternative function bases: (1) eigenvectors of the dense FC matrix (Data S1, page 2), (2) eigenvectors of the adjacency matrix (Data S1, page 3), (3) principal components (PCA) (Data S1, page 4), and (4) independent components (ICA) (Data S1, page 55). The first two of these bases were chosen in order to test the effect of the two major processing steps that are necessary when transforming the dense FC matrix into the Laplacian: (1) the adjacency matrix is obtained from the dense FC matrix by binarizing (using k = 300 nearest neighbors), and (2) the graph Laplacian is obtained from the adjacency matrix by normalizing, i.e., taking into consideration the degree of each vertex (). Furthermore, to relate the performance of functional harmonics to other well-known function bases, we also performed the analysis with the basis functions of (3) PCA and (4) ICA. As shown in Figure 4A, although all approaches exhibited moderate to high modified silhouette values, indicating that they all reflected the functional organization well, we found that functional harmonics had significantly higher modified silhouette values than three of the four function bases; only eigenvectors of the adjacency matrix had equally high modified silhouette values.

Figure 4

Functional harmonics capture specialized brain areas

The modified silhouette values (y-axes in all panels) quantifies the degree to which gradients described by the functional harmonics as well as control basis sets are flat within the parcels defined by the HCP parcellation. A modified silhouette value close to 1 indicates homogeneous values within HCP parcels.

(A) A comparison between the basis sets indicates that functional harmonics and adjacency eigenvectors have significantly higher modified silhouette values than the other three (Wilcoxon rank-sum test, black bars indicate significant differences at pcorr < 0.05 after Bonferroni correction). Each data point is computed from a matrix that is an average over 812 subjects.

(B–F) Modified silhouette values of the first 11 non-constant components of each basis set (colored circles; each data point is computed from a matrix which is an average over 812 subjects) compared to their rotations (gray crosses). The stars above each column indicate significant silhouette values (pcorr < 0.05 after Bonferroni correction, Monte Carlo test with 220 spherical rotations per component). Functional harmonics had the highest number of significant modified silhouette values (10 out of 11).

We established significance of the modified silhouette values obtained from each functional harmonic by computing the same measure also for n = 220 spherical rotations of the functional harmonics (Alexander-Bloch et al., 2018) (stars above each column in Figures 4B–4F). By doing so, we ensured that homogeneous gradients within parcels are not due to the properties of the harmonic decomposition alone, but are specific to their location on the cortical surface. As shown in Figure 4B, we found that the modified silhouette values of 10 of the first 11 functional harmonics were larger than those of the rotated harmonic basis (pcon < 0.05 after Bonferroni correction, Monte Carlo tests; see STAR Methods for details). The only exception to this finding was functional harmonic 4 (ψ4), which captures the retinotopic organization of early visual regions (Figure 3B; see above for a discussion of retinotopic organization of functional harmonics), and which is relatively flat over the rest of the cortex. When applied to the four alternative function bases, the same analysis showed that for eigenvectors of the dense FC matrix, only 7 out of 11 basis functions had a significantly higher modified silhouette value than their spherical rotations; for eigenvectors of the adjacency matrix, this number was 8 out of 11, for PCA, only 3 PCs reached significance, and for ICA, this number was 4 out of 11 (Figures 4C–4F). For qualitative evaluation, the overlap between parcels and functional harmonics as well as other bases are shown in Data S1, pages 1–5.

Functional harmonics yield graph frequency-specific brain networks

In Margulies et al. (2016), the authors investigated the principal cortical gradients captured by the first two functional harmonics. Although it should always be considered that harmonic modes describe the structure of the underlying graph in a multi-dimensional fashion (as is illustrated for two dimensions in Figure 3A), in the following, we provide some insight into the functional significance of each of the functional harmonics shown in Figure 2. Functional harmonics 1 (ψ1) and 2 (ψ2) correspond to previously identified large-scale gradients (Margulies et al., 2016) that delineate the separation between the major sensory and the uni- versus multimodal cortices in the brain, respectively (see Figure S1D). Figures 2A and 2B demonstrate the overlap between the visual and sensorimotor networks as defined in Yeo et al. (2011) and the gradiental patterns of the first and second functional harmonics. We observed that functional harmonic 3 (ψ3) reveals a finer subdivision of the somatosensory/motor system (Zeharia et al., 2012, 2015; Kuehn et al., 2017). The overlay of the borders of the five somatotopic areas defined by the HCP (Glasser et al., 2016; Barch et al., 2013) on the third functional harmonic are shown in Figure 2C. Similarly, in functional harmonic 4 (ψ4), we found a finer segregation of the visual system, following a retinotopic eccentricity gradient (Benson et al., 2018). The overlay of the borders of early visual areas (V1–V4) on functional harmonic 4 (ψ4) are shown in Figure 2D (for further details on retinotopic and somatotopic mapping, see above).

Qualitative evaluation of higher frequency functional harmonics systematically revealed their link to more specialized complex brain function. The pattern observed in functional harmonic 5 (ψ5) (Figure 2E) is consistent with a functional network in which action and perception interact. In the negative polarity, we found primary visual, auditory, and somatosensory cortices, whereas the regions in the positive polarity closely resemble the sensory-motor pathway, which has been shown to mediate selective interactions between resting state networks along the visual hierarchy (Yeo et al., 2011). Parts of this pathway are known to be modulated by visuospatial attention (Corbetta and Shulman, 2002). The fact that we also found auditory and somatosensory regions is in line with the idea that the interplay between action and perception circuits—also known as active inference (Friston et al., 2011; Adams et al., 2013)—is a multimodal process (Keysers et al., 2010; Hauk and Pulvermüller, 2004; Pulvermüller and Fadiga, 2010). In functional harmonic 6 (ψ6), auditory and visual areas were both localized in the positive polarity, forming a network related to audiovisual object (including faces) recognition (Beauchamp et al., 2004; Levy et al., 2001; Hasson et al., 2003) (i.e., recognition of the “outer world”). The negative polarity of functional harmonic 6 (ψ6) segregates the somatotopic face area as well as parts of the default mode network (DMN), a network of regions whose activity has been related to self-referential tasks (Gusnard et al., 2001). Thus, the negative polarity of functional harmonic 6 (ψ6) forms a self-referential processing stream (Gusnard et al., 2001; Haxby et al., 2000; Leslie et al., 2004). Functional harmonic 7 (ψ7) provides a finer somatotopic gradient, including a higher hand area, 24 dd, in the medial cortex (Zeharia et al., 2015) (see Figure 2G and annotations in Figure 2C). Functional harmonics 8–10 (ψ8, ψ9, and ψ10) correspond to different subdivisions of higher order networks such as the frontoparietal network and DMN (see Figures S3A–S3C). In particular, the DMN (Raichle et al., 2001) is delineated in the positive polarity of functional harmonic 9 (ψ9) (borders of the DMN as defined by Yeo et al. [2011] are overlaid on functional harmonic 9 (ψ9) in Figure 2I). In functional harmonic 10 (ψ10), we found a significant correlation (r = −0.63, p = 4 ⋅ 10−21) with the degree of auditory involvement of the functional areas (Figure S3D). Functional harmonic 11 (ψ11), the first clearly asymmetric harmonic between the two hemispheres, yields the separation between the right and left somatotopic hand areas (Pool et al., 2014).

Overall, these results demonstrate that functional harmonics provide a multitude of functionally relevant macroscopic networks, each associated with a unique graph frequency.

Relating rest and task with functional harmonics

Functional harmonics, by definition, yield the extension of the Fourier basis to the functional connectivity of the human brain. Laplace eigenfunctions on a one-dimensional domain with cyclic boundary conditions (i.e., a circle) yield sine and cosine waves with different frequencies, which constitute the well-known Fourier basis. In a similar manner, we estimated the functional harmonics by computing the Laplace eigenfunctions on a graph encoding the human brain’s functional connectivity (illustrated in Figure 1D). As such, the functional harmonics provide a frequency-specific function basis, in which any pattern of brain activity on the graph can be represented as a weighted combination of functional harmonics. Given the experimental evidence showing that resting state functional connectivity reflects connectivity during task (Biswal et al., 1995; Deco et al., 2011; Buckner et al., 2013; Cole et al., 2016), we tested how efficiently the functional harmonics derived from the resting state dense FC can represent task-induced cortical activity.

To this end, we reconstructed 47 group-level task maps provided by the HCP (Barch et al., 2013) from the superposition of functional harmonics (see STAR Methods for details). The 47 maps consist of activation maps as well as contrasts derived from 7 groups of tasks (working memory, motor, gambling, language, social, emotional, and relational; see STAR Methods for summaries). This reconstruction yields a coefficient (weight) for each functional harmonic, quantifying how much it contributes to a certain task map. The set of all coefficients forms a spectrum equivalent to the power spectrum obtained from a Fourier transform, in this case the power spectrum of the functional harmonic basis. We quantified the goodness of fit by measuring the distance between the original and the reconstructed task maps.

We determined how well the first 11 non-constant functional harmonics shown in Figure 2 were able to approximate task maps and compared their performance to the four alternative function bases mentioned above (eigenvectors of the FC, eigenvectors of the adjacency matrix, principal components, and independent components), as well as spherical rotations of functional harmonics (Alexander-Bloch et al., 2018). Note that the first 11 functional harmonics constitute ∼0.02% of the total functional harmonic spectrum. Due to the intrinsic ordering of functional harmonics by graph frequency, using only the first few functional harmonics omits more localized details of brain activity and thus, the omitted information results in a reconstruction error.

Figures 5A–5G show the average normalized reconstruction errors for all groups of tasks and for all compared function bases. First focusing on the functional harmonic basis (red line), the error drops from between 1.00 (for emotion) and 1.40 (for language) to between 0.65 (for emotion) and 0.78 (for language). This corresponds to a level of correlation between the original task maps and the task maps reconstructed from the first 11 non-constant functional harmonics of between 0.78 (for emotion) and 0.69 (for language; see Figure S4A). Figure 5H illustrates the reconstruction procedure for one specific task (working memory: body) (see Data S1, page 6 for all tasks).

Figure 5

A small set of functional harmonics suffices to reconstruct diverse task activity maps

(A–G) Mean reconstruction errors for each of the 7 task groups and all 6 basis function sets when only the first 11 non-constant components are used (see Figure S4 for results when using more components).

(H) One example for a reconstruction using a working memory task. The bottom panel is the original task activation map (working memory—body; see also Data S1, page 6, panel l), and top panels use the number of harmonics indicated on the left to reconstruct it.

(I) Results of significance tests comparing normalized reconstruction errors of functional harmonics to other function bases (significant differences at pcorr < 0.05 after Bonferroni correction, Monte Carlo tests with 1,000 permutations). Top: When using only the first 11 components, functional harmonics outperformed each of the other function bases. Bottom: When using the first 100 components (Figure S4), eigenvectors of the dense FC outperformed functional harmonics in 3 task groups. All basis sets were derived from matrices which are averages over the same 812 subjects. The task activity maps (Cohen’s D activation contrast maps) are based on 997 subjects.

We next compared the reconstruction performance of the functional harmonics to each alternative function basis, employing a Monte-Carlo analysis. Again, using only the first 11 non-constant components, we found that reconstruction errors of functional harmonics were significantly lower than those of their rotations for each of the task groups (all pcorr < 0.035, Monte-Carlo tests with 1,000 permutations, Bonferroni corrected for multiple comparisons), and significantly lower than those of the adjacency eigenvectors in six out of seven task groups (all pcorr < 0.035, Monte-Carlo tests with 1,000 permutation, Bonferroni corrected for multiple comparisons, except language, where p = 0.18 before correction for multiple comparisons, n.s.). In comparison to FC eigenvectors, we found that functional harmonics performed significantly better in the reconstruction of motor tasks (pcorr < 0.035, Monte-Carlo tests with 1,000 permutations, Bonferroni corrected for multiple comparisons), whereas there was no significant difference in other task groups (all p > 0.01 before correction for multiple comparisons, n.s.). Compared to PCA and ICA, the reconstruction errors of functional harmonics were significantly lower for motor and working memory task groups (all pcorr < 0.035, Monte-Carlo tests with 1,000 permutations, Bonferroni corrected for multiple comparisons), whereas for all other task groups, there were no significant differences (all p > 0.01 before correction for multiple comparisons, n.s.). We also tested the performance of the alternative function bases when the first 100 components were used (Data S1, page 6). We found several changes in contrast to the results with only the first 11 components. First, FC eigenvectors outperformed functional harmonics in working memory, language, and social tasks, whereas there was no significant difference in the other task groups. Second, functional harmonics outperformed PCA in all seven task groups instead of only in motor and working memory. Third, for ICA, only 15 ICs were extracted, and when all of them were used, functional harmonics outperformed ICA in social and emotional tasks in addition to working memory and motor tasks as before. The results remained the same for adjacency eigenvectors.

Figure 5I provides a summary of the results. The first 11 non-constant functional harmonics performed at least as well as the tested alternative function bases in all task groups (i.e., none of the alternative function bases outperformed functional harmonics in any of the task groups). When 100 components were used, FC eigenvectors outperformed functional harmonics in three task groups.

Next, we wanted to confirm that the decrease in reconstruction error did not merely reflect the fact that functional harmonics capture global features of brain activity maps common to all tasks. We designed an analysis that allowed us to test whether functional harmonics are able to capture brain activity that is specific to a task or group of tasks. Figures 6A, 6D, and 6G and Figures 6B, 6E, and 6H show three examples of task activation maps and the corresponding normalized power of the first 11 non-constant functional harmonics, respectively, revealing how strongly each of the 11 functional harmonics shown in Figure 2 contributes to these particular task maps. In other words, Figures 6B, 6E, and 6H show spectral representations of the task maps in Figures 6A, 6D, and 6G. For qualitative evaluation, we display the task activation maps reconstructed by superimposing functional harmonics in the order of their contribution strength for varying numbers of functional harmonics in Figures 6C and 6F (see Data S1, page 6 for all tasks). Across all 47 task maps that were evaluated, the functional harmonic that was the strongest contributor was always either the constant functional harmonic or one of the first 11 non-constant harmonics shown in Figure 2.

Figure 6

Functional harmonics provide a characterization of task activity maps

(A) Original task map of the contrast between working memory (face) and average working memory from the HCP task dataset (Barch et al., 2013).

(B) Spectral representation of the task map shown in (A) (i.e., the normalized coefficients of the graph Fourier transform quantify the contribution of the first 11 non-constant functional harmonics to the task map). The color indicates the task group (see legend in L).

(C) Reconstruction of the task map in (A) when using the functional harmonic with the strongest contribution (highest coefficient) only, the four functional harmonics with the strongest contributions, and the forty functional harmonics with the strongest contributions.

(D–F) The same as (A)–(C) using the map of the contrast between motor (right hand) and average motor.

(G–I) The same as (A)–(C) using the map of the contrast between motor (trunk) and average motor.

(J–L) Confusion matrices. Black entries mark the task map-reconstruction-pair that has the lowest reconstruction error; colored squares indicate the task group.

(M) Proportion of reconstructions, for each number of harmonics, which have the minimum reconstruction error with their exact original task map (thick line) and a task map belonging to the same group of tasks as the original map (thin line). Functional harmonics were derived from the HCP’s dense functional connectivity matrix, which is an average of over 812 subjects. The task activity maps (Cohen’s D activation contrast maps) are based on 997 subjects.

In order to confirm that each spectral representation was specific to the task map from which it was obtained, we computed the distance between a given reconstructed map and all original task maps, resulting in a confusion matrix for each number of harmonics with maximum contribution. If a spectral representation is indeed specific to a task map, the error should be minimal between a reconstruction and its corresponding task map compared to the error of the reconstruction of the other 46 task maps. The confusion matrices in Figures 6J–6L show the pairs of the original and reconstructed task activation maps with the minimum distance when using 1, 4, and 40 functional harmonics with maximum contribution. Colored squares mark the 7 task groups as in Figure 5. The proportion of unambiguously identified tasks depending on the number of functional harmonics is shown in Figure 6M. We found that sparse representations using the 4 functional harmonics with the largest power for each task are sufficient to unambiguously characterize the seven task groups with the exception of one working memory task (Figure 6K), and 70% of all individual tasks. When the 40 functional harmonics with maximum contribution were used, which corresponds to 0.1% of the complete spectrum of functional harmonics, 44 out of 47 task maps were correctly identified from their reconstructions (Figure 6L).

These results demonstrate that functional harmonics provide a functionally relevant representation, where the brain activity accompanying different tasks can be unambiguously identified from the activation profiles of a small range of functional harmonics.

Discussion

In this paper, we show that harmonic modes of the vertex-level functional connectivity of the human cortex, termed functional harmonics, link observations spanning several spatial scales, from sub-areal topographic mappings (somatotopy, retinotopy), to delineation of specialized functional brain areas, to the level of large-scale networks and combinations thereof. In this view, functional integration and segregation on all these levels can be accounted for by the same mathematical framework (i.e., harmonic modes).

Although it is abundantly clear from the literature that modular and gradiental organization are complementary principles in the human brain, a mathematical framework that can account for both has thus far been missing. Within the functional harmonic framework, specialized regions as well as networks of functionally related regions correspond to patches of the cortical surface within which the gradiental patterns of functional harmonics are relatively constant. Such constant regions are then increasingly segregated into finer regions through the interaction of functional harmonics across multiple dimensions, as was recently shown to be the case for the subcortex by Tian et al. (2020).

We have interpreted the first 11 functional harmonics in relation to existing characterizations of functional anatomy, like the concordance between functional harmonic 5 and the sensory-motor pathway (Yeo et al., 2011). Thus, in this framework, functional networks observed during rest and task emerge as a graph frequency-specific function basis derived from the human resting state functional connectivity matrix.

Harmonic modes are initially a dimensionality reduction technique: they are estimated as the eigenvectors of the graph Laplacian. As such, they are ordered by graph frequency and orthogonal to one another. By definition, any pattern of cortical activity can be expressed in this function basis as a superposition of functional harmonics. In this new function basis, brain activation patterns measured during tasks are expressed in a convenient and compact manner as task-specific spectra that quantify the contribution of each functional harmonic.

We tested how well functional harmonics are able to capture task activity patterns (reconstruction errors, as shown in Figure 5), as well as specialized brain regions (modified silhouette values, as shown in Figure 4), and compared them to other function bases on these two spatial scales. In both regards, functional harmonics outperformed PCA and ICA, two very popular dimensionality reduction techniques, suggesting that the non-linear nature of harmonic modes captures important properties of cortical organization. The remaining two function bases were chosen because they constitute steps in our processing pipeline. Eigenvectors of the dense FC performed equally well or better than functional harmonics in task reconstruction but performed significantly worse in capturing brain regions. This is interesting because it suggests that although details of the FC that are not preserved in the k = 300 nearest neighbor graph shape task activity, the same does not seem to be true for parcel boundaries. In fact, with regard to the parcellation, eigenvectors of the adjacency matrix performed almost as well as functional harmonics, and significantly better than FC eigenvectors, suggesting that keeping only the nearest neighbors is actually beneficial for delineating parcel borders. At the same time, eigenvectors of the adjacency matrix were outperformed by functional harmonics in task reconstruction. It should be noted, however, that the differences between functional harmonics, FC eigenvectors, and adjacency eigenvectors were slight in reconstruction performance. We speculate that the reason for these contrasting results could be that the parcellation captures average functional regions, whereas during task, parcel borders can shift and adapt (Salehi et al., 2020), consistent with the idea that transitions between parcels may be gradiental (Bailey and Von Bonin, 1951; Bajada et al., 2020).

The approach of functional harmonics allows interpreting cortical activity patterns as different combinations of the same functional networks that are activated in order to fulfill task demands. In addition, the continuous nature of functional harmonics enables them to combine brain regions into macroscopic networks in a manner that allows the shape of the regions to vary depending on context (Salehi et al., 2020). Future research should evaluate whether the exact shapes of (combinations of) functional harmonics is in line with activation studies using different tasks.

Beyond dimensionality reduction, the interpretability of functional harmonics demonstrated in this paper, their multi-scale properties, and their orthogonality make them a candidate organizational principle that could explain how the brain flexibly switches between overlapping functional networks: functional harmonics are, by definition (Belkin and Niyogi, 2003), the patterns with the least overall variation on the cortex that respect the constraints posed by the functional relationships given by the FC. This implies that the average difference between neighboring nodes in a graph representation is minimized. Intriguingly, theoretical work has shown that activation patterns on graphs in which neighboring nodes co-activate lead to patterns with minimum free energy or entropy (Tomasi et al., 2013; Gu et al., 2018; Friston, 2010), and the transition between such patterns requires minimal energy (Gu et al., 2015). This means that transitioning between the functional networks instantiated by the functional harmonics in response to changing task demands is optimally efficient. In addition, functional harmonics and the networks to which they correspond emerge in a wholly self-organizing fashion from the functional connectivity.

Another important conceptual aspect of harmonic modes is their grounding in physical theories linking them to standing waves. This suggests that each spatial pattern on the graph might be linked to a unique set of temporal frequencies. Such a link would make it possible to interpret functional harmonic modes as “communication channels” along the lines of the ideas of multiplexing (Akam and Kullmann, 2014). In Haak et al. (2018), the authors demonstrated that harmonic modes applied to the visual system were able to “tease apart” retinotopic gradients in separate harmonic modes. In parallel to this idea, we assert that functional harmonics are able to “tease apart” overlapping functional networks subserved by specialized brain regions with multiple functions. For example, area MT, annotated in Figures 2E and 2F, is part of the sensory-motor pathway (Figure 2E) as well as audiovisual processing (Figure 2F). The idea of “communication channels” suggests that depending on which function area MT is fulfilling at any given moment, the network that it is engaged in would synchronize using different temporal frequencies.

Taken together, we conclude that although functional harmonics perform very well in both capturing the HCP parcellation and task map reconstruction, they offer specific benefits in terms of interpretability. This makes them an excellent choice not only as a tool for dimensionality reduction, but for exploring cortical organization.

Despite the relative simplicity of the approach presented here as compared to, e.g., ICA, there are several technical issues that need to be addressed in the future in order to ensure robustness of the eigenmaps discussed in this work. In particular, the choice of how to set the weight in the adjacency matrix, the choice of Laplacian, and the choice of k (number of nearest neighbors) are free parameters, the impact of which needs to be systematically explored in the future. Here, we obtained highly interpretable cortical patterns using the technique of Laplacian eigenmaps (Belkin and Niyogi, 2003) with a binary adjacency matrix with k = 300 nearest neighbors. Although this can be seen as a standard approach, connection weights in the adjacency matrix could be set in a graph distance-dependent manner, or the original correlation values could be used (Belkin and Niyogi, 2003); instead of the combinatorial Laplacian, one could use the symmetric normalized Laplacian. In this study, we did not find any significant difference when using the normalized Lapalcian. The choice of k was previously discussed regarding a related technique, the Isomap algorithm (Tenenbaum et al., 2000; Balasubramanian et al., 2002): k should be chosen such that a balance between local and global features is achieved. The amount of local and global connections has been suggested to have a major impact on the performance of the harmonic decomposition (Haak et al., 2018; Naze et al., 2021). At the same time, k might vary in a subject- or state-specific manner, revealing further interesting properties of cortical organization.

It should be noted that our approach differs from that used in Margulies et al. (2016) and Haak et al. (2018), in that it uses an adjacency matrix derived directly from the FC instead of quantifying the similarity of connectivity profiles in an intermediate step. The fact that the first two functional harmonics look virtually identical to what was presented in Margulies et al. (2016) speaks to the robustness of these principal cortical gradients toward the specific technique (they also used a slightly different definition of the Laplacian, and k was around 6,000, i.e., 10% of the connections were retained).

Considering that the principle of harmonic modes, when applied to the structural connectivity of the human brain—the human connectome—has been shown to reveal functional networks (Atasoy et al., 2016), our results point to the existence of the same fundamental principle in multiple aspects of human brain function, including functional integration and segregation. Beyond the results presented here, functional harmonics suggest novel ways to understand the dynamics of the human brain in health and in pathology, as well as to explore individual differences within this multi-dimensional harmonic representation.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Katharina Glomb (katharina.glomb@gmail.com).

Materials availability

This study did not generate new unique reagents.

Experimental models and subject details

General information

The data used in this study was acquired and made publicly available by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. All study protocols were approved by the Washington University institutional review board, and informed consent was obtained in all cases (Glasser et al., 2013; Van Essen et al., 2013).

For visualization purposes, we used the surfaces provided with the functional data (see Key resources table for surface files).

While all results presented here can be reproduced from the freely available data described above with the help of our published code, the functional harmonics and adjacency matrix referred to throughout this paper have been shared on Zenodo. The DOI is provided in the Key resources table.

Dense functional connectivity matrix

In this study, we used the dense functional connectivity matrix (dense FC), which is part of the Human Connectome Project’s 900 subjects data release (Glasser et al., 2013; Van Essen et al., 2012, 2013; Moeller et al., 2010; Feinberg et al., 2010; Setsompop et al., 2012; Xu et al., 2012; Jenkinson et al., 2002). It is available under https://db.humanconnectome.org/data/projects/HCP_1200. Clicking on “812 Subjects, recon r227, Dense Connectome” will download the appropriate .zip-archive (user login necessary). The list of names of all the files used in this study is shown in the Key resources table. Note that in this release, many of the subjects are related to at least one other subject of the group. The group average functional connectivity matrix was obtained by correlating group-PCA eigenmaps from 812 out of the 900 subjects included in this release, which are the subjects that completed all four sessions of 15-minute resting state fMRI.

Data are encoded in CIFTI file format (Glasser et al., 2013), which means that coordinates are defined on the cortical surface (“grayordinates”), i.e., using n vertices rather than voxels (Van Essen et al., 2012). The file was read using connectome workbench functions (Marcus et al., 2011) and converted to a single precision vector of length (due to its symmetry) using MATLAB (The MathWorks, Natick, MA). We also excluded the medial wall. This reduced the size of the FC matrix in memory from 33 GB to approximately 6 GB, greatly easing subsequent computations.

Retinotopy data

For the analyses involving retinotopic maps, we used data available on https://osf.io/bw9ec/ and described in Benson et al. (2018). The relevant file is named “prfresults.mat” and contains a variable “allresults” of dimensionality 91282 (grayordinates) × 6 (quantities) × 184 (181 subjects plus 3 different group averages) × 3 (model fits). We used only the quantities ‘ang’ and ‘ecc’, the first model fit, of the group average across all available subjects, which uses all available time points. See https://osf.io/bw9ec/wiki/home/ for details.

Task maps

For task reconstructions, we used data contained in the S1200 group average data release, which is available on https://www.humanconnectome.org/study/hcp-young-adult/document/extensively-processed-fmri-data-documentation, as “HCP_S1200_GroupAvg_v1 Dataset” (see Key resources table). Here we provide a summary of the tasks that form part of the HCP task battery (Barch et al., 2013). There are 7 groups of tasks: working memory, motor, gambling, language, social, emotional, relational. Subjects performed all tasks in two separate sessions (working memory, gambling, and motor in the first session, language, social cognition, relational processing, and emotion processing in the second).

Working memory. Four different stimulus types were used, presented in separate blocks: pictures of faces, places, tools and body parts. Two different task types were used: a 2-back working memory task, where subjects had to respond if a stimulus matched that two trials back, and a 0-back working memory task, where subjects had to respond whenever a single stimulus returned that was presented at the beginning of the block. This resulted in a total of 19 different working memory task maps, consisting of 14 activation maps (such as 0-back, 2-back, face, body, etc.) and 5 contrasts (between the two task types, between each stimulus type and the average across all stimuli, etc.).

Motor. Visual cues indicated whether participants should move their left or right fingers, left or right toes, or move their tongue. The goal was to identify the motor areas that correspond to these five body parts. This resulted in 26 different task maps (7 activation maps for 5 body parts plus visual cue plus average, and 6 contrast maps).

Gambling (Incentive processing). Subjects played a game in which they could win or lose money. The game was to guess whether the number on a “mystery card” that could range between 1 and 9 would be less or more than 5. The numbers were given after subjects made their guess and were chosen according to the trial type: “win” - the number would correspond to their guess and they would win 1$; “neutral” - the number would equal 5 and they would neither win nor lose any money; “loss” - the number would not correspond to the guess and participants would lose $0.50. Separate blocks are used in which trials are either mostly win or mostly lose, resulting in two conditions, punish and reward. This resulted in 3 different task maps (2 activation maps, i.e., one for each condition, and 1 contrast).

Language. Two different task types were used, “story” and “math.” “Story” consisted of participants listening to 5-9 sentences of a story, and answering a 2-alternative forced choice question thereafter. “Math” required participants to solve simple addition and subtraction problems. The two task types are similar in terms of auditory input and attentional load, but different in terms of semantic and numerosity related processing. As for gambling, the two task types resulted in 3 task maps (2 activation, 1 contrast).

Social (Theory of Mind, TOM). Subjects viewed videos of objects (squares, circles, triangles) that moved around in one of two ways: “Random” - there was no interaction between the objects, or “TOM” - the objects moved as if they were reacting to the other objects’ “thoughts and feelings.” The subjects then had to judge whether the objects were interacting or not, or respond with “not sure.” As with gambling and language, the two task types resulted in 3 task maps (2 activation, 1 contrast).

Emotional. Subjects viewed one of two types of stimuli, “faces” or “shapes,” and had to decide which of two stimuli presented at the bottom of the screen matched the stimulus at the top of the screen. The faces included emotional stimuli, i.e., angry or fearful expressions. Again, the two task types resulted in 3 task maps (2 activation, 1 contrast).

Relational. There were two conditions, “match” and “relational.” In all cases, stimuli could have one of six shapes combined with one of six textures. In the “match” condition, which served as a control condition, two shapes were presented at the top and one at the bottom of the screen. A word (“shape” or “texture”) that appeared in the middle of the screen instructed subjects to decide whether the bottom stimulus matched either of the top stimuli in the dimension indicated by the word. In the “relational” condition, two stimuli were presented each at the top and at the bottom of the screen, with no word in the middle. Instead, participants had to determine themselves across which dimension the top pair differed, and, subsequently, indicate whether the bottom pair differed over the same dimension. Again, the two task types resulted in 3 task maps (2 activation, 1 contrast).

Task maps were computed using FSL’s FEAT and FLAME (Woolrich et al., 2004; Jenkinson et al., 2012), conducting a between-subject (“level 2”) analysis, resulting in effect sizes (Cohen’s d). We used the task maps with minimal smoothing (2mm total smoothing); see 1200 subjects data release reference manual (https://www.humanconnectome.org/storage/app/media/documentation/s1200/HCP_S1200_Release_Reference_Manual.pdf), pp. 45-54 and 100-104.

Method details

Software

All data analysis was performed using MATLAB versions 2014b and 2017b. All custom code is available on Zenodo, the DOI is provided in the Key resources table. We also used scripts and functions from the following freely available software packages: Fieldtrip (version 20180903) (Oostenveld et al., 2011); Connectome workbench (Marcus et al., 2011). Links to these toolboxes are provided in the Key resources table.

Background: Functional Harmonics

The approach presented here relies on representing the human brain’s communication structure (dense FC) as a graph and estimating the eigenfunctions of the graph Laplacian applied to this structure. The graph representation of the dense FC is created by using the vertices sampled from the gray matter cortical surface as the nodes with n being the total number of nodes ( in this study) and the connections between the vertices as the edges ; the connections are taken as the correlations between vertices in the dense FC matrix. We represent this graph structure by its adjacency matrix that is formed by connecting each node i to its k-nearest neighbors ( in this study) according to its correlations in the dense FC matrix, i.e.:where is the set of the k largest values in row i in the dense FC matrix. In order to ensure A is symmetric, we also set , if . Defining in this way results in a symmetrical sparse binary matrix. Note that this symmetrization is necessary in order to ensure that an eigenvalue decomposition remains feasible.

Then we estimate the graph Laplacian defined aswhere is the adjacency matrix as defined above, and is the degree matrix, which is defined as a diagonal matrix with diagonal elementsAs such, the degree matrix contains each node’s degree in its diagonal. Due to the symmetrization of , it is possible that vertices have slightly more than nearest neighbors, as functional connectivity in the brain is not fully symmetrical (i.e., if j is among the k nearest neighbors of i, the reverse is not necessarily the case). By subtracting the adjacency from the degree matrix, a normalization is effectively performed, as the resulting graph Laplacian’s rows and columns sum to 0.

Finally, we estimate the functional harmonics as the eigenfunctions of by solving:where are the eigenvectors and are the corresponding eigenvalues.

Control function bases

Spherical rotations: We performed comparisons against spherical rotations of surface maps. We followed Alexander-Bloch et al. (2018), adapting freely available code (https://github.com/spin-test/spin-test) to be used with HCP surfaces. In this approach, surface maps are projected to a spherical surface and rotated by a random angle. Values are then mapped back to the nearest vertex, and the map is symmetrized in order to preserve this property. Parts of the corpus callosum that are rotated to the cortical surface are labeled as missing data (NaNs) and are ignored in any subsequent calculations (e.g., within- and across area distances, see below). Since we used multi-dimensional function bases, we rotated the surface maps corresponding to each dimension by the same angle. Note that, however, the resulting rotated function basis is no longer orthonormal due to the symmetry preserving step.

Eigenvectors of the dense FC: An intuitive basis is to take the eigenvectors of the dense FC without applying a threshold as done for obtaining the adjacency matrix. These eigenvectors have been shown to contain valuable information about dynamical FC (Cabral et al., 2017). The first 20 eigenvectors of the dense FC are shown in Data S1, page 2.

Eigenvectors of the adjacency matrix: In order to test the effect of thresholding/binarizing on the one hand and the effect of using the graph Laplacian instead of the adjacency matrix itself on the other, we also compared to the eigenvectors of the adjacency matrix, i.e., the dense FC thresholded such that only the 300 nearest neighbors of each vertex are retained and set to 1. The first 20 eigenvectors of the adjacency are shown in Data S1, page 3.

Principal components (PCs): PCA (principal component analysis) is a popular dimensionality reduction technique which preserves the maximum amount of variance in the data. It consists of taking the eigenvectors of the covariance matrix of the time series. These principal components are provided by the HCP via Connectome DB (see Key resources table). The first 20 PCs are shown in Data S1, page 4.

Independent components (ICs): A popular dimensionality reduction technique in resting state fMRI (Beckmann et al., 2005), independent component analysis is the foremost method for obtaining resting state networks. It consists of analyzing the time series of the data and finding those spatial patterns that are maximally independent from each other. We tested all sets of ICs that are provided by the HCP (see Key resources table), and found that the set with the lowest number of components, i.e., , performs best. Therefore, we restricted our comparisons to this set of ICs. Note that ICs are not orthonormal and therefore do not form a basis in the strictly mathematical sense. The 15 ICs used in our comparisons are shown in Data S1, page 5.

Modified silhouette values

To test whether functional harmonics are flat within the parcels as defined in the HCP parcellation (Glasser et al., 2016), we modified the silhouette value which is a commonly used measure to evaluate the degree to which a cluster structure captures the distances in the dataset (de Amorim and Hennig, 2015). We computed this modified silhouette value of each functional harmonic as:where is the average Euclidean distance on the functional harmonic between vertices belonging to a parcel i and vertices belonging to all other parcels, while is the average distance between vertices within the parcel i. Note that while the commonly used silhouette coefficient takes on values between −1 and 1, this modified version lies between 0 and 1: If all vertices belonging to a parcel i have the same value, and at least some vertices outside the parcel i have different values, then , and , indicating a perfectly flat gradient within the parcel i. If, on the other hand, there is no difference between and , . By averaging over the modified silhouette values of all parcels, one obtains a measure of how well the data fit the parcellation. Note that we replaced the somatosensory/motor core areas 1, 2, 3a, 3b, and 4 with the somatotopic sub-areas given by the HCP (Glasser et al., 2016) for a more detailed evaluation.

To evaluate the somatotopic organization of the functional harmonics, we used a measure that is similar to the silhouette value, but adapted to measure the separation from the rest of the cortex and from other somatotopic areas:where is the average Euclidean distance between vertices belonging to a somatotopic area and all vertices belonging to other somatotopic areas. The first term of the equation is between 1 and 2 and is close to 2 if both and are equal. Multiplying by ensures that is not large if both and are small.

Reconstructing task maps

The spatial pattern of each task map on on the cortex was decomposed into and reconstructed from the functional harmonics as:where the coefficient of each functional harmonic was estimated by projecting the task map onto that particular harmonic . As such, are estimated as:Then, each task map is reconstructed using Equation 7. In this study, we limit our reconstructions to using a maximum of 100 non-constant functional harmonics .

In order to make reconstruction errors comparable across tasks, and to allow us to average across tasks, we normalized the original and the reconstructed task maps to zero mean and unit variance across the cortical surface. For a reconstruction , where m indicates a binary vector of dimensionality which contains ones for harmonic basis functions that are used in the reconstruction and zeros otherwise, we then computed the normalized reconstruction error as:We normalized by the sum of squares of the normalized original task map in order to obtain values of RE that are comparable across tasks.

We also computed the Pearson correlations between and (reported in Figure S4).

Somatotopic areas

In the visual and somatosensory/motor cortices, some functional harmonics are rather determined by retinotopy and somatotopy than by anatomical or microstructural features. For the former, somatotopic areas occupy exactly the same surface area as the sensorimotor core areas, 1, 2, 3a, 3b, and 4. We therefore replaced, where appropriate, the borders of the HCP parcellation by the borders of the five somatotopic regions.

Parcel borders for visualization

In order to discuss the meaning of the functional harmonics, we show borders of certain parcels on the cortical surfaces (Figure 2). We used three different methods to select which borders to show. First, for some functional harmonics, it was feasible to select these areas manually (for example, early visual areas in functional harmonic 4, somatotopic areas in functional harmonics 3 and 4). In the supplementary neuroanatomical results of Glasser et al. (2016), the authors introduce a functional grouping of many regions that we often used as a guideline, for instance to distinguish between early and association auditory cortex. Second, for some functional harmonics (for instance, functional harmonics 1 and 2), we show the borders of parcels that belong to resting state networks as defined by Yeo et al. (2011). The 7-network parcellation is provided by the HCP, but does not perfectly overlap with the HCP parcellation. We adjusted the network borders slightly to align the network borders to follow those of the parcels defined in HCP. Thereby we assigned each parcel to the RSN with which it had the most overlap. Third, some functional harmonics are too complex to manually select areas or networks (namely, functional harmonics 5, 6, 8, and 10). Here we employed simple k-means clustering on the functional harmonic, using k = 2 (functional harmonics 5, 6, and 8) or k = 3 (functional harmonic 10); the clustering results are provided together with the code on Zenodo, see the Key resources table for the DOI. To obtain meaningful clusters in the somatosensory/motor cortex, we again replaced the sensorimotor core regions 1, 2, 3a, 3b and 4 with the somatotopic areas. For this purpose, we used vertices within the core regions and re-assigned them to the somatotopic areas based on their distances to the sub-area borders.

Third-party software

We used the following software packages to aide us in visualization.

plot_mesh by Gabriel Peyr, https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/5355/versions/5/previews/toolbox_graph/plot_mesh.m/index.html

(also includes getoptions.m and check_face_vertex.m)

freezeColors by John Iversen, https://www.mathworks.com/matlabcentral/fileexchange/7943-freezecolors-unfreezecolors

rgb by Ben Mitch, https://www.mathworks.com/matlabcentral/fileexchange/1805-rgb-m

subtightplot by Felipe G. Nievinski, https://www.mathworks.com/matlabcentral/fileexchange/39664-subtightplot

MarkerTransparency toolbox by Free Software Foundation, https://www.mathworks.com/matlabcentral/fileexchange/65194-peterrochford-markertransparency

Quantification and statistical analysis

Monte Carlo simulations

We used a Monte-Carlo approach for statistical validation.

For the silhouette values, we followed Alexander-Bloch et al. (2018), where permutations consist of rotated surface maps of the functional harmonics as well as principal components, independent components, eigenvectors of the dense FC, and eigenvectors of the adjacency matrix. Silhouette values were then computed for the original, non-rotated map as well as for rotated maps, and p values were computed based on the number P of rotations that performed better than the original map:We performed Bonferroni correction by multiplying the resulting p value by 11, i.e., the number of dimensions that was tested.

We used the same approach for the somatotopy index, but only applied to the functional harmonics and their rotations. Since in this case, we had five somatotopic areas (we averaged over the two hemispheres) and tested three of the 11 functional harmonics (, , and ), we required rotations in order to achieve a significance level of with 15 comparisons.

We also applied a Monte-Carlo permutation test to the mean reconstruction errors by permuting the labels of the basis 1000 times for each control basis. Here, we pooled the reconstruction errors over the first 11 non-constant components. For the overall reconstruction performance, we also pooled all 47 task maps; for ad hoc tests of each task category, we pooled only over the tasks in each category.

Auditory hierarchy

We observed that functional harmonic 10 captures the hierarchical organization of the auditory system (Figure S3D). To quantify this agreement, we measured the correlation between the spatial pattern of functional harmonic 10 and the extent to which each area is associated with the auditory network in the resting state (degree of auditory involvement) (Glasser et al., 2016). We found a significant correlation (, ) between functional harmonic 10 and the degree of auditory involvement of the functional areas.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
functional harmonics (HCP_S900_CORR_manifold_knn300.mat)	this manuscript	https://doi.org/10.5281/zenodo.5046434
k-means cluster memberships (cluster_kmeans.mat)	this manuscript	https://doi.org/10.5281/zenodo.5046434
adjacency matrix with 300 nearest neighbors (adjacency_nn300.mat)	this manuscript	https://doi.org/10.5281/zenodo.5046434
custom formatted version of RSN (Yeo 7 networks) surface labels	Human connectome project/this manuscript (converted for convenience from .nii to .gii using connectome workbench)	https://doi.org/10.5281/zenodo.5046434
Dense functional connectivity matrix (HCP_S1200_812_rfMRI_MSMAll_groupPCA_d4500ROW_zcorr_recon2.dconn.nii)	Human connectome project	https://db.humanconnectome.org/data/projects/HCP_1200 ; click on “812 Subjects, recon r227, Dense Connectome”
Medial wall index file (Human.MedialWall_Conte69.32k_fs_LR.dlabel.nii)	Human connectome project	BALSA database, https://balsa.wustl.edu/file/show/J2Vn
Cortical surfaces, inflated (S900. < hemisphere > .inflated_MSMAll.32k_fs_LR.surf)	Human connectome project	BALSA database, https://balsa.wustl.edu/file/show/QlL6
Cortical surfaces, flat (S900. < hemisphere > .flat.32k_fs_LR.surf.gii,)	Human connectome project	BALSA database, https://balsa.wustl.edu/file/show/2VVz4
Cortical surfaces, spherical (< hemisphere > .sphere.32k_fs_LR.surf.gii)	Human connectome project	BALSA database, https://balsa.wustl.edu/file/show/55xB
Cortical surface labels (Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii)	Human connectome project	BALSA database: https://balsa.wustl.edu/file/show/3VLx
Cortical areal borders (Q1-Q6_RelatedParcellation210. < hemisphere > .CorticalAreas.32k_fs_LR.border)	Human connectome project	BALSA database, https://balsa.wustl.edu/file/show/gmXm
HCP atlas colors	Human connectome project	https://osf.io/bw9ec/
HCP task maps (HCP_S1200_997_tfMRI_ALLTASKS_level2_cohensd_hp200_s2_MSMAll.dscalar.nii)	Human connectome project	https://db.humanconnectome.org/data/projects/HCP_1200; click on “Download group average dataset”
Retinotopy maps	Human connectome project	prfresults.mat from https://osf.io/bw9ec/; using ”group subject” (ID 999999) and full model fit
Principal components of dense FC	Human connectome project	https://db.humanconnectome.org/data/projects/HCP_1200; click on “812 Subjects, recon r227, Group-PCA Eigenmaps”
Independent components of HCP resting state data	Human connectome project	https://db.humanconnectome.org/data/projects/HCP_1200; click on “812 Subjects, recon r227 only, PTN Release”
Software and algorithms
custom code (MATLAB) used to produce analysis results and plot figures	this manuscript	https://doi.org/10.5281/zenodo.5046434
Fieldtrip (version 20180903)	Oostenveld et al., 2011	https://www.fieldtriptoolbox.org/download/ (registration required)
Connectome workbench	Human connectome project	https://www.humanconnectome.org/software/connectome-workbench