Cortical surface parcellation via dMRI using mutual nearest neighbor condition

In this paper, we present a method that aims at parcellating the cortical surface from individual anatomy. The parcellation is obtained using the Mutual Nearest Neighbor (MNN) criterion to obtain regions with similar structural connectivity. The structural connectivity is obtained by applying a probabilistic tractography on the diffusion MRI (dMRI), a non-invasive modality allowing access to the structural information of the white matter. The results are compared to some of the atlases that can be found in the literature. We show that these atlases have lower similarity of structural connectivity than the proposed algorithm implying that the regions of the atlases may have lower functional homogeneity.


INTRODUCTION
Understanding the human cortex organization is a very challenging research field.One of the challenges is to find a relation between anatomical and functional regions.Physiologically, the brain is organized in functional areas, and structural connectivity of the cortical region is the primary indicator of its function [1].In [2], the author described cortical maps suggesting a specific organization of the brain functions in relation to the cortical anatomy.Diffusion Magnetic Resonance Imaging (dMRI) is a non-invasive imaging modality that allows the understanding and exploration of the underlying tissue structure of the human brain.
There have been several methods to label the cortical surface [3,4,5,6,7] and usually researchers register their subjects to some atlases that can be found in the literature; Destrieux (DX) [3], Desikin-Killiany (DK) [4], and Mindboggle (ME) [5], see Fig. 1.The registration can be misleading due to variations in the cortical surface from one subject to another.Previous methods for dMRI based parcellation use atlases as a pre-parcellation step to reduce the computational cost and to obtain sub-regions that have similar structural connectivity [8,9].The method presented in this paper uses the underlying white matter pathways to separate regions.It is similar to the one presented in [10], but diffusion MRI is used instead of functional MRI information.This allows us to parcellate the whole cortical surface and not only an active region of the brain.We also compare connectivity profiles of the resulting regions to both atlases and subdivided atlases using k-means.
To the best of our knowledge, no one has yet investigated the connectivity profiles of the atlases' regions.

Image acquisition
Structural and diffusion MRI data were taken from 6 healthy subjects [11], 4 males and 2 females with an average age of 26.83 years.The T1 weighted images of size 256×256×192 were acquired by Siemens 3T Trio with GRAPPA 3D MPRA-GE sequence (TR = 2250 ms; TE = 2.99 ms; flip-angle = 9 • ; acceleration factor = 2) at 1 mm isotropic resolution.The diffusion weighted images of size 96×96×68 were collected by the same scanner at 2 mm isotropic resolution (64 gradient directions and b-value =1000 s/mm 2 ), with one b0 image.

Image processing
The white/gray matter interfaces (cortical surface) were extracted using Freesurfer [12] from T1 images and remeshed to 10 4 vertices, then projected from the anatomical to the diffusion space.The projection was obtained by registering linearly the brain between the two spaces using FSL [13].

Cortex parcellation
Connectivity profile of a seed: Each node in the cortical mesh is considered as a seed.For each seed, we compute the connectivity profile, i.e. its connectivity with every voxel in the diffusion image space.Connectivity profiles are obtained by running probabilistic tractography for all seeds.We assume a model of multiple fiber orientations for each voxel developed in FSL [14] with 10 4 samples.Clustering the cortical surface: The goal of parcellation algorithm is to divide the cortical surface into a set of non overlapping and connected regions that have similar connectivity profiles.In the first step, regions are singletons of seeds and are possible candidates for merging.The algorithm is explained in Algorithm 1. CP represents the matrix of the seeds' connectivity profiles and M C is the mesh connectivity Initialize: for Each region do 5: end for 8: for Each region do end for return Label 20: end procedure For every iteration, we obtain the set of neighboring regions N r to the region r (line 5), and compute the similarity measure, SM , between region r and N r (line 6).The vector B of merging regions candidates contains the labels of the possible candidates: B(r) is the neighbor of region r that has the highest SM value.We define the SM between two regions R and S as; where corr is the P earson s correlation and |R| is the number of seeds in R. Region r and B(r) are merged if they satisfy the Mutual Nearest Neighbor (MNN) condition:

