Abstract A fundamental assumption in neuroscience is that brain function is constrained by its structural properties. This motivates the idea that the brain can be parcellated into functionally coherent regions based on anatomical connectivity patterns that capture how different areas are interconnected. Several studies have successfully implemented this idea in humans using diffusion weighted MRI, allowing parcellation to be conducted in vivo. Two distinct approaches to connectivity-based parcellation can be identified. The first uses the connection profiles of brain regions as a feature vector, and groups brain regions with similar connection profiles together. Alternatively, one may adopt a network perspective that aims to identify clusters of brain regions that show dense within-cluster and sparse between-cluster connectivity. In this paper, we introduce a probabilistic model for connectivity-based parcellation that unifies both approaches. Using the model we are able to obtain a parcellation of the human brain whose clusters may adhere to either interpretation. We find that parts of the connectome consistently cluster as densely connected components, while other parts consistently result in clusters with similar connections. Interestingly, the densely connected components consist predominantly of major cortical areas, while the clusters with similar connection profiles consist of regions that have previously been identified as the ‘rich club’; regions known for their integrative role in connectivity. Furthermore, the probabilistic model allows quantification of the uncertainty in cluster assignments. We show that, while most clusters are clearly delineated, some regions are more difficult to assign. These results indicate that care should be taken when interpreting connectivity-based parcellations obtained using alternative deterministic procedures.

Citation: Hinne M, Ekman M, Janssen RJ, Heskes T, van Gerven MAJ (2015) Probabilistic Clustering of the Human Connectome Identifies Communities and Hubs. PLoS ONE 10(1): e0117179. https://doi.org/10.1371/journal.pone.0117179 Academic Editor: Gaolang Gong, Beijing Normal University, CHINA Received: September 16, 2014; Accepted: December 14, 2014; Published: January 30, 2015 Copyright: © 2015 Hinne et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited Data Availability: All data are available from the DANS database under persistent identifier urn:nbn:nl:ui:13-j9q9-o2. Funding: The research was funded by the BrainGain Smart Mix Programme of the Netherlands Ministry of Economic Affairs and the Netherlands Ministry of Education, Culture and Science (http://www.senternovem.nl/smartmix). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction The brain can be described as a vast network of interconnected neurons. The acquisition and subsequent analysis of this network has spawned a discipline known as connectomics [1, 2]. Part of connectomics focusses on structural connectivity, which concerns the layout of physical tracts consisting of axonal fibers. At a macroscopic scale, structural connectivity is defined in terms of brain regions that consist of neuronal populations which are inter-connected via white matter fiber bundles. As neuronal activity is constrained by neuroanatomy, correctly identifying these structural connections aids in understanding how spatially remote regions of the brain cooperate [3, 4]. Both functional and structural connectivity have been shown to be relevant in clinical applications, for instance by characterizing the connectomes that correspond to neurological and/or psychological disorders [5–7]. Similarly, connectivity studies have become a useful aid in understanding cognition by elucidating which networks are related to particular functions [8–11]. An important application of structural connectivity is the delineation of functionally specialized clusters of brain regions based on their structural connectivity patterns. Importantly, connectivity-based parcellation can be obtained in vivo, contrasting it with other methods such as histological analysis of cytoarchitecture [12, 13]. Other non-invasive methods that delineate structural cortical boundaries exist, based on identification of major anatomical landmarks, but these approaches have been shown to be susceptible to large inter-subject variability [14]. Intuitively, a cluster consists of a set of regions that are more similar to other regions in the cluster than to regions outside of it, but this idea can be operationalized in different ways. Under the first interpretation, which we refer to as profile-based clustering, regions may be clustered based on similarity of their connectivity profiles (e.g. [15–19], or see [20] for a literature review). In other words, it assumes that regions in the same cluster connect with the same areas. Several studies have used this approach to parcellate regions of interest such as the frontal pole [21], posteromedial cortex [22], occipital lobes [23], cingulate cortex [24] and thalamus [25–27]. Note that in this interpretation of connectivity-based clusters, regions within a cluster are not necessarily mutually connected. Implicitly, this approach aims to find regions that similarly integrate information from other parts of the brain. Under the second interpretation, which we refer to as community-based clustering, parcellations are taken to consist of densely connected clusters that are only sparsely connected to regions outside the cluster. This approach is typically used for whole-brain parcellation [28–32]. Here, one implicitly assumes that clusters are structurally (and thus, indirectly, functionally) specialized and mostly interact with the other regions in their respective cluster. Either operationalization of connectivity-based clustering is applicable to the human brain. As a consequence, choosing either perspective may hinder finding meaningful clusters that only adhere to the other definition of a cluster. Furthermore, clustering behavior may not be uniform across the brain, and instead may be a mixture of profile-based and community-based clusters. In this study we introduce a probabilistic model that is able to parcellate structural connectivity and is sufficiently flexible to incorporate the different cluster interpretations. The model reveals community-based clusters in the human connectome, but also a set of nodes that are assigned to small profile-based clusters that function as hubs. These small clusters all contain regions that have been labeled the ‘rich club’ [33]. In addition, we find that the whole-brain parcellations which the model provides, are up to par with algorithms that have previously been used to parcellate connectivity. Moreover, in contrast to parcellations obtained with the other approaches, parcellations estimated by the probabilistic model can be obtained without the need to choose the number of clusters in advance. Finally, because the model is probabilistic, it is able to explicitly represent the uncertainty in the obtained parcellations. These visualizations show that the cluster assignment for particular regions is uncertain, indicating that care should be taken when interpreting connectivity-based parcellations obtained using alternative deterministic procedures.

