bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678.

The copyright holder for this preprint (which

was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

1 Functional harmonics reveal multi-dimensional

2 basis functions underlying cortical organization
3 Katharina Glomb1,* , Morten L. Kringelbach2,3 , Gustavo Deco4,5,6,7 , Patric Hagmann1 , Joel
4 Pearson8 , and Selen Atasoy2,3

5
1 Department of Radiology, Lausanne University Hospital and University of Lausanne (CHUV-UNIL), Lausanne,
6 Switzerland
7
2 Department of Psychiatry, University of Oxford, Oxford, United Kingdom

8
3 Center of Music in the Brain (MIB), Clinical Medicine, Aarhus University

9
4 Center of Brain and Cognition, Universitat Pompeu Fabra, Barcelona, Spain

10
5 ICREA, Institució Catalana de Recerca i Estudis Avancats (ICREA), Spain

11
6 Department of Neuropsychology, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany

12
7 School of Psychological Sciences, Monash University, Melbourne, Australia

13
8 School of Psychology, University of New South Wales, Sydney Australia

14
* [email protected]

15 ABSTRACT

The human brain consists of functionally specialized areas, which flexibly interact and integrate forming a multitude of complex
functional networks. The principles underlying this functional differentiation and integration remain unknown. Here, we
demonstrate that a fundamental principle ubiquitous in nature - harmonic modes - explains the orchestration of brain’s functional
organization. Applied to the functional connectivity in resting state averaged across 812 participants, harmonic modes give rise
to functional harmonics revealing the communication channels of the human brain. Each associated with a different spatial
16
frequency, the functional harmonics provide the frequency-ordered building blocks, which can reconstruct any pattern of brain
activity. We show that 47 brain activation patterns elicited by 7 different task categories in the Human Connectome Project
task battery can be reconstructed from a very small subset of functional harmonics, uncovering a parsimonious description of
the previously unknown relationship between task and resting state brain activity. Our findings evidence that the functional
specialization of the human cortex occurs in a gradiental manner across multiple dimensions in the functional harmonic basis,
where isolines of gradients follow the borders of functional areas as defined in a state-of-the-art parcellation of the human
cortex. Hence, the multi-dimensional functional harmonic representation unifies two competing views of the brain’s functional
organization, i.e. modular vs gradiental perspective.

17 Introduction
18 The human brain is topographically organized into functionally specific brain areas1 . Integration of these areas in various
19 different constellations allows for the immense complexity of human brain function2 . Despite remarkable progress in mapping
20 the brain into functionally meaningful subdivisions, known as cortical areas3, 4 , and in identifying functionally relevant
21 combinations of these areas forming the functional networks of the brain5 , the principles governing this functional segregation
22 and integration in the human brain have remained unknown. Here we demonstrate that a fundamental principle ubiquitous in
23 nature, i.e. the harmonic modes, when applied to functional connectivity data in humans, reveals both, the brain’s functional
24 networks as well as its topographic organization.
25 The topographic organization of the brain into functionally specific areas is one of its fundamental properties, observed
26 in evolution as early as the last common ancestor of vertebrates4, 6 . The individuality of each brain area is determined by its
27 functional specification, its microstructure (cyto- and myeloarchitecture)4 , and its inter- and intra-area connectivity3 . Significant
28 effort in neuroscience has been directed towards subdividing the brain into adjoining parcels, which are assumed to have
29 uniform functional specification and homogeneous connectivity3, 4 . A multitude of functionally distinct brain areas coordinate
30 through synchronous fluctuations in their activity7 . Coherent oscillations among distinct brain areas have been shown to be
31 another evolutionarily conserved aspect of brain activity8 . The overlap of the networks formed through these spontaneous
32 system oscillations, termed the functional connectivity patterns, with the functional networks of the human brain identified by
33 various sensory, motor, and cognitive task paradigms9–12 , strongly indicates their relevance for the brain’s functionality.
34 However, this modular view of brain organization, where separate, adjoining brain areas with uniform functionality and

1
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

35 homogeneous structural connectivity integrate into functional networks through coherent oscillations, has been challenged by
36 the presence of gradually varying boundaries between brain areas suggesting a degree of transition instead of sharply separated
37 brain areas13 , as well as by the existence of topographic mappings, which characterize the differences within a functionally
38 specific brain area14–16 . Topographic mappings including retinotopy14 , somatotopy15 , tonotopy16 , show that representation
39 of our visual field, body and auditory frequency spectrum are spatially continuously represented across the areas of the
40 primary visual, somatomotor and auditory cortices, respectively, challenging the assumption of uniform functionality within
41 the determined brain areas and demonstrating a smoothly varying functionality13 . As an alternative, theoretical work17, 18 and
42 recent experimental findings13 suggested a ”gradiental perspective”, where the functional organization of the cortex is argued to
43 be continuous, interactive and emergent as opposed to mosaic, modular and prededicated17 . Similar to the smoothly varying
44 functionality of primary sensory and motor areas, association cortices functioning as integration centres for more complex or
45 elaborated mental processes, are hypothesized to emerge from the convergence of information across sensory modalities18
46 with increasing spatial distance on the cortex from the highly functionally specialized primary cortices19 . Supporting this
47 hypothesis, a principal connectivity gradient of cortical organization in the human connectome has been identified, where the
48 functional networks of the human brain are located according to a functional spectrum from perception and action to more
49 abstract cognitive functions13 . Although converging evidence13, 20, 21 supports the continuous and emergent view of cortical
50 organization, the principles underlying the functional organization in the brain remain largely unknown.
51 Here, we demonstrate that the functional segregation and integration in the brain are governed by the same natural principle
52 of harmonic modes that underlies a multitude of physical and biological phenomena including the emergence of harmonic waves
53 (modes) encountered in acoustics22 , optics23 , electron orbits24, 25 , electro-magnetism26, 27 and morphogenesis28, 29 . By solving
54 the time-independent (standing) wave equation30, 31 on the communication structure of the human brain given by functional
55 magnetic brain imaging (fMRI) connectivity, we uncover the shape of the harmonic modes emerging from synchronous
56 oscillations in large scale brain activity. These harmonic modes decompose the functional connectivity into a hierarchical set of
57 frequency-specific communication channels, which naturally emerge from coherent, spontaneous brain activity, and unveil
58 both, the principal connectivity gradient13 , as well as cortical parcellations3 . Our results indicate that the functional segregation
59 and integration in the brain are governed by a multidimensional harmonic representation that we call ”functional harmonics”.
60 Finally, the decomposition of the brain activity maps elicited by various cognitive tasks into the set of functional harmonics
61 reveals that each task primarily involves activation of a very small subset of functional harmonics, suggesting that the functional
62 harmonics reveal fundamental building blocks of not only resting state activity, but also various cognitive functions.

63 Estimation of functional harmonics

64 Mathematically, the patterns of harmonic modes of a dynamical system are estimated by the eigendecomposition of the Laplace
65 operator, which lies at the heart of theories of heat, light, sound, electricity, magnetism, gravitation and fluid mechanics32 . In
66 vibrating systems, eigenfunctions of the Laplacian constitute standing waves, which also have been proposed as the mechanism
67 underlying cortical communication observed in electroencephalogram (EEG) data33 . Theoretical studies as well as experimental
68 findings have shown that spherical harmonics, i.e. Laplace eigenfunctions on a sphere, underlie cortical activation patterns
69 in fMRI34 . Harmonic modes of the structural connectivity of the human brain, i.e. Laplace eigenfunctions on the human
70 connectome, have been found to predict the collective dynamics of cortical activity at the macroscopic scale, and reveal resting
71 state networks31 .
72 In this work, we hypothesized that the harmonic modes of the brain’s communication structure given by its functional
73 connectivity underlie its functional integration and segregation. There are several crucial properties of harmonic modes that
74 lead us to form our hypothesis:
75 1) The dense functional connectivity (dense FC) matrix, in our study estimated from the pairwise temporal correlations
76 between all pairs of vertices on the cortical surface (59.412 vertices in total), encodes the communication structure of the
77 human brain. In order to find a multidimensional representation that best preserves this functional communication structure, we
78 utilized the discrete counterpart of the harmonic modes defined on a graph, i.e. the eigenvectors of the graph Laplacian, which
79 have been shown to optimally preserve the local graph structure while embedding it into a lower-dimensional space35 . Hence,
80 the functional harmonics estimated by the harmonic modes of the dense FC in this work, reveal the optimal multidimensional
81 mapping between the communication structure of the brain given by the dense FC and the cortical surface in the sense that the
82 strongest functional relationships given by the largest correlation values are optimally preserved.
83 2) Functional harmonics are the smoothest patterns that respect the constraints posed by the functional relationships given
84 by the FC35 . This implies that the average difference between neighbouring nodes in a graph representation is minimized.
85 Intriguingly, theoretical work has shown that activation patterns on graphs in which neighbouring nodes co-activate lead to
86 patterns with minimum free energy or entropy36, 37 , and that the transition between such patterns requires minimal energy38 .
87 3) Due to their orthogonality, Laplace eigenfunctions provide a new function basis. When applied to a one-dimensional
88 domain with cyclic boundary conditions, i.e. to a circle, Laplace eigenfunctions constitute the well-known Fourier basis,

