Parcellating connectivity in spatial maps

Probabilistic model

Intuitively, we wish to find a parcellation z which identifies local regions, such that all elements in a region have the same connectivity “fingerprint.” Specifically, for any two parcels m and n, all pairwise connectivities between an element in parcel m and an element in parcel n should have a similar value. Our method uses the full distribution of all pairwise connectivities between two parcels, and finds a clustering for which this distribution is highly peaked. This makes our method much more robust than approaches which greedily merge similar clusters (Thirion et al., 2014; Blumensath et al., 2013) or define parcel edges where neighboring voxels differ (Thirion et al., 2006; Wig, Laumann & Petersen, 2014; Gordon et al., 2014). The goal of identifying modules with similar connectivity properties is conceptually similar to weighted stochastic block models (Aicher, Jacobs & Clauset, 2014), but it is unclear how these models could be extended to incorporate the spatial-connectivity constraint.

We would like to learn the number of regions automatically from data, and additionally impose the requirement that all regions must be spatially-connected. We can accomplish both goals more efficiently in a single framework, by using an infinite clustering prior on our parcellation z which simultaneously constrains regions to be spatially coherent and does not limit the number of possible clusters. Specifically, since the mere existence of a element (even with unknown connectivity properties) changes the spatial connectivity and thus affects the most likely clustering, we must employ a nonparametric prior which is not marginally invariant. Other Bayesian nonparametric models allow for spatial dependencies between datapoints, but the only class of CRPs which is not marginally invariant is the distance-dependent Chinese Restaurant Process (dd-CRP) (Blei & Frazier, 2011). Instead of directly sampling a label for each element, the dd-CRP prior assigns each element i a link to a neighboring element c_i. The actual parcel labels z(c) are then defined implicitly as the undirected connected components of the link graph. Intuitively, this allows for changes in the labels of many elements when a single connection c_i is modified, since it may break apart or merge together two large connected sets of elements. Additionally, this construction allows the model to search freely in the space of parcel links c, since every possible setting of the parcel links corresponds to a parcellation satisfying the spatial-coherence constraint.

Given a parcellation, we must then specify a generative model for the data matrix D. Analogous to the approach taken in stochastic block modeling (Aicher, Jacobs & Clauset, 2014), we model the connectivity between each pair of parcels as a separate distribution with latent parameters. To allow efficient collapsed sampling (see below), we utilize a Normal distribution for each set of connectivities between parcels, and the conjugate prior for the latent parameters.

Mathematically, our generative clustering model is: ${\displaystyle \mathbf{c}\sim \text{dd-CRP}\left(\alpha ,f\right)}$ ${\displaystyle A_{mn},{\sigma }_{mn}^2\sim \text{Normal-Inverse-}{\chi }^2\left({\mu }_0,{\kappa }_0,{\sigma }_0^2,{\nu }_0\right)}$ ${\displaystyle D_{ij}\sim \text{Normal}\left(A_{\mathbf{z}{\left(\mathbf{c}\right)}_i\mathbf{z}{\left(\mathbf{c}\right)}_j},{\sigma }_{\mathbf{z}{\left(\mathbf{c}\right)}_i\mathbf{z}{\left(\mathbf{c}\right)}_j}^2\right).}$

For N elements and K parcels: c is a vector of length N which defines the cluster links for all elements (producing a region labeling vector z(c) of length N, taking values from 1 to K); α and f are the scalar hyperparameter and N × N distance function defining the dd-CRP; A and σ² are the K × K connectivity strength and variance between regions; μ₀ and κ₀ are the scalar prior mean and precision for the connectivity strength; ${\sigma }_0^2$ and ν₀ are the scalar prior mean and precision for the connectivity variance; and D is the N × N observed connectivity between individual elements.

