Direct segmentation of cortical cytoarchitectonic domains using ultra-high-resolution whole-brain diffusion MRI

We acquired whole-brain MAP-MRI data with 200μm resolution in a fixed macaque monkey brain and processed the cortical voxels with a local anisotropic Gaussian denoising filter determined voxel-wise by the plane tangent to the cortical surface. We computed DTI-derived parameters, (i.e., fractional anisotropy – FA; mean, axial, and radial diffusivities - MD, AD, and RD, respectively) and MAP-MRI-derived microstructural parameters (i.e., propagator anisotropy – PA, non-Gaussianity – NG, return-to-origin probability – RTOP, return-to-axis probability – RTAP, and return-to-plane probability – RTPP). We directly clustered all cortical voxels using only the MAP-derived microstructural imaging biomarkers using a Gaussian mixture model (GMM) analysis. MAP-based 3D cytoarchitectonic segmentation revealed laminar patterns similar to those observed in the corresponding histological images. Moreover, transition regions between these laminar patterns agreed more accurately with histology than the borders between cortical areas estimated using conventional atlas/template-warping cortical parcellation. By cross-tabulating all cortical labels in the atlas- and MAP-based segmentations, we automatically matched the corresponding MAP-derived clusters (i.e., cytoarchitectonic domains) across the left and right hemispheres. High-resolution MAP-MRI biomarkers can effectively delineate three-dimensional cortical cytoarchitectonic domains in single individuals.

Methods
We acquired diffusion MRI data brain with 200μm isotropic resolution in a fixed macaque monkey brain using a 3D diffusion-weighted spin-echo (SE) echo-planar imaging (EPI) sequence with 50 ms echo time (TE) and 650 ms repetition time (TR), imaging matrix 375x320x230 on a 7.5x6.4x4.6 cm field-of-view (FOV) with 17 segments per k-z plane and 1.5 partial Fourier acceleration. The diffusion encoding table contained multiple b-value shells: 0.1, 1.0, 2.5, 4.5, 7.0, and 10.0 ms/μm2 with diffusion-encoding gradient orientations (3, 9, 15, 21, 28, and 36, respectively) uniformly sampling the unit sphere on each shell and across shells. The diffusion gradient pulse duration and separation were 6 ms and 28 ms, respectively.

We processed the diffusion MRI data using the TORTOISE software package (https://tortoise.nibib.nih.gov/). From a structural scan, we derived the brain tissue segmentation using FSL-FAST (https://fsl.fmrib.ox.ac.uk) and computed the directions parallel and perpendicular with respect to the cortical surface at each voxel (i.e., the cortical reference frame) using LAYNII (https://layerfmri.com/category/laynii/). Next, we denoised the signals in gray matter voxels using a spatially varying anisotropic Gaussian filter whose orientation was determined by the local cortical reference frame. In each voxel we estimated the diffusion propagator using a MAP-MRI series expansion truncated at order 4 and computed DTI (i.e., fractional anisotropy – FA; mean, axial, and radial diffusivities - MD, AD, and RD, respectively) and MAP-MRI microstructural parameters (i.e., propagator anisotropy – PA, non-Gaussianity – NG, return-to-origin probability – RTOP, return-to-axis probability – RTAP, and return-to-plane probability – RTPP).

We directly clustered individual cortical voxels in each hemisphere based on their values of scalar-valued MAP/DTI parameters (FA, RD, AD, PA, NG, RTOP, RTAP, and RTPP) using a Gaussian mixture model (GMM) clustering algorithm. The best segmentation consistency across the left and right hemispheres was obtained when the MAP parameters PA, NG, RTAP, and RTPP were used as features and the number of clusters was set to 14.

We used 3D morphological processing to separate disjoint (i.e., disconnected) components larger than 100 voxels initially assigned to the same GMM cluster and relabeled them as new clusters. Isolated cluster components with volumes smaller than 100 voxels were likely due to noise and were merged with the neighboring cluster with the largest shared boundary. This step reduced the dependence of the final segmentation on the number of clusters chosen in the initial GMM clustering analysis. 

We automatically matched and labeled the cytoarchitectonic domains segmented using MAP-based voxelwise GMM clustering in the left and right hemispheres. First, we derived a canonical layer segmentation of each label in the symmetric D99 atlas using six layers with equivolumetric spacing obtained with LayNii and computed the cross-tabulation matrix between the MAP- and D99 layer segmentations and applied the Kuhn-Munkres assignment algorithm. Because the D99 layer labels are symmetric across the left and right hemispheres, we directly match the MAP labels by multiplying the MAP-D99 cross-tabulation matrices of the left and right hemispheres. 

Finally, we compared the segmented MAP cytoarchitectonic regions with cortical area estimates from the warped D99 cortical atlas using matched histological tissue sections obtained from the same brain.