2
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

89 whereas on a sphere, they yield the spherical harmonics. Each eigenfunction corresponds to a unique eigenvalue related
90 to its spatial frequency, and the set of all eigenfunctions forms a function basis, in which any signal can be represented
91 in the frequency domain. Considering this particular aspect of harmonic patterns, the functional harmonics provide a new
92 frequency-specific function basis driven by the brain’s communication structure (dense FC), where each dimension provides a
93 frequency-specific communication channel on the cortex.
94 4) The eigenfunctions of the Laplacian explain self-organizing patterns in many dynamical systems24–26, 28, 39 , ranging from
95 relatively simple physical phenomena like vibrating strings and metal plates22 to complex biological processes such as pattern
96 formation and morphogenesis28, 29 .
97 Considering that functional harmonics provide an optimal, frequency-specific mapping of the brain’s communication
98 structure to the cortex; that they represent the most energy-efficient activation patterns, which respect the constraints posed by
99 this communication structure; and given the ubiquity of harmonics in nature, we hypothesized that functional harmonics provide
100 the ideal candidate to explain functional segregation and integration in the brain. In order to test this hypothesis, we used the
101 dense FC computed from resting state fMRI (rfMRI) data averaged across 812 subjects, provided by the Human Connectome
102 Project (HCP), 900 subjects data release40–47 . We obtained the functional harmonics by estimating the eigenvectors of the
103 graph Laplacian computed on the graph representation of the FC (Figure 1).

104 Functional harmonics reveal functionally relevant communication channels

105 We first investigated whether functional harmonics yield functionally meaningful communication channels, i.e. patterns of
106 correlated activity. Figure 2 shows the first 11 non-constant functional harmonics (referred to as ψ1 -ψ11 ), ordered starting from
107 the lowest spatial frequency, illustrating that each harmonic is a smoothly varying pattern on the cortex between a positive and a
108 negative polarity; i.e., a gradient. Functional harmonics 1 and 2 correspond to previously identified large-scale gradients13 that
109 delineate the separation between the major sensory and the uni- vs. multimodal cortices in the brain, respectively (see SI Figure
110 1a). Figure 2a and b demonstrate the overlap between the visual and sensorimotor networks as defined in Yeo et al. (2011)48
111 and the gradiental patterns of the first and second functional harmonics. We observed that functional harmonic 3 reveals a
112 finer subdivision of the somatosensory/motor system49–51 . The overlay of the borders of the five somatotopic areas defined by
113 the HCP3, 52 on the third functional harmonic are shown in Figure 2c. Similarly, in functional harmonic 4, we found a finer
114 segregation of the visual system, following a retinotopic eccentricity gradient (for further details on retinotopic mapping see the
115 following section)53 . The overlay of the borders of early visual areas (V1-V4) on functional harmonic 4 are shown in Figure 2d.
116 The regions found in the positive polarity of functional harmonic 5 (borders shown on the left 3 panels in Figure 2e)
117 closely resemble the sensory-motor pathway48 (Figure 2e), which are known to be modulated by visuospatial attention54 . In
118 the negative polarity, we found the auditory cortex and parts of the somatosensory/motor network (Figure 2e). In contrast, in
119 functional harmonic 6, auditory and visual areas were both localized in the positive polarity, forming a network related to
120 audiovisual object (including faces) recognition55–57 , i.e. recognition of the ”outer world”. The negative polarity of harmonic
121 6 segregates the somatotopic face area as well as parts of the default mode network (DMN), a network of regions whose
122 activity has been related to self-referential tasks58 . Thus, the negative polarity of functional harmonic 6 forms a self-referential
123 processing stream58–60 . Functional harmonic 7 provides a further somatotopic gradient, including a higher hand area, 24dd,
124 in the medial cortex49 (see Figure 2g and annotations in Figure 2c). Functional harmonics 8 to 10 correspond to different
125 subdivisions of higher order networks such as the frontoparietal network and DMN (see SI Figure 2). In particular, the DMN61
126 is delineated in the positive polarity of functional harmonic 9 (borders of the DMN as defined by Yeo et al. (2011)48 are
127 overlaid on functional harmonic 9 in Figure 2i). Functional harmonic 11, the first asymmetric harmonic between the two
128 hemispheres, yields the separation between the right and left somatotopic hand areas62 . Overall, these results demonstrate that
129 functional harmonics provide a multitude of functionally relevant communication channels, each associated with a unique
130 spatial frequency, and enable a set of parallel processing streams in the human brain.

131 Functional harmonics reveal brain areas and topographic mappings

132 The fact that functional harmonics display both, well-delineated specialized regions (as for example in functional harmonics
133 3, 7, and 11) as well as gradients that integrate brain regions from different functional systems (as for example in functional
134 harmonics 5 and 6) led us to hypothesize that functional harmonics provide a unifying explanation for two seemingly opposing
135 perspectives of cortical organization: the gradiental perspective arguing that cortical organization is continuous on the one hand
136 and the modular perspective stating that brain function emerges from modular organization of specialized brain regions on the
137 other. We therefore explicitly tested whether functional harmonics fulfill the constraints posed by both of these views.
138 We first quantified the overall overlap between functional harmonic patterns and the brain regions (parcels) delineated
139 by the HCP parcellation3 . To this end, for each of the functional harmonics shown in Figure 2, we compared the within-
140 and between-area-variability of each cortical region (see SI Figures 3-5), where a large difference between the within- and