The probability of choosing a particular c_i in the dd-CRP is defined by a distance function f; we use f_ij = 1 if i and j are neighbors, and 0 otherwise, which guarantees that all clusters will be spatially connected. A hyperparameter α controls the probability that a voxel will choose to link to itself. Note that, due to our choice of distance function f, a random partition drawn from the dd-CRP can have many clusters even for α = 0, since elements are only locally connected.

The connectivity strength A_mn and variance ${\sigma }_{mn}^2$ between each pair of clusters m and n is given by a Normal-Inverse-χ² (NIχ²) distribution, and the connectivity D_ij between every element i in one region and j in the other is sampled based on this strength and variance. The conjugacy of the Normal-Inverse-χ² and Normal distributions allows us to collapse over A_mn and ${\sigma }_{mn}^2$ and sample only the clustering variables c_i. Empirically, we find that the only critical hyperparameter is the expected variance ${\sigma }_0^2$, with lower values encouraging parcels to be smaller (we set α = 10, μ₀ = 0, κ₀ = 0.0001, ν₀ = 1 for all experiments).

To allow the comparison of hyperparameter values between problems with the same number of elements (e.g., the functional and structural datasets), we normalize the input matrix D to have zero mean and unit variance. We then initialize the model using the Ward clustering (see below) with the most likely number of clusters under our model, and setting the links c to form a random spanning tree within each cluster.

In summary, we have introduced a novel connectivity clustering model which (a) uses the full distribution of connectivity properties to define the parcellation likelihood, and (b) employs an infinite clustering model which automatically chooses the number of parcels and enforces that parcels be spatially-connected.

Derivation of Gibbs sampling equations

To infer a maximum a posteriori (MAP) parcellation z based on the dd-CRP prior, we perform collapsed Gibbs sampling on the element links c. A link c_i for element i is drawn from (1)${\displaystyle p\left(c_i^{\left(new\right)}|{\mathbf{c}}_{-\mathbf{i}},D\right)\propto p\left(c_i^{\left(new\right)}\right)p\left(D|\mathbf{z}\left({\mathbf{c}}_{-\mathbf{i}}\cup c_i^{\left(new\right)}\right)\right)=p\left(c_i^{\left(new\right)}\right)p\left(D|{\mathbf{z}}^{\left(new\right)}\right)\propto \left\{\begin{array}{ll}\alpha \hfill & \text{if }{c}_i^{\left(new\right)}=i\hfill \\ 1\hfill & \text{else}\hfill \end{array}\right\}\prod _{k_1,k_2=1}^{\left|{\mathbf{z}}^{\left(new\right)}\right|}p\left(D_{\underset{k_1}{\overset{\left(new\right)}{z}},\underset{k_2}{\overset{\left(new\right)}{z}}}\right).}$ To compare the likelihood term for different choices of $c_i^{\left(new\right)}$, we first remove the current link c_i, giving the induced partition z(c_−i) (which may split a region). If we resample c_i to a self-loop or to a neighbor j that does not join two regions, the likelihood term is based on the partition z(c_−i) = z. Alternatively, c_i can be resampled to a neighbor j such that two regions K′ and K″ in z(c_−i) are merged into one region K in $\mathbf{z}\left({\mathbf{c}}_{-\mathbf{i}}\cup c_i^{\left(new\right)}\right)=\hat{\mathbf{z}}$. Numbering the regions so that z_i∈{1⋯(K−1), K′, K″} and ${\hat{z}}_i\in \left\{1\cdots \left(K-1\right),K\right\}$ gives (2)${\displaystyle \frac{p\left(D|\hat{\mathbf{z}}\right)}{p\left(D|\mathbf{z}\right)}=\frac{\prod _{k=1}^{K}p\left(D_{{\hat{z}}_k,{\hat{z}}_K}\right)\prod _{k=1}^{K-1}p\left(D_{{\hat{z}}_K,{\hat{z}}_k}\right)}{\prod _{k=1}^{K^{\prime }}p\left(D_{z_k,z_{K^{\prime }}}\right)\prod _{k=1}^{K^{″}}p\left(D_{z_k,z_{K^{″}}}\right)\prod _{k=1}^{K-1}p\left(D_{z_{K^{\prime }},z_k}\right)\prod _{k=1}^{K^{\prime }}p\left(D_{z_{K^{″}},z_k}\right)}.}$

