Diffusion kurtosis imaging, MAP-MRI and NODDI selectively track gray matter myelin density in the primate cerebral cortex

Diffusion magnetic resonance imaging (dMRI) is increasingly used to study the development (Ouyang et al., 2019; Paydar et al., 2014), areal divisions (Avram et al., 2022; Howard et al., 2022; Liu et al., 2018), and pathology (Andica et al., 2020; McKavanagh et al., 2019; Torso et al., 2021) of the human cerebral cortex gray matter in-vivo. However, most investigation into the tissue basis of the dMRI signal has focused on the densely packed axonal projections of the white matter (Johansen-Berg & Behrens, 2013; Song et al., 2002, 2003, 2005), where the cell membrane is coated with myelin, a dense composite of lipids and proteins (Kandel et al., 2013; Morell & Quarles, 1999) that presents a strong barrier to the diffusion of water molecules. Much less is known about how dMRI parameters behave in gray matter, which is primarily composed of neuropil, a conglomeration of neurites, glial processes, cell bodies, vasculature, and other unmyelinated structures whose influence on diffusion has yet to be fully established (Jelescu et al., 2020; Novikov et al., 2019). Gray matter also contains a proportion of thin myelinated axons that are enmeshed within the neuropil, and whose structural organization varies by anatomical region and cortical layer (Bock et al., 2011; Braitenberg, 1962; Lewis & Van Essen, 2000; Vogt & Vogt, 1919). These diverse axon arrangements tend to be more diffuse than fiber bundles in the white matter (Lewis & Van Essen, 2000; Vogt & Vogt, 1919), and they may affect the dynamics of diffusing water molecules differently. Moreover, the aggregate dMRI signal in a cortical voxel simultaneously reflects the influence of both neuropil and myelinated axons. Ideally, we would like to distinguish the effects of the myelinated and non-myelinated compartments of gray matter tissue on diffusion, since these two tissue components are affected differently by disease.

Early work showed that the presence of myelin is not necessary for neurites to exhibit diffusion anisotropy (Beaulieu, 2002; Ono et al., 1995) and this was confirmed by dMRI experiments in the developing gray matter, where axons have not yet acquired their myelin sheath (Jespersen et al., 2007, 2012). Anisotropy occurs partly because the cell membrane reduces diffusivity across the neurites, and partly because there may be greater intrinsic diffusivity in the cytoplasm of neurites compared to the perikaryal cytoplasm or extracellular medium, so that diffusion is inherently faster along axons or dendrites (Beaulieu, 2002; Song et al., 2003). At the same time, foundational studies revealed that the presence of the myelin sheath around axons elevates anisotropy by further reducing diffusivity across their membranes (Jito et al., 2008; Song et al., 2002, 2005). These early results were obtained in isolated axons, or in selected regions of white matter tissue where large-caliber, myelinated neurites are uniformly aligned. However, myelinated axon fiber bundles traverse much of the human forebrain in all three spatial axes, so that diffusion anisotropies within them induce a heterogenous white matter vector field. This complicates the anatomical interpretation of the diffusion tensor because bundles often intersect within a white matter voxel (Stikov et al., 2011; Wheeler-Kingshott & Cercignani, 2009). More complex models of the dMRI signal often approach this problem by introducing assumptions that are specific to the white matter (Behrens et al., 2007; Zhang et al., 2012).

The mixture of neuropil and myelinated axons in cortical gray matter is a different and highly complex anatomical regime. It is a varied, yet structured environment for which the appropriate modeling assumptions have yet to be fully established (Jelescu et al., 2020; Novikov et al., 2019). For example, myelinated axons do not constitute the bulk of gray matter tissue (Schüz et al., 2002), nor do they form tightly packed bundles that intersect at unpredictable angles; nevertheless, the largest myelinated axons and non-myelinated dendrites often form a variety of interleaved columnar structures that run vertically from the white matter to the pial surface (Peters, 2010; Peters & Sethares, 1996; Rockland & Ichinohe, 2004), along which cell bodies tend to cluster (Buxhoeveden & Casanova, 2002; Peters, 2010). In contrast to the heterogenous vector field in the white matter, diffusion anisotropies in the gray matter tend to produce a vector field whose orientations are uniformly vertical, which may be connected to this columnar anatomy (Jespersen et al., 2010, 2012; McNab et al., 2013; Reveley et al., 2015).

In a recent study, we examined the relationship between the diffusion tensor and several gray matter histological parameters (Reveley et al., 2022). We found that the fractional anisotropy (FA) showed a very weak correlation to histological myelin density; however, it did track anisotropy in the cellular anatomy of the cortex, suggesting that neuropil influences gray matter diffusion. We found that radial diffusivity (RD) showed some correlation to cortical myelin density as might be expected (Song et al., 2002, 2005); however, RD was also strongly affected by the arrangement of myelinated axons in the tissue (Reveley et al., 2022). No tensor parameter tracked myelin density specifically, and separately from the arrangement of axons in the gray matter tissue. A number of non-diffusion MRI strategies have been developed to map myelin density in the cerebral cortex. These include the T1w signal (Bock et al., 2011; Glasser et al., 2014), myelin water fraction (Heath et al., 2018; Lazari & Lipp, 2021; Mancini et al., 2020), and the magnetization transfer ratio (MTR) (Grossman et al., 1994; Heath et al., 2018; Koenig, 1991; Lazari & Lipp, 2021; Mancini et al., 2020). Magnetization transfer imaging estimates the density of protons bound to large macromolecules based on the magnetization transfer from the free pool of protons to the bound (Grossman et al., 1994; Underhill et al., 2009; Wolff & Balaban, 1989). Previously, we found that the MTR signal was strongly correlated with myelin density in co-registered histological sections of the gray matter (Reveley et al., 2022). It would be of both theoretical value and practical benefit to find a diffusion parameter that also has this property.

