Mapping the individual human cortex using multidimensional MRI and unsupervised learning

Abstract

Human evolution has seen the development of higher-order cognitive and social capabilities in conjunction with the unique laminar cytoarchitecture of the human cortex. Moreover, early-life cortical maldevelopment has been associated with various neurodevelopmental diseases. Despite these connections, there is currently no noninvasive technique available for imaging the detailed cortical laminar structure. This study aims to address this scientific and clinical gap by introducing an approach for imaging human cortical lamina. This method combines diffusion–relaxation multidimensional MRI with a tailored unsupervised machine learning approach that introduces enhanced microstructural sensitivity. This new imaging method simultaneously encodes the microstructure, the local chemical composition and importantly their correlation within complex and heterogenous tissue. To validate our approach, we compared the intra-cortical layers obtained using our ex vivo MRI-based method with those derived from Nissl staining of postmortem human brain specimens. The integration of unsupervised learning with diffusion–relaxation correlation MRI generated maps that demonstrate sensitivity to areal differences in cytoarchitectonic features observed in histology. Significantly, our observations revealed layer-specific diffusion–relaxation signatures, showing reductions in both relaxation times and diffusivities at the deeper cortical levels. These findings suggest a radial decrease in myelin content and changes in cell size and anisotropy, reflecting variations in both cytoarchitecture and myeloarchitecture. Additionally, we demonstrated that 1D relaxation and high-order diffusion MRI scalar indices, even when aggregated and used jointly in a multimodal fashion, cannot disentangle the cortical layers. Looking ahead, our technique holds the potential to open new avenues of research in human neurodevelopment and the vast array of disorders caused by disruptions in neurodevelopment.

Graphical Abstract

multidimensional MRI, microstructure imaging, machine learning, cortical parcellation, diffusion–relaxation

Introduction

The development of higher-order cognitive and social capacities that are distinctively human emerged in concert with changes in the layered cytoarchitecture of the human brain cortex.^1-3 These cortical layers define grey matter (GM) architecture and neuroanatomical regions, which play a critical role in development and connectivity, as well as function and pathology.^4,5 Furthermore, region-specific cortical changes are observed in a range of neurodevelopmental and neuropsychiatric diseases,⁶ which are modified by a combination of genetic and molecular factors, as well as aberrant early cell migration.⁷ Although these present exciting new avenues for targeted therapeutics in the treatment of neurodevelopmental diseases, a significant scientific and clinical gap exists owing to the lack of in vivo noninvasive imaging techniques to characterize the structure of the human cortex. To date, defining cortical lamina is only possible using laborious cytological methods ex vivo.

Cortical imaging using MRI has been the subject of several studies, generally using two different approaches: atlas and contrast based. Atlas-based cortical parcellation relies on the registration of structural MRI data to a ‘one-size-fits-all’ standard template or atlas, in which the cortical labels have been determined in advance.⁸ While robust, convenient and requires very little data, this approach bases the cortical segmentation only on the sulcal and gyral geometry of the subject, without considering other, potentially important drivers of MR contrast, such as diffusion and relaxation. MR contrast-based cortical segmentation methods have also shown sensitivity towards laminar architectonic features. These can be roughly divided into relaxation- and diffusion-based approaches. The former includes T₂-weighted,⁹ T₂*-weighted,¹⁰ T₁-weighted,¹¹ MR phase contrast¹² and combinations of T₁, T₂ and T₂* contrasts.¹³ While using these multiple relaxation contrasts provided segmentation of specific regions such as the primary somatosensory and motor areas, other cortical areas could not be detected and elucidated.^9,14,15 Diffusion MRI (dMRI) provides complementary information and is particularly sensitive to micro- and meso-scale architectural features in biological tissue.^16-20 Although mostly used to visualize white matter (WM) orientational and microstructural features,^21-30 dMRI has been repeatedly shown to provide meaningful information about the microstructural properties of GM, capturing cortical-depth variations in diffusion characteristics, especially when examining diffusion anisotropy.^14,23,31-41

Despite these significant advances, current noninvasive techniques for characterizing GM laminar composition remain inadequate.¹⁴ So far, GM characterization has been done using either relaxation or dMRI-based methods, separately. While MRI may not achieve the same level of microscopic resolution as optogenetic methods like Brainbow,⁴² recent progress in acquisition, pulse sequence design^43-46 and data processing strategies^47-52 has advanced significantly. These developments have now made it possible to simultaneously encode both diffusion and relaxation imaging data and analyse them jointly. This innovative approach, known as diffusion–relaxation multidimensional MRI, enables the concurrent extraction of chemical and microstructural information through the use of relaxation (T₁ and T₂) and diffusion mechanisms, respectively. As a result, it offers a plethora of sub-voxel morphological and compositional features, making it applicable in various biological contexts.^52-61 Recently, this technique has been combined with unsupervised machine learning, resulting in discrete clustering of data, further enhancing its potential and usefulness.⁶² Multidimensional MRI yields high-dimensional images, in which each voxel may contain a 2D or higher dimensional probability distribution encoding the correlations among the investigated MR ‘dimensions’, i.e. T₁, T₂ and diffusion.