3
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

Figure 1. Workflow for the estimation of functional harmonics. a: Brain activity measured with functional magnetic resonance
imaging (fMRI) in resting state for 820 subjects provided by the Human Connectome Project (HCP, 900 subjects data
release)40–47 . b: Illustration of brain time series activity of three representative vertices on the cortex (x1 , x2 , · · · , xn ). 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 location of the vertices are colour-coded3 . 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 64984 × 1 dimensional vectors (middle) and on the cortical surface (bottom). For illustrative
purposes, the graph representations are shown for a parcellated version of the FC matrix using the HCP parcellation3 in d and e.
We note that the computation of the functional harmonics have been performed on the dense FC using 64984 × 64984 without
using any parcellation.

141 between-area variability, indexed by a large silhouette value, indicates that that particular region is well-separated from the rest
142 of the cortex63 . We found an alignment between the isolines of the functional harmonics and parcel borders, as verified by
143 significantly larger silhouette values for functional harmonics compared to harmonics computed on randomized connectivity
144 graphs (values for functional harmonics between 0.67 and 0.87, maximum random value 0.20, significantly different with
145 pcorr < 0.05 after Bonferroni correction, Monte Carlo tests; see Online Methods for details). For qualitative evaluation, the
146 overlap between parcels and functional harmonics is shown in SI Figure 6.
147 In addition to the parcels delineated in the HCP parcellation3 , we investigated whether functional harmonics also capture
148 somatotopy15 and retinotopy14 , two major topographic mappings found in the brain. Topographic mappings represent sensory
149 input on the cortical surface such that the relative positions of the receptors, which receive these inputs, are preserved. Five
150 somatotopic sub-areas (in each hemisphere) as defined by the HCP3 form a topographic map of the surface of the body on the
151 cortex, i.e., the face, hands, eyes, feet, and trunk.
152 To quantify the degree to which each somatotopic sub-area is delineated within a particular functional harmonic, we again

4
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

Figure 2. a-k: The first 11 non-constant harmonics plotted on the cortical surface. White lines show borders of HCP parcels.
V1-V4: visual areas 1 to 4; MT: middle temporal visual area; 24dd: an area that contains a higher order representation of the
hand; fusiform face complex: an area that responds specifically to images of human faces.

153 utilized the within- and between-area-variability as above, but applied specifically to the 10 somatotopic sub-areas (see SI Figure
154 7). We measured their separation both from the rest of the brain as well as from other somatotopic areas. We found that for all
155 functional harmonics shown in Figure 2 except for the one with the lowest frequency that delineates the complete somatomotor
156 cortex, at least one somatotopic region is significantly separated (pcorr < 0.05 after Bonferroni correction, Monte Carlo tests).
157 This finding indicates that functional harmonics capture somatotopic organization in the cortex. Figure 3a illustrates the
158 two-dimensional subspace formed by functional harmonics 3 and 11, which strikingly accounts for the precise mapping of the
159 human body onto the somatotopic regions of the cortex (see SI Figure 1b-d for further examples).
160 We next investigated the presence of retinotopic mapping of early visual regions (V1-V4), where cortical representations of
161 the visual field reflect the positions of the receptors such that each vertex within the patterns of functional harmonics is assigned
162 an eccentricity (distance from the fovea) and an angle (top, bottom, left, right)53 . To investigate the degree of agreement
163 between the functional harmonics and the retinotopic mappings, we measured the correlation between eccentricity as well
164 as polar angle maps and the functional harmonic patterns in V1-V4. We found significant correlations (pcorr < 0.05 after
165 Bonferroni correction) between the retinotopic eccentricity map and all harmonics except functional harmonic 9; and between
166 the retinotopic angular map and harmonics 1-4, 7-9, and 11. Examples of polar plots of the retinotopic gradients are shown in

5
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

Figure 3. a: Functional harmonics 3 (ψ3 ) and 11 (ψ11 ) in their own space. The location of 4 somatotopic areas in this space is
annotated. b: Correlation between the degree to which areas are related to auditory regions3 and the value of functional
harmonic 10 (ψ10 ), averaged within each of the 360 parcels. The color code is taken from the parcellation in Glasser et al.
(2016)3 , see also figure 1d. c and d: Retinotopies of functional harmonics 4 (c; ψ4 ) and 8 (d; ψ8 ). Each panel shows, on the left,
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 right, the respective functional
harmonic is shown on a flat map of early visual cortex (left hemisphere). V1, V2, V3, V4: visual areas 1, 2, 3, 4.

167 Figure 3c, d (all polar plots available in SI Figure 8).

168 Besides the two major topographic mappings of the cortex, we observed that functional harmonic 10 captures the hierarchical
169 organization of the auditory system. To quantify this agreement, we measured the correlation between the spatial pattern of
170 functional harmonic 10 and the extent to which each area is associated with the auditory network in the resting state (degree of
171 auditory involvement)3 . We found a significant correlation (r = −0.63, p = 4 · 10−21 ) between functional harmonic 10 and the
172 degree of auditory involvement of the functional areas (Figure 3b).

173 Functional harmonics are basis functions of human cognition