The tensor model assumes that the distribution of diffusing molecule displacements is Gaussian. However, the physical barriers imposed by the tissue might yield different molecule displacement distributions. Since myelin forms a very strong barrier to diffusion, it may be that higher-order diffusion parameters beyond the diffusion tensor are strongly correlated to myelin levels in gray matter. Multiple modeling approaches have been developed to assess non-Gaussian diffusion profiles for dMRI data acquired with multiple b-values. The most conceptually straightforward of these is diffusion kurtosis imaging (Jensen et al., 2005), which offers scalar parameter maps of mean, radial, and axial kurtoses, summarizing deviation from a Gaussian displacement distribution. Previous studies of the white matter have reported a link between these parameter maps and myelin density (Falangola et al., 2014; Guglielmetti et al., 2016; Kelm et al., 2016). Kurtosis measures have also been linked to tissue “complexity” (Wu & Cheung, 2010) as well as to changing tissue constituents during cortical development (Jespersen et al., 2007; Ouyang et al., 2019) and a range of pathological conditions, including astrogliosis, glioma, stroke, and neurodegeneration (Zhuo & Gullapalli, 2020) that are not directly related to myelin levels. It is not yet fully established what tissue properties diffusion kurtosis, and other higher-order diffusion parameters, track in healthy adult gray matter.

Here, we investigate the relationship of diffusion parameters derived from a multishelled dMRI acquisition to the myelin distribution within the cortical gray matter, focusing on the mean diffusion kurtosis (MK) but also including the NG parameter map from the MAP-MRI (Özarslan et al., 2013) model and the NDI map from NODDI (Zhang et al., 2012). We acquired multi-shelled dMRI scans of the fixed brains of two common marmosets (Callithrix jacchus) at isotropic spatial resolutions of 80–150 µm. Tissue myelin density was estimated with MTR scans in the same brains at 75–80 µm. We subsequently sectioned, stained for myelin, and non-linearly co-registered the tissue of one of the brains using a careful, manually guided process, which allowed for a pointwise comparison of the dMRI parameter maps and tissue myelin levels across the gray matter. We found that MK, as well as the NODDI NDI map and the MAP-MRI NG parameter, closely tracked the myelin content in the cerebral cortex, as assessed both by direct histological comparison and comparison to precisely registered MTR data, used as a surrogate for myelin.

The postmortem brains of two adult marmosets “case M” and “case P”, both used in previous studies (Liu et al., 2018, 2020; Reveley et al., 2022), were employed in this work. Both animals were scanned for multi-shelled diffusion MRI (dMRI) and for Magnetization Transfer Ratio (MTR) MRI on the same 7T Bruker preclinical scanner with a 30 cm bore and 450 mT/m gradients. Case M was scanned and then carefully co-registered to an atlas (Paxinos et al., 2013) of the marmoset cerebral cortex as part of the marmoset brain mapping project (Liu et al., 2018, 2020). The dMRI was obtained at a resolution of 80 µm isotropic in a whole brain scan lasting 15 days, conducted for that project (Liu et al., 2020). The left hemisphere of “case P” was scanned specifically for the present study at a resolution of 150 µm isotropic. As part of a previous study (Reveley et al., 2022), case P was sectioned for Gallyas myelin histology (Gallyas, 1979), which was then non-linearly registered to the dMRI data. All the procedures and provision of materials for this study were in full compliance with the Guidelines for the Care and Use of Laboratory Animals by National Institute of Health and approved by the Animal Care and Use Committee of the National Institute of Mental Health.

Detailed scan information can be found in Liu et al. (2020). Briefly, the brain of case M (male common marmoset, 4.5 years old “Case M”) was soaked with 0.2% gadolinium (1 mmol ml^–1 Gadavist, Bayer) in 1× PBS for 2 weeks and then with 0.2% gadolinium in pure water for 1 day. Case M was then scanned using a coil specially tailored for ex-vivo scanning of marmoset brains. 204 dMRI volumes were acquired at 80 µm isotropic resolution in four shells (8 b = 0, 6 b = 30, 64 b = 2,400, and 126 b = 4,800). The pulse width for the diffusion weighting gradient was 6.4 ms and pulse separation was 14 ms, for a diffusion time of approximately 14 ms. The scan time was 15 days.

MTR data were acquired was acquired for case M, also at 80 µm isotropic resolution using a 3D FLASH sequence. Five MTR scans were acquired, each comprising two volumes with (M_sat) and without (M_off) an offset magnetization transfer (±2,000 Hz off-resonance, Gaussian-shaped), which were averaged before calculating the MTR. The MTR value is calculated as 100(1-M_sat/M_off). All five MTR scans were then averaged. The total acquisition time was 50 hours.

