Structural MRI Based Parcellation of The Human Brain Using Spatial Hierarchical Clustering


 Parcellating the human brain into neurobiologically meaningful regions serves two (key) purposes , i.e. to gain neurobiological insight on the modules of the brain and to achieve dimensionality reduction for large -scale studies . Data- driven parcellation approaches based on magnetic resonance imaging ( MRI ) mainly differ from each other by the modalities that they are derived from and the parcellation algorithm. Here we extensively evaluate parcellation approaches which are based on one specific modality , i.e. grey matter volume deduced from T1- weighted structural MRI scans . The newly proposed spatial hierarchical variable clustering method SPARTACUS is employed as parcellation algorithm and its performance is compared with other spatial hierarchical agglomerative clustering ( SHAC ) algorithms as well as with spatial spectral clustering on both simulated data and a single- site structural MRI data set of older adults using two metrics , i.e. clustering quality and clustering stability across subsamples . Our analysis reveals that SPARTACUS and Ward's SHAC outperform the other competing methods . Using spatial ensemble clustering methods we could further improve the quality of the parcellations . Performing parameter tuning of the granularity ( considering 2-1000 brain regions ) we identify multiple interesting numbers of brain regions , where the corresponding parcellations may reflect different levels of brain organization . Finally , we observe that our parcellations converge well above chance with other structural MRI based parcellations but hardly better or even worse than chance with parcellations from another modality .


Introduction 4
The human brain is clearly the most complex organ of the human body. In 5 order to get a better understanding on how the brain works, a crucial and 6 widely accepted concept is that the brain is organized into spatially contigu-7 ous, specialized brain regions (cortical areas and subcortical nuclei), which are 8 inter-connected by large-scale networks (Eickhoff et al, 2018a). These brain 9 regions should be of large within homogeneity and between heterogeneity with 10 respect to different neurobiological modalities, where the boundaries should 11 be consistent among different modalities (Eickhoff et al, 2018a). While the 12 first modalities are histological-based, e.g., investigating cyto-and myeloar-13 chitecture in postmortem brains, the development of high-quality magnetic 14 applying clustering methods to structural MRI data, where an extensive inves-68 tigation of the performance of SHAC methods in application to structural MRI 69 data is missing in the literature up till now. 70 Therefore, in this article, our goal is to derive brain parcellations using 71 SHAC methods based on one specific modality, namely T1-weighted structural 72 MRI data. For this, since the objects to be clustered are voxels, i.e. variables, 73 we propose a new spatially constrained hierarchical agglomerative variable 74 clustering procedure, referred to as SPARTACUS, introducing contiguity con- stable clusters for resting-state fMRI data. 89 A critical challenge for any parcellation procedure is to select the number 90 of brain regions. Since the human brain is organized into multiple levels, the 91 correct number of brain regions might not exist. Instead, it is more likely that This article is structured as follows. First, the relevant spatial clustering 104 methods and the methods to identify interesting numbers of brain regions 105 are presented. Subsequently, we evaluate and compare the performance of 106 these methods in a simulation study. Moreover, we apply these methods to 107 the 1000BRAINS data set including structural brain scans of older subjects.

108
As further quality feature, the convergence of the final brain parcellations 109 with existing anatomical atlases as well as alternative atlases generated by 110 (semi-)algorithmic approaches based on MRI data is analyzed. all 54 smaller clusters are of size 6 ⇥ 6 ⇥ 3. In the unbalanced setting, 9 larger 158 clusters are of size 3 ⇥ 6 ⇥ 6, 9 larger clusters are of size 6 ⇥ 6 ⇥ 6 and 9 larger 159 clusters are of size 9 ⇥ 6 ⇥ 6 which entails that 18 smaller clusters are of size 160 3 ⇥ 6 ⇥ 3, 18 smaller clusters are of size 6 ⇥ 6 ⇥ 3 and 18 smaller clusters are 161 of size 9 ⇥ 6 ⇥ 3. A visualization of the true parcellations is shown in Figure 1. Hereby, σ S is the correlation between voxels from the same smaller cluster 171 and σ L is the correlation between voxels that are in the same larger cluster 172 but not in the same smaller cluster. Moreover, the correlation between voxels 173 from different larger clusters is chosen to be zero.
where t 2{B, U} indicates the setting. In Section A.1 of the Online Resource In the beginning of any SHAC algorithm, each voxel forms its own cluster.