174 Considering the parallel between functional harmonics and the well-known Fourier basis, i.e. the fact that they both are defined
175 as Laplacian eigenfunctions, the former applied to a one-dimensional domain with cyclic boundaries (a circle) and the latter
176 to the communication structure of the human brain (dense FC matrix), the functional harmonics provide an extension of the
177 Fourier basis to the communication structure of the human brain. As such, they provide per definition a frequency-specific
178 function basis, in which any pattern of brain activity can be represented as a weighted combination of functional harmonics.
179 Given the experimental evidence showing that resting state functional connectivity reflects connectivity during task9–12 , we
180 hypothesized that functional harmonics provide building blocks of task activity measured on the cortex. In order to test this
181 hypothesis, we reconstructed 47 group-level task maps provided by the HCP52 from the functional harmonics (see Online
182 Methods). The 47 maps are activation maps as well as contrasts derived from 7 groups of tasks (working memory, motor,
183 gambling, language, social, emotional, relational - see Online Methods for summaries). The functional harmonic reconstruction
184 yields a coefficient (weight) for each functional harmonic, quantifying how much it contributes to a certain task map. The set of
185 all coefficients forms a spectrum equivalent to the power spectrum obtained from a Fourier transform, in this case the power
186 spectrum of the functional harmonic basis.
187 We first tested whether it is possible to approximate task maps as superpositions of subsets of functional harmonics, linearly
188 combining them in the order of their eigenvalues. We quantified the goodness of fit by measuring the distance between the
189 original and the reconstructed task maps. Figure 4a shows the normalized reconstruction errors for all groups of tasks: The
190 error drops from about 1.00 to about 0.40 when using the 11 exemplary harmonics of Figure 2, corresponding to 0.02% of the
191 total functional harmonic spectrum. Figure 4b illustrates this for one specific task (working memory - body; see also SI Figures

6
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

Figure 4. a: Reconstruction errors for all 7 task groups. Thick line: mean, the shaded area contains all values within that task
group (all curves shown in SI Figure 9a); b: One example for a reconstruction using a working memory task. The top panel is the
original task activation map (working memory - body), and subsequent panels use the number of harmonics indicated on the left
to reconstruct it. c: Like a, but showing correlation between reconstructed and original task maps (all curves shown in SI Figure
9b).

192 10-16).
193 As a control, we performed 418 Monte Carlo simulations, in which we randomized the dense FC matrix while preserving
194 the degree of each node. We computed the functional harmonics of each randomized connectivity matrix. This control still
195 yields harmonic patterns emerging on the cortex computed as the eigenfunctions of the Laplacians of the randomized FC graphs,
196 but these random harmonics are estimated from a random connectivity instead of the communication structure of the human
197 brain. The black line in Figure 4a shows the lowest reconstruction error across all randomizations (0.65), for any number
198 of functional harmonics between 0 (only the constant harmonic) and 100. When using 30 functional harmonics (0.05% of
199 the functional harmonic spectrum), the reconstruction with the maximum reconstruction error across all 47 task maps still
200 outperforms the best reconstruction using the randomized harmonic basis. As a complementary measure to the reconstruction
201 error, we computed the correlation between the reconstructed and the original task activation maps, shown in Figure 4c. When
202 using 9 or more functional harmonics, all task reconstructions significantly correlate with their original task map (pcorr < 0.05
203 after Bonferroni correction; maximum random correlation: 0.05, shown in black in Figure 4c; see Online Methods for details).
204 Given that functional harmonics constitute functionally relevant communication channels, we hypothesized that the task
205 activation maps can be characterized by their power spectrum. Figure 5a, d and Figure 5b, e, show two examples of task
206 activation maps and the corresponding normalized power of the first 11 non-constant functional harmonics, respectively,
207 revealing how strongly each of the 11 functional harmonics shown in Figure 2 contributes to these particular task maps. For
208 qualitative evaluation, we display the task activation maps reconstructed by superimposing functional harmonics in the order of
209 their contribution strength for varying numbers of functional harmonics in Figure 5c, f (see also SI Figures 10-16). Across all
210 47 task maps that were evaluated, the functional harmonic which was the strongest contributor was always either the constant
211 functional harmonic or one of the first 11 non-constant harmonics shown in Figure 2.
212 In order to evaluate the uniqueness of the functional harmonic power spectrum of each task activation map, we computed
213 the distance between a given reconstructed map and all original task maps, resulting in a confusion matrix for each number of
214 harmonics with maximum contribution. If task maps can indeed be characterized by their functional harmonics power spectra,
215 the error should be minimal between a reconstruction and its corresponding task map compared to the error of the reconstruction
216 of the other 46 task maps. The confusion matrices in Figure 5g-i show the pairs of the original and reconstructed task activation
217 maps with the minimum distance when using 1, 4, and 40 functional harmonics with maximum contribution. Coloured squares
218 mark the 7 task groups as in Figure 4. The proportion of unambiguously identified tasks in relation to the number of functional
219 harmonics is shown in Figure 5j. We found that sparse representations using the 4 functional harmonics with the largest power
220 for each task are sufficient to unambiguously characterize the seven task groups (Figure 5h), and 70% of all individual tasks.

7
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

Figure 5. a: Map of the contrast between working memory (face) and average working memory from the HCP task dataset52 ,
b: Contributions (normalized coefficients of the graph Fourier transform) for the first 11 non-constant functional harmonics, c:
Reconstruction of the task map in panel 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: Confusion
matrices. Black entries mark the task map-reconstruction-pair which has the lowest reconstruction error, colored squares
indicate the task group as in Figure 4. j: 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).

221 When the 40 functional harmonics with maximum contribution are used, which corresponds to 0.1% of the complete spectrum
222 of functional harmonics, all of the 47 task maps are correctly identified from their reconstructions (Figure 5i). The maximum
223 number of correctly classified tasks was significantly larger for functional harmonics than for randomized harmonics (47 vs 4
224 correctly classified tasks, p < 0.005, Monte Carlo simulations with 418 trials).
225 Overall, our results demonstrate that that functional harmonics provide a novel functionally relevant representation, where
226 the brain activity accompanying different tasks can be uniquely identified from the activation profiles of a small range of
227 functional harmonics.

228 Discussion
229 We reveal a previously unknown principle of cortical organization by applying a fundamental principle ubiquitous in nature
230 - harmonic modes - to the communication structure of the human brain. The resulting modes termed functional harmonics
231 reveal a data-driven, frequency-specific function basis derived from the human resting state functional connectivity matrix and
232 constitute the optimal mapping of the communication structure encoded in this matrix onto the cortex. Here, we demonstrate
233 the meaning of the first 11 functional harmonics as brain’s functional communication channels, which are distinctly separated
234 and ordered in the frequency domain. In this view, functional harmonics provide a multidimensional representation of brain
235 activity that is not confined to a certain spatial scale. In particular, the mechanism by which a brain region is able to fulfill a
236 multitude of functions depending on the context is enabled by its simultaneous membership in several communication channels.
237 Thus, functional harmonics unify the competing views that brain activity arises either from smoothly varying gradients or from
238 the modular and specialized regions: specialized regions emerge from the interaction of functional harmonics across multiple
239 dimensions. Hence our findings provide, to our knowledge, the first principle that unifies the gradiental and modular aspects
240 and reveals the multi-dimensional nature of cortical organization.

8
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

241 Furthermore, by definition, functional harmonics are the extension of the well-known Fourier basis to the functional
242 connectivity of the human brain. As such they provide a function basis to reconstruct any pattern of brain activity as
243 superpositions of these harmonic patterns. We explicitly show that functional harmonics are building blocks of cognitive
244 activity in the brain by characterizing a multitude of task activation maps from their functional harmonic reconstructions.
245 Considering that the principle of harmonic modes when applied to the structural connectivity of the human brain - the
246 human connectome - have been shown to reveal the functional networks31 , our results point to the emergence of the same
247 fundamental principle in multiple aspects of human brain function. Beyond the results presented here, functional harmonics
248 suggest novel ways to understand the dynamics of the human brain in health and in pathology as well to explore individual
249 differences within this multi-dimensional harmonic representation.