Probabilistic model for connectivity-based clustering First, we describe a probabilistic model for clustering based on structural connectivity. Next, we extend this model to allow for direct estimation of cluster structure from probabilistic tractography data. The complete model is shown in Fig. 1A. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 1. Probabilistic model for connectivity-based parcellation. A. The streamline infinite relational model combines a forward model for streamline data S with an infinite relational model that allows estimation of the cluster assignment matrix Z as well as the connection probability matrix ρ. Hyperparameters {ξ, α, β, δ T , δ F } complete the model. B. The probabilistic model supports both profile-based clustering as well as community-based clustering. The top row shows simulated connection probability matrices that correspond to profile-based clustering (left) and community-based clustering (right). Example networks that correspond to these probabilities are shown at the bottom. https://doi.org/10.1371/journal.pone.0117179.g001 Stochastic block model Let us assume that each cluster has a distinct connectivity profile. In other words, each cluster has a set of parameters that govern its connectivity behavior. This is known as a stochastic block model (SBM) [34]. Within a cluster, all nodes (i.e. brain regions) use the same connectivity parameters. As a consequence, regions within a cluster are stochastically exchangeable. Formally, each node 1, …, N is assigned a cluster label, using the unobserved (latent) cluster assignment variable Z nk = 1 if node n is in cluster k, and 0 otherwise. Each node is assigned to exactly one cluster. Structural connectivity is described by a symmetric and binary adjacency matrix A. Because of the stochastic exchangeability assumption, the probability of a connection a ij depends solely on the clusters to which nodes i and j are assigned. The cluster connection probabilities between clusters a and b are given by ρ ab . The set of cluster connection probabilities is collected in the connection probability matrix ρ. In practice, ρ is unknown so we use a prior on ρ to reflect our assumptions about it. Here, this is a Beta prior that depends on two hyperparameters α and β, which model the probability of a connection or non-connection between different clusters, respectively. We assume an uninformative prior with α = β = 1. After observing the data, the posterior expectation for ρ will reflect cluster connection probabilities and therefore informs about the contributions of profile-based and community-based clusters [35]. As an example of how block models may represent different connectivity patterns, Fig. 1B shows the posterior expectations for the connection probability matrix using a toy network of twelve nodes distributed evenly over three clusters. Here, the network on the left shows how profile-based clusters may be captured by the off-diagonal weights of ρ, while the network on the right shows how strong weighs on the diagonal of ρ generate a traditional community-based network. SBM have seen widespread application in literature, ranging from discovery of roles in social networks [34] to identification of protein-protein interactions [36, 37]. In addition, several model variants have been introduced, such as approaches that deal with overlapping clusters [38] or SBM tailored to weighted networks [39]. However, these approaches assume that the number of clusters K is known a priori, which is frequently not the case with empirical data. A nonparametric extension was introduced to learn the number of clusters from data as well [40]. This is achieved by placing a prior distribution on the cluster assignment matrix Z. Doing so allows the model to accommodate a (potentially) infinitely large number of clusters, rather than needing to specify the number of clusters K beforehand. Specifically, we draw Z from a Chinese restaurant process (CRP) [41]. This distribution over partitions can be used to generate samples from in the following way. Consider nodes that are assigned to clusters one by one, as customers entering a restaurant and choosing a table to sit at. Each customer is assigned to a non-empty table k with probability m k N − 1 + ξ , with m k the number of customers already assigned to table k and with probability ξ N − 1 + ξ to an empty (new) table. Its concentration parameter ξ determines how likely it is for a customer to sit at an empty table, which affects the total number of tables with customers, i.e. the number of clusters. Using the CRP, the generative model is then given by (1) where we introduce the notation m i to indicate the ith row of a matrix M. The model is known as the infinite relational model (IRM) [35, 40, 42, 43]. The infinite relational model is easily generalized to encompass multiple (conditionally independent) subjects that share a parcellation [42] by changing the model into (2) where the superscript m indicates the subject index. Note that although the clustering Z is shared across subjects, the cluster-to-cluster connection probabilities are subject-specific. When the IRM is applied to binary connectivity data, we will refer to the model as the bIRM. Forward model for streamlines To infer cluster assignments, structural connectivity data must be provided in the form of a binary adjacency matrix A. This matrix can be obtained using probabilistic tractography [44]. Specifically, probabilistic tractography proceeds by drawing streamlines between brain regions based on local estimates of anisotropic diffusion. These streamlines are collected in a streamline count matrix S, reflecting the number of streamlines between pairs of brain regions. Thresholding of the streamline count matrix produces a binary matrix that reflects structural connectivity. However, as a threshold results in a point estimate, we make use of a probabilistic model that describes how a structural network generates the distributions of probabilistic streamlines that are obtained through tractography [45]. Using this forward model, a streamline threshold is no longer required. Ideally, probabilistic streamlines show a distribution that perfectly reflects the underlying structural connectivity. However, tractography is prone to noise and errors, in particular in the presence of kissing, splaying and crossing fibres [46]. Hence, we distinguish between true connections along which we expect to observe streamlines, and false connections that may occasionally display streamlines, but do not correspond to actual anatomical pathways. The probability of a streamline between a pair of regions is represented by the matrix X. Formally, the streamline probability vector x i for a particular region i is determined by a Dirichlet distribution with parameter δ T for true connections and δ F for false connections. According to these probabilities, streamlines are distributed amongst the target regions using a multinomial distribution. By integrating out the streamline probability vectors, we obtain the following forward model: (3) where DirMul(α) stands for the Dirichlet compound multinomial distribution with hyper-parameters α. Note that this formulation assumes undirected structural connectivity, but allows streamline counts s ij and s ji to be different. To estimate connectivity from streamline data, (3) is supplied with a prior on A. In the most straightforward case, this prior is uniform, i.e. P(A) ∝ 1 [45]. Alternatively, this uninformative prior may be replaced by the IRM (see Fig. 1A), expressing our assumption of clustering in the network. The interpretation of the integrated model is that most observed streamlines indicate a structural connection (although there is some noise in the tractography process), and regions in the same cluster share their connectivity preferences (although some exceptions are allowed). We refer to the combined model that operates on streamline data as sIRM. Details of the approximate inference algorithm used to compute posterior estimates of interest are provided in SI1. Visualizing clustering uncertainty In addition to the display of individual parcellations, such as the maximum a posteriori (MAP) estimate, the proposed method may be used to visualize uncertainty about a parcellation. However, as the number of clusters may vary within a clustering distribution, and cluster labels are arbitrary for each sample, such a visualization is not trivial. A representation that allows different samples to be compared is the cluster co-assignment matrix M = ZT Z. The expectation of M, i.e. E [ M ] = 1 S ∑ s = 1 S M s , with S the number of obtained samples (see SI1), describes the posterior probability that two regions are assigned to the same cluster. The probability of co-assignment is used to color a region i with a weighted color coding given by c i = ∑ j m i j c ^ j / ∑ j m i j , with c ^ j the color representation in the MAP estimate for region j. For example, suppose that the MAP estimate consists of two clusters, one colored red and one colored yellow. If a region i is co-assigned to a region in the red cluster in half the samples, and to a region in the yellow cluster in the other half of the samples, it will be colored orange. Evaluating parcellation quality Since a clustering ground truth is unavailable, we used the reproducibility of the parcellations as a indicator of parcellation quality [14, 47, 48]. We quantify the reproducibility as the similarity between the parcellations of different participants (or groups of participants), expressed using adjusted mutual information (AMI) [49]. The AMI measure differs from the more traditional normalized mutual information measure in that it compensates for possible bias as a result of different numbers of clusters per parcellation. The measure is defined as (4) with MI(Z 1 , Z 2 ) the mutual information between two clusterings, H(Z) the entropy of a clustering and E [ MI ( Z 1 , Z 2 ) ] the expected mutual information between two clusterings. For comparison, we also obtain parcellations using K-means, the canonical algorithm for profile-based clustering, and the Infomap algorithm [50], which relies on community-based clustering (see SI3). Both methods have been used before in the context of brain parcellation [18, 32]. Each of the comparison algorithms was applied to the empirical streamline distributions. For K-means and Infomap the streamline matrices were made symmetric, i.e. S′ = S + ST. The number of clusters for K-means and Infomap is fixed to be the same as for the MAP estimate of the sIRM approach. To apply the bIRM, we first obtained the MAP estimate of connectivity with a flat prior, i.e. P(A) ∝ 1, which was subsequently clustered using the IRM approach. For the group-level analysis, a group-level streamline matrix was created by summing the streamline counts of all ten participants per split. This matrix was provided as input for K-means and Infomap. The bIRM method uses the ten individual MAP estimates of the forward model as input, while the sIRM method uses the ten streamline matrices.