We hypothesize that multidimensional T₁–T₂–diffusion distributions vary across the human cortex in a voxelwise manner, and that these changes correspond to the known laminar structure. To investigate this hypothesis, we used a unique distance measure between multidimensional spectra based on the physics of multidimensional MRI employing a modified Wasserstein distance derived from a linearized version of the optimal transport distance (LOT).⁶³ This study makes a unique contribution by combining multidimensional MRI and physics-based unsupervised machine learning to detect cortical lamina without human input. We include validation studies on postmortem human brain specimens with comparison to histologic segmentation, as well as computational validation experiments of the machine learning technique. We demonstrate that various imaging techniques such as T₁ and T₂ mapping, diffusion tensor imaging (DTI),⁶⁴ neurite orientation dispersion and density imaging (NODDI)⁶⁵ and the mean apparent propagator (MAP)²⁴ models, while previously shown to be sensitive to microstructural variations in human cerebral cortical GM,^9,15,66-69 are unable to resolve intra-cortical layers effectively. We show that even when used jointly in a multimodal fashion with the proposed unsupervised framework, scalar maps derived from relaxometry, DTI, NODDI and MAP-MRI significantly underperform when compared with diffusion–relaxation correlation MRI in terms of their ability to resolve the layers within the cortex. In the future, this technique might be extended to in vivo studies that reveal the structure of the human cortex, potentially sparking new research discoveries in normal human neurodevelopment and neurodevelopmental diseases, as well as targeted therapeutics.

Materials and methods

Subject cohort

We investigated fixed ex vivo cerebral cortical tissue specimens (see Table 1 for location) from three control brain donors (ages 48,30 and 83 years) obtained from the University of Washington BioRepository and Integrated Neuropathology (BRaIN) laboratory. Neuropathological changes were limited to those associated with age, including low Alzheimer’s disease neuropathologic change in the 83-year-old donor, which by consensus are not considered sufficient pathology for significant neurodegeneration or cognitive impairment.^70,71 Brain donation occurred through both the Pacific Northwest Brain Donor Project Study, approved by the University of Washington School of Medicine, under the auspices of the anatomical gift act, and under the Seattle Longitudinal Study. Participants in the Seattle Longitudinal Study are invited to participate in brain donation upon death. The Seattle Longitudinal Study was approved by the University of Washington Institutional Review Board. Table 1 summarizes demographic details, the manner of death, postmortem interval and anatomical region for each subject. Briefly, every donor brain undergoes fixation in 10% neutral-buffered formalin for at least 2 weeks, followed by agarose embedding, conventional ex vivo MRI, coronal slicing at precisely 4 mm intervals and sampling (see Latimer et al.⁷² for complete details). Fixed tissue sample blocks of ∼20 × 20 × 10 mm³ were obtained from cortical regions. The specimens were sent to the National Institute on Aging specifically for the current study and subsequently scanned using multidimensional MRI acquisitions.

Table 1

Main demographic and histopathological findings in subjects included in this study

Donor	Age	Sex	Postmortem interval (h)	Manner of death	Anatomic region
1	48	M	25.3	Suicide	Left middle frontal cortex
2	30	M	7.5	Suicide	Right superior frontal gyrus
3	83	F	2.7	Natural	Left middle temporal gyrus

MRI acquisition

Upon receipt at the National Institute on Aging and prior to MRI scanning, each formalin-fixed brain specimen was transferred to a phosphate-buffered saline–filled container for 12 days to ensure that any residual fixative was removed from the tissue. Each specimen was then placed in a 20 mm tube and immersed in perfluoropolyether (Fomblin LC/8, Solvay Solexis, Italy), a proton-free fluid void of a proton MRI signal. Specimens were imaged using a 7 T Bruker vertical bore MRI scanner equipped with a microimaging probe and a 20 mm quadrupole radio frequency (RF) coil.

Multidimensional data were acquired using a 3D inversion recovery diffusion-weighted sequence with a repetition time of 1000 ms and isotropic resolution of 200 µm. We used a previously published hierarchical sampling scheme^56,60 to encode the multidimensional MR space spanned by T₁ and T₂ (i.e. T₁–T₂), by T₁ and mean diffusivity (i.e. T₁–MD) and by T₂ and MD (i.e. T₂–MD). Briefly, T₁ was encoded by varying the inversion time between 14.6 and 980 ms; T₂ was encoded by varying the echo time between 12.3 and 125 ms; and diffusion was encoded by varying the diffusion weighting b-value between 2400 and 14 400 s/mm², at 13 independent directions,⁷³ with gradient duration of 4 ms and diffusion time of 15 ms. The total number of imaging volumes to generate voxelwise T₂–MD, T₁–MD and T₁–T₂ distributions was 660. Full details regarding the multidimensional MRI acquisition can be found in the Supplementary material.

