Functional Parcellation of fMRI data using multistage k-means clustering

Purpose: Functional Magnetic Resonance Imaging (fMRI) data acquired through resting-state studies have been used to obtain information about the spontaneous activations inside the brain. One of the approaches for analysis and interpretation of resting-state fMRI data require spatially and functionally homogenous parcellation of the whole brain based on underlying temporal fluctuations. Clustering is often used to generate functional parcellation. However, major clustering algorithms, when used for fMRI data, have their limitations. Among commonly used parcellation schemes, a tradeoff exists between intra-cluster functional similarity and alignment with anatomical regions. Approach: In this work, we present a clustering algorithm for resting state and task fMRI data which is developed to obtain brain parcellations that show high structural and functional homogeneity. The clustering is performed by multistage binary k-means clustering algorithm designed specifically for the 4D fMRI data. The results from this multistage k-means algorithm show that by modifying and combining different algorithms, we can take advantage of the strengths of different techniques while overcoming their limitations. Results: The clustering output for resting state fMRI data using the multistage k-means approach is shown to be better than simple k-means or functional atlas in terms of spatial and functional homogeneity. The clusters also correspond to commonly identifiable brain networks. For task fMRI, the clustering output can identify primary and secondary activation regions and provide information about the varying hemodynamic response across different brain regions. Conclusion: The multistage k-means approach can provide functional parcellations of the brain using resting state fMRI data. The method is model-free and is data driven which can be applied to both resting state and task fMRI.