Each term p(D_zm,zn) is a marginal likelihood of the NIχ² distribution, which can be computed in closed form (Murphy, 2007): ${\displaystyle p\left(D_{z_m,z_n}\right)=\frac{\Gamma \left({\nu }_{mn}/2\right)}{\Gamma \left({\nu }_0/2\right)}{\left(\frac{{\kappa }_0}{{\kappa }_{mn}}\right)}^{\frac{1}{2}}\frac{{\left({\nu }_0{\sigma }_0^2\right)}^{{\nu }_0/2}}{{\left({\nu }_{mn}{\sigma }_{mn}^2\right)}^{{\nu }_{mn}/2}}{\left(\pi \right)}^{-n/2}}$ ${\displaystyle L=\left|z_m\right|\left|z_n\right|{\kappa }_{mn}={\kappa }_0+L;\hspace{0.17em}{\nu }_{mn}={\nu }_0+L{\mu }_{mn}=\frac{{\kappa }_0{\mu }_0+L\bar{d}}{{\kappa }_{mn}}}$ ${\displaystyle \bar{d}=\frac{1}{L}\sum _{\genfrac{}{}{0ex}{}{i\in \underset{m}{z}}{j\in z_n}}{D}_{ij}s=\sum _{\genfrac{}{}{0ex}{}{i\in \underset{m}{z}}{j\in z_n}}{\left(D_{ij}-\bar{d}\right)}^2{\sigma }_{mn}^2=\frac{1}{{\nu }_{mn}}\left({\nu }_0{\sigma }_0^2+s+\frac{L{\kappa }_0}{{\kappa }_0+L}{\left({\mu }_0-\bar{d}\right)}^2\right).}$ Intuitively, Eq. (2) computes the probability of merging or splitting two regions at each step based on whether the connectivities between these regions’ elements and the rest of the regions are better fit by one distribution or two.

In practice, the time-consuming portion of each sampling iteration is computing the sum of squared deviations s. This can be made more efficient by computing the s values for the merged $\hat{\mathbf{z}}$ in closed form. Given that the connectivities D_K′ = {D_iK′}_i∈k between parcel k and K′ have sum of squares deviations s_K′ and mean ${\bar{d}}_{K^{\prime }}$, and similarly for K″, then the sum of squares s_K for the connectivities between parcel k and the merged parcel K (merging K′ and K″) is: ${\displaystyle s_K=\sum _{d\in \underset{\stackrel{\prime }{K}}{D}\cup D_{K^{″}}}{\left(d-\bar{d}\right)}^2=\left(\sum _{d\in \underset{\stackrel{\prime }{K}}{D}\cup D_{K^{″}}}{d}^2\right)-\left(\left|D_{K^{\prime }}\right|+\left|D_{K^{″}}\right|\right)\cdot {\left(\frac{\left|D_{K^{\prime }}\right|\cdot {\bar{d}}_{K^{\prime }}+\left|D_{K^{″}}\right|\cdot {\bar{d}}_{K^{″}}}{\left|D_{K^{\prime }}\right|+\left|D_{K^{″}}\right|}\right)}^2=\left(\sum _{d\in \underset{\stackrel{\prime }{K}}{D}\cup D_{K^{″}}}{d}^2\right)-\frac{|D_{K^{\prime }}{|}^2}{\left|D_{K^{\prime }}\right|+\left|D_{K^{″}}\right|}{\bar{d}}_{K^{\prime }}^2-\frac{|D_{K^{″}}{|}^2}{\left|D_{K^{\prime }}\right|+\left|D_{K^{″}}\right|}{\bar{d}}_{K^{″}}^2-2\frac{\left|D_{K^{\prime }}\right|\left|D_{K^{″}}\right|}{\left|D_{K^{\prime }}\right|+\left|D_{K^{″}}\right|}{\bar{d}}_{K^{\prime }}{\bar{d}}_{K^{″}}=\left(\sum _{d\in \underset{\stackrel{\prime }{K}}{D}}{d}^2-\left|D_{K^{\prime }}\right|{\bar{d}}_{K^{\prime }}^2\right)+\left(\sum _{d\in \underset{\stackrel{″}{K}}{D}}{d}^2-\left|D_{K^{″}}\right|{\bar{d}}_{K^{″}}^2\right)+\frac{\left|D_{K^{\prime }}\right|\left|D_{K^{″}}\right|}{\left|D_{K^{\prime }}\right|+\left|D_{K^{″}}\right|}\left({\bar{d}}_{K^{\prime }}^2+{\bar{d}}_{K^{″}}^2-2{\bar{d}}_{K^{\prime }}{\bar{d}}_{K^{″}}\right)=s_{K^{\prime }}+s_{K^{″}}+\frac{\left|D_{K^{\prime }}\right|\left|D_{K^{″}}\right|}{\left|D_{K^{\prime }}\right|+\left|D_{K^{″}}\right|}{\left({\bar{d}}_{K^{\prime }}-{\bar{d}}_{K^{″}}\right)}^2.}$