Importantly, the MD variable used in this study was derived from an orientationally averaged diffusion-weighted signal, calculated for every inversion time–echo time–b-value combination. As a result, it can be considered a rotationally invariant diffusion scalar parameter that is not anticipated to vary based on cortical orientation or the selected laboratory coordinate system. It is important to highlight that due to this orientational averaging, the MD in this study is influenced by concomitant alterations in cellular dimensions, shapes and orientations.

A high-resolution MRI scan with an isotropic voxel dimensions of 100 µm was acquired using a fast low angle shot sequence⁷⁴ with a flip angle of 49.6° to serve as a high-resolution reference image and facilitate co-registration of histological and MR images.

A standard DTI protocol was applied with the same imaging parameters as the multidimensional data and using 21 diffusion gradient directions and four b-values ranging from 0 to 1400 s/mm².

T₁ and T₂ map processing

Conventional quantitative relaxation maps were computed by fitting the signal decay to monoexponential functions. The T₁ value was computed by fitting a subset of the multidimensional data that included 20 images with inversion times in the range of 14.6–980 ms and minimal echo time (12.3 ms). The T₂ value was computed by fitting a subset of the multidimensional data that included 20 images with echo times in the range of 12.3–125 ms and no inversion recovery pulse. In both cases, no diffusion weighting was applied.

DTI processing

DTI parameters,⁶⁴ axial diffusivity, radial diffusivity, MD and fractional anisotropy, were calculated using in-house MATLAB (The Mathworks, Natick, MA) code based on previous work.⁷⁵

Higher-order dMRI model processing

Two higher-order diffusion models, NODDI and the MAP-MRI, were used to allow computation of scalar maps using the full range of b-values in the multidimensional data (i.e. ∼14 400 s/mm²). These models were evaluated by fitting to a subset of the multidimensional data that included 16 b-values between 2400 and 14 400 s/mm², at 13 independent directions,⁷³ with gradient duration of 4 ms, diffusion time of 15 ms and with an echo time of 26.5 ms.

For the NODDI model, the NODDI toolbox for MATLAB was used for fitting the data.⁶⁵ The software default value of intrinsic diffusivity for ex vivo fixed tissue was 0.6 µm²/ms; however, similarly to previous ex vivo studies,⁷⁶ we established this parameter at 1.3 µm²/ms to allow convergence of the model. The isotropic diffusivity was set to 1.9 µm²/ms.⁷⁷ For ex vivo NODDI model, a dot-compartment (the diffusion is restricted in all directions) was included for the fitting.^18,78 The parameters derived from the NODDI model were orientation dispersion index, neurite density index and isotropically restricted volume fraction.

For the MAP model, we estimated the voxelwise diffusion propagator (i.e. the probability density function of 3D net displacements of diffusing water molecules) using a MAP-MRI series expansion truncated at Order 6.²⁴ Using the MAP-MRI series expansion, we computed images of the MAP parameters: the zero-displacement parameter maps, return-to-origin probability, return-to-axis probability and the return-to-plane probability, which are all inversely proportional to the spatial dimensions within the microstructure were obtained. We also mapped the non-Gaussianity, which quantifies the dissimilarity between the full propagator and its Gaussian component and reflects the deviation from a simple tensor model, and the propagator anisotropy, which quantifies the directional dependence of the diffusion process.

Multidimensional MRI processing

Prior to processing, multidimensional MRI data were denoised using an adaptive nonlocal multispectral filter (i.e. NESMA⁷⁹), which was shown to reduce noise and improve the accuracy of the resulting injury MRI biomarker maps.⁸⁰ The filtered data were then processed using marginally constrained, ${\ell }_2$-regularized, nonnegative least square optimization to compute the multidimensional distribution in each voxel, as previously described.^47,51 It is a well-tested approach that had been proven robust and reliable,^81-84 which in this study had resulted in three types of distributions in each voxel: T₁–T₂, T₂–MD and T₁–MD.

Histology and MRI co-registration

After MRI scanning, each tissue block was transferred for histological processing. Tissue blocks from each brain specimen were processed using an automated tissue processor (ASP 6025, Leica Biosystems, Nussloch, Germany). After tissue processing, each tissue block was embedded in paraffin, and the entire block was cut into 10 µm thick sections. A Nissl stain was performed at every 100 µm to visualize cellular morphology. All stained sections were digitally scanned using an Aperio scanner system (Aperio AT2—High Volume, Digital whole slide scanning scanner, Leica Biosystems, Inc., Richmond, IL) for further assessment and analyses. A Zeiss Imager A2 (ImagerA2 microscope, Zeiss, Munich, Germany) bright-field microscope with ×40 and ×63 magnification lenses was used to identify histologic and pathologic details, as needed.

