Cerebral cortical areas are fundamental units of brain function and historically defined by local transitions in microarchitecture^1,2. The idea of local transitions was adopted by resting-state fMRI (rs-fMRI) parcellations using functional connectivity gradients³. Consequently, local approaches appeared to capture cytoarchitectonic boundaries⁴. These local approaches contrasted with global approaches^5-8 that directly or indirectly maximized the functional homogeneity within parcels. Here we combined both local and global approaches in a single fusion model yielding a cortical parcellation with excellent parcel homogeneity and alignment with cytoarchitectonic and visuotopic areal boundaries.

We utilized rs-fMRI data from the Genomics Superstruct Project (GSP⁹). The data underwent rs-fMRI preprocessing previously reported in Holmes et al.⁹ and was projected onto the fsaverage6 surface space with roughly 82000 vertices. The data were divided into training (N = 744) and test (N = 744) sets. To assess whether our parcellation would generalize well to data from a different scanner and preprocessing pipeline, we considered a group-averaged connectivity matrix (n=468) from the Human Connectome Project (HCP¹⁰).

For each subject, the fMRI timecourses at each vertex was normalized to have zero mean and 1 standard deviation. At each vertex $$$n$$$, the normalized timecourses of all subjects were concatenated into a long vector $$$y_n$$$ and normalized to be unit length. The resulting data were modeled as a mixture model with a Markov random field (MRF) prior. The mixture model assumed the data $$$y_n$$$ followed a von Mises-Fisher distribution conditioned on the parcellation label $$$l_n$$$:

$$p(y_n| l_n, \mu_{1:L},\kappa_{1:L})=p(y_n| l_n, \mu_{l_n},\kappa_{l_n})=z(\kappa_{l_n}) e^{\kappa_{l_n} y_{n}^T \mu_{l_n}},$$

where $$$z(\kappa_{l_n})$$$ is a normalization constant, $$$L$$$ is the number of parcels, $$$\mu_{l_n}$$$ is the mean direction and $$$\kappa_{l_n}$$$ is the concentration parameter of the $$$l_n$$$-th von Mises-Fisher distribution.

The MRF prior encouraged neighboring vertices to have the same parcellation labels. The prior was weighted so that there was lower (higher) penalty in the presence (absence) of local functional connectivity gradients:

$$p(l_1, \dots, l_n) = \sum_{u \in V} \sum_{v \in N_G(u)} \delta(l_u,l_v) (e^{-kGrad(u,v)}-e^{-k}),$$

where V are all brain locations and $$$N_G(u)$$$ are the neighboring vertices of vertex $$$u$$$; $$$\delta(l_u,l_v)$$$ is equal to zero if $$$l_u=l_v$$$ and one otherwise. $$$Grad(u,v)$$$ is the local connectivity gradient between brain locations $$$u$$$ and $$$v$$$, and $$$k$$$ is a tunable parameter set to 15. The local connectivity gradients were obtained from Gordon et al.¹¹.

An additional simple geometric prior ensured each parcel $$$l$$$ remained spatially connected. Let $$$s_n$$$ be the $$$x,y,z$$$ coordinates of vertex $$$n$$$ on the fsaverage6 spherical meshes. The geometric prior assumed that $$$s_n$$$ conditioned on the label $$$l_n$$$ follows a von Mises-Fisher distribution

$$p(s_n| l_n, \nu_{1:L},\gamma_{1:L})=p(s_n| l_n, \nu_{l_n},\gamma_{l_n})=z(\gamma_{l_n}) e^{\gamma_{l_n} s_{n}^T \nu_{l_n}},$$

where $$$z(\gamma_{l_n})$$$ is a normalization constant, $$$\nu_{l_n}$$$ is the mean direction and $$$\gamma_{l_n}$$$ is the concentration parameter of the $$$l_n$$$-th von Mises-Fisher distribution.

Inference proceeded by setting $$$\gamma_1, \dots ,\gamma_L = 5 \times 10^6$$$ and iterating between estimating the labels $$$l_{1:N}=\{l_1, \dots, l_N\}$$$ and the von Mises-Fisher parameters $$$\mu_{1:L}$$$ and $$$\kappa_{1:L}$$$. Once the algorithm has converged, the geometric prior was gradually weakened by iteratively reducing $$$\gamma_{1:L}$$$ (by a factor of ten) and re-estimating the parcellation while ensuring each parcel remained spatially connected. This procedure was initialized with 100 random initializations, and the solution with the highest log likelihood (excluding the geometric prior) was selected.

The model was applied to the GSP training set to create a cerebral cortex parcellation. The number of parcels ($$$L = 333$$$) was selected to match a recent rs-fMRI gradient-based cortical parcellation by Gordon et al.¹¹. The two parcellations were compared in terms of alignment with cytoarchitectonic/visuotopic borders^12,13 and functional homogeneity on the GSP test set and external HCP dataset. Functional homogeneity was computed by averaging pairwise correlations within each parcel.

Fig. 1 shows our cortical parcellation estimated from the GSP training set. Compared with Gordon et al.¹¹ the average homogeneity of our parcellation was 9% higher on the GSP test set ($$$p \approx 0$$$) and 6% higher on the HCP data. The two parcellations have comparable alignment with cyto-architectonic and visuotopic areal boundaries (Table 1). It is worth noting that our parcellation appeared to further fractionate BA3 and BA4 along the homunculus (consistent with Gordon et al.¹¹) and fractionate V1 along visual eccentricity (consistent with Yeo et al.⁷).We have created a cerebral cortex parcellation by fusing local and global connectivity features. Previous results suggest that parcellations using local features give better cytoarchitectonic alignment, while parcellations using global features achieve better parcel homogeneity^4,14. By fusing local and global features, our parcellation enjoys both excellent areal alignment and parcel homogeneity.