1 Introduction

Cortical parcellation, the division of the human cortex into distinct functional regions, is a rapidly evolving sub-field of neuroscience. The hierarchical organization of the brain, as well as the wide variety of imaging modalities and algorithms used for cortical parcellation render agreement on a ‘gold-standard’ parcellation impossible. There are however, features of parcellation schemes that are particularly desirable to clinicians and researchers: (i) Cortical parcellation can vary across individuals due to noise and therefore parcellation schemes that preserve individual differences while maintaining group-level correspondence are advantageous and desirable. (ii) Parcellation schemes that do not incorporate anatomical priors, i.e., data-driven schemes, are less likely to be biased by confounding information. (iii) Connectivity-driven approaches are expected to generate more fine-grained parcellations as they are directly driven by connectivity data.

In parallel, many previous studies have used measures of long-range white matter connections to parcellate the cortex [1]. However most of these have relied on coarsely discretizing the cortex using manually selected prior anatomical information. In contrast, the current work described here does not employ any neuroanatomical priors and discretizes the cortex at an unprecedented level of detail, resulting in unbiased, potentially more accurate, and more detailed parcellations. Further, this work demonstrates a simple multi-layer graph formulation that provides simultaneously subject-specific and group-consistent parcellations at varying degrees of granularity. Importantly, qualitative similarity between these parcellations and existing fMRI-based parcellations supports a strong link between structural and functional organizations of the human brain.

In this paper, we model the underlying structural geometry of the human cortex as a graph and harness the neuroanatomically meaningful spatial patterns of the modes or harmonics of the graph for cortical parcellation. Atasoy et al. [2] reported that the harmonics of the human structural connectome are predictive of the functional connectome. Note that several important differences in the present model, including but not limited to superior sampling density of structural connectivity, the weighting strategy used to construct the graph, and the multi-layer graph formulation, all render our model distinct from that of [2].

The orthonormal bases yielded by eigendecomposition of the Laplacian imposed on complex domains are in many instances particularly convenient for representing physical dynamics [3]. Here, we demonstrate that this is extensible to the human cortex by imposing the discrete Laplacian on a domain representing the structure of the human brain, the structural connectome graph (SCG). It has previously been shown that the eigenvectors of the discrete Laplacian (i.e., Laplacian matrix) imposed on a polygonal mesh representing the interface between white matter and the cortex (i.e., white matter surface) can be used to segment the cortex into its major distinct lobes [4]. This result shows that the harmonic modes of a fully surface-based graph model of the cortex are predictive of its large-scale functional organization.

2 Methods

Our parcellation workflow can be broadly divided into four major steps (Fig. 1). First, estimates of cortical geometry and structural connectivity are computed from the high resolution T1-weighted and diffusion MR scans of each subject. Next, the adjacency matrices of the group-average, multi-layer, and individual subject structural connectome graphs are populated. The Laplacian matrix of each graph is then calculated, and the first \(N_{\text {harm}}\) non-trivial eigenvectors of each Laplacian matrix are computed. Finally, hierarchical clustering of vertices based on their eigenprofile, or position in spectral space, is conducted at varying resolutions, i.e., number of clusters, \(N_{\text {clust}}\), for all graphs and with varying inter-layer connection weight \(\omega \) for the multi-layer graph. Resulting parcellations are then comprehensively evaluated, and numerical comparison of the cross-subject variability of clustering results from both the fully-individual and multi-layer formulations is conducted.

2.1 Data and Preprocessing

We used the preprocessed structural and diffusion data of the HCP minimal preprocessing pipeline [5] for 21 subjects from the HCP test-retest dataset [6], including 1.25 mm isotropic 3-shell 3T diffusion scans and the native-space 32k LR white matter surfaces.

Fig. 1.
figure 1

Overview of the proposed workflow for one subject. (a) Tractogram of 5 million long-range white matter fiber streamlines is generated from dMRI; (b) White matter surface mesh is generated from high resolution T1-weighted MRI; (c) Adjacency matrix of structural connectome graph (SCG) is constructed from (a) and (b); (d) \(N_{\text {harm}}\) harmonic modes, \(\{\lambda _k,\psi _k\}\), \(k \in \{0,...,N_{\text {harm}}\}\), of SCG are computed from Laplacian matrix of (c); (e) Cortical parcellation is computed from hierarchical clustering of white matter surface vertices based on spectral coordinates.