Introduction
Functional Magnetic Resonance Imaging (fMRI) is a non-invasive technique used to study the neural activation dynamics of the brain. A localized neuronal activation causes an increase in blood flow and blood oxygenation concentration in that region. The difference in the blood oxygenation causes a change in the local magnetization that can be detected in the MR s ignal [Bandettini et al., 1994;Ogawa et al., 1990]. Since fMRI measures blood oxygenation (BOLD: Blood Oxygenation Level Dependent) changes following a neuronal activation, it is an indirect measure of neuronal activity. fMRI studies can be broadly classified into two main types, task fMRI and resting-state fMRI. In task fMRI, the main aim is to identify regions of the brain corresponding to different physical and mental activities [Bandettini et al., 1992]. In rest fMRI, the main aim is to study the underlying spontaneous activation and communication between different brain regions [Biswal et al., 1995]. One of the methods to determine the communication among regions is by using the temporal correlation between the fMRI time series of different brain regions [Van Den Heuvel & Hulshoff, 2010]. Because of the large number of voxels in the brain acquired by fMRI, voxel-wise correlation calculation is highly computationally expensive and not practically feasible. To overcome the limitation of voxel wise correlation computation, earlier analyses used a seed-based approach, [Biswal et al., 1995] where a random seed-voxel is selected within the brain, and a region is formed by combining all the voxels that have a high temporal correlation with the seed voxel time series. Because the seed selection is arbitrary, the regions obtained in such a fashion are highly dependent on the initialization of the seed voxel. As an alternative, researchers perform clustering to obtain different brain regions. A region-wise analysis is then performed with predefined regions. One main challenge is to obtain homogenous brain parcellation from the resting state data. Several brain parcellation schemes have been discussed in the literature [Thirion et al., 2014], but each has its advantages and disadvantages. When resting-state fMRI was being developed, anatomical atlases were commonly used. It was believed that regions having anatomical similarity might also have functional similarity [Nieto-Castanon et al., 2003]. In this technique, the time series of all the voxels within a given anatomical region are averaged to obtain a representative time series for that region. In reality, some anatomical regions are small and may be functionally similar, while certain anatomically large regions, such as Brodmann regions 6, 11, 37 and 48, may not be functionally homogenous and have subdivided or overlapping functional regions. For those large anatomical regions, such a representative time series is not always an appropriate descriptor since the region itself is not functionally homogenous, and temporal variations within different subregions are lost due to averaging, thereby yielding inaccurate analysis. Some researchers also use data-driven techniques like Independent Component Analysis (ICA) [Beckmann et al., 2005;Calhoun et al., 2001a;Calhoun et al., 2001b], mixture models [Golland et al., 2007], graph-based approach [Shen et al., 2010], dictionary learning [Wang et al., 2016], etc. to obtain functionally homogenous parcellation of the brain. Sometimes, clusters obtained by such an approach may lack spatial homogeneity as the clusters may be scattered across the brain regions. To have spatially uniform parcels, functional atlases based on spatially constrained spectral clustering [Craddock et al.,2012], and seed-based region growing [Blumensath et al., 2013], have been developed. The main limitation of a global atlas is the inability to account for subject-level differences. Moreover, it is important to note that atlas-based approaches may have alignment and boundary issues with improper coregistration scheme.
To address these shortcomings, we propose a novel multistage k-means clustering algorithm that can be used to generate spatially and functionally homogenous brain regions from fMRI data. The multistage k-means approach is an unsupervised clustering approach. Because no mapping of regions from a separate atlas is required, the proposed approach avoids the problem caused by inappropriate coregistration. Moreover, the algorithm takes advantage of both spatial and functional constraints to obtain spatially and functionally homogenous parcels. The two of the most commonly used clustering algorithms are k-means clustering [Hartigan & Wong, 1979] and hierarchical clustering [Johnson, 1967]. Both k-means and hierarchical are unsupervised clustering algorithms, but when applied to fMRI data, they have their limitations. In k-means, we start by selecting random centroids within the dataset. Each data point is assigned to the cluster with the shortest distance between the data point and cluster centroid. After all the data points have been assigned to a cluster, new corresponding cluster centroids are computed, and the process is repeated until the change in cluster centroid position for all clusters is less than some predefined threshold. For fMRI data, the clustering results can easily become biased due to arbitrary cluster centroid initialization. In hierarchical clustering, we start by computing the pairwise distance matrix between every pair of data points. Initially, each data point is treated as a separate cluster.
Based on the pairwise distances, the closest data point pairs are merged until all data points finally belong to a single large cluster. The entire merging process is represented by a hierarchy tree known as the dendrogram. Different clusters are obtained by splitting the dendrogram at a different distance threshold. Because of the large number of voxels in the fMRI data, typically ranging from 50,000 to 200,000 voxels per brain volume, the calculation of the pairwise distance matrix for hierarchical clustering is computationally impractical.
The proposed clustering approach combines the advantage of both hierarchical and k-means clustering approach for fMRI dataset while reducing the limitations. Our method supports the notion of a hierarchical structure. Within each stage or level of the hierarchy, the clustering is done based on the k-means algorithm. The level-wise clustering avoids the computation of a full distance matrix, and the hierarchical structure gives the flexibility to get any desired number of clusters. Each k-means clustering step is performed to obtain at most 2 clusters. A small number of clusters reduces the bias of centroid initialization and yield stable clusters. The multistage kmeans approach is also not limited to resting-state fMRI. It can be used with task fMRI data as well to get information about activation regions and the shape of the Hemodynamic Response Function (HRF). The advantage of a model-free approach in estimating the shape of the hemodynamic response eliminates the bias introduced because of an HRF model. Moreover, different clusters in task fMRI data can be used to investigate more complicated issues like the variability of the hemodynamic response [Handwerker et al., 2004;Handwerker et al., 2012] and the whole-brain activation [Gonzalez-Castillo et al., 2012].