To allow 3D volume co-registration of the MRI and histological data sets, a reconstruction pipeline was assembled to render 3D Nissl volumes. For each case, the series of 2D Nissl images were first zero padded to ensure uniformity in in-plane dimensions and stored as a 3D matrix. Individual slices were aligned consecutively with respect to the centre slice using an affine transformation. Once completed, the matrix dimension and resolution of the 3D histological volume were matched to the high-resolution MR images, which were used as anatomical reference for co-registration using rigid geometric transformation (Image Processing Toolbox, MATLAB, The Mathworks, Natick, MA). The transformed histology images were overlaid on the MR images to assess the quality of the co-registration.

A representative slice from the MR-Nissl co-registered data sets was selected from each case. Histological tissue classification into cortical regions based on the cellular patterns of the different layers: Layer I (molecular layer), Layer II (external granular layer), Layer III (external pyramidal layer), Layer IV (internal granular layer), Layer V (internal pyramidal) and Layer VI (multiform layer) were performed using a deep learning convolutional neural network (DenseNet, HALO AI, New Mexico, USA), which uses DenseNet-121 architecture.⁸⁵ To train the network, Nissl sections from four additional subjects, each included 10 large representative areas from each layer and five transitional regions between layers, were used. When the tissue segmentation was not possible due to an arrangement of the cells different to a coronal orientation, a manual correction of the annotation was necessary. As such, analyses based on the cortical laminations are based on deep learning using DenseNet with oversight and adjustments as needed by an expert neuroanatomist.

Linear optimal transport

Optimal mass transport can provide a distance measure between two probability distributions by computing the cost of transforming one distribution into another while minimizing the total transportation cost. Given the multidimensional distributions T₁–T₂, T₂–MD and T₁–MD in each voxel, we measured the pairwise distance, or similarity, between the distributions of two voxels. In contrast to information-theoretic dissimilarity measures based on the concept of entropy (i.e. f-divergences such as Jensen–Shannon), the Wasserstein distance enables probability densities to be jointly modelled within their diffusion–relaxation domains. Furthermore, even for nonoverlapping domains, the Wasserstein distance can provide a useful gradient and distance measure with respect to the underlying geometric space.⁸⁶ Specifically, we utilize a LOT, which provides a computationally efficient way to compare sets of probability distributions using a modified two-Wasserstein distance.⁶³ The system diagram in Fig. 1 summarizes the approach used in this paper.

Figure 1

System diagram. Voxelwise diffusion–relaxation probability distributions in the image domain were transformed to the transport domain, where each spectrum is represented as a point on a high-dimensional Riemannian manifold. The shortest distance between a point on the manifold and the mean distribution I₀ can be computed using the modified two-Wasserstein distance. The linear optimal transport (LOT) embedding provides a linearized version of this distance and enables geodesic distances to be approximated as Euclidean distances in the transport domain. Unsupervised learning was performed in the transport domain using the diffusion–relaxation distributions alone, without their spatial coordinates, and the cluster assignments were visualized in the image domain as segmented cortical images. Note that when k = 3 clusters are used, WM and two GM components are distinguished.

First, a small constant was added to all distributions to ensure that they were strictly positive, followed by normalization of the T₁–T₂, T₂–MD and T₁–MD distributions to ensure the same total mass. Next, within-subject reference distributions were calculated for each subject; the reference distribution is the mean distribution obtained by computing the Euclidean mean of all voxelwise spectra belonging to T₁–T₂, T₂–MD and T₁–MD distributions within a subject. Subsequently, all pairwise distances were calculated between this reference distribution and each individual voxel distribution.

Mathematically, for a given mean probability distribution $I_0(\mathit{x})$ as in Fig. 1, the LOT distance⁶³ was computed for each distribution with respect to this common reference, $I_0(\mathit{x})$ for all $n=1,\dots ,N$ voxels across the 3D volume within each individual high-dimensional brain image, such that $I_0(\mathit{x})=(1/N){\sum }_{n=1}^N{I}_n(\mathit{x})$⁠. For each voxelwise distribution $I_n(\mathit{x})$⁠, if $f_n:\mathrm{\Omega }\to \mathrm{\Omega }$ is a mass preserving mapping from $I_0$ to $I_n$⁠, then Eq. 1 summarizes the analysis equation to transform distributions in the image domain to their representation in the transport domain.