Discussion Human brain connectivity is shown to exhibit clustering according to two different principles. Multiple brain regions show community-based clustering where clusters are both spatially contiguous and densely intraconnected. These clusters are tied together by brain regions that reveal connectivity-based clustering. The latter contain only one or two nodes and do not ‘fit in’ with other clusters due to substantially different connection profiles. These clusters consist predominantly of the superior frontal gyri, the superior parietal gyri, the precunei, the thalami and the putamina. All these regions have previously been identified as members of the ‘rich club’ [33]. In addition, all regions except for the thalamus have been pointed out as being the most vulnerable and central regions in structural brain networks [51]. Since the thalamus can be considered to act as a relay station, it stands to reason that it connects different cortical clusters. The rich-club regions stand apart from the rest of the connectome, as they are not part of the community-like clusters. Rather, they are each assigned their own cluster or group together with a similar region. Intuitively, this is an appealing observation. The community-like clusters contain regions dedicated to specialized processing. Here, an abundance of local connectivity, required for extensive within-cluster communication, results in dense intra-cluster connectivity. Simultaneously, the signal from these clusters must be integrated and disseminated, which is presumably achieved via the rich club regions [52]. Note that although small clusters consisting of one or a few regions will not be community-based clusters by definition, because there are simply not enough possible internal edges, it is the finding that these regions are not assigned to bigger clusters containing other nodes that is interesting. This means that a substantial fraction of the connections of these regions is to several different clusters, which emphasizes their integrative role. Validating the empirical results from connectivity-based parcellation remains a difficult task [20]. Since ground truth is not available, reproducibility across subjects is often used as a proxy for parcellation quality [14, 47, 48]. Based on this criterion, we have shown that our approach performs well. Yet, by visual inspection of the parcellations (see S2 Fig., S4–S6) some pairs of regions appear grouped together that are not immediately obvious. For instance, in many parcellations we find that thalamus and putamen are assigned to the same cluster, as well as amygdala and hippocampus and bilateral superior frontal gyri. One may argue that in particular for the subcortical areas, these regions should appear in separate clusters as they have specialized functionality. However, since these parcellations are based on connectivity and these regions project similarly to the cortex, they are put together in a cluster. This is inherent to connectivity-based clustering and occurs for K-means and Infomap as well. In fact, of the algorithms we considered, sIRM appears to be the only method that consistently assigns many of the rich club regions to small clusters instead of agglomerating them into large clusters. Finally, we note that anatomical constraints such as enforcing that subcortical areas constitute singleton clusters may easily be added into the prior distribution. Previously used methods are limited to one particular clustering behavior. A notable exception is a recent study that reveals that the connectome of the C. elegans roundworm also consists of a number of densely intraconnected clusters that are integrated via a core cluster that strongly connects to each community [53]. As we described above, the human connectome reveals both densely connected clusters as well as disconnected clusters with very similar connectivity. Consequently, picking one perspective inevitably neglects part of the available structure, which is exemplified by the parcellations found using K-means and Infomap that assign the rich club regions to larger clusters. Although, in the case of Infomap, this results in more consistent parcellations, the more detailed picture of communities and integrating regions is lost. Estimating parcellation with the sIRM approach provides a number of additional interesting quantities apart from the parcellation itself. The cluster connection probabilities ρ act as a cluster-level estimation of connectivity, expressed in terms of the rich-club regions and major cortices. Region-to-region connectivity A is estimated alongside clustering, and is potentially less prone to noise due to the prior that encourages regions in the same cluster to connect similarly, unless the data provides strong evidence to the contrary. Furthermore, other connectivity-based parcellation strategies, such as K-means, which groups regions with similar connectivity, and Infomap, which optimizes for densely connected components, provide point estimates of a parcellation. Although these approaches have provided valuable insight in the organization of both structural and functional connectivity [32, 54, 55], they do not quantify the uncertainty in their results (this issue is also discussed in [48]). Instead, our proposed method characterizes the full posterior distribution of all variables involved and thus provides a richer representation of parcellation. We find that a number of regions show substantial uncertainty in their assignments to a cluster. This illustrates that point estimates should be used with care, as a number of regions could easily be assigned erroneously. Integrating the forward model for structural connectivity with the IRM as prior leads to qualitatively different parcellations than when the bIRM is applied post-hoc to the connectivity estimates from the forward model. This is visible from the different number of clusters that both pipelines provide (the combined approach results in roughly 50% more clusters). In addition, the clear distinction between the two kinds of clusters is only marginally visible in the bIRM parcellations. Note that regardless of these differences, the reproducibility of the parcellations is similar for both methods. The IRM is a nonparametric method, which is helpful since the number of clusters is not known a priori. Still, parameters remain that affect the resulting parcellation. In particular the parameters δ T and δ F that govern the estimated connectivity and concentration parameter ξ that affects the number of identified clusters (although we observe that its influence is drowned out by the contribution of the likelihood) may be estimated from the data using empirical Bayes [56]. However, this will incur substantial additional computation costs that must be overcome to arrive at an efficient model. There are a number of directions one could take the proposed approach. To start with, alternative generative models may be used instead. For instance, the forward model that now consists of a Dirichlet compound multinomial distribution, modeling the probabilities of streamlines, may be replaced by a Poisson model [57], modeling streamline counts instead. Additional analyses will be needed to identify which model best captures the underlying connectivity. Furthermore, alternative stochastic block models may be used, such as described in a recent study [35] that applies the IRM and two variants to resting-state functional MRI. Lastly, an interesting avenue to pursue is identifying what causes the uncertainty in the posterior distribution of parcellation. Presumably, this can be attributed to a large extent to noise in data acquisition and tractography. But aside from methodological reasons, there may be other causes that could provide insight in the functional organization of the brain. For instance, uncertainty may in some cases be a result of overlapping clusters [58, 59]. It is likely that regions that are part of multiple, overlapping clusters show higher uncertainty in their cluster assignments. This may be clarified by embedding the infinite relational model within a larger framework to infer clusters at different levels of a nested hierarchy [30, 60]. One would expect that those clusters that show overlap and uncertainty become merged at a higher level of the hierarchy.