Multistage k-means Algorithm
The multistage k-means clustering algorithm is used for functional brain parcellation using fMRI data. The clustering is performed on the preprocessed 4D fMRI time series data. Before applying the clustering algorithm to fMRI data, the preprocessed 4D fMRI data structure is modified. The preprocessed data is vectorized into a 2D matrix whose rows correspond to voxels and whose columns correspond to time points. The original 4D fMRI volume is converted to an N x T matrix, where N is the total number of voxels inside the brain region, and T is the number of time points.
Thus, each of the N voxel's time series is treated as a T dimensional feature vector.
All feature vectors start as a single parent cluster. Before clustering is begun, the desired number of hierarchical stages and a correlation threshold is specified, denoted by 'NS' and 'CT', respectively. The parent cluster is then split into two children clusters using k-means. The k-means split also yields cluster centroids. For fMRI data, the centroids are T-dimensional and are considered as representative time series for the clusters. Next, a correlation is computed between the representative time series of the children clusters. If the correlation between the representative time series of the children clusters is greater than the threshold specified, it is assumed that the child clusters are a part of the same cluster and further clustering of the parent cluster is stopped.
The parent cluster is considered as a converged cluster and saved for further use. Such correlation stopping condition prevents the subdivision of a single cluster into multiple sub-clusters. However, if the correlation between children clusters is less than the threshold, then each of the children clusters is considered as an independent parent cluster in the next level of the hierarchy. This process is repeated until all the parent clusters have converged or the last stage of the hierarchy is reached. Figure 1 shows the pseudocode for the multistage k-means algorithm. A dummy example for multistage k-means is shown in Appendix A1.
In k-means clustering, the most common distance metric is the Euclidean distance. In this approach, however, a correlation-based distance metric has been used. Because for functional parcellation, the aim is to combine voxels having similar temporal fluctuation, a correlation distance metric is preferred over Euclidean. Moreover, the maximum number of iterations for each k-means split is set to 100 with 5 replicates. The iteration corresponds to a single centroid update in the k-means algorithm while replicates correspond to repeating the entire clustering processing.
The centroid initialization is done using the k-means++ algorithm of MATLAB. At the end of the last stage of the hierarchy, pairwise temporal correlation is computed between centroids of all the converged clusters. Any cluster pair for which the correlation coefficient is higher than the correlation threshold specified earlier is merged.

Similarity Index Measure
A measurable metric to test the similarity between two different spatial clusters was used by Mezer et al. [Mezer et al., 2009]. This metric is called the similarity index analysis. To compare the similarity between spatial clusters obtained by the multistage k-means approach to existing functional and anatomical atlases, we used a slightly modified form of the similarity index analysis. The mathematical form to compute the cluster wise similarity index ( ) and the maximum similarity index (mSI) between two different cluster sets α, and β is given in equations (1) and (2): Here, i = 1, 2, …, Nα and j = 1, 2, …, Nβ, with Nα and Nβ are the maximum number of clusters in α and β, respectively. αi and βj are the number of elements in the i th and j th clusters of α and β, respectively. αi ∩ βj is the number of elements in αi that overlap with elements of βj. For two identical datasets, the mSI value is 1. Moreover, for two datasets with a different number of clusters, the mSI value is always less than 1. Thus, if Nα ≠ Nβ, mSI < 1.
Because the ordering of the cluster number may be different in different sets, the first step is to find the closest cluster in set β corresponding to each cluster in set α. For each cluster in α, similarity index is computed with every other cluster in β. The cluster in β having the highest similarity index measure for cluster αi is considered the closest cluster, and the similarity index measure is assigned to that particular cluster. This process is represented by equation (1). For every cluster in α, the closest cluster is detected, and the similarity index measure is computed. To compute the similarity index between different atlases and spatial maps, it is necessary that the physical scaling of both the spatial maps be the same. To make sure that the spatial maps being compared are of the same size, a 3D affine transformation is applied using the 'Coregister