Comparison methods

In order to evaluate the performance of our model, we compared our results to those of four existing methods. All of them require computing a dissimilarity measure between the connectivity patterns of elements i and j. For a connectivity matrix D, (3)${\displaystyle W_{i,j}=\sqrt{\sum _{a\ne i,j}{\left(D_{i,a}-D_{j,a}\right)}^2+\sum _{a\ne i,j}{\left(D_{a,i}-D_{a,j}\right)}^2}.}$

“Local similarity” computes the edge dissimilarity W_i,j between each pair of neighboring elements, and then removes all edges above a given threshold. Here we set the threshold in order to obtain a desired number of clusters. This type of edge-finding approach has been used extensively for neuroimaging parcellation (Cohen et al., 2008; Wig, Laumann & Petersen, 2014; Gordon et al., 2014). Additionally, this is equivalent to using a spectral clustering approach (Thirion et al., 2006) if clustering in the embedding space is performed using single-linkage hierarchical clustering.

“Normalized cut” computes the edge similarity S_i,j = 1/W_i,j between each pair of neighboring elements, then runs the normalized cut algorithm of Shi & Malik (2000). This draws partitions between elements a and b when their edge similarity S_a,b is low relative to their similarities with other neighbors. Although computing the globally optimal normalized cut is NP-complete, an approximate solution can be found quickly by solving a generalized eigenvalue problem. This approach has been specifically applied to neuroimaging data (Craddock et al., 2012).

“Region growing” is based on the approach described in Blumensath et al. (2013). First, a set of seed points is selected which have high similarity to all their neighbors, since they are likely to be near the center of parcels. Seeds are then grown by iteratively adding neighboring elements with high similarity to the seed. Once every element has been assigned to a region, Ward clustering (see below) was used to cluster adjacent regions until the desired number of regions is reached.

“Ward clustering” requires computing W_i,j between all pairs of elements (not just neighboring elements). Elements are each initialized as a separate cluster, and neighboring clusters are merged based on Ward’s variance-minimizing linkage rule (Ward, 1963). This approach has been previously applied to neuroimaging data (Thirion et al., 2014; Eickhoff et al., 2011).

We also compared to random clusterings. Starting with each element in its own cluster, we iteratively picked a cluster uniformly at random and then merged it with a neighboring cluster (also picked uniformly at random from all neighbors). The process continued until the desired number of clusters remained.

Synthetic data