RESULTS AND DISCUSSION
The mean of the similarity values (MSV) between all the connectivity profile pairs inside the regions is used to measure the homogeneity of the parcellation's regions.High MSV means that the region has close connectivity profiles.Fig. 2 (a) shows the MSV of the different atlases and subjects.We applied the cortical surface parcellation algorithm explained in Section 2.3 to the dMRI data [11].Fig. 2 (b) shows the MSV of the resulting parcellation for all subjects and T ∈ {1000, 800, 600, 400, 200}.To compare the MSV of the atlases, Fig. 2 (a), and the ones of the proposed parcellation, we chose T ∈ {148, 62}, Fig. 2 (c).These values give parcellations with close number of regions to the ones of the atlases (see Fig. 3).Fig. 2 (d) shows the MSV and variance of the similarity measure values of the regions obtained by subdividing the atlas' regions using k-means.
Several papers used the atlases as a pre-parcellation step to reduce the computation time [8,9].Each atlas' region is parcellated using k-means.This approach is assumed to group seeds that have similar CP inside the atlas' region.We use the eigenvalues of the CP s' cross-correlation to identify the number of clusters in each region.Based on the eigengap (difference between consecutive eigenvalues), we sort the eigenvalues, α i for i ∈ {1, 2, ..., |r|}, of the cross-correlation matrix in decreasing order and pick the number of clusters, k, as: k = argmax i (α i − α i+1 ).The resulting number of regions for the six subjects is shown in Table 1.
As shown in Fig. 2 (a), the MSV of the DX atlas is higher than the others and is followed by DK and ME atlases for all subjects.This follows the decreasing order of the number of regions of the atlases, see Fig. 1.
In average and as expected, the MSV is proportional to T due to including more dissimilar connectivity profiles for low T values.We observe also that the order is homogeneous across subjects, see Fig. 2 (b).For all subjects in Fig. 2 S1   (c), the parcellations of the proposed method provided regions with higher MSV when compared to the atlases (Fig. 2 (a)).This is due to merging only regions with the highest similarity measure due to the MNN condition.In Fig. 3, which represents the result of the parcellations for T = 62, we can distinguish several classical regions like Fusiform gyrus, Superior Frontal, Occipital and others.
Even if the MSV is increased after applying the k-means on the atlases, Fig. 2 (d), but it still less than the MSV of the proposed parcellation with bigger regions (T ∈ {148, 62}), Fig. 2 (c).This implies that the atlas' regions have low homogeneity of the structural connectivity because they do not include the structural information, which implies that they have lower functional homogeneity.
Table 2 shows the mean and the standard deviation of the Dice coefficients vector (D).The Dice coefficient measures the extent of spatial overlap between regions.Let A i be a region in the parcellation A (n regions) and B j be a region in the parcellation B (m regions).D is a vector of n elements that compares the regions of A to B: |Ai|+|Bj | for j = 1, .., m.The mean of D between the atlases and the proposed parcellation is low for all subjects, see the first four rows in Table 2. Atlases and parcellations based on MNN condition do not only differ in the MSV but also have low spatial overlapping.The MSV for ME and DK are close, Fig. 2 (a), this can be explained by their high spatial overlapping (90 %), last two rows in Table 2.This is why we do not compare the Dice coefficient between the parcellation with T = 68 and the DK   [8].Other similarity measures like T animoto measure can be used and should be tested in this framework to see the effect of the similarity measure on the resulting parcellation.The implementation of the algorithm will be available for comparison studies.

CONCLUSION
We have presented algorithm to subdivide the human cerebral cortex into neuro-anatomical parcels that have similar connectivity profiles.The results show that the algorithm provides regions with higher similarity between the CP with compared to atlases even after subdividing them.These atlases have low structural homogeneity because they do not include any structural information.The algorithm may provide valuable tool for research studies involving group and functional studies.Even though and as expected the parcellation using MNN gives regions with higher MSV, it needs to be validated using functional data.

Fig. 1 :Algorithm 1
Fig. 1: The four views of the atlases obtained from FreeSurfer for Subject 1. From left to right; DX (148 regions), DK (68 regions) and ME (62 regions).The Thalamus is not included in the labeling.

Fig. 2 :
Fig. 2: The MSV and variance of the cross-correlation of the CP over regions of (a) different atlases and subjects, (b) the resulting parcellation for different subjects and T values, (c) parcellations with T ∈ {148, 62} and (d) the regions obtained by subdividing the atlases using k-means.

Fig. 3 :
Fig.3: The Resulting parcellation of six subjects on the inflated cortical surfaces with T = 62.Each panel correspond to four views of one subject. .The Thalamus is not included in the parcellation.We obtain number of regions for the different subjects between 60 and 66.

Table 1 :
Number of parcels for different subjects and atlases obtained from k-means parcellation on the atlases.

Table 2 :
The Dice coefficient (mean/standard deviation in %) between the different parcellations (T ∈ {148, 62}) and atlases.atlas.DX has low spatial overlapping with DK and ME, see row five to eight in Table2, this is because DX has more regions than DK and ME.We used P earson's correlation in Equation.1 because it was used successfully in different papers like in Anwander et al.