Full size image

2.2 Diffusion Tractography

All processing of diffusion data subsequent to the HCP minimal preprocessing pipeline was performed using tools contained in the open-source software package MRtrix3Footnote 1 [7]. For each subject, a tissue segmentation for use in anatomically constrained response estimation, FOD computation, and tractography was derived from the registered, high resolution T1 images using FSL’s FAST and FIRST algorithms. Next, the fiber orientation distribution was calculated from the diffusion signal and an estimated response function using the multi-shell, multi-tissue constrained spherical deconvolution algorithm with the default parameters. 50 million streamlines were then generated from the fiber orientation distribution using the Second-order Integration over Fiber Orientation Distributions (iFOD2) probabilistic tractography algorithm with 0.06 for the FOD cutoff amplitude and 250 mm as the maximum streamline length. The seeding and termination of tracks was constrained to the white matter/gray matter boundary, defined by FreeSurfer segmentation, using the “anatomically constrained tractography” option in the track generation algorithm. The tractogram of 50 million streamlines and the FOD were then inputted to the Spherical-deconvolution Informed Filtering of Tractograms (SIFT) algorithm to reduce the number of streamlines to 5 million. The spatial coordinates of the endpoints of each streamline were extracted.

2.3 Graphs

A graph G is defined by a collection of vertices V and the edges E between those vertices, \(E=\{e_{ij}|(v_i,v_j)\in (V\times V) \}\). Although the vertices V may be embedded in space in a particular way, the connectivity of G and all of its topological properties can be described by its adjacency matrix A, which is populated by