To generate synthetic connectivity data, we created three different parcellation patterns on an 18 by 18 grid (see Fig. 2), with the number of regions K = 5, 6, 9. Each element of the K × K connectivity matrix A was sampled from a standard normal distribution. For a given noise level σ, the connectivity value D_i,j between element i in cluster z_i and element j in cluster z_j was sampled from a normal distribution with mean A_zi,zj and standard deviation σ. This data matrix was then input to our method with ${\sigma }_0^2=0.01$, which returned the MAP solution after 30 passes through the elements (approximately 10,000 steps). Both our method and all comparison methods were run for 20 different synthetic datasets for each noise level σ and the results were averaged.

We also performed a supplementary experiment using a more challenging three-spiral dataset (Chang & Yeung, 2008). We generated the connectivity matrix as above, and defined elements to be spatially adjacent if they were consecutive along a spiral or adjacent between neighboring spirals. In addition to our standard initialization scheme using the Ward clustering with highest probability according to our model, we also considered initializations with fixed numbers of clusters derived from Ward clustering (K = 2, 10) or initializations in which the links c were chosen are random. The ${\sigma }_0^2$ hyperparameter was set to 0.01 as above, and the MAP solution was returned after 100 passes (or 1,000 passes for the random initialization).

Parcellations were evaluated by calculating their normalized mutual information (NMI) with the ground truth labeling. We calculate NMI as in Strehl & Ghosh (2002). This measure ranges from 0 to 1, and does not require any explicit “matching” between parcels. For N total elements, if z assigns n_h elements to cluster h, z_gt assigns $n_l^{gt}$ elements to cluster l, and n_h,l elements are assigned to cluster h by z and cluster l by z_gt, this is given by (4)${\displaystyle \text{NMI}\left(\mathbf{z},{\mathbf{z}}_{\mathbf{gt}}\right)=\frac{I\left(\mathbf{z},{\mathbf{z}}_{\mathbf{gt}}\right)}{\sqrt{H\left(\mathbf{z}\right)H\left({\mathbf{z}}_{\mathbf{gt}}\right)}}=\frac{\sum _h\sum _l{n}_{h,l}log\left(N{n}_{h,l}/\left(n_h{n}_l^{gt}\right)\right)}{\sqrt{\left(\sum _h{n}_{h}log\left(n_h/N\right)\right)\left(\sum _l{n}_l^{gt}log\left(n_l^{gt}/N\right)\right)}}.}$

Human brain functional data

We utilized group-averaged resting-state functional MRI correlation data from 468 subjects, provided by the Human Connectome Project’s 500 Subjects release (Van Essen et al., 2013). Using a specialized Siemens 3T “Connectome Skyra” scanner (Siemens AG, Berlin, Germany), data was collected during four 15-min runs, during which subjects fixated with their eyes open on a small cross-hair. A multiband sequence was used, allowing for acquisition of 2.0 mm isotropic voxels at a rate of 720 ms. Data for each subject was cleaned using motion regression and ICA + FIX denoising (Smith et al., 2013; Salimi-Khorshidi et al., 2014) and then combined across subjects using an approximate group-PCA method yielding the strongest 4,500 spatial eigenvectors (Smith et al., 2014). The symmetric 59,412 by 59,412 functional connectivity matrix D_a,b was computed as the correlation between the 4,500-dimensional eigenmaps of voxels a and b. For each of ${\sigma }_0^2=2\text{,}000,3\text{,}000,4\text{,}000,5\text{,}000$, we ran Gibbs Sampling for 10 passes (approximately 600,000 steps) to find the MAP solution. For comparison with individual subjects, we also computed functional connectivity matrices for the first 20 subjects with resting-state data in the 500 Subjects release.

The map of retinotopic regions in visual cortex was created by mapping the volume-based atlas from (Wang et al., 2014) onto the Human Connectome group-averaged surface.

Human brain structural data