$$\begin{array}{rl}& f_n^{*}(\mathit{x})={\displaystyle \underset{f_n\in \mathrm{MP}}{min}}{\displaystyle \underset{\mathrm{\Omega }}{\overset{}{\int }}}(f(\mathit{x})-\mathit{x}{)}^2{I}_0(\mathit{x})d\mathit{x}\\ & \mathrm{s}.\mathrm{t}.det(D{f}_n(\mathit{x}))I_n(f_n(\mathit{x}))=I_0(\mathit{x})\mathrm{\forall }\mathit{x}\in \mathrm{\Omega }\end{array}$$

(1)

Here, $D{f}_n$ is the Jacobian of the mapping $f_n$ and MP is the family of all mass preserving mappings from $I_0$ to $I_n$⁠. Each transport map $f_n^{*}(\mathit{x})$ is a vector field that transforms a sample distribution $I_n(\mathit{x})$ to the reference distribution $I_o(\mathit{x})$⁠, requiring the minimum optimal transport distance. Equation 1 has a unique minimizer, as previously shown.⁸² The isometric linear embedding for the distribution is then provided by⁶³  $\widehat{I_n}=(f_n^{*}(\mathit{x})-\mathit{x})\sqrt{I_0(\mathit{x})}.$

Within each subject, this process was repeated for T₁–T₂, T₂–MD and T₁–MD, and the resultant transport maps were concatenated into a data matrix $X\in {\mathbb{R}}^{n\times D}$⁠, where n is the number of voxels in the image and D is the dimension of the transport maps when concatenated across all distributions in each voxel T₁–T₂, T₂–MD and T₁–MD.

Unsupervised machine learning

Unsupervised learning using k-means clustering was performed on the data matrix in the transport space. The process, as illustrated in Fig. 1, was done in two steps: first, segmenting GM and WM, and background noise using k = 3 clusters. Next, based on the k = 3 step results, the entire GM was further segmented into cortical lamina. The histological segmentation showed that Layer IV in all specimens was thinner than the single MRI voxel size (Supplementary Fig. 1). Based on that, we decided to use k = 5 GM clusters to reflect the number of feasibly detectable cortical layers in the human cortex samples we examined. Each k-means clustering step was iterated 100 times with random initializations according to the k-means++ algorithm.⁸⁷ All computations were performed in MATLAB (Mathworks, USA).

Statistical analysis

To assess the stability of k-means cluster assignments over multiple iterations with random initializations, we computed the adjusted Rand index⁸⁸ for each clustering iteration compared to the optimal clustering. The adjusted Rand index corrects for assignments by random chance.

Results

Histologically derived cortical architectonic features

Cellular morphology is demonstrated in representative Nissl-stained sections from the three cases in Fig. 2A–C, corresponding to Subjects 1–3, respectively. Region-specific laminar variations in the density of neuronal cells and their architectural patterns, which provide the biological basis for subdividing the cortex into structurally distinct areas, are seen in the magnified regions below. The cortical laminae were delineated using the convolutional neural network as detailed above and are shown in Fig. 2 as overlayed lines on the magnified regions from each subject.

Figure 2

Histological findings from the three examined cases. Cellular morphology is demonstrated in representative Nissl-stained sections from Subjects 1–3 in A–C, respectively. Entire cortical cross-sections are shown in the top panel and magnified regions in the bottom panel. Delineated cortical layers were overlayed on the magnified regions from each subject, depicting region-specific laminar variations in the density of neuronal cells and their architectural patterns.

Cortical depth dependency of diffusion–relaxation distributions

Distributions from transcortical voxels are shown in Fig. 3 to illustrate our main hypothesis that microstructural alterations between cortical layers lead to notable changes in the diffusion–relaxation correlation. A transcortical region of interest that included 12 voxels was taken from Subject 1, and voxelwise T₁–T₂, T₁–MD and T₂–MD spectra are shown from the subcortical WM margin out to the cortical surface. We observed that the peak T₂ component shifted from 57.8 ms at the cortical surface to 30.7 ms at the cortical depth. Compared to T₂, the effect on the main T₁ peak was less pronounced, shifting from 339.3 ms at the cortical surface to 281.2 ms at the cortical depth; however, a clear short T₁ component appears at about 1 mm into the cortex and intensifies as a function of depth. The main diffusion peak shifted from 0.85 to 0.54 µm²/ms as a function of depth, while the most significant change was observed with the appearance of a slow diffusion component that intensifies as a function of depth from the cortical surface.

Figure 3

Intra-cortical diffusion–relaxation signature. A transcortical region of interest from Subject 1, and single voxel depth-specific T₁–T₂, T₁–MD and T₂–MD spectra are shown with logarithmic scale from the subcortical WM margin out to the cortical surface. The main diffusion peak shifted from 0.85 to 0.54 µm²/ms as a function of depth, along with the appearance of a slow diffusion component that intensifies as a function of depth from the cortical surface. In addition, the main T₁ and T₂ components shifted from 339.3 and 57.8 ms at the cortical surface to 281.2 and 30.7 ms, respectively, at the cortical depth. Finally, a clear short T₁ component appears at about 1 mm into the cortex and intensifies as a function of depth. To facilitate visualization of the logarithmic scale, an asterisk symbol marks the location of the main signal components at the cortical surface. MD, mean diffusivity.