Dataset and Preprocessing
The performance of the multistage k-means was first validated on a synthetic time series data. The synthetic ground truth data consists of 8 x 8 2D grid with six distinct regions. For each region, a representative time series was obtained by sampling a sinusoid with 50-time points. The sampling frequency was 10 Hz, and the frequencies for each region were 0.65 Hz, 0.7 Hz, 0.7 Hz, 0.8 Hz, 0.5 Hz, and 0.9 Hz with a random phase. Entire synthetic time series data was created by assigning the corresponding representative time series to each 8 x 8 pixel and then adding Gaussian random noise with different variance. Figure 2A shows the 2-D ground truth regions a sample visualization of the synthetic time series data. Four different synthetic datasets were created using different noise variance to obtain a Signal to Noise Ratio (SNR) value of 2, 1, 0.5, and 0.25. The structure of the dataset was such to have spatiotemporal data similar to 4D fMRI.
Next, a simulated resting-state fMRI data was generated using f-Sim fMRI simulator [Parmar et al., 2018]. The simulated resting-state fMRI data was generated using a multivariate covariance matching scheme with Craddock atlas as the functional reference parcellation. The spatial and temporal properties of the simulated data were similar to an ordinary resting-state fMRI scan. The data were simulated for a 240 s long resting-state experiment with a Repetition Time (TR) of 2 s yielding 120 brain volumes. Each brain volume consists of 48 axial slices of thickness 3.5 mm with an in-plane matrix size of 64 x 64 and a resolution of 3 mm x 3 mm.
The simulated data validation acts as a bridge between the synthetic and the real data. The simulated data has the structure of real fMRI data with an added advantage of having the ground truth to compare with. For the simulated data, the correlation threshold was set to 0.7 as the stopping criteria with a maximum of 7 stages in the hierarchy. The selection of correlation threshold was arbitrary based on the results from the synthetic data output. The number of stages was set to 7, which may yield a maximum of 128 clusters. The simulated data was generated using the Craddock atlas with 50 parcellations. Thus, the number of stages were set to a power of 2, which can produce more than 50 clusters but not the immediate next. Hence NS was set to 7 and not 6.  subjects, the performance of the multistage k-means clustering output was compared with simple k-means, ICA clustering and Craddock functional atlas using student's t-test.
The task fMRI dataset was obtained from the 100 functional runs per subject dataset collected at the National Institute of Health (NIH) by Dr Javier Gonzalez-Castillo fluctuations are less than 2%-5% the absolute signal intensity and thus global signal regression is an essential preprocessing step. Lastly, the voxels outside the brain region were thresholded and masked out. Masking reduces the effective number of voxels and hence the complexity. Moreover, for both task and resting-state fMRI data results shown in the paper, the multistage clustering was applied with CT = 0.7 and NS = 7 unless specified otherwise. The value for CT and NS were selected based on the results of simulated fMRI data.

Synthetic Data
First, the multistage k-means clustering algorithm was validated using the synthetic data.
Clustering output was computed for each of the four cases, different SNR values. The clustering output was compared with the ground truth using the mean similarity index measure. For the multistage k-means, the correlation threshold (CT) stopping criteria was set to 0.7. The ground truth matrix has 6 clusters, so the number of stages (NS) was set to 3, the next highest power of 2 closest to 6. The clustering output of the multistage k-means was also compared with the clustering output of simple k-means. For comparison to be uniform, the same number of clusters were used in simple and multistage k-means. To have a uniformity in the cluster numbers, the multistage approach was run first, and the number of clusters were determined. Then simple k-means was run to get the same number of clusters in the output. Figure 2C shows the clustering output for both the approach for different values of SNR and different combination of CT and NS. Different colors in the 2D matrix represent different clusters. The accuracy of the clustering was computed using the mean of the similarity index measure for all clusters. As expected, the clustering accuracy reduces with reducing SNR value. Also note, the mean similarity index for multistage k-means is always slightly better than simple k-means. Along with that, the number of clusters detected by the multistage approach increases with reducing SNR. With increasing noise, the intra-cluster correlation reduces and hence some of the clusters that converge for low noise variance fail to converge for a higher noise variance.
Next, the effect of NS was investigated by repeating the clustering with 4 stages in the hierarchy instead of 3. The clustering results for that is shown in Figure 2D. For lower noise variance, SNR = 2, as all the clusters converge before reaching the maximum number of stages, the clustering output remains the same. However, for larger noise variances, the clusters do not converge until the last stage of the hierarchy and thus an increase in the total number of clusters in observed. The increase in the total number of clusters also reduces the similarity index as compared to lower NS.
Nevertheless, for this case also, the multistage approach performs slightly better than the simple k-means.
Lastly, the effect of CT was computed using synthetic data. Clustering was performed for a higher, 0.85, and a lower, 0.4, value of CT keeping the NS as 3 is shown in Figure 2E and 2F. For a lower value of CT, there is a higher probability that the clusters will converge because of the low threshold. Higher convergence probability results in a lower number of total clusters obtained. It can be noted from the results that even for a lower SNR, 0.25, clusters converge, and the total number of clusters obtained is less than the maximum number of cluster possible. One more thing to note is that for a lower CT, simple k-means performs better than the multistage counterpart. On the other hand, a higher value of CT means fewer clusters will converge because of a tighter threshold. Thus, even for high SNR, 2, clusters will fail to converge and will result in a higher number of clusters.