$$\begin{aligned} {A}_{i,j}= {A}_{j,i}= {\left\{ \begin{array}{ll} w_{i,j} &{} \text {if} \, (i,j) \in E,\\ 0, &{} \text {otherwise,} \end{array}\right. } \end{aligned}$$
(1)

where \(w_{i,j}\) is the weight of the connection between vertices i and j. Importantly, it has been shown that for a graph whose vertices are a sampling of a Riemannian manifold, in the limit of \(N_{\text {vert}} \rightarrow \infty \), the properties of the graph approach the analogous properties of the continuous manifold represented by the graph [8].

2.4 Structural Connectome Graph

By virtue of the cross-subject registration of the white matter surfaces and the invariance of graphs under spatial embedding of their vertices, only one adjacency matrix \({A}_{\text {mesh}} \) is necessary to represent the surface mesh of all subjects. Here, the spatial resolution of \(N_{\text {vert}}\) = 64,984 represents an unprecedented level of resolution for graph models of structural connectivity. \({A}_{\text {mesh}}\) is populated according to (1), where the edges are between vertices connected by the polygonal connections that define the mesh, and the edge weights \(w^{s}_{i,j}\) are all equal to 1.

White matter surfaces are constructed separately for each hemisphere; when both are represented in a single graph, there exists a physiologically erroneous disconnection between them. We address this problem by introducing a unique weighting strategy where an inter-hemisphere connection adjacency matrix \({A}_{\text {IC}}\) is constructed to connect all vertices lying on the midline structures of the surface of each hemisphere to their nearest neighboring vertex on the opposite hemisphere with edge weight \(w^{\text {IC}}_{i,j}\) = 1. The full surface adjacency matrix \({A}_{\text {surf}}\) is then calculated as \({A}_{\text {surf}} = {A}_{\text {mesh}}+{A}_{\text {IC}}\). Notably, the only other study [2] seeking to calculate connectome harmonics detailed no attempt to correct for the false disconnection between the hemispheres of white matter surface meshes, and thus inaccurately represented the geometry of the cortex.

In order to encode the connectivity information contained in each subject’s tractogram of 5 million white matter fiber streamlines in the structural connectome graph, each surface vertex is associated with its 60 nearest streamline endpoints within 2 mm, and the other endpoint of each of these identified streamlines is associated with its nearest neighboring surface vertex. Nearest-neighbor computation was conducted using SciKit-Learn’s kd-tree nearest neighbors algorithm [9].

The long-range fiber adjacency matrix \({A}^m_{\text {fiber}}\) for each subject m is formed according to (1), with edges given by \(w^f_{i,j} = N^{\text {conn}}_{i,j}\times \frac{N^{\text {nz}}_{\text {surf}}}{N^{\text {nz}}_{\text {fiber}}}\), where \(N^{\text {conn}}_{i,j}\) is the number of streamlines connecting surface vertices i and j, and \(N^{\text {nz}}_{\text {surf}}\) and \(N^{\text {nz}}_{\text {fiber}}\) are the number of nonzero elements in \({A}_{\text {surf}}\) and \({A}^m_{\text {fiber}}\) respectively. The magnitude of \(w^f_{i,j}\)—which is on average 0.3—relative to \(w^{s}_{i,j}\) and \(w^{\text {IC}}_{i,j}\) is justified by the longstanding observations that local coupling of the cortex (through short fibers and other short-range circuits) is generally stronger, denser, and more efficient than long-range coupling of the cortex [10]. The structural connectome graph of an individual subject m is described by its adjacency matrix \({A}^m\), given by \({A}^m={A}^m_{\text {fiber}}+{A}_{\text {surf}}\).

The Laplacian matrix L of graph G(A) and its \(N_{\text {harm}}\) harmonics, \(\psi _k\) (\(\lambda _1< \lambda _k < \lambda _{N_{\text {harm}}}\)) are calculated by

$$\begin{aligned} {L}=L({A}) = {D}^{1/2}{A}{D}^{1/2}, \quad {D}_{i,i}= \sum _{j=0}^{N_{\text {vert}}} {A}_{i,j}, \end{aligned}$$
(2)

and

$$\begin{aligned} {L} \psi _k=\lambda _k \psi _k,\, k \in \{0,...,N_{\text {harm}}\}. \end{aligned}$$
(3)

Here, \(N_{\text {harm}}=99\), as we observed that the first \(\approx \)100 harmonics are generally meaningful and hypothesize that this low-frequency sub-space has strong overlap with the space spanned by global functional organization. Individual structural connectome normal modes, \(\psi ^m_k\), were calculated from (2) and (3), with \({L}^m=L({A}^m)\), \(m \in \{1,...,N_{\text {sub}}\}\). Group average structural connectome harmonics, \(\overline{\psi }_k\), were calculated from (2) and (3), with \(\overline{{L}}=L(\overline{{A}})\), where \(\overline{{A}}= \frac{1}{N_{\text {sub}}} \sum _{m=1}^{N_{\text {sub}}} {A}^m\).

Multi-layer approaches have been increasingly employed in neuroimaging studies, and have been shown to be efficacious in producing subject-specific results with group-level correspondence [4, 11]. The multi-layer structural connectome graph is represented by its supra-adjacency matrix \(\varvec{{A}}\), a size (\(N_{\text {sub}}N_{\text {vert}}\,\times \,N_{\text {sub}}N_{\text {vert}}\)) block matrix whose diagonal blocks of size (\(N_{\text {vert}}\,\times \,N_{\text {vert}}\)) are the individual adjacency matrices \({A}^m\), and all other blocks are the identity matrix multiplied by an inter-layer weighting strength \(\omega \). As it is difficult to determine the optimal \(\omega \), we used 5 different values of \(\omega \) ranging from 0.5 (less group constraint) to 2.0 (larger group consistency). Multi-layer eigenvectors \(\varvec{\psi }_k\) of length \(N_{\text {vert}}\,\times \,N_{\text {sub}}\) (1,364,664) were computed from (2) and (3), with \(\varvec{{L}}=L(\varvec{{A}})\). Individual multi-layer eigenvectors \(\tilde{\psi }^m_k\), \(m \in \{1,...,N_{\text {sub}}\}\), are given by the length \(N_{\text {vert}}\) section of \(\varvec{\psi }_k\) corresponding to the position of \({A}^m\) within \(\varvec{{A}}\).

3 Parcellation via Clustering

We define the eigenprofile (spectral coordinates) of vertex i of the white matter surface mesh, \(v_i\), as \(x_i\), where \(x_i=[\psi _1(v_i),\psi _2(v_i),\ldots , \psi _{N_{\text {harm}}}(v_i)]\). Hierarchical clustering of vertices based on their spectral coordinates was performed using SciKit-Learn’s agglomerative clustering algorithm [9]. Clustering was carried out for each set of individual eigenvectors, \(\psi ^{m}\), for the group-average eigenvectors, \(\overline{\psi }\), and for the multi-layer eigenvectors \(\varvec{\psi }\) calculated using 5 different \(\omega \)’s, all at 4 different resolutions \(N_{\text {clust}}\)’s. In each case, vertices lying within the midbrain structures are colored black.

Fig. 2.
figure 2

Cortical parcellations with \(N_{\text {clust}}=80\), \(N_{\text {clust}}=150\), \(N_{\text {clust}}=250\), and \(N_{\text {clust}}=500\) obtained by hierarchical agglomerative clustering of the spectral coordinates defined by (a–d) a typical subject’s individual eigenvectors, \(\{\psi ^1_k|k=1,...,N_{\text {harm}}\}\), and (e–h) Group average harmonics, \(\{\overline{\psi }_k|k=1,...,N_{\text {harm}}\}\), color-mapped onto the group-average white matter surface from medial and lateral views.

Full size image

4 Results and Discussion

When \(\omega \) is smaller, inter-subject consistency of multi-layer parcellation results is lower, as expected [11]. The parcellations are more consistent in the center of each region, but less consistent on the border. Our results based on the harmonic modes of our SCG formulation are very similar to previous functional parcellation studies [12,13,14,15]. This is especially true for the large-scale brain functional networks, as Fig. 3 shows many regions that are similar in location and shape to widely accepted parcellations derived using functional connectivity [12, 13], cortical folding patterns [14], or combinations of modalities [15]. For example, it is established that the default mode network has two large components in the medial part of the brain (the medial prefrontal cortex and the precuneus/posterior cingulate cortex), which are both reflected in our parcellation results. Another interesting feature of our parcellation results is that the medial part of the cortex is generally parcellated into fewer, larger regions, whereas the lateral cortex is parcellated into many smaller regions. This could be due to the fact that the lateral cortex contains more complex and more folded gyri/sulci, while the medial cortex is relatively flat and less folded.

Fig. 3.
figure 3

Left: Cortical parcellations with \(N_{\text {clust}}=150\) for (a,c) \(\omega =2\) and (b,d) \(\omega =0.5\) color-mapped onto group-average white matter surface from different views for (a,b) Group-consensus multi-layer parcellations (vertices that do not share the same label across all subjects appear in gray) and (c,d) A typical subject’s mutli-layer individual parcellation computed from multi-layer eigenvectors \(\tilde{\psi }^m_k\)’s. Right: Plots of group-average parcel-wise dice score between individual multi-layer parcellations and the consensus parcellation, and between individual parcellations and group-average parcellation (labelled ‘Ind’) at varying resolutions.

Full size image
Fig. 4.
figure 4

Average between-subject normalized mutual information (MI, left) and adjusted rand index (ARI, right) for individual multi-layer parcellations (\(\omega =2.0,\ldots ,0.5\)), and fully individual parcellations (Ind) at varying \(N_{\text {clust}}\).

Full size image

At coarse resolutions, average (Fig. 2 a–d), individual (Fig. 2 e–h), and multi-layer parcellations (Fig. 3. left) are reflective of large scale functional networks. The finer parcellation of the medial prefrontal coretex, insula, and the precuneus are similar to those obtain via functional connectivity based parcellation [16]. At higher resolution, subdivisions of large areas resemble more detailed atlases in the literature. Our results (Fig. 3, left and Fig. 4) also demonstrate that the multi-layer formulation achieves group-wise consistency superior to the individual formulation. From (a) and (b) of Fig. 3, it is evident that high-level cognitive function-related areas (such as medial prefrontal cortex, lateral frontal-parietal association regions) show less consensus across subjects, indicating high inter-subject variability in fiber connections associated with these regions. Subject-specific, group-consistent parcellation is critical to future efforts seeking to conduct group comparisons of the functional or structural connectivity between two or more populations. With a predefined unified atlas or parcellation scheme, individual variability can lead to reduced statistical power and false negatives.

5 Conclusion

This work demonstrates that hierarchical clustering of the harmonic modes of our multi-layer graph model of the structural connectome produces meaningful parcellations of the cortex that preserve subject-specific differences while achieving high group-level consistency. Our highly scalable model can be used to address longstanding problems in network neuroscience on the individual and group levels.