250 References
251 1. Felleman, D. J. & Van, D. E. Distributed hierarchical processing in the primate cerebral cortex. Cereb. cortex (New York,
252 NY: 1991) 1, 1–47 (1991).
253 2. Tononi, G., Sporns, O. & Edelman, G. M. A measure for brain complexity: relating functional segregation and integration
254 in the nervous system. Proc. Natl. Acad. Sci. 91, 5033–5037 (1994).
255 3. Glasser, M. F. et al. A multi-modal parcellation of human cerebral cortex. Nat. 536, 171–178 (2016).
256 4. Eickhoff, S. B., Constable, R. T. & Yeo, B. T. Topographic organization of the cerebral cortex and brain cartography.
257 Neuroimage 170, 332–347 (2018).
258 5. Damoiseaux, J. et al. Consistent resting-state networks across healthy subjects. Proc. national academy sciences 103,
259 13848–13853 (2006).
260 6. Krubitzer, L. & Kaas, J. The evolution of the neocortex in mammals: how is phenotypic diversity generated? Curr. opinion
261 neurobiology 15, 444–453 (2005).
262 7. Varela, F., Lachaux, J.-P., Rodriguez, E. & Martinerie, J. The brainweb: phase synchronization and large-scale integration.
263 Nat. reviews neuroscience 2, 229 (2001).
264 8. Vincent, J. et al. Intrinsic functional architecture in the anaesthetized monkey brain. Nat. 447, 83–86 (2007).
265 9. Biswal, B., Zerrin Yetkin, F., Haughton, V. & Hyde, J. Functional connectivity in the motor cortex of resting human brain
266 using echo-planar MRI. Magn. resonance medicine 34, 537–541 (1995).
267 10. Deco, G., Jirsa, V. K. & McIntosh, A. R. Emerging concepts for the dynamical organization of resting-state activity in the
268 brain. Nat. Rev. Neurosci. 12, 43 (2011).
269 11. Buckner, R. L., Krienen, F. M. & Yeo, B. T. Opportunities and limitations of intrinsic functional connectivity mri. Nat.
270 neuroscience 16, 832 (2013).
271 12. Cole, M. W., Ito, T., Bassett, D. S. & Schultz, D. H. Activity flow over resting-state networks shapes cognitive task
272 activations. Nat. neuroscience 19, 1718 (2016).
273 13. Margulies, D. S. et al. Situating the default-mode network along a principal gradient of macroscale cortical organization.
274 Proc. Natl. Acad. Sci. 113, 12574–12579 (2016).
275 14. Sereno, M., Pitzalis, S. & Martinez, A. Mapping of contralateral space in retinotopic coordinates by a parietal cortical area
276 in humans. Sci. 294, 1350–1354 (2001).
277 15. Penfield, W. & Rasmussen, T. The cerebral cortex of man; a clinical study of localization of function. (Macmillan, 1950).
278 16. Perrone-Capano, C., Volpicelli, F. & di Porzio, U. Biological bases of human musicality. Rev. Neurosci. 28, 235–245
279 (2017).
280 17. Goldberg, E. Gradiental approach to neocortical functional organization. J. Clin. Exp. Neuropsychol. 11, 489–517 (1989).
281 18. Mesulam, M.-M. From sensation to cognition. Brain: a journal neurology 121, 1013–1052 (1998).
282 19. Buckner, R. L. & Krienen, F. M. The evolution of distributed association networks in the human brain. Trends cognitive
283 sciences 17, 648–665 (2013).
284 20. Huntenburg, J. M., Bazin, P.-L. & Margulies, D. S. Large-scale gradients in human cortical organization. Trends cognitive
285 sciences 22, 21–31 (2018).
286 21. Haak, K. V., Marquand, A. F. & Beckmann, C. F. Connectopic mapping with resting-state fmri. Neuroimage 170, 83–94
287 (2018).

9
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

288 22. Chladni, E. F. F. Die Akustik (Breitkopf & Härtel, 1802).

289 23. Bedzyk, M. J., Bilderback, D. H., Bommarito, G. M., Caffrey, M. & Schildkraut, J. S. X-ray Standing Waves: A Molecular
290 Yardstick for Biological Membranes. Sci. 1–4 (1988).
291 24. Schrödinger, E. An undulatory theory of the mechanics of atoms and molecules. Phys. Rev. 28, 1049 (1926).
292 25. Moon, C. R. et al. Quantum Phase Extraction in Isospectral Electronic Nanostructures. Sci. 319, 782–787 (2008).
293 26. Roos, C. Quantum physics: Simulating magnetism. Nat. 484, 461–462 (2012).
294 27. Britton, J. W. et al. Engineered two-dimensional Ising interactions in a trapped-ion quantum simulator with hundreds of
295 spins. Nat. 484, 489–492 (2012).
296 28. Murray, J. D. How the leopard gets its spots. Sci. Am. 258, 80–87 (1988).
297 29. Xu, Y., Vest, C. M. & Murray, J. D. Holographic interferometry used to demonstrate a theory of pattern formation in
298 animal coats. Appl. optics 22, 3479–3483 (1983).
299 30. Levy, B. Laplace-beltrami eigenfunctions towards an algorithm that” understands” geometry. In IEEE International
300 Conference on Shape Modeling and Applications 2006 (SMI’06), 13–13 (IEEE, 2006).
301 31. Atasoy, S., Donnelly, I. & Pearson, J. Human brain networks function in connectome-specific harmonic waves. Nat.
302 communications 7 (2016).
303 32. Stewart, I. Mathematics: Holes and hot spots. Nat. 401, 863–865 (1999).
304 33. Nunez, P. L. & Srinivasan, R. A theoretical basis for standing and traveling brain waves measured with human eeg with
305 implications for an integrated consciousness. Clin. neurophysiology 117, 2424–2435 (2006).
306 34. Robinson, P. A. et al. Eigenmodes of brain activity: Neural field theory predictions and comparison with experiment.
307 NeuroImage 142, 79–98 (2016).
308 35. Belkin, M. & Niyogi, P. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation
309 15, 1373–1396 (2003).
310 36. Tomasi, D., Wang, G.-J. & Volkow, N. D. Energetic cost of brain functional connectivity. Proc. Natl. Acad. Sci. 110,
311 13642–13647 (2013).
312 37. Gu, S. et al. The energy landscape of neurophysiological activity implicit in brain network structure. Sci. reports 8, 2507
313 (2018).
314 38. Gu, S. et al. Controllability of structural brain networks. Nat. communications 6, 8414 (2015).
315 39. Britton, J. W. et al. Engineered two-dimensional ising interactions in a trapped-ion quantum simulator with hundreds of
316 spins. Nat. 484, 489 (2012).
317 40. Glasser, M. F. et al. The minimal preprocessing pipelines for the human connectome project. Neuroimage 80, 105–124
318 (2013).
319 41. Van Essen, D. C. et al. The WU-Minn human connectome project: an overview. Neuroimage 80, 62–79 (2013).
320 42. Moeller, S. et al. Multiband multislice ge-epi at 7 tesla, with 16-fold acceleration using partial parallel imaging with
321 application to high spatial and temporal whole-brain fmri. Magn. Reson. Medicine 63, 1144–1153 (2010).
322 43. Feinberg, D. A. et al. Multiplexed echo planar imaging for sub-second whole brain fmri and fast diffusion imaging. PloS
323 one 5, e15710 (2010).
324 44. Setsompop, K. et al. Blipped-controlled aliasing in parallel imaging for simultaneous multislice echo planar imaging with
325 reduced g-factor penalty. Magn. resonance medicine 67, 1210–1224 (2012).
326 45. Xu, J. et al. Highly accelerated whole brain imaging using aligned-blipped-controlled-aliasing multiband epi. In
327 Proceedings of the 20th Annual Meeting of ISMRM, vol. 2306 (2012).
328 46. Jenkinson, M., Bannister, P., Brady, M. & Smith, S. Improved optimization for the robust and accurate linear registration
329 and motion correction of brain images. Neuroimage 17, 825–841 (2002).
330 47. Van Essen, D. C., Glasser, M. F., Dierker, D. L., Harwell, J. & Coalson, T. Parcellations and hemispheric asymmetries of
331 human cerebral cortex analyzed on surface-based atlases. Cereb. cortex 22, 2241–2262 (2011).
332 48. Yeo, B. T. et al. The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J.
333 neurophysiology 106, 1125–1165 (2011).

10
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