We obtained diffusion MRI data for 10 subjects from the Human Connectome Project’s Q3 release (Van Essen et al., 2013). This data was collected on the specialized Skyra described above, using a multi-shell acquisition over 6 runs. Probabilistic tractgraphy was performed using FSL (Jenkinson et al., 2012), by estimating up to 3 crossing fibers with bedpostx (using gradient nonlinearities and a rician noise model) and then running probtrackx2 using the default parameters and distance correction. 2000 fibers were generated for each of the 1.7⋅10⁶ white-matter voxels, yielding 3.4⋅10⁹ total sampled tracks per subject (approximately 34 billion tracks in total). We assigned each of the endpoints to gray-matter voxels using the 32 k/hemisphere Conte69 registered standard mesh distributed for each subject, discarding the small number of tracks that did not have both endpoints in gray matter (e.g., cerebellar or spinal cord tracks). Since we are using distance correction, the weight of a track is set equal to its length. In order to account for imprecise tracking near the gray matter border, the weight of a track whose two endpoints are closest to voxels a and b is spread evenly across the connection between a and b, the connections between a and b’s neighbors, and the connections between a’s neighbors and b. Since the gray-matter mesh has a correspondence between subjects, we can compute the group-average number of tracks between every pair of voxels. Finally, since connectivity strengths are known to have a lognormal distribution (Markov et al., 2014), we define the symmetric 59,412 by 59,412 structural connectivity matrix D_a,b as the log group-averaged weight between voxels a and b. The hyperparameter ${\sigma }_0^2$ was set to 3,000, and Gibbs Sampling was run for 10 passes (approximately 600,000 steps) to find the MAP solution.

Human migration data

We used the February 2014 release of the 2007–2011 county-to-county U.S. migration flows from the U.S. Census Bureau American Community Survey (ACS). This dataset includes estimates of the number of annual movers from every county to every other county, as well as population estimates for each county. We restricted our analysis to the continential U.S. To reduce the influence of noisy measurements from small counties, we preprocessed the dataset by iteratively merging the lowest-population county with its lowest-population neighbor (within the same state) until all regions contained at least 10,000 residents. This process produced 2,594 regions which we continue to refer to as “counties” for simplicity, though 306 cover multiple low-population counties. For visualization of counties and states, we utilized the KML Cartographic Boundary Files provided by the U.S. Census Bureau (KML).

One major issue with analyzing this migration data is that counties have widely varying populations (even after the preprocessing above), making it difficult to compare the absolute number of movers between counties. We correct for this by normalizing the migration flows relative to chance flows driven purely by population. If we assume a chance distribution in which a random mover is found to be moving from county a to county b based purely on population, then the normalized flow matrix is (5)${\displaystyle D_{a,b}=\frac{M_{a,b}}{\left(\sum _{i,j}{M}_{i,j}\right)\cdot \frac{P_a{P}_b}{{\left(\sum _i{P}_i\right)}^2}}}$ where M_i,j is the absolute number of movers from county i to county j, and P_i is the population of county i. This migration connectivity matrix D is therefore a nonnegative, asymmetric matrix in which values less than 1 indicate below-chance migration, and values greater than 1 indicate above-chance migration. Setting ${\sigma }_0^2=10$, we ran Gibbs Sampling for 50 passes (approximately 130,000 steps) to find the MAP solution.

Comparison on synthetic data

In order to understand the properties of our model and quantitatively compare it to alternatives on a dataset with a known ground truth, we performed several experiments with synthetic datasets. We compared against random parcellations (in which elements were randomly merged together) as well as four existing methods: local similarity, which simply thresholds the similarities between pairwise elements (similar to (Thirion et al., 2006; Cohen et al., 2008; Wig, Laumann & Petersen, 2014; Gordon et al., 2014)); normalized cut (Craddock et al., 2012) which finds parcels maximizing within-cluster similarity and between-cluster difference; region growing (Blumensath et al., 2013), an agglomerative clustering method which selects stable points and iteratively merges similar elements; and Ward clustering (Thirion et al., 2014), an agglomerative clustering method which iteratively merges elements to minimize the total variance. Since these methods cannot automatically discover the number of clusters, they (and the random clustering) are set to use the same number of clusters as inferred by our method. We varied the noise level of the synthetic connectivity matrix from low to high, and evaluated the learned clusters using the normalized mutual information with the ground truth, which ranges from 0 to 1 (with 1 indicating perfect recovery).