210
A sparse V ⇥ V dissimilarity matrix including all pairwise distances between 211 neighboring voxels is determined using a distance metric. We consider the 212 Euclidean distance and the squared correlation distance, which, for two voxels 213 x j and x`, is calculated as More precisely, let the columns of the data matrix X be centered, and under the constraints c T k c k =1 (Vigneau and Qannari, 2003) and that all vox-262 els belonging to cluster C k ,k =1,...,K, are spatially contiguous. Cov (x j , c k ) 263 is the covariance between x j and c k . Let further X k 2 R N ⇥C k be the matrix 264 which columns consist of the voxels x j 2 C k ,i.e.X k is the data matrix of clus- and Qannari (2003) show that the maximization criterion can be rewritten as is the first eigenvalue of cluster C k , i.e. the first eigenvalue of Furthermore, Vigneau and Qannari (2003) show that Hence, any cluster aggregation causes the criterion T to decrease. The idea 277 behind the SPARTACUS method is to merge at each iteration the two neighbor 278 clusters that cause the smallest decrease in T . Consequently, the SPARTACUS 279 method uses as distance measure between neighbor clusters C k and C m in a SHAC 281 algorithm. such as K-means. Moreover, (spatial) constraints can be easily incorporated.

289
The SSPEC algorithm we employ consists of three steps. First, a spatially 290 constrained adjacency matrix W = w j` j,`=1,...,V is calculated by parcellation. Combining these K smallest eigenvectors in a matrix F 2 R V ⇥K , 300 the K-means algorithm (Lloyd, 1982;MacQueen, 1967) is applied to parcellate 301 the rows of F into K clusters C 1 ,...,C K . Finally, voxel x j is assigned to 302 cluster C k , if and only if the j-th row of F is assigned to cluster C k .

303
Note that in order to make clustering decisions, SSPEC only uses the 304 distances between voxels that are neighbors or have a common neighbor. In contrast, the SHAC algorithms consider all pairwise distances between vox-306 els at least once over the run of the algorithm. I.e. the SHAC algorithms use 307 much more information in order to make clustering decisions than the SSPEC 308 algorithms.

323
Next, as base clustering method, a SHAC algorithm is applied to all subsam-324 ples to generate B base parcellations C (1) is fixed). This yields a homogeneous cluster ensemble where diversity is solely explained by the subsampling approach. Note that 327 multiple cluster ensembles for multiple numbers of clusters can be obtained SPARTACUS in short time, since, once the dendrograms are calculated for the subsamples, 329 they can be cut to give any number of clusters with low computational cost.

330
A pairwise similarity based approach is used as consensus function, i.e.

331
to obtain the final ensemble parcellation based on the cluster ensemble. Let Then, the ensemble distance between any two voxels is Using d E as distance measure between voxels, the final ensemble parcellation 336 with K clusters can be obtained by employing a SHAC algorithm with, e.g.,

337
single or average linkage.  between-cluster separation is the silhouette coefficient (SC) (Rousseeuw, 1987; where a j is the average distance of x j to all other voxels in C k and b j is the 350 minimum average distance of x j to all voxels from another cluster. As distance 351 measure we choose the absolute correlation based distance d absCorr ( The SC for C K is then given by 2.6.3 Spatial adaptation of (simplified) silhouette coefficient 368 Both, the SC and SSC ignore the spatial information provided by the data.

369
However, it is, e.g., known that brain regions in one hemisphere may interact  x j 2 C k to the centroids of the neighbor clusters of C k .  The usual approach to identify interesting numbers of clusters using internal is a dummy variable coding for SC, SSC, SC spatial or SSC spatial . Again, local 446 maxima of IVS may indicate interesting numbers of brain regions.

463
We determine the convergence of the 1000BRAINS based parcellations 464 derived in this manuscript with multiple anatomical and alternative brain SPARTACUS Ta b l e 1 Overview over brain atlases considered for comparison.   Table 1.

479
In order to compare the performance of SPARTACUS with the competing 480 methods, i.e. SHAC and SSPEC, and in order to investigate whether using 481 SEC as additional layer improves clustering quality, these methods are applied 482 to the simulation study as well as to the 1000BRAINS data set.   Thus, in total we consider eight different clustering methods, where we indicate 499 clustering methods that are applied to standardized data sets by method S , 500 e.g., SPARTACUS S .

501
In a first step, we apply all eight spatial clustering methods to all six 502 simulation scenarios (SimBS, SimBM, SimBW, SimUS, SimUM and SimUW) SPARTACUS Ta b l e 2 The mean over 25 ARI and SSC scores for each of eight clustering methods based on SimBW and SimUW, where each ARI value compares a predicted parcellation with K =27orK =54clusterswiththerespectiv etrueparcellationandeac hSSCscore evaluates the quality of a predicted parcellation on the training data. approach, the quality of each predicted parcellation is evaluated on the same 509 data set this parcellation is created on using the SSC. In both approaches, the 510 mean is taken over the 25 scores corresponding to the same simulation sce-511 nario, the same clustering method and the same number of clusters. The results presented in Table 3 (for SimBW and SimUW) and in Section     clustering quality approach we consider SSC spatial .

638
By looking at the results generated by the subsampling based clustering 639 quality approach in Figure 4 we observe that the SSC spatial curves progress  Interestingly, it is observed from these figures that there is some agree-651 ment between the @ @K SSC spatial and ARI curves. Consistently, maxima in 652 the ARI curves are accompanied by changes, e.g., elbow points, in the 653 @ @K SSC spatial curves. E.g., at K = 154 or K = 146 there is an elbow point 654 in the @ @K SSC spatial curve and a local maximum in the ARI curve based on 655 SPARTACUS S or SHAC S AL, corr , respectively. Therefore, at least three interesting numbers of brain regions can be deduced, i.e. K ⇡ 70, K ⇡ 150 32 SPARTACUS Fig. 5 The first derivative of the SSC spatial curve (blue) generated by the subsampling based clustering quality approach and the ARI curve (green) generated by the subsampling based clustering stability approach, both with K =2,...,1000 and using SPARTACUS S or SHAC S AL, corr as spatial clustering algorithm and K ⇡ 600. When only looking at the curves due to SHAC S AL, corr and 658 SHAC S AL, Eucl , K ⇡ 300 is another interesting number of brain regions.  to produce a few large and multiple smaller brain regions.  the simulation study, that SSPEC S tends to produce spherical shaped brain 703 regions of equal size. As expected, geometric clustering produces equally sized 704 brain regions that reflect the spatial structure of the brain. atlases with similar numbers of brain regions is quantified by the ARI.

719
However, some issues occur when quantifying the convergence between 720 atlases, i.e. the brain atlases neither are registered to the same brain template 721 in MNI space, nor have the same voxel resolution, nor have the same brain 722 coverage, nor have the same numbers of brain regions. Thus, we employ some 723 preprocessing steps in order to match atlases from two different atlas families.

724
These preprocessing steps are described in Section B.5 of the Online Resource.

725
After pairwise matching, those atlases from two matched families are paired 726 that have the most similar numbers of brain regions, but only, if the absolute 727 difference between their numbers of brain regions is less than or equal to 20.   modalities seem to capture different aspects of brain organization.

839
As the human brain is assumed to be organized in a multi level fashion, a 840 single true number of brain regions is unlikely to exist. In our analysis based ponents. Thus, our estimation K ⇡ 300 is associated with these estimations.

861
Finally, K ⇡ 600 is not established as interesting granularity in the literature.

862
These observations illustrates a major issue of granularity estimation, i.e.

863
interesting numbers of brain regions are variable with respect to, e.g., data 864 sets, modalities or resolutions. Thus we suggest to make granularity decisions 865 individually adapted to the respective data scenario. for the final prediction represented by the weights. 892 We showed that standardizing the voxels prior to clustering has a positive 893 effect on clustering quality. For future work it is interesting to investigate 894 whether standardizing subjects instead of voxels has a similar effect.