Segmenting cortical laminae using multidimensional MRI

First, k-means clustering was performed with three clusters to identify gross anatomical features and to separate WM and GM regions. For all subjects, the WM regions were robustly identified without human guidance or knowledge of spatial contiguity.

Our unsupervised segmentation framework was applied separately on multidimensional and scalar MRI data from the three subjects we investigated. Scalar maps of T₁, T₂, fractional anisotropy, axial diffusivity, radial diffusivity, MD, neurite density index, orientation dispersion index, isotropically restricted volume fraction, return-to-axis probability, return-to-origin probability, return-to-plane probability, non-Gaussianity and propagator anisotropy were aggregated together and fed into the same unsupervised segmentation framework instead of the multidimensional T₁–T₂–MD signatures to enable comparisons and to assess the relative advantage of one method over the other. Figure 4 shows the 3D clustered images of the three subjects resulting from the automatic cortical layer parcellation. A representative slice (x–y plane) is shown for each subject, along with four horizontal and four vertical cross-sections in the z direction. The top and bottom panels in Fig. 4 show the multidimensional spectra- and scalar-based clustering results.

Figure 4

Segmentation spatial contiguity and comparison of 1D and multidimensional MRI. Segmented cortical volumes from Subjects 1 to 3 are shown in A–C, respectively. A representative slice (x-y plane) is shown for each subject, along with four horizontal and four vertical cross-sections in the z direction. The top and bottom panels show the multidimensional spectra- and scalar-based clustering results. Diffusion–relaxation MRI data result in five distinct and spatially consistent clusters, labelled GM 1–GM 5 from the outermost to innermost cortical layer, respectively. Layered and contiguous laminar structures were demonstrated across the volume in all samples. Conversely, scalar MRI data result in roughly three consistently distinguishable regions, and the cluster maps appear noisy and spatially inconsistent. GM, grey matter.

When the diffusion–relaxation correlation spectra are used as inputs to the cortical segmentation framework, five distinct and spatially consistent clusters can be seen, labelled GM 1–GM 5 from the outermost to innermost cortical layer, respectively. Layered structures were demonstrated, with contiguous laminae across the volume in all samples. Contiguity and consistency can be observed in all directions. Conversely, when the scalar maps are used as inputs to the cortical parcellation framework, the cluster maps do not capture contiguous cortical lamina and appear noisy. In addition, although five clusters were specified, the scalar-based segmentation was only able to robustly distinguish roughly three regions in each subject.

Figure 5A–C show representative slices of MRI (top)- and histology (bottom)-based cortical lamina segmentation from Subjects 1 to 3, respectively. In all cases, the MRI-based maps delineated contiguous cortical lamina based on T₁–T₂–MD spectral signatures, despite the clustering approach being agnostic to location in space. Visual inspection confirms that although not perfectly aligned due to limitations in inter-modality registration between the MRI and histological data sets, the unsupervised MRI- and histology-based laminar segmentations agreed. Nevertheless, the differences in the slice plane of the MRI and histology could contribute to the discrepancy in the apparent layer thicknesses.

Figure 5

Multimodal comparison of MRI- and histology-based cortical parcellation. Representative slices of MRI (top)- and histology (bottom)-based cortical lamina segmentation from Subjects 1 to 3 are shown in A–C, respectively. In all cases, the MRI-based maps enabled delineation of contiguous cortical lamina based on T₁–T₂–MD spectral signatures, although our clustering approach was agnostic to location in space. Qualitative agreement is evident from visual inspection. Limitations in inter-modality co-registration of MRI and histological data sets hinder optimal image alignment, which may result in slice plane differences, leading to discrepancy in the apparent layer thicknesses.

To assess the stability of cluster assignments over each iteration, the adjusted Rand index was calculated for each clustering compared to the optimal clustering, resulting in mean (range) values of 0.85 (0.39, 1.00), 0.99 (0.47, 1.00) and 0.9 (0.40, 1.00), for Subjects 1, 2 and 3, respectively, across iterations of k-means for the LOT distance. When using scalar values instead of distributions, the adjusted rand indices were 0.98 (0.46, 1.00), 0.82 (0.51, 1.00) and 0.91 (0.50, 1.00), respectively. The high mean adjusted Rand indices across all subjects point to the stability of the identified lamina clusters using the proposed unsupervised learning approach.

Layer-specific intra-cortical microstructure and architecture