Before MRI scanning, the formalin-fixed marmoset brain of case P was soaked with 0.15% gadopentetate dimeglumine (Magnevist, Bayer Leverkusen, Germany) for 3 weeks to reduce the T1 relaxation time as described in a previous study (Liu et al., 2018; Reveley et al., 2022). In order to improve SNR, brain P was cut through the midline. Only the left hemisphere was used in the present study. The left hemisphere was scanned with a 25-mm birdcage volume coil (Bruker). The imaging parameters were: TR = 470 ms, TE = 34 ms, flip angle = 90°, FOV = 38.4 x 24 x 21 mm, matrix size = 256 x 160 x 140, and resolution = 0.15 mm isotropic. The pulse width for the diffusion weighting gradient was 8 ms and pulse separation was 18 ms, for a diffusion time of approximately 18 ms. In order to keep the echo time short, 10 segments were used in the EPI acquisition. The DWI data were acquired using 2 averages, and the acquisition time for each 3D volume was about 22 minutes. In total, 142 DWI volumes were collected to sample the q-space on 5-shells. The 5 shells were defined by b-values: 419, 1,677, 3,773, 6,708, and 10,481 s/mm², obtained by setting the gradient magnitude at 80, 160, 240, 320, and 400 mT/m. 3, 21, 21, 37, and 58 volumes were collected in each shell respectively. Two volumes were collected at b = 0.

The left hemisphere of case P underwent scanning for MTR at 75 µm as described in a previous study (Reveley et al., 2022). Briefly, the MTR images were collected with a 3D FLASH sequence, FOV = 34.8 x 20.4 x 20.1 mm, and matrix size = 464 x 272 x 268. As with case M, the magnetization transfer was achieved in M_sat scans by a Gaussian-shaped pulse at +2,000 Hz offset frequency in one scan and at -2,000 Hz in the other scan of a paired acquisition. This magnetization transfer pulse was turned off in M_off scans. The MTR value was calculated as 100(1-M_sat/M_off). The number of averages was 5. The total MTR acquisition time for the left hemisphere of case P was about 60 hours.

Here, we describe MRI pre-processing steps and computation of diffusion tensor, kurtosis, NODDI, and MAP-MRI parameters in each ex-vivo brain.

We used the publicly available 80 µm MTR and dMRI data from case M, which was previously preprocessed (Liu et al., 2020) using the DIFFPREP of TORTOISE (Pierpaoli et al., 2010). We removed the left hemisphere and closely cropped the right hemisphere using fslroi, part of FSL 6.0.5.2 (Woolrich et al., 2009), to save computation, since only the right hemisphere was aligned to the atlas labels. We calculated the mean, axial, and radial kurtoses of the 80 µm data using dtifit part of FSL 6.0.5.2 (Woolrich et al., 2009). For the remaining parameter maps, we downsampled the diffusion data from 80 µm to 160 µm using the “fslmaths -subsamp2offc” command in FSL 6.0.5.2 (Woolrich et al., 2009). We computed the mean kurtosis a second time with “RobustDKIFitting” MATLAB script (Henriques et al., 2021) and compared this to the DTIFIT results, finding good agreement. We computed the NODDI model using AMICO (Daducci et al., 2015) version 1.5.4 with default diffusivity parameters and the “ex-vivo” parameter set to 1. We estimated the MAP model order 4 using TORTOISE (Pierpaoli et al., 2010) version 3.2.0, EstimateMAPMRI with default parameters.

For case P, we pre-processed the data using the DIFFPREP of TORTOISE (Pierpaoli et al., 2010). We registered MTR data to the histology slice plane and as described in a previous publication (Reveley et al., 2022), and performed a rigid body (6 parameter) registration of the dMRI data to the MTR using FSL 6.0.5.2 FLIRT (Woolrich et al., 2009), thus upsampling the dMRI from 150 to 75 µm with spline interpolation. We estimated the MAP model at order 4 using TORTOISE (Pierpaoli et al., 2010) version 3.2.0, EstimateMAPMRI with default parameters. We computed the NODDI model using the NODDI toolbox 1.4 (Zhang et al., 2012), with parameters disoIdx 2.0E-9, diIdx 0.4E-9. We calculated the mean, axial, and radial kurtoses of the upsampled 75 µm data using DTIFIT part of FSL 6.0.5.2 (Woolrich et al., 2009), and compared these results to those of the “RobustDKIFitting” MATLAB script (Henriques et al., 2021), finding good agreement.

The cortex of case M was carefully labeled and aligned to the Paxinos atlas (Paxinos et al., 2013) in previous work (Liu et al., 2018, 2020). We removed the left hemisphere of the atlas MTR template and atlas labels since these were a mirror image of the right. We fine-tuned the alignment of the atlas labels and MRI data using a shift/translate only FLIRT transform with no interpolation. We used the medium granularity labels of the marmoset brain mapping labels (Liu et al., 2018) in this work. We manually trimmed the anatomical labels at the white/gray and pial boundaries to avoid partial volume effects and non-brain voxels. The white matter voxels were defined as part of the marmoset brain project white matter atlas (Liu et al., 2020). We downsampled each of the three label volumes released as part of the project to 160 µm using the “fslmaths -subsamp2offc” command in FSL 6.0.5.2 (Woolrich et al., 2009). We manually eroded each label slightly to avoid partial volume effects and combined the three volumes into a single mask file.