Simulated fMRI Data
Next, the performance of multistage clustering was validated using simulated resting-state fMRI data. The total number of clusters detected by the multistage approach was 54. Thus, the number of clusters obtained was closer to the actual ground truth (50) than the maximum number of clusters possible (128). The number of clusters being closer to 50 rather than 128 indicates that the majority of the clusters converge before reaching the final stage of the hierarchy suggesting stable clustering. Figure 3 shows the spatial cluster comparison for the output of multistage clustering and the Craddock atlas with 50 parcels. Each cluster is represented by a different color. The color- bar represents the corresponding cluster number. The maximum similarity index measure, mSI, is also computed between the clustering output and the Craddock atlas. A mSI value of 0.771 is obtained that indicates high spatial similarity, which is also observed in the figure. The results obtained on synthetic and simulated time series data show the potential of the multistage k-means algorithm for functional parcellation.

Resting-State fMRI Data
The lack of ground truth for resting-state fMRI data makes a quantitative analysis challenging.
Here, the performance of multistage clustering is compared with the output of data driven Functional ROIs and representative time series for each ROI were obtained for each of the Independent Components (ICs). In the case of ICA analysis, the time series of each voxel is assumed to be obtained by a weighted combination of different ICs. Thus, a given voxel can belong to mode than one functional ROI. For simplicity, to have a comparable clustering output, each voxel was assigned to the IC with the maximum contribution to that voxel, hard clustering. All the spatial analysis of the ICA analysis were performed using the hard clustered ICA spatial maps.
The spatial comparison was performed using the maximum similarity index (mSI) measure. The mSI values, mean ± standard deviation across all 16 subjects, between different data driven and atlases is shown in Figure 4. Along with the mSI values, Figure 4 also shows the t-value for all the cases. A student's t test was performed to test if there was a significant difference in the mSI values when using the multistage k-means clustering algorithm as opposed to the other. Cases for which the mSI value for multistage k-means is significantly higher (p < 0.05) than the other atlas is highlighted in blue. The comparison with both functional and anatomical atlas gives an insight into the structural homogeneity of the parcels. If the clusters have very low similarity with both atlases, then it means the clusters are scattered without any uniform spatial structure. If the similarity of the clusters with anatomical atlas is much larger as compared to a functional atlas, then it might indicate that the parcellation is performed based on anatomical differences and not functional differences. A good functional parcellation scheme should have high similarity with both functional and anatomical atlas with a slightly higher similarity with the functional atlas. At first glance the similarity index values might appear very low but different number of clusters in both the clustering output being compared reduces the mSI values as discussed earlier.
Moreover, for reference, the maximum similarity index was computed between two synthetic clusters of size 16 x 16 x 16 with 8 randomly distributed spatial clusters each. The mean mSI value obtained for 100 random simulation was 0.025. Compared to the random simulation, the mSI values obtained for clustering outputs is much higher. The number of output clusters for both kmeans approach was 58, the number of ICs were 50 and thus, to have a better comparison, Craddock atlas with 60 parcels was used. The limitation with the anatomical atlas is the fixed number of parcels (100 ROIs inside the brain region); thus, a slightly lower mSI is expected for comparison with anatomical atlas as there exists a difference in the total number of clusters. The multistage approach performs better than simple k-means and ICA in terms of mSI values with both anatomical and functional atlases. Even when comparing with the data driven clusters (ICA output), the multistage approach performs better than simple k-means. Craddock atlas' similarity is higher than the similarity of the multistage clustering output for anatomical atlas. The lowest mSI value is obtained between ICA clusters and anatomical atlas. As discussed earlier, data driven techniques like ICA yields cluster that have high functional similarity compared to structural similarity. The t-test result also suggests that using multistage approach, a significant improvement is observed in structural similarity of the output clusters as compared to other clustering approach and data driven techniques.
The spatial clustering output for simple and multistage k-means and ICA is shown in Figure 5.
Each cluster is represented by a different color. The color-bar indicated the color corresponding to different cluster number. Visually, the clusters for both the sets of k-means look spatially uniform while the ICA clusters look little scattered. One thing to observe for multistage k-means clusters is a visually apparent difference between the gray matter and the white matter which is less visible in the simple k-means clusters or ICA. To check if the clusters belong strictly to gray matter or overlap between gray matter and white matter, percentage overlap was computed for all the clusters to gray matter and white matter segmentation mask. Figure 6 shows the percentage overlap for each cluster to gray matter and white matter in sorted order for multistage & simple k-means and ICA. For multistage k-means, nearly half of the clusters (30 out of 62) have more than 70% overlap with gray matter whereas for simple k-means it is only one third (20 out of 62) and for ICA only one fourth (12 out of 50).
Moreover, the clusters that overlap significantly with both gray matter and white matter, i.e. gray matter overlap between 40% and 60%, are almost twice for both simple k-means (15 clusters) and ICA (14 clusters) as compared to multistage (8 clusters).
The overlap results show that for multistage k-means, the parcellations are made based on functional information. For simple k-means, as large number of centroids are initialized at once, more weight is given to spatial similarity. As opposed to that, because parcellations are performed in a hierarchical fashion based on temporal correlation between centroids, the clusters obtained by multistage k-means are both functionally and spatially homogenous. For ICA, however, because a hard clustered spatial map was used, clusters are scattered across the brain as seen in Figure 5. The spatial homogeneity can be quantified using the mSI. To quantify the functional similarity, mean intra-cluster correlation histogram was computed. For each cluster, the correlation coefficient is computed between the time series of each voxel with the representative time series.
The mean of all the correlation value is computed. A mean correlation value of close to 1 suggests high functional homogeneity within the cluster, while that of close to 0 suggests poor functional homogeneity. In reality, no region of the brain is totally homogenous at the cluster scale in discussion, thus a mean correlation value of 1 is practically not possible. The mean intra-cluster correlation computation is repeated for all the clusters, and a histogram is plotted using all the mean intra-cluster correlation values. Figure 7 shows the mean correlation histogram for multistage k-means, simple k-means, ICA clusters, and Craddock atlas for a single subject. The maximum functional similarity is observed in the ICA clusters, which is expected. One more observation is the wider range of histogram for ICA as compared to others. Next is the multistage approach which performs better than simple k-means and atlas-based approach. Finally, the simple k-means which shows better functional similarity than the atlas based parcellation. The atlas-based approach does not account for subject-level differences, which might be the reason for lower intracluster functional homogeneity. Figure 8 shows the spatial visualization for some of the clusters that correspond to some of the standard resting-state brain networks. The visualization shows the clusters overlaid on the anatomical scan for one of the randomly selected subject. The most commonly identifiable networks from the figure include visual, default mode, primary motor, auditory, frontal, executive, sensorimotor and cerebellar network. Figure 8: Extracted resting state brain networks using multistage k-means approach. The brain networks (red) are overlaid on the anatomical image (grayscale). The different brain networks shown here are (left to right, top to bottom) visual, default mode, primary motor, auditory, frontal, executive, sensorimotor and cerebellum.