Examining the cluster-specific T₁–T₂, T₂–MD and T₁–MD distributions averaged across the three subjects could facilitate biophysical interpretation of layer-specific intra-cortical microstructure and architecture. The mean distributions from each cluster, averaged across the entire sample volume from all three subjects, are visualized in Fig. 6 to investigate the spectral signatures of each cluster. From the outermost GM layer to WM, a gradual shift towards shorter T₁ and T₂ values, as well as the emergence of slow diffusivity and short T₁ components, is evident. Specifically, the macromolecular content-related (e.g. myelin water⁸⁹) short T₁ peak⁹⁰ was observed for all subjects in the T₁–T₂ and T₁–MD WM spectral signatures, consistent with previous multidimensional MRI findings in WM.^47,48,56 This short T₁ component gradually disappears as a function of radial distance (T₁–T₂ and T₁–MD spectral signatures in Fig. 6). Similarly, the characteristic slow diffusivity peak in WM^47,48 is gradually diminished as we move further into the cortex. The effect on T₂ can be seen in the T₁–T₂ and T₂–MD spectral signatures, in which longer values are captured further away from WM.

Figure 6

Layer-specific multidimensional MRI signatures. The mean distributions from each cluster, averaged across the entire sample volume from all three subjects in the image domain, are shown for (A) T₁–T₂, (B) T₁–MD and (C) T₂–MD. The average pairwise LOT distance between each cluster centroid is summarized in distance matrices, shown in the rightmost column. A gradual shift towards shorter T₁ and T₂ values, as well as the emergence of slow diffusivity and short T₁ components, is evident from the outermost GM layer to WM. From the distance matrices, it is evident that radially contiguous clusters have centroids that are more similar in the LOT metric space, and that T₁–T₂ and T₂–MD yield consistently higher inter-layer distances compared with T₁–MD. LOT, linear optimal transport; MD, mean diffusivity; WM, white matter; GM, grey matter.

The average pairwise LOT distance between each cluster centroid is summarized in distance matrices, shown in the rightmost column of Fig. 6. We observe that clusters that are radially contiguous have centroids that are more similar in the LOT metric space. For example, a gradual increase as a function of cortical depth in the LOT distance between the outermost cluster (GM 1) and the neighbouring clusters was observed in all subjects and for all the 2D spectral sets. We can also observe from these distance matrices that T₁–T₂ and T₂–MD yield consistently higher inter-layer distances compared with T₁–MD. This finding indicates that the T₁–T₂ and T₂–MD distributions contain more layer-specific intra-cortical information than T₁–MD, which supports the importance of T₂ in disentangling cortical architecture.

Discussion

As many neurological disorders result in perturbations of the cortical laminar architecture, this study sought to develop a noninvasive technique to probe the microstructural and architectural organization of the human cortex. By harnessing the rich information content of multidimensional diffusion–relaxation MRI and combining it with optimal mass transport-based unsupervised learning, we developed and demonstrated mapping of GM laminar composition in individual human brains ex vivo. We found that these MRI-derived segmented maps delineate contiguous cortical lamina based on multidimensional T₁–T₂–MD MRI spectral signatures, when trained using an unsupervised learning technique agnostic to location in space. We showed that our approach has the spatial sensitivity to map the cortical layers at the individual level in an unsupervised manner, and that the lamina maps were in good agreement with gold standard histology-based segmentation. These novel maps, together with histological validation, suggest that each cortical layer contains a unique microstructural and compositional signature that can be measured using multidimensional MRI.

Our goal in this study was to investigate whether cortical layers in the human brain produce distinct multidimensional MRI signatures and, if so, to design an unsupervised and automated proof-of-concept framework to distinguish and map them. We believe that unsupervised learning is most appropriate in this context as there is typically no ground truth available, especially in vivo when histological segmentation is impossible. Our proposed unsupervised learning approach using a modified version of the Wasserstein distance yielded stable cluster assignments across all subjects. Furthermore, despite lacking spatial coordinates, our approach revealed spatially contiguous lamina without human guidance. By moving through the cortex from subcortical WM to Layer I, our results show the gradual microstructural and compositional changes, reflected in the layer-specific T₁–T₂–MD spectral signatures (Fig. 6). As the myelin content diminishes from WM to the outermost GM layer, reductions in the relative fractions of the short T₁ and slow MD components were observed in all cases, along with a gradual shift in T₁ and T₂ towards longer values. The MD dimension measures orientationally averaged tissue water mobilities that is therefore rotationally invariant and, presumably, is not expected to vary with respect to cortical orientation or the choice of the laboratory coordinate system. While the MD in this study is modulated by coupled changes in cellular dimensions, shapes and orientations, the gradual disappearance of the slow MD component is likely indicative of increased cellular size, decreased anisotropy and reduced orientational coherence, from subcortical WM towards the cortical surface. In line with previous studies,^91-93 the decrease in both relaxation times and diffusivities at the deeper cortical levels likely indicates the increase in myelin content and microstructural changes, reflecting both cytoarchitecture and myeloarchitecture. We also found that the T₁–T₂ and T₂–MD distributions contain more layer-specific information than T₁–MD (Fig. 6), which supports the importance of T₂ in disentangling intra-cortical architecture.