For estimation of cytoarchitectural boundaries, we non-linearly registered the right hemisphere of case M MTR to the left hemisphere of case P MTR using ANTS 2.4.1 (Avants et al., 2011), as described in a previous study (Reveley et al., 2022). We then manually entered the coordinates of the region boundaries to the analysis code described below.

Myelin histology was obtained and processed from case P only as described in a previous study (Reveley et al., 2022). Briefly, the left hemisphere of case P was sectioned at 50 µm on a freezing microtome. We used the Gallyas silver stain (Gallyas, 1979) which stains individual axons at high resolution. Each section was digitally scanned at 0.88 μm resolution (×10 magnification) using a Zeiss Axioscan microscope slide scanner. Seven Gallyas stained sections were selected for nonlinear registration based on criteria of even staining and minimal damage from a previous study (Reveley et al., 2022).

To estimate the anisotropy content of the myelin histology we employed the structure tensor (ST) method (Budde & Frank, 2012), as described in a previous study (Reveley et al., 2022). Briefly, the 0.88 μm per pixel resolution histology images were processed using the structure tensor as implemented in the “OrientationJ” ImageJ package using a local Gaussian window σ = 85 pixels and image gradients computed with a cubic spline. Since the dMRI data contain diffusion information gathered over a 150 μm cubic area, the histology structure tensor at each pixel was a distance-weighted function of the structure tensors within a radius of 75 μm. For each pixel, the structure tensor coherence was computed from the ST eigenvalues.

The 75 µm MTR was manually rotated to the histology slice plane and the dMRI data registered to this using FLIRT as described above. Histology parameter maps of stain density and anisotropy at 75 µm resolution were then non-linearly registered to the MRI. We performed careful nonlinear alignment of seven local frontal and parietal ROI from the histology to the MRI. Stain intensity maps of the 75 µm maps were first binarized using a threshold adjustment in imageJ (Schindelin et al., 2012) and then registered to their corresponding MRI slices to correct for shrinkage and deformation, first using manual 2D alignment in imageJ and then MIND (Heinrich et al., 2012) MATLAB code (α = 0.1) as described in a previous study (Reveley et al., 2022).

The data in Figure 2 were plotted together in MATLAB 2023b at a resolution of 80 µm by concatenating a mask of all the gray matter regions obtained above (“Anatomical Mapping”), while those in Figure 4 used 160 µm labels. For Figure 3, each region label was evaluated separately.

Once rotated into the histology slice plane, each coronal slice of the dMRI parameter maps acquired for this study was converted to quantitative tif format and integrated into a software system written for a previous study on MRI histology comparisons (Reveley et al., 2022). Briefly, this system took advantage of the fact that cortex consists of vertically oriented tissue components to produce a matrix representation of the non-linearly registered dMRI and histology data whose columns were derived from data sampled along the vertical dMRI tractography and whose rows spanned the horizontal direction. Tractography lines were integrated through the 2 d projected principal eigenvectors in each gray matter ROI. A vector of unique pixels whose coordinates intersected one streamlines was acquired, each of which was then resampled to a fixed length using bilinear interpolation, to yield a simple matrix representation of the gray matter. This operation was invertible such that data in matrix form could be displayed in image space.

Our statistical analysis proceeded by parcellating these matrix representations of the nonlinearly registered histology and MRI data. The parcellation strategy served to separate out regions of potentially differing stain intensity, and to assess relationships between histology and MRI variables in different anatomical regions in a statistically independent manner. For each parcel we took the Pearson correlation between sample variables, in, for example, Figures 5C and 6A. To assess the average laminar properties of histology or dMRI variables, we averaged the rows of the matrix representation and plotted the resulting vector, for example, in Figures 5B and 6B.

Our principal goal was to determine which diffusion MRI (dMRI) parameters most closely track the density of myelin in the gray matter tissue of the cerebral cortex. To this end, we analyzed high-resolution, ex-vivo dMRI data from the brains of two marmoset monkeys which were collectively scanned for over 3 weeks at 7T (see Section 2). We then compared this dMRI data to the myelin content of the tissue, first as assessed by a non-diffusion MRI sequence, and second as assessed by myelin-stained histology gathered from one brain, and carefully co-registered to the MRI (Reveley et al., 2022). We employed the magnetization transfer ratio as our non-diffusion MRI measure, since this sequence is known to be sensitive to large macromolecules such as those comprising myelin (Grossman et al., 1994; Heath et al., 2018; Henkelman et al., 2001). In previous work, we showed that MTR intensity closely matched the fine-scale distribution of myelin in the cortical sheet, using histology data gathered from the scanned tissue (Reveley et al., 2022).