334 49. Zeharia, N., Hertz, U., Flash, T. & Amedi, A. New whole-body sensory-motor gradients revealed using phase-locked
335 analysis and verified using multivoxel pattern analysis and functional connectivity. J. Neurosci. 35, 2845–2859 (2015).
336 50. Zeharia, N., Hertz, U., Flash, T. & Amedi, A. Negative blood oxygenation level dependent homunculus and somatotopic
337 information in primary motor cortex and supplementary motor area. Proc. Natl. Acad. Sci. 201119125 (2012).
338 51. Kuehn, E. et al. Body topography parcellates human sensory and motor cortex. Cereb. Cortex 27, 3790–3805 (2017).
339 52. Barch, D. M. et al. Function in the human connectome: task-fmri and individual differences in behavior. Neuroimage 80,
340 169–189 (2013).
341 53. Benson, N. C. et al. The human connectome project 7 tesla retinotopy dataset: Description and population receptive field
342 analysis. J. vision 18, 23–23 (2018).
343 54. Corbetta, M. & Shulman, G. L. Control of goal-directed and stimulus-driven attention in the brain. Nat. reviews
344 neuroscience 3, 201 (2002).
345 55. Beauchamp, M. S., Lee, K. E., Argall, B. D. & Martin, A. Integration of auditory and visual information about objects in
346 superior temporal sulcus. Neuron 41, 809–823 (2004).
347 56. Levy, I., Hasson, U., Avidan, G., Hendler, T. & Malach, R. Center–periphery organization of human object areas. Nat.
348 neuroscience 4, 533 (2001).
349 57. Hasson, U., Harel, M., Levy, I. & Malach, R. Large-scale mirror-symmetry organization of human occipito-temporal
350 object areas. Neuron 37, 1027–1041 (2003).
351 58. Gusnard, D. A., Akbudak, E., Shulman, G. L. & Raichle, M. E. Medial prefrontal cortex and self-referential mental activity:
352 relation to a default mode of brain function. Proc. Natl. Acad. Sci. 98, 4259–4264 (2001).
353 59. Haxby, J. V., Hoffman, E. A. & Gobbini, M. I. The distributed human neural system for face perception. Trends cognitive
354 sciences 4, 223–233 (2000).
355 60. Leslie, K. R., Johnson-Frey, S. H. & Grafton, S. T. Functional imaging of face and hand imitation: towards a motor theory
356 of empathy. Neuroimage 21, 601–607 (2004).
357 61. Raichle, M. et al. A default mode of brain function. Proc. Natl. Acad. Sci. 98, 676–682 (2001).
358 62. Pool, E.-M., Rehme, A. K., Fink, G. R., Eickhoff, S. B. & Grefkes, C. Handedness and effective connectivity of the motor
359 system. Neuroimage 99, 451–460 (2014).
360 63. de Amorim, R. C. & Hennig, C. Recovering the number of clusters in data sets with noise features using feature rescaling
361 factors. Inf. Sci. 324, 126–145 (2015).
362 64. Marcus, D. et al. Informatics and data mining tools and strategies for the human connectome project. Front. neuroinfor-
363 matics 5, 4 (2011).
364 65. MATLAB. version R2014b (The MathWorks Inc., Natick, Massachusetts, 2014).
365 66. Molloy, M. & Reed, B. A critical point for random graphs with a given degree sequence. Random structures & algorithms
366 6, 161–180 (1995).
367 67. Woolrich, M. W., Behrens, T. E., Beckmann, C. F., Jenkinson, M. & Smith, S. M. Multilevel linear modelling for fmri
368 group analysis using bayesian inference. Neuroimage 21, 1732–1747 (2004).
369 68. Jenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W. & Smith, S. M. Fsl. Neuroimage 62, 782–790 (2012).

11
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

370 Online Methods

371 Data
372 The data used in this study was acquired and made publicaly available by the Human Connectome Project, WU-Minn
373 Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes
374 and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience
375 at Washington University. All study protocols were approved by the Washington University institutional review board, and
376 informed consent was obtained in all cases40, 41 .
377 In this study, we used the dense functional connectivity (FC) matrix, which is part of the Human Connectome Project’s 900
378 subjects data release40–47 . It is available under db.humanconnectome.org/data/projects/HCP 120064 . Clicking on ”812 Subjects,
379 recon r227, Dense Connectome” will download the appropriate .zip-archive (user login necessary). The list of names of all
380 the files used in this study is shown in Table 1. Note that in this release, many of the subjects are related to at least one other
381 subject of the group. The group average functional connectivity matrix was obtained by correlating group-PCA eigenmaps
382 from 812 out of the 900 subjects included in this release, which are the subjects that having completed all four sessions of
383 15-minute resting state fMRI.
384 For task reconstructions, we used data contained in the S1200 group average data release, which is available on
385 www.humanconnectome.org/study/hcp-young-adult/document/extensively-processed-fmri-data-documentation, as
386 ”HCP S1200 GroupAvg v1 Dataset”.
387 For the analyses involving retinotopic maps, we used data available on osf.io/bw9ec/ and described in Benson et al. (2018)53 .
388 The relevant file is named ”prfresults.mat” and contains a variable ”allresults” of dimensionality 91282 (grayordinates) ×
389 6 (quantities) × 184 (181 subjects plus 3 different group averages) × 3 (model fits). We used only the quantities ’ang’
390 and ’ecc’, the first model fit, of the group average across all available subjects, which uses all available time points. See
391 osf.io/bw9ec/wiki/home/ for details.
392 Data are encoded in CIFTI file format40 , which means that coordinates are defined on the cortical surface (”grayordinates”),
393 i.e. using n vertices rather than voxels47 . The file was read using connectome workbench functions64 and converted to a single
394 precision vector of length (n · n − n)/2 (due to its symmetry) using Matlab65 . We also excluded the medial wall. This reduced
395 the size of the FC matrix in memory from 33 GB to approximately 6 GB, greatly easing subsequent computations. The loss
396 in precision is negligible compared to the accuracy with which pairwise correlation can be estimated from noisy fMRI time
397 courses.
398 For visualization purposes, we used the surfaces provided with the functional data.

399 Software
400 All data analysis was performed using MATLAB 2014b or 2017b, using also scripts and functions from the following freely
401 available software packages:

402 • Fieldtrip version 20180903