As shown in Fig. 2, our method identifies parcels that best match the ground truth, across all three datasets and all noise levels. The naive local similarity approach performs very poorly under even mild noise conditions, and becomes worse than chance for high noise levels (for which most parcellations consist of single noisy voxels). Normalized cut is competitive only when the ground-truth parcels are equally sized (matching results from (Blumensath et al., 2013)), and is near-chance in the other cases. Region growing is more consistent across datasets, but does not reach the performance of Ward clustering, which outperforms all methods other than ours. Our model correctly infers the number of clusters with moderate amounts of noise (using the same hyperparameters in all experiments), and finds near-perfect parcellations even at very high noise levels (see Fig. 2C).

We also evaluated our model on a three-spiral dataset previously used in clustering work (Chang & Yeung, 2008), showing that we outperform other methods regardless of initialization scheme (see Figure S1).

Functional connectivity in the human brain

To investigate the spatial structure of functional connectivity in the human brain, we applied our model to data from the Human Connectome Project (Van Essen et al., 2013). Combining data from 468 subjects, this symmetric 59,412 by 59,412 matrix gives the correlation between fMRI timecourses of every pair of vertices on the surface of the brain (at 2 mm resolution) during a resting-state scan (in which subjects fixated on a blank screen). Using the anatomical surface models provided with the data, we defined vertices to be spatially adjacent if they were neighbors along the cortical surface.

Evaluating cortical parcellations is challenging since there is no clear ground truth for comparison, and different applications could require parcellations with different types of properties (e.g., optimizing for fitting individual subjects or for stability across subjects (Thirion et al., 2014)). One simple measure of an effective clustering is the fraction of variance in the full 3.5 billion element matrix which is captured by the connectivity between parcels (consisting of only tens of thousands of connections). As shown in Fig. 3A, our parcellation explains more variance for a given number of clusters than greedy Ward clustering; in order to achieve the same level of performance as our model, the simpler approach would require approximately 30 additional clusters. We can also measure how well this group-level parcellation (using data averaged from hundreds of subjects) fits the data from 20 individual subjects. Although the variance explained is substantially smaller for individual subjects, due both to higher noise levels and inter-subject connectivity differences, our model explains significantly more variance than Ward clustering with 140 clusters (t₁₉ = 2.97, p < 0.01 one-tailed t-test), 155 clusters (t₁₉ = 3.67, p < 0.01), or 172 clusters (t₁₉ = 1.77, p < 0.05). The 220-cluster solutions from our model and Ward clustering generalize equally well, suggesting that our method’s largest gains over greedy approximation occur in the more challenging regime of small numbers of clusters.

One part of the brain in which we do have prior knowledge about cortical organization is in visual cortex, which is segmented into well-known retinotopic field maps (Wang et al., 2014). We can qualitatively examine the match between our 172-cluster parcellation (Fig. 3C) and these retinotopic maps on an inflated cortical surface, shown in Fig. 3D. First, we observe a wide variety in the size and shape of the learned parcels, since the model places no explicit constraints on the clusters except that they must be spatially connected. We also see that we correctly infer very similar parcellations between hemispheres, despite the fact that bilateral symmetry is not enforced by the model. The earliest visual field maps (V1, V2, V3, hV4, LO1, LO2) all radiate out from a common representation of the fovea (Brewer & Barton, 2012), and in this region, our model generates ring parcellations which divide the visual field based on distance from the fovea. The parcellation also draws a sharp border between peripheral V1 and V2. In the dorsal V3A/V3B cluster, V3A and V3B are divided into separate parcels. In medial temporal regions, parcel borders show an approximate correspondence with known VO and PHC borders, with an especially close match along the PHC1-PHC2 border. Overall, we therefore see a transition from an eccentricity-based parcellation in the early visual cluster to a parcellation corresponding to known field maps in the later dorsal and ventral visual areas.