We began by re-establishing that MTR is a valid surrogate measure for myelin intensity in the cerebral cortex. To do this, we directly compared the spatial distribution of MTR intensity to seven non-linearly registered myelin-stained sections obtained from the scanned brain. This comparison is shown for one segment of the adult marmoset cerebral cortex in Figure 1A, where the myelin-stained histological section (top) has been resampled to match the 75 µm isotropic resolution of the MTR scan (bottom). The myelin in the histological sections and the MTR signal in the MRI data followed a similar gradient, showing highest intensity in deep cortical layers and declining toward the pial surface. In addition to this laminar correspondence, the two modalities also matched in their tangential cortical variation. An obvious example of this is the higher signal intensity in the middle temporal (MT) area, which is known for its high myelin content. For this registered section, Figure 1B shows a pixel-by-pixel comparison of the two modalities, demonstrating a strong correlation (ρ = 0.94). Applying this analysis to multiple sections, we found a similarly strong correlation between MTR and histological data (Supplementary Fig. 1), confirming that the high-resolution MTR acquisition from the gray matter of fixed brains could serve as a reliable surrogate for cortical myelin distribution.

With this tool in hand, we focused on our target question of which diffusion parameters might accurately reflect cortical myelin content, a potentially important use of dMRI for clinical neurology. Given that neither FA nor RD are effective predictors of gray matter myelin density (Reveley et al., 2022), we investigated the hypothesis that dMRI model parameters reflecting non-Gaussian diffusion may track it more closely. The idea that myelin restricts, rather than merely slows diffusion across axons, has been discussed in previous diffusion studies (Guglielmetti et al., 2016; Jelescu et al., 2016; Jespersen et al., 2010; Kelm et al., 2016). In this study, we compare cortical myelin density captured by MTR and by finely co-registered histology to higher-order diffusion parameters, using an ultra-high spatial resolution that captures the detailed laminar and tangential variation of myelin content across large areas of the primate cortical sheet.

We first compared the MK mapped across the whole cerebral cortex with the myelin content as estimated by the MTR scan. Since both maps are MRI measures from the same ex-vivo sample, the spatial registration was very precise. Importantly, the MK and MTR measures capture fundamentally different physical processes (Henkelman et al., 2001). Nonetheless, we found the two measures exhibited strong similarity in the monkey cortex.

Figure 2A and B illustrate this similarity in representative coronal, sagittal, and axial slices at 80 μm isotropic resolution. Within these sections, the close correspondence in both the laminar and tangential features is evident. In the laminar direction, the MK and MTR measures both exhibited the expected myelin gradient, showing high values in the lower layers near the white matter and decreasing values toward the pia. In the tangential direction, the changes across cortical areas were also similar. For example, the elevated MTR in MT (middle temporal cortex) and M1 (primary motor cortex) were reflected by elevated MK in the same areas, perhaps due to putative restricted diffusion within the relatively more numerous axonal myelin sheaths in these cortical regions.

We quantified the pixelwise relationship between MK and MTR (Fig. 2C), which revealed a strong correlation in the gray matter (ρ = 0.81, black Fig. 2C). The mean kurtosis measure, therefore, appeared sensitive to variation in the vertical and horizontal cortical myelin distribution, in a similar way to the MTR. Interestingly, when this analysis was extended to the white matter, with its higher overall myelin content due to dense composition of myelinated axons, the correlation remained strong (ρ = 0.82 red inset Fig. 2C). White matter pixels with higher MTR values, reflecting a higher spin magnetization transfer from the bound pool of protons to the free pool water protons, also exhibited higher MK, reflecting more restricted water diffusion. This robust agreement between diffusion and large-molecule MRI measures in the precisely registered maps of the same brain suggests that mean diffusion kurtosis is a strong marker for myelin content. At the same time, the strength of the gray matter correlation was slightly weaker than that between myelin histology and MTR, shown in Figure 1B. To investigate this discrepancy further, we examined the regional variation in the relationship between MK and MTR.

We divided the cortex into anatomical regions based on the areal definitions of the Paxinos atlas (Paxinos et al., 2013) and marmoset brain mapping project (Liu et al., 2018, 2020), using the registered maps provided by a previous study (Liu et al., 2018, 2020). We found that the correlation coefficient between MK and MTR varied across cortical regions (Fig. 3A). In contrast to a surface-based analysis of values averaged through the cortical depth, our voxel-based methods at ultra-high-resolution captured rich laminar detail of the cortex (Fig. 3B). Nevertheless, the variation in regional correlation coefficients was related to each region’s total myelin content. For example, in cortical regions where the mean MTR value was high, and which were known to have high myelin content, for example M1 or MT (Fig. 3B) correlation levels approached 1, whereas regions with low MTR, such as the piriform cortex (Pir), exhibited lower correlations, closer to 0.6. The lower correlations may be linked to a decline in contrast relative to noise where myelin content is lower, rather than a failure of either measure to track myelin content (Supplementary Fig. 2).

We next expanded our investigation of dMRI and cortical myelin to include two further popular diffusion models, namely the non-Gaussian (NG) component of the mean apparent propagator (MAP-MRI) MRI model (Özarslan et al., 2013) and the neurite density index (NDI) calculated as part of neurite orientation dispersion and density imaging (NODDI) (Zhang et al., 2012). The NG parameter computed from MAP-MRI is conceptually similar to mean kurtosis. However, NDI is not a measure of non-Gaussian diffusion per se, but is instead part of NODDI, a biophysically motivated tissue model of Gaussian diffusion that seeks to estimate the variation in neurite density (Zhang et al., 2012). In NODDI, contributing gray matter neurites can potentially include both myelinated and unmyelinated axons, as well as dendrites, all of which are modeled as sticks with negligible radial diffusivity (Zhang et al., 2012).