403 • Connectome workbench (https://www.humanconnectome.org/software/connectome-workbench)
404 • CNM toolbox (https://www.mathworks.com/matlabcentral/fileexchange/45734-cnm)
405 • gifti toolbox (https://www.artefact.tk/software/matlab/gifti/)

406 Background: Functional Harmonics

407 The approach presented here relies on representing the human brain’s communication structure (dFC) as a graph and estimating
408 the eigenfunctions of graph Laplacian applied to this structure. The graph representation of the brain’s communication
409 structure G = (V , E ) is created by representing the vertices sampled form the gray matter cortical surface as the nodes
410 V = {vi |i ∈ 1, · · · , n} with n being the total number of nodes (n = 59.412 in this study) and by representing the connections
411 between the vertices as the edges E = {ei j |(vi , v j ) ∈ V × V }, which come from the connections in the dFC matrix. We
412 represent this graph structure G by its n × n adjacency matrix A = [ai j ] that is formed by connecting each node i to its k-nearest
413 neighbours (k = 300 in this study) according to its correlations in the dFC matrix, i.e.:

(
1 ci j ∈ κi , ∀ j : 1 ≤ j ≤ n, j 6= i
ai j = (1)
0 ci j ∈
/ κi , ∀ j : 1 ≤ j ≤ n, j 6= i ,

414 where κi is the set of the k largest values in row i in the dFC matrix. In order to ensure A is symmetric, we also set a ji = 1, if
415 ai j = 1. Defining A as such results in a symmetrical sparse binary matrix.

12
bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

Table 1. Files used in our computations. All data was downloaded from the human connectome project database
(db.humanconnectome.org/data/projects/HCP 1200) unless otherwise specified.

Purpose file name comment

Dense functional con- HCP S900 820 rfMRI MSMAll groupPCA d4500ROW...
nectivity matrix ... zcorr.dconn.nii

Medial wall index file Human.MedialWall Conte69.32k fs LR.dlabel.nii

Cortical surfaces S900.<hemisphere>.inflated MSMAll.32k fs LR.surf, replace

S900.<hemisphere>.flat.32k fs LR.surf.gii <hemisphere>
with ”L” for left
hemisphere, ”R” for
right

Cortical surface la- Q1-Q6 RelatedValidation210.<hemisphere>. ... replace

bels ...CorticalAreas Final Final Areas Group Colors.32k fs LR.label.gii <hemisphere>
with ”L” for left
hemisphere, ”R” for
right

Borders Parcellation: Q1-Q6 RelatedParcellation210. ... replace

...<hemisphere>.CorticalAreas.32k fs LR.border, <hemisphere>
Somatotopy: Q1-Q6 RelatedParcellation210. ... with ”L” for left
...<hemisphere>.SubAreas.32k fs LR.border hemisphere, ”R” for
right

HCP atlas colours atlas.mat from osf.io/bw9ec/

Task maps HCP S1200 997 tfMRI ALLTASKS level2 cohensd... 2 mm smoothing ker-
... hp200 s2 MSMAll.dscalar.nii nel
www.humanconnectome.org/study/hcp-young-
adult/document/extensively-processed-fmri-data-documentation

Retinotopy maps prfresults.mat from osf.io/bw9ec/ using ”group subject”

(ID 999999) and full
model fit

Yeo et al. 7-networks RSN-networks.<hemisphere>.32k fs LR.label replace

parcellation <hemisphere>
with ”L” for left
hemisphere, ”R” for
right

13
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

Then we estimate the graph Laplacian defined as

LG = D − A , (2)

where A is the adjacency matrix as defined above, and D is the degree matrix, which is defined as a diagonal matrix with
diagonal elements
n
dii = ∑ ai j . (3)
j=1

As such, the degree matrix D contains each node’s degree in its diagonal. Finally, we estimate the functional harmonics as the
eigenfunctions Ψ = {ψ1 , ψ2 , · · · , ψn } by solving:

LG ψi = λi ψi , i ∈ {0, 1, · · · , n} , (4)

416 where ψi are the n × 1 eigenvectors and λi are the corresponding eigenvalues.

417 Monte Carlo simulations

418 For statistical tests, we perform comparisons against a randomised connectivity graph. For randomization, we use an algorithm
419 described in Molloy et al. (1995)66 and adapt an implementation contained in the CNM toolbox (https://ch.mathworks.
420 com/matlabcentral/fileexchange/45734-cnm). The randomization is performed at the level of the sparse nearest
421 neighbor graph. First, all edges are ”cut”, producing vertices that have as many ”stubs” as they had connections. These
422 stubs are then randomly reconnected, effectively preserving the degree of each node. Random harmonic basis functions were
423 then obtained from the resulting nearest neighbour graphs as for the original FC matrix, as described above. We produced
424 (1/α − 1) · H randomizations, where α = 0.05 and H = 22 = 11 · 2 since we tested the first 11 functional harmonics shown in
425 Figure 2 and performed two separate sets of tests (average silhouette value over the whole cortex as well as maximum modified
426 silhouette value over all somatotopic regions; see below).

427 Silhouette values

To test whether isolines of the functional harmonics follow the boundaries of the parcels as defined in the HCP parcellation3 ,
we compute the silhouette value63 of each functional harmonic as:
1
S= (Mbetween (i) − Mwithin (i))/ max(Mbetween (i) − Mwithin (i)) , (5)
N∑i

428 where Mbetween (i) is the average Euclidean distance between vertices belonging to a parcel i and vertices belonging to all other
429 parcels, while Mwithin (i) is the average distance between vertices within the parcel i. If all vertices belonging to a parcel i have
430 the same value, and at least some vertices outside the parcel i have different values, then Mbetween (i) > 0, Mwithin (i) = 0 and
431 S(i) = 1. By averaging over the silhouette values of all parcels, one obtains a measure of how well the data fit the parcellation.
432 Note that we replaced the somatosensory/motor core areas 1, 2, 3a, 3b, and 4 with the somatotopic sub-areas given by the HCP3
433 for a more detailed evaluation.
434 To evaluate the somatotopic organization of the functional harmonics, we use a measure that was similar to the silhouette
435 value, but adapted to measure the separation from the rest of the cortex and from other somatotopic areas.

Ssom = maxsom ((Mbetween,som + Mbetween )/ max(Mbetween,som , Mbetween ) · Mbetween,som ) , (6)

436 where Mbetween,som is the average Euclidean distance between vertices belonging to a somatotopic area and all other vertices
437 belonging to all other somatotopic areas. The maximum is taken over the somatotopic areas. The first term of the equation is
438 between 1 and 2 and is close to 2 if both Mbetween,som and Mbetween are equal.

439 Task maps

440 We used group-averaged task activation maps provided with the S1200 group average data release of the HCP (see Table 1,
441 www.humanconnectome.org/study/hcp-young-adult/document/extensively-processed-fmri-data-documentation). Here we
442 provide a summary of the tasks that form part of the HCP task battery52 . There are 7 groups of tasks: working memory, motor,
443 gambling, language, social, emotional, relational. Subjects performed all tasks in two separate sessions (working memory,
444 gambling, and motor in the first session, language, social cognition, relational processing, and emotion processing in the
445 second).

14
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

446 Working memory. Four different stimulus types were used, presented in separate blocks: pictures of faces, places, tools and
447 body parts. Two different task types were used: a 2-back working memory task, where subjects had to respond if a stimulus
448 matched that two trials back, and a 0-back working memory task, where subjects had to respond whenever a single stimulus
449 returned that was presented at the beginning of the block. This results in a total of 19 different working memory task maps,
450 consisting of 14 activation maps (such as 0-back, 2-back, face, body, etc.) and 5 contrasts (between the two task types, between
451 each stimulus type and the average across all stimuli, etc.).
452 Motor. Visual cues indicated whether participants should move their left or right fingers, left or right toes, or move their
453 tongue. The goal was to identify the motor areas that correspond to these five body parts. This results in 26 different task maps
454 (7 activation maps for 5 body parts plus visual cue plus average, and 6 contrast maps).
455 Gambling. (Incentive processing.) Subjects played a game in which they could win or lose money. The game was to guess
456 whether the number on a ”mystery card” that could range between 1 and 9 would be less or more than 5. The numbers were
457 given after subjects made their guess and were chosen according to the trial type: ”win” - the number would correspond to their
458 guess and they would win 1$; ”neutral” - the number would equal 5 and they would neither win nor lose any money; ”loss”
459 - the number would not correspond to the guess and participants would lose $0.50. Separate blocks are used in which trials
460 are either mostly win or mostly lose, resulting in two conditions, punish and reward. This results in 3 different task maps (2
461 activation maps, i.e. one for each condition, and 1 contrast).
462 Language. Two different task types were used, ”story” and ”math”. ”Story” consisted of participants listening to 5-9
463 sentences of a story, and answering a 2-alternative forced choice question thereafter. ”Math” required participants to solve
464 simple addition and subtraction problems. The two task types are similar in terms of auditory input and attentional load, but
465 different in terms of semantic and numerosity related processing. As for gambling, the two task types result in 3 task maps (2
466 activation, 1 contrast).
467 Social. (Theory of Mind, TOM.) Subjects viewed videos of objects (squares, circles, triangles) that moved around in one of
468 two ways: ”Random” - there was no interaction between the objects, or ”TOM” - the objects moved as if they were reacting to
469 the other objects’ ”thoughts and feelings”. They then had to judge whether the objects were interacting or not, or respond with
470 ”not sure”. As with gambling and language, the two task types result in 3 task maps (2 activation, 1 contrast).
471 Emotional. Subjects viewed one of two types of stimuli, ”faces” or ”shapes”, and had to decide which of two stimuli
472 presented at the bottom of the screen matched the stimulus at the top of the screen. The faces included emotional stimuli, i.e.
473 angry or fearful expressions. Again, the two task types result in 3 task maps (2 activation, 1 contrast).
474 Relational. There were two conditions, ”match” and ”relational”. In all cases, stimuli can have one of six shapes combined
475 with one of six textures. In the ”match” condition, which served as a control condition, two shapes were presented at the top
476 and one at the bottom of the screen. A word (”shape” or ”texture”) that appears in the middle of the screen instructs subjects to
477 decide whether the bottom stimulus matches either of the top stimuli in the dimension indicated by the word. In the ”relational”
478 condition, two stimuli are presented each at the top and at the bottom of the screen, with no word in the middle. Instead,
479 participants have to determine themselves across which dimension the top pair differs, and, subsequently, indicate whether the
480 bottom pair differs over the same dimension. Again, the two task types result in 3 task maps (2 activation, 1 contrast).
481 Task maps were computed using FSL’s FEAT and FLAME67, 68 and conducting a between-subject (”level 2”) analysis,
482 resulting in effect sizes (Cohen’s d). We used the task maps with minimal smoothing (2mm total smoothing); see 1200 subjects
483 data release reference manual, pp. 45-54 and 100-104.

484 Reconstructing the task maps from functional harmonics

The spatial pattern of each task map on on the cortex s(v) was decomposed into and reconstructed from the functional harmonics
Ψ = {ψk }nk=1 as:
n
ŝ = α1 ψ1 + α2 ψ2 + · · · + αn ψn = ∑ αk ψk (v), (7)
k=1

where the coefficient αk of each functional harmonic ψk was estimated by projecting the task map ŝ(v) onto that particular
harmonic ψk . As such αk are estimated as:

αk = hŝ, ψk i . (8)

485 Then, each task map is reconstructed using Eq. 7. In this study, we limit our reconstructions to using a maximum of 100
486 non-constant functional harmonics (n = 101).

15
 bioRxiv preprint first posted online Jul. 11, 2019; doi: http://dx.doi.org/10.1101/699678. The copyright holder for this preprint (which
was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
It is made available under a CC-BY-NC-ND 4.0 International license.

For a reconstruction s∗ (m) , where m indicates a binary vector of dimensionality 101 × 1 which contains ones for harmonic
basis functions that are used in the reconstruction and zeros otherwise, we then compute the reconstruction error as:
r
RE(m) = ∑(si − s∗(m),i )2 / ∑ s2i (9)
i i

487 We also computed the Pearson correlations between s and s∗ (m) . For comparing the correlations between task maps and
488 reconstructions obtained from real functional harmonics versus randomized connectivity harmonics, we considered the number
489 of comparisons to be nC = nTasks · nLevels, where the number of tasks equals 47 and the number of levels refers to the different
490 numbers of harmonics used in the reconstructions, i.e. 0, 1, 2, 3, ..., 20, 30, 40, ..., 100, in 29 levels. From this we obtained a
491 corrected alpha level of αcorr = 0.05/nC, and we computed the critical value as Fisher’s z-transform of the correlation which a
492 sample has to exceed in order to be significantly higher than the random correlation:
q
zα · ( N11−3 + N21−3 )
zcrit = (10)
zrand
493 We obtain zcrit = 0.44, which corresponds to a minimum required empirical correlation of 0.41, with N1 = N2 = 59.412
494 (the number of vertices that contribute to the correlation values), zα = 0.438 (the inverse Student’s t distribution with
495 N1 = N2 = 59.412 degrees of freedom evaluated at 1 − αcorr ), and zrand = atanh(0.05) (Fisher’s z-transform of the maximal
496 random correlation between any reconstruction - with any number of functional harmonics - and any task).

497 Visualization
498 Somatotopic areas. In the visual and somatosensory/motor cortices, functional harmonics are rather determined by retinotopy
499 and somatotopy than by anatomical or microstructural features. For the former, somatotopic areas occupy exactly the same
500 surface area as the sensorimotor core areas, 1, 2, 3a, 3b, and 4. We therefore replaced, where appropriate, the borders of the
501 HCP parcellation by the borders of the five somatotopic regions.
502 Parcel borders for visualization. In order to discuss the meaning of the functional harmonics, we show borders of certain
503 parcels on the cortical surfaces (Figure 2). We used three different methods to select which borders to show. First, for some
504 functional harmonics, it was feasible to select these areas manually (for example, early visual areas in functional harmonic 4,
505 somatotopic areas in functional harmonics 3 and 4). The anatomical supplementary information from Glasser et al. (2016)3
506 uses a functional grouping of many regions that we often used as a guideline, for instance to distinguish between early and
507 association auditory cortex. Second, for some functional harmonics (for instance, functional harmonics 1 and 2), we show
508 the borders of parcels that belong to resting state networks as defined by Yeo et al. (2011)48 . The 7-network parcellation is
509 provided by the HCP, which does not perfectly overlap with the HCP parcellation. We adjusted the network borders slightly to
510 align the network borders to follow those of the parcels defined in HCP. Thereby we assigned each parcel to the RSN with
511 which it had the most overlap. Third, some functional harmonics are too complex to manually select areas or networks (namely,
512 functional harmonics 5, 6, 8, and 10). Here we employed simple k-means clustering on the functional harmonic, using k=2
513 (functional harmonics 5, 6, and 8) or k=3 (functional harmonic 10). To obtain meaningful clusters in the somatosensory/motor
514 cortex, we again replaced the sensorimotor core regions 1, 2, 3a, 3b and 4 with the somatotopic areas. For this purpose, we used
515 vertices within the core regions and re-assigned them to the somatotopic areas based on their distances to the sub-area borders.

516 Author contributions statement

517 S. A. and K. G. designed the methodology and the analysis. J. P. and G. D. contributed to the design of the study. M. L. K.
518 contributed to design of the methodology and the statistical analysis. M. L. K. and P. H. aided in the interpretation of the results.
519 S. A., K. G. and M. L. K. wrote the manuscript. All authors reviewed the manuscript.

520 Data availability

521 All data generated in this study are available from the corresponding author upon reasonable request.

522 Code availability

523 All custom scripts used in this study are available from the corresponding author upon reasonable request.

524 Additional information

525 Competing financial interests The authors declare no competing financial interests.

16