Conclusion In this paper, we have described an approach for connectivity-based parcellation that encompasses both community-like clusters as well as clusters that have similar connections, without being densely intraconnected. We find that both kinds of clusters are represented in the human connectome, and that the division into these two types corresponds to previous knowledge of structural connectivity. The model is able to quantify which regions are difficult to assign to a cluster, and it learns the number of clusters from the data. Finally, it does not depend on thresholded connectivity, but derives connectivity simultaneously with the parcellation. We hope that connectivity-based parcellation based on probabilistic models such as the one presented here will help to better understand the structural organization principles of human brain networks.

Materials and Methods Ethics statement Twenty healthy volunteers were scanned after giving informed written consent in accordance with the guidelines of the local reviewing committee CMO Arnhem-Nijmegen. This study was approved by CMO Arnhem-Nijmegen (CMO 2001/095 and amendment “Imaging Human Cognition”). Data acquisition The data used here was previously described in [61]. A T1 structural scan and diffusion-weighted images were obtained using a Siemens Magnetom Trio 3T system at the Donders Centre for Cognitive Neuroimaging, Radboud University Nijmegen, The Netherlands. An optimized acquisition order described by [62] was used in the DWI protocol (voxel size 2.0 mm isotropic, matrix size 110×110, TR = 13000 ms, TE = 101 ms, 70 slices, 256 directions at b = 1500 s/mm2 and 24 images at b = 0). Preprocessing DWI data were preprocessed using FSL FDT [63], which consisted of motion correction, correction for eddy currents and estimation of the diffusion parameters. To obtain a measure of white-matter connectivity, we used FDT Probtrackx 2.0 [63, 64] using seed voxel to target voxel tracking. Structural scans were segmented using FAST [65] and FIRST [66] to generate seed and target voxels. Seed voxels were those voxels in the cortical gray matter mask with a non-zero white-matter partial volume estimate and the outermost voxels of the subcortical masks. The remainder of the cortical and subcortical voxels served as target voxels. All FSL preprocessing was performed using version 5.0 and default settings unless otherwise specified. In addition, streamlines were terminated once they hit the target mask. This prevents transitive connections being erroneously interpreted as direct connections. The cortical surface was reconstructed using Freesurfer version 5.0 [67] and network nodes were defined according to the atlas provided in [68] (see S1 Fig.) in combination with subcortical segmentation using FIRST [66]. This resulted in a total of 162 nodes. Seed and target voxels were assigned to the nearest vertex in the reconstructed surface and by extension to the corresponding region. All streamline counts between two regions were summed to construct the final streamline matrix.

Acknowledgments The authors gratefully acknowledge the support of the BrainGain Smart Mix Programme of the Netherlands Ministry of Economic Affairs and the Netherlands Ministry of Education, Culture and Science. The authors thank Erik van Oort and David Norris for the acquisition of the DWI data. Finally, the authors express their gratitude towards Morten Mørup for providing the IRM sampling software and his valuable comments on the manuscript.

Author Contributions Conceived and designed the experiments: MvG TH ME MH. Performed the experiments: MH. Analyzed the data: MH RJ. Contributed reagents/materials/analysis tools: MvG TH ME MH RJ. Wrote the paper: MH MvG.