We found that all three of the dMRI parameter maps we considered (MK, NG, and NDI) resembled the MTR and one another in the cortical gray matter but differed substantially from the FA maps of the diffusion tensor (Fig. 4A). The correlation between the three dMRI measures of interest was extremely high, with the correlation between MK and NDI at 0.86 and between MK and NG at 0.92. The correlations with each of these parameters to the MTR maps were also high, albeit with the correlation of NDI to MTR taking a slightly lower value (0.73) than those of MK and NG (0.81 and 0.8 respectively). By contrast, the relationship of FA to the three putatively myelin-sensitive dMRI measures, and to MTR, were all much lower, between 0.26 and 0.3 (Fig. 4B). When we examined the spatial relationships of dMRI parameters to MTR, we found that MK and NG exhibited a nearly identical, roughly linear voxel-wise relationship to MTR, and NODDI showed a slightly different but still similar relationship. The relationship of FA to MTR, however, was much weaker (Fig. 4C).

Having established the close relationship between MTR and the diffusion parameters MK, NG, and NDI over the whole brain, we next focused on a fine-grained comparison of the dMRI parameters with myelin histology obtained from the same brain. We compared seven non-linearly registered sections of myelin histology to matching dMRI and MTR MRI. In the example scatter plots from one section shown in Figure 5A, the MRI parameters of interest all showed a clear linear relationship with histological myelin, with MTR, MK, and NG showing Pearson correlations of 0.85, 0.82, and 0.84 respectively. NDI showed a slightly weaker correlation of 0.75. FA showed a lower overall correlation and a more complex, piecewise relationship to myelin density, which we previously reported (Reveley et al., 2022).

Continuing to investigate the histological basis of the dMRI parameters, we next examined their laminar variation, their variation by cortical area, and how both of these tracked histological myelin content. To assess the similarity in laminar variation, we sampled the cortex along columnar lines derived from 2D tractography and resampled the resulting data into a matrix formulation (see Section 2). This allowed us to calculate and compare the mean laminar profile of the various MRI measures to the co-registered myelin histology (Fig. 5B). We found that NG and MK exhibited a laminar profile that was nearly identical to that of myelin intensity, declining monotonically from deep to superficial layers. However, as expected, the laminar profile of FA differed substantially from that of myelin, and from the other diffusion parameters, with minima at the white/gray and pial boundaries and a peak in middle cortical layers. This laminar profile reflects the influence of neurite organization as well as density, and non-myelinated components as well as myelinated (Reveley et al., 2022).

To provide a more granular analysis, we next sampled the matrix representation of the cortex into columnar parcels (Reveley et al., 2022) running vertically from the white matter to the pia and with a horizontal width of about 300 µm (Fig. 5C). Within each of these parcels, we measured the correlation of the MRI parameters with histological myelin. This analysis showed a high spatial correlation of voxels within most columnar parcels, with the highest median values observed for MTR (0.95) and MK (0.94), slightly lower values for NG (0.92) and NDI (0.89), but much lower values observed for FA (0.42).

Together, these analyses reveal that the diffusion-based measures—obtained either by modeling the spatial distribution of diffusion (MK, NG) or by biophysical forward modeling of the neurite distribution (NDI)—were all good measures of the cortical myelin distribution. This result is especially clear when the higher-order parameters are contrasted to FA, which also reflects the influence of non-myelinated components in gray matter.

The above results show that several higher-order models obtained from multi-shelled diffusion data closely track myelin content in the cortex. Although FA in the cortex is not specific to myelinated structures (Beaulieu, 2002; Reveley et al., 2022), radial diffusivity (RD), another measure derived from the diffusion tensor, has shown more specificity to myelin in both gray (Reveley et al., 2022) and white matter (Song et al., 2002). In some studies, RD has been shown to correlate with histological measures of myelin content and has been used as a marker for assessing myelin integrity and alterations in white matter microstructure (Lazari & Lipp, 2021; Mancini et al., 2020). Previously, we found that RD tended to track myelin organization rather than myelin intensity in gray matter (Reveley et al., 2022). We next examined the relationships between RD and myelin, using the same approach as above (Fig. 5).

We found that the laminar profile of RD diverged significantly from that of myelin histology (Fig. 6A, B), but more closely matched that of myelin histological anisotropy. Although RD did show a correlation to myelin intensity, this correlation is notably weaker than MK’s, exhibiting a much wider spread of correlation coefficients between columnar parcels of cortex (Supplementary Fig. 3). We found that while the spatial pattern of RD was very sensitive to the directional content (Reveley et al., 2022) (Supplementary Fig. 3; Supplementary Text 1) of the myelinated axons as well as their density, it did not serve as a good quantitative measure of myelin content in the same way as mean kurtosis.

There is wide regional variation in cytoarchitecture, myeloarchitecture, and neurite structure over the expanse of the cortical sheet. We were interested whether MK could assess the laminar structure of the primary visual cortex (V1), a well-studied cortical area with a unique laminar patterning of cell density and myelin (Balaram & Kaas, 2014; Casagrande & Kaas, 1994). Figure 7 demonstrates the capacity of MTR, MK, NDI, and FA to reveal aspects of V1’s laminar architecture in the postmortem sample. Importantly, the laminae of highest intensity did not match among all the registered MRI measures. MTR, MK, and NDI matched well and showed their highest levels in layer 4B, the so-called Gennari band, reflecting its known high myelin content (Balaram & Kaas, 2014; Casagrande & Kaas, 1994) (Fig. 7A).