The approach we propose here provides a direct link between intra-cortical organization and the MRI signal response and, importantly, showed that 1D T₁, T₂, DTI, NODDI and MAP-MRI measurements, even when combined and aggregated in a multimodal fashion, cannot disentangle the cortical layers (Fig. 4). And indeed, previous studies that aimed at cortical layer segmentation using either relaxation or dMRI were only partially successful, capturing changes in MR parameters as a function of cortical depth, but unable to fully disentangle all layers.^{9,11,69,93-96} The current study appears to be the first one to use multidimensional MRI, simultaneously encoding both relaxation and diffusion jointly, to address this problem.

We used cortical brain samples from three donors to enhance the generality of this proof-of-concept methodological study. Assessing the cytoarchitectonic classification of the cortical layers was done independently of MRI modality, based on the analysis of shape and organization of neuronal cells as observed using Nissl-stained histological sections. Including a gold standard cortical layer segmentation based on an independent modality established that the layers were preserved faithfully and also provided a qualitative reference to the proposed method. While visual inspection of the Nissl- and multidimensional MRI-based cortical segmentations suggests qualitative agreement (Fig. 5), a quantitative comparison and analysis was not possible due to a few reasons. First, cross-modality comparison relies on the accuracy of the co-registration between the microscopy and MRI data sets. In this case, we are interested in investigating the intra-cortical structures, and the cortical thickness precludes the use of diffeomorphic registration⁹⁷ that would be required for a more accurate co-registration. We used affine registration because diffeomorphic registration heavily distorted the relatively thin cortex and, in it, the cortical layers, which would significantly hinder our study. A second important confound in comparing the histological- and MRI-based results is that the sectioning planes may be different, which would lead to differences in the apparent layer thicknesses (see Supplementary Fig. 2). The last important difference relates to the type of classification; while Nissl stain provides cytoarchitectonic information, diffusion–relaxation MRI signal is modulated by both the organization of neuronal cells (cytoarchitecture) and the myelin content (myeloarchitecture).

This study represents the first step towards a more comprehensive noninvasive cortical parcellation of the human brain, and as such, it was an ex vivo study by design. Along with the crucial benefit of complementary and independent histological data, ex vivo human studies are affected by postmortem degeneration, fixation and resulting dehydration. While a direct comparison with in vivo data is not possible because the fixation process by itself changes the tissue properties (which affect the measured MRI parameters^98,99), we believe that the MRI signal trends we have uncovered here are of great value. Future ex vivo studies, in which the whole brain is imaged and processed, may yield additional information regarding how the cortical laminae structure change in different brain regions. Acquiring in vivo multidimensional MRI data in a clinical setting has become feasible following recent developments of clinical acquisition protocols by multiple groups.^52,58 Specifically in our context, in vivo clinical application of our method will be determined by spatial resolution, which in most cases is limited to 1–3 mm³. With an overall average human cortical thickness of ∼2.5 mm,¹⁰⁰ the spatial resolution presents a current barrier for widespread clinical translation. However, 7 T MRI scanners, which are becoming more widely available, offer ultra-high resolution^101,102 on the order of 0.5 mm³, which would allow to parse cortex laminar arrangement in vivo.

Diffusion–relaxation multidimensional MRI has been gaining momentum and been recently used to address challenging biological questions in cancer,^58,103,104 placenta dysfunction,⁵³ spinal cord injury^48,55 and traumatic brain injury.^56,60 While multiple ‘top-down’ studies of advanced MRI strategies to segment the human cortex suggest the potential clinical value of these methods, fewer ‘bottom-up’ MRI microscopy studies to directly associate MRI metrics with cortical layers in humans have been performed. By combining multidimensional diffusion–relaxation MRI, optimal mass transport–based unsupervised learning and histological methods, our ‘bottom-up’ study provides novel specificity for noninvasive human brain cortical parcellation. In the future, this technique may allow to follow the cortical laminar arrangement in normal development and aging, elucidate function–structure relationships in the brain and probe pathologies perturbing laminar architecture such as focal cortical dysplasia,¹⁰⁵ which can be clinically silent over decades.^106,107

Supplementary material

Supplementary material is available at Brain Communications online.

Acknowledgements

The authors thank C. Harker Rhodes for useful initial discussions about cortical parcellation.

Funding

This work utilized the computational resources of the National Institutes of Health (NIH) Biowulf cluster (http://hpc.nih.gov). This research was supported by the Intramural Research Program of the National institute on Aging and by the Intramural Research Program of the Eunice Kennedy Shriver National Institute of Child Health and Human Development.

Competing interests

The authors report no competing interests.

Data availability

The data and code that support the findings of this study are available from the corresponding author upon reasonable request.