Task fMRI Data
Task fMRI analysis is mainly done for identification of brain regions responsible for different tasks. Multistage k-means clustering can be applied to task fMRI data to obtain the activation region through clustering. For task fMRI, the voxels inside the activation region would have similar temporal fluctuation and hence would be grouped in a single cluster. It has been shown that functional clustering of noise reduced task fMRI data can be used to obtain information about the region having time-locked temporal fluctuation with the experimental paradigm [Gonzalez-Castillo et al., 2012]. Apart from the primary activation region, many brain regions have temporal changes like transients and negative response, which are influenced by the experimental paradigm.
Such regions can be defined as secondary activation region. Multistage k-means clustering can be used to identify primary and secondary activation regions and the corresponding hemodynamic response.
The multistage approach identified a total of 32 clusters for the task fMRI data. Figure 9 shows the spatial clusters obtained from the multistage k-means approach. The spatial clusters for task and resting state are different. A possible reason for that is the difference in the temporal BOLD response for task and resting state fMRI. In resting state fMRI, the BOLD response mainly consists of the spontaneous fluctuation inside different brain regions. In task fMRI, the regions that respond to the task have a task specific BOLD response. representative time series of each cluster is used as the target time series to fit the regression model. 25 The weights corresponding to each regressor correspond to the amplitude of the HR at that time point.      Figure 11: Spatial visualization of task fMRI clusters with corresponding hemodynamic response. The hemodynamic response for the brain regions involved in the task show a positive correlation with the experiment paradigm while the hemodynamic response in the default mode network, which is usually considered a task negative network, shows a negative correlation with the experiment paradigm.