The abundant myelinated axons in 4B are visible in the corresponding histology to have a matted structure (Fig. 7B) and this lack of spatial coherence explains the low FA in this layer, and also in the deepest layer of cortex adjacent to the white matter. By contrast, layer 4°C was revealed by a high-intensity band of FA, which likely reflects an abundance of vertically oriented, non-myelinated neurites in this layer (Fig. 7A, right) (Reveley et al., 2022). Additionally, we found that radial diffusivity was relatively lower in layer 4°C where the myelinated axons are coherently organized, than in 4B where they are matted, despite the higher myelin content of the Gennari band in 4B (Supplementary Fig. 4; Supplementary Text 2). This likely reflects the sensitivity of RD to the organization of the myelinated fibers as shown in Figure 6B. Together, these results indicate that models based on multi-shelled diffusion acquisitions can distinguish and quantify the diffusion profiles of myelinated and non-myelinated neurites in the gray matter at a laminar level.

In this study, we assessed the relationship between cortical myelin content and higher-order diffusion MRI model parameters based on multi-shelled acquisitions. We analyzed ultra-high-resolution dMRI and MTR scans of whole primate brains, as well as carefully co-registered histological data covering large areas of the frontal and parietal lobes. We examined the differences between cortical regions, as well as the laminar properties of the parameters within those regions. Our guiding hypothesis was that water restricted to the intra-axonal space by the myelin sheath would be reflected in a non-Gaussian molecule displacement profile, so we initially focused our investigation on mean kurtosis (MK). MK strongly correlated with MTR—a separate MRI contrast unrelated to diffusion, but closely linked to myelin levels—globally over the whole cortical sheet. Regions with the highest myelin content, and therefore the highest contrast-to-noise ratio, exhibited the strongest correlation with MK. Broadly similar results were obtained with the non-Gaussian (NG) parameter map obtained from MAP-MRI (Özarslan et al., 2013), and neurite density index (NDI) parameter map obtained from the NODDI biophysical model (Zhang et al., 2012). Each of these high-order diffusion measures closely tracked myelin, as well as one another. Specifically, they matched the monotonic decay of histological myelin levels from the white gray boundary to the pial surface, as well as the tangential variation across the cortical sheet. In the primary visual cortex, the MTR and MK measures captured the known distribution of myelin fibers in layer 4B. Of the three higher-order diffusion parameters we tested, NDI was unlike MK and NG since it does not explicitly reflect non-Gaussian diffusion processes. Potentially, both myelinated and unmyelinated neurites are intended to drive the NDI contrast (Zhang et al., 2012), but our analysis shows that NDI behaves similarly to MK, a measure of non-Gaussianity, in the environment of fixed primate gray matter.

The behavior of MK, NG, and NDI may be contrasted with the behavior of the Gaussian diffusion tensor in gray matter. Radial diffusivity (RD) exhibited sensitivity to the organization of the myelinated fibers as well as their density, rather than being attuned to myelin density alone. Fractional anisotropy (FA) showed only weak correspondence to myelin parameters, and in previous work we showed that FA is strongly driven by non-myelinated tissue components. Employing both the diffusion tensor and higher-order diffusion terms may offer a way to separate out the specific influences of myelinated and non-myelinated tissue components on the diffusion signal. Moreover, we found that MTR and MK were both myelin sensitive, despite obtaining their contrasts via different mechanisms. Thus, although they have correlated values in healthy tissue, MK and MTR might diverge for pathological cases in specific ways that have diagnostic value (Edwards et al., 2018; Kamiya et al., 2020).

Our results cohere with and draw together several lines of converging evidence (Fukutomi et al., 2018; Guglielmetti et al., 2016; Jelescu et al., 2016; Kelm et al., 2016; Yoshida et al., 2013) at different spatial scales, to suggest that mean kurtosis has both sensitivity and specificity to myelin levels in healthy gray matter.

Mean kurtosis values have previously been shown to track changing myelination levels in vivo. For example, Guglielmetti et al. (2016) studied dMRI and immunohistochemical material from mice which had been treated with cuprizone, a drug that reversibly induces demyelination through toxicity to oligodendrocytes. The authors found that mean kurtosis decreased in the gray matter of motor and somatosensory cortex of the mouse following cuprizone administration and regained their values on recovery from the treatment. This was true except in areas of the cortex where the histology showed remyelination to be incomplete. Likewise, in a developmental study, Praet et al. (2018) assessed the correlation of mean kurtosis based on the optical density of co-registered myelin basic protein immunohistochemistry. They found that as myelin levels in the motor cortex increased during development, cortical MK tended to increase, and FA tended to decrease. These rodent studies emphasized longitudinal changes in experimental cohorts, rather than detail within the cortex of particular subjects.