Structural connectivity in the human brain

Based on diffusion MRI data from the Human Connectome Project (Van Essen et al., 2013), we used probabilistic tractography (Behrens et al., 2007) to generate estimates of the strength of the structural fiber connections between each pair of 2 mm gray-matter voxels. Approximately 34 billion tracts were sampled across 10 subjects, yielding a symmetric 59,412 by 59,412 matrix in which about two-thirds of the elements are non-zero. Applying our method to this matrix parcellates the brain into groups of voxels that all had the same distribution of incident fibers. This problem is even more challenging than in the functional case, since this matrix is much less spatially smooth.

Figure 4A shows a 190-region parcellation. Our clustering outperforms greedy clustering by an even larger margin than with the functional data, explaining as much variance as a greedy parcellation with 55 additional clusters. Figure 4B also shows how the model fit evolves over many rounds of Gibbs sampling, when initialized with the greedy solution. Since our method can flexibly explore different numbers of clusters, it is able (unlike a greedy method) to perform complex splitting and merging operations on the parcels. Qualitatively evaluating our parcellation is even more challenging than in the previous functional experiment, but we find that our parcels match the endpoints of major known tracts. For example, Fig. 4C shows 35,000 probabilistically-sampled tracts intersecting with a parcel in the left lateral occipital sulcus, which (in addition to many short-range fibers) connects to the temporal lobe through the inferior longitudinal fasciculus, to the frontal lobe through the inferior fronto-occipital fasciculus, and to homologous regions in the right hemisphere through the corpus callosum (Wakana et al., 2004). Note that the full connectivity matrix was constructed from a million times as many tracks as shown in this figure, in order to estimate the pairwise connectivity between every pair of gray-matter voxels.

Human migration in the United States

Given our successful results on neuroimaging data, we then applied our method to an entirely distinct dataset: internal migration within the United States. Using our probabilistic model, we sought to summarize the (asymmetric) matrix of migration between US counties as flows between a smaller number of contiguous regions. The model is essentially searching for a parcellation such that all counties within a parcel have similar (in- and out-) migration patterns. Note that this is a challenging dataset for clustering analyses since the county-level migration matrix is extremely noisy and sparse, with only 3.8% of flows having a nonzero value.

As shown in Fig. 5A, we identify 83 regions defined by their migration properties. There are a number of interesting properties of this parcellation of the United States. Many clusters share borders with state borders, even though no information about the state membership of different counties was used during the parcellation. This alignment was substantially more prominent than when generating random 83-cluster parcellations, as shown in Fig. 5B. As described in the Discussion, this is consistent with previous work showing behavioral differences caused by state borders, providing the first evidence that state membership also has an impact on intranational migration patterns. Greedy clustering performs very poorly on this sparse, noisy matrix, producing many clusters containing only one or a small number of counties, and has a lower NMI with state borders than even the random parcellations.

The 10 most populous clusters (Fig. 5C) cover 18 of the 20 largest cities in the US, with the two largest parcels covering the Northeast and the west coast. Some clusters roughly align with states or groups of states, while other divide states (e.g., the urban centers of east Texas) or cut across multiple states (e.g., the “urban midwest” cluster consisting of Columbus, Detroit, and Chicago). As shown in Fig. 5D, our method succeeds in reordering the migration matrix to be composed of approximately piecewise constant blocks. In this case (and in many applications) the blocks along the main diagonal are most prominent, but this assortative structure is not enforced by the model. Though largely symmetric, some flows do show large asymmetries. For example, the two most asymmetrical flows by absolute difference are between the urban midwest and Illinois (out of Illinois = 1.3, into Illinois = 2.0), and Florida and Georgia (out of Georgia = 1.3, into Georgia = 2.0).