Discussion
Spatial parcellation of brain volume based on functional similarity is of prime importance in fMRI data, especially for resting-state studies. Numerous techniques ranging from atlas-based to datadriven techniques, are used for this task. Each has its limitations and advantages. If using the atlasbased techniques, appropriate mapping and registration algorithms are of prime importance,  Apart from the hyperparameters, the effect of spatial smoothing was also studied for resting state data. Structural and functional similarity was computed on fMRI data preprocessed with a smoothing kernel with different FWHM. Figure 12 shows  For task fMRI data, this clustering approach can also be used to explore the aspect of HR variability and whole-brain activation. The idea of whole-brain activation was introduced by Gonzales-Castillo et al. in 2012. In their paper, they showed that if the noise is suppressed by a significant level, which was achieved by averaging over 100 functional runs of the same experiment, activation is observed throughout the whole brain. The shape of the hemodynamic response, however, varies from region to region. Despite varying shape, all hemodynamic response shows a time-locked behavior to the experimental paradigm.
The results shown for the task fMRI analysis were from the same dataset. However, as opposed to averaging over 100 functional runs, the results shown were for a single randomly selected functional run. Figure 13 shows

Conclusions
The effectiveness of multistage k-means approach can be inferred from the results discussed above. Functional clusters can be obtained which show high spatial and functional homogeneity. time series corresponding to different spatial clusters in task fMRI data are used to obtain hemodynamic response for different regions. Using this technique, we were also able to replicate the hemodynamic response corresponding to whole brain activation. As a future work, we wish to expand the utility of multistage clustering algorithm for group level studies. For task fMRI data, we plan to obtain hemodynamic response for different subjects, and investigate inter-subject hemodynamic response variability along with variability in the spatial region. Moreover, the multistage k-means technique itself is not limited to fMRI data. With some minor modifications this technique can be applied to any unsupervised clustering problem. Figure 14 also shows the block diagram of multistage k-means clustering applied to a simple example. Here, the clustering approach is used with 3 stages of hierarchy, i.e. NS = 3. Consider that A is the collection of all feature vectors. In the first stage, k-means is used to split A into 2 clusters, B and C. Once that is done, the similarity is checked between B and C. The similarity is less than the stopping criteria, so we proceed to the next stage. In the second stage, again k-means is used to split both B and C into 2 clusters each. B is split into D and E while C is split into F and G. A high similarity is observed between D and E, shown in green, and hence it is likely that they are the subdivision of the same cluster, B in this case. Hence, B is considered as one independent cluster, shown in yellow. For F and G, a high similarity is not observed; thus, we move to the third stage. In the third stage F is split into H and I, and G is split into J and K. Again, a high similarity is observed between J and K, so G is considered as an independent cluster. Finally, as the maximum number of stages in the hierarchy is reached, H and I are considered as independent clusters. Thus, we end with 4 clusters from A, which are B, G, H and I.

Appendix A2: Timing and Memory requirement comparison
The timing and memory requirement performance of the multistage algorithm was analyzed on synthetic data and compared to both hierarchical and simple k-means clustering approach. First, the timing complexity is analyzed analytically. The operation that is repeated the greatest number of times during a clustering algorithm is the distance computation between pair of data points.
During each iteration of simple k-means, distance is computed between data points and each to 200,000 voxels thus the above assumption holds true and the multistage approach requires lesser computation time than the hierarchical approach. To have lower computation time than simple k-means, multistage approach needs to satisfy the condition: 2 NS<k. For the configuration mentioned above, for any k > 14, the multistage approach will perform faster than simple k-means which is reasonable assumption and most of the time the number of functional parcels is greater than twice the number of stages in hierarchy. Moreover, the space (memory) complexity for k-means is O (n) while for hierarchical clustering is O (n 2 ).           The hemodynamic response for the brain regions involved in the task show a positive correlation with the experiment paradigm while the hemodynamic response in the default mode network, which is usually considered a task negative network, shows a negative correlation with the experiment paradigm.