A few studies have examined the distribution of higher-order diffusion metrics in groups of human subjects at clinical resolutions over the whole brain. Bester et al. (2015) examined diffusion kurtosis in a cohort of subjects suffering from multiple sclerosis (MS), a demyelinating disease. They found that MK was significantly lower in the gray matter of MS patients compared to controls, and that MK was inversely correlated with executive function where poor function presumably reflects cortical demyelination. Fukutomi et al. (2018) found that NDI was strongly correlated to myelin density as estimated by the T1w/T2w ratio, in a group-based cortical surface-based study of 505 healthy subjects from the human connectome project. These studies both suggest that the present findings in an ex-vivo primate model are strongly relevant to human cortical tissue scanned in-vivo, at much lower spatial resolutions.

Zhu et al. (2021) compared MK in a cohort of six postmortem macaque brains at lower spatial resolutions of around 600 µm in plane, gathered in 2 mm slices, to publicly available neurofilament (SMI-32) immunohistochemistry from a different animal. They found qualitative correspondence between the two measures over the cortex as a whole; however, they were not able to examine the histology/dMRI relationship at a fine spatial scale in the same brains. SMI-32 is reactive in a subset of large pyramidal cells in the cortical gray matter and is closely associated with myelination (Kirkcaldie et al., 2002), because those same cells give rise to the largest myelinated axons in cortex. The broad correspondence Zhu et al. (2021) noted at comparatively low spatial resolution may be a proxy for average cortical myelin levels, rather than indicating a direct link between MK and the neurofilaments themselves. More detailed work is required to address this question.

A handful of studies (Lazari & Lipp, 2021) have examined the relationship of diffusion kurtosis to the detailed properties of myelinated axons in murine white matter using electron microscopy. Kelm et al. (2016) studied the relation of myelin to diffusion kurtosis in an ex-vivo cohort of knockout and control mice with varying degrees of hypomyelination. They found that MK correlated to the fraction of myelin as estimated from the electron microscopy data in four white matter ROI across their subject pool. Jelescu et al. (2016) also found that MK exhibited sensitivity to cuprizone-induced demyelination in an analysis of electron microscopy in the mouse corpus callosum.

Connecting MRI data to underlying tissue parameters requires certain assumptions, so it is important to identify where such assumptions might break down. Importantly, while the focus of the present study was on myelin, both MTR and MK may also have sensitivity to other tissue properties. For example, MTR is known to vary with water content, as brought about by edema, inflammation, and other immune responses (Heath et al., 2018; Vavasour et al., 2011). Likewise, beyond marking axonal myelination, MK varies with other factors that increase tissue “complexity”, such as gliosis or gliomas (Wu & Cheung, 2010). These caveats may be particularly important in the translation of the ex-vivo data presented here to the living brain. In this study, measures of diffusion kurtosis or non-Gaussian diffusion were found to correlate to the tissue myelin content. We presume this is due to restricted water diffusion in the myelinated axons, but it is worth mentioning that the relationship is empirical and has limitations in theory. For example, in a given region of brain, there could be more than one way for increased myelin content to occur—thicker myelin on the same number of axons or more myelinated axons per volume. Our results do not demonstrate whether MK can distinguish these scenarios; with that said, we did find empirically that MK closely tracked both MTR and histological myelin density in the gray matter, which contains axons of a wide variety of calibers. Isolating any myelin-unique component of MRI signals is particularly difficult in the cortical gray matter, because multiple tissue parameters correlate spatially with one another. For example, in the cortex, the largest myelinated and non-myelinated neurites both have a strong vertical orientation bias and are spatially enmeshed (Buxhoeveden & Casanova, 2002; Peters et al., 1997; Peters & Sethares, 1996). Nonetheless, our demonstration of a laminar gradient of MK in the cortex does suggest that myelin is the determining influence, since this is a signature of cortical myelin levels (Braitenberg, 1962), whereas the density of unmyelinated neurites is typically high in superficial cortical layers, near the pia (Peters, 2010; Peters et al., 1997).

Diffusion parameters are sensitive to the state of the tissue. In the present study, we exclusively addressed fixed, ex-vivo tissue scanned at diffusion times in the range 10–20 ms (see Section 2). Acquisition in living tissue might give contrasting results, since both diffusivities and higher-order diffusion measures differ in living tissue. Likewise, while the MTR method is suitable for ex-vivo studies, where T1 times are equalized by Gd doping, T1-weighted imaging is likely more useful for in-vivo studies, due to its higher contrast-to-noise ratio and lower specific absorption rate (SAR) exposure to patients. We note that Fukutomi et al. (2018) found good agreement between NODDI NDI and myelin as estimated by the T1w/T2w ratio over the whole human brain in-vivo; so far, the ex-vivo and in-vivo findings are supportive to each other. Work to assess these issues and develop sophisticated, empirically grounded models of complex diffusion processes in gray matter is ongoing (Jelescu et al., 2022; Lee et al., 2020).

The data employed here were used in two previous studies (Liu et al., 2018, 2020; Reveley et al., 2022).

C.R. conceived the project, analyzed the data, and wrote the paper. F.Q.Y. acquired the data, advised on analyses, and edited the paper. D.A.L. advised on analyses, edited the paper.

The authors declare no competing interests.

The research was supported by the Intramural Research Program of the National Institute of Mental Health (ZIAMH002786 and ZICMH002899).

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00368.