1. Introduction

By quantifying the microscopic motions of water molecules diffusion MRI (dMRI) provides a sensitive clinical tool to non-invasively probe the tissue structures at length scales (≈5μm) much smaller than the voxel size. In isotropic and anisotropic tissues, the dMRI signal at low diffusion sensitizations (b-values) can be described phenomenologically using diffusion tensor imaging (DTI) (Basser et al., 1994a,b). In DTI, the diffusion signal attenuation in each voxel is modeled using a diffusion tensor, D, which has 6 degrees of freedom. The diffusion tensor can be decomposed or diagonalized in an orthogonal reference frame whose principal coordinate axes are characterized by the eigenvectors ϵ1, ϵ2, and ϵ3. The normalized orthogonal unit vectors along the principal tensor axes represent 3 degrees of freedom of D that define its orientation with respect to the laboratory reference frame. The scalar principal diffusivities λ1, λ2, λ3 corresponding to these directions represent the other 3 degrees of freedom of D and determine the mean diffusivity and diffusion anisotropy. In general, D can be written as:

D=λ1ϵ1ϵ1T+λ2ϵ2ϵ2T+λ3ϵ3ϵ3T,    (1)

where ϵ1ϵ1T,ϵ2ϵ2T,ϵ3ϵ3T are the principal coordinate axes dyads (or rank-1 tensors) derived from the eigenvectors of the diffusion tensor while the positivity of the principal diffusivities (i.e., eigenvalues of D) guarantees that D is positive definite.

However, at b-values larger than 1, 500s/mm2 the dMRI tissue signal is more sensitive to the intravoxel variation of water diffusion properties, and the DTI approximation may no longer hold. To quantify the intravoxel diffusion heterogeneity many approaches have been proposed, including using signal representations with higher-order terms, such as diffusion kurtosis imaging (DKI) (Jensen et al., 2005), generalized diffusion tensor imaging (GDTI) (Özarslan and Mareci, 2003; Liu et al., 2004), mean apparent propagator (MAP) MRI (Özarslan et al., 2013; Avram et al., 2016), as well as multi-exponential, multi-tensor, or multi-compartment tissue diffusion models (Stanisz et al., 1997; Mulkern et al., 1999; Assaf and Basser, 2005; Zhang et al., 2012).

Jian et al., extended the multi-tensor signal representations to describe intravoxel diffusion heterogeneity using a Wishart distribution of microscopic diffusion tensors (Jian et al., 2007). Even though this parametric distribution is limited in its ability to accurately quantify the range of diffusion heterogeneity in healthy and diseased tissues, it nonetheless inspired great interest in measuring subvoxel distributions of microscopic diffusion tensors (DTDs). In general, however, to disentangle microscopic processes with arbitrary diffusivities, diffusion anisotropies, and orientations, it is necessary to sensitize the measurement to diffusion-diffusion correlations (Cory et al., 1990; Mitra, 1995; Callaghan and Komlosh, 2002) by preparing the signal with multiple pulsed-field gradients (mPFG), or multiple diffusion encodings (MDE). Historically, biological and clinical applications of mPFG or MDE methods (Komlosh et al., 2007) have focused on estimating microstructural parameters such as the average axon diameters (Koch and Finsterbusch, 2008; Avram et al., 2013a,b; Komlosh et al., 2018) or pore size distributions (Benjamini et al., 2016). More recently, MDE-prepared MRI measurements were described using tensor-valued diffusion encoding (Westin et al., 2016; Topgaard, 2017) in the context of probing subvoxel diffusion heterogeneity described using an ensemble of non-exchanging Gaussian diffusion tensor processes whose corresponding ellipsoids have distinct sizes, shapes, and orientations, i.e., the DTD.

While, at least in principle, one can reconstruct DTDs from a very large number of measurements with encodings sampling the 6D space of b-tensors, in practice, the limited signal-to-noise ratio (SNR) and long scan duration make such clinical or biological experiments very challenging (Topgaard, 2017; Song et al., 2022). To reduce the requirements for the high SNR level and a large number of measurements some have made simplifying assumptions such as cylindrical symmetry of microscopic tensors (Topgaard, 2017) which reduce the dimensionality of non-parametric DTD reconstructions from six to four degrees of freedom. Alternatively, one can use parametric models (e.g., analytical functions) to estimate features of the DTDs (Jian et al., 2007; Lasic et al., 2014; Szczepankiewicz et al., 2016; Westin et al., 2016; Magdoom et al., 2021) from data acquired using MDE and conventional single diffusion encoding (SDE) (Stejskal and Tanner, 1965).

Meanwhile, numerous studies using dMRI and other modalities provide converging evidence that, at a sufficiently small (i.e., mesoscopic) length scale, neuronal tissues, including cortical gray matter (GM) are organized preferentially along local orthogonal frames of reference. Ever since the earliest observations of cortical cyto- and myeloarchitecture (Brodmann, 1909; Cajal, 1909; Vogt, 1910), histochemistry and immunohistochemistry studies have consistently shown that cellular and subcellular structures at the microscopic scale are oriented predominantly along orthogonal, i.e., radial and tangential, orientations with respect to the cortical surface. This orthogonal reference frame persists at larger, mesoscopic scales of tens and hundreds of micrometers, and can be clearly seen in the arrangements of cells with various sizes, shapes, and densities forming tissue architectural patterns along the same radial and tangential orientations such as cortical columns and laminae, respectively (Amunts and Zilles, 2015; Rubenstein and Rakic, 2020). Most recently, studies using state-of-the-art electron microscopy (EM) in cortical GM (Lichtman and Denk, 2011; Shapson-Coe et al., 2021) have mapped the 3D organization of neuronal cells in gray matter with nanometer resolution over fields-of-view (FOVs) of hundreds of micrometers. These studies revealed in unprecedented detail anisotropic tissue structures, such as the microvasculature (Zhang et al., 2015), branching dendrites, neurofilaments, and other cell processes in various neuronal and non-neuronal cells (pyramidal neurons, intrinsic neurons, glial cells, etc.) roughly aligned along a local orthogonal frame of reference.

At mesoscopic length scales of a few 100 μm, diffusion processes in neural tissues align closely with the dominant orientations in the local tissue microstructure. Histological validation studies using ultra high-resolution dMRI have consistently found a good correspondence between the orientations of the underlying tissue microstructure and the orthogonal DTI reference frame (Budde and Annese, 2013; Seehaus et al., 2013, 2015) defined by ϵ1ϵ1 T,ϵ2ϵ2 T,ϵ3ϵ3 T, or the fiber orientation distribution functions (FOD) (Tournier et al., 2004) measured with high-angular resolution diffusion MRI (HARDI) (Tuch et al., 2002) in the brain (Leergaard et al., 2010). Numerous dMRI studies of cortical microstructure in fixed tissues (McNab et al., 2009, 2013; Dyrby et al., 2011; Miller et al., 2011; Kleinnijenhuis et al., 2013; Leuze et al., 2014; Aggarwal et al., 2015; Avram et al., 2022) and in vivo (Jaermann et al., 2008; Heidemann et al., 2010; McNab et al., 2013; Kleinnijenhuis et al., 2015; Gulban et al., 2018; Wang et al., 2021), for review see Assaf (2019), suggest that at submillimeter spatial resolution diffusion in the cortex is anisotropic and varies with the cortical folding geometry (McNab et al., 2013; Cottaar et al., 2018; Avram et al., 2022), in good agreement with the cortical cyto- and myeloarchitectonic features observed with histology and other modalities (Nieuwenhuys, 2013). Moreover, HARDI-derived FODs show preferentially radial and tangential components (Kleinnijenhuis et al., 2013; Leuze et al., 2014; Aggarwal et al., 2015) which evoke cortical columns (Petersen, 2007; Yacoub et al., 2008) and layers (Nagy et al., 2013; Bastiani et al., 2016), respectively, that can be observed with post-mortem histological staining. In addition, studies of laminar-specific intra-cortical connectivity measured with diffusion fiber microtractography (Leuze et al., 2014) of cortical FODs (Aggarwal et al., 2015; Gulban et al., 2018) suggest a similar orthogonal (radial and tangential) organization.

Increasing the spatial resolution in dMRI reduces the intravoxel angular dispersion of subvoxel diffusion processes and implicitly the orientational variance of the DTD. At submillimeter spatial resolution, dMRI is sensitive to cortical diffusion anisotropy and allows us to identify the radial and tangential orientations along which diffusion processes align. Recently, a careful survey of the high-resolution dMRI literature (Assaf, 2019) suggests that when different contrast preparations are used to vary the relative contributions of microscopic tissue water pools to net voxel dMRI signal in the cortex, the dominant diffusion orientations, as measured using the DTI eigenvectors or the directions of FOD peaks, remain unaffected even though the relative diffusivities or FOD amplitudes along these orientations may change. At mesoscopic spatial resolutions of a few 100 μm, the orientational characteristics of the dMRI signal remain remarkably consistent across experiments with fixed and live cortical tissues using different T1- and/or T2-weightings, i.e., different echo time (TE), repetition time (TR), or inversion time (TI), diffusion sensitizations (b-values) or diffusion/mixing times. These findings imply that at mesoscopic spatial resolutions, subvoxel cortical diffusion tensors from microscopic water pools are coincident along the same dominant (radial and tangential) orientations but may have potentially different diffusion anisotropies and diffusivities. Implicitly, at the mesoscopic length scale, the DTD is predominantly determined by the variations in the shapes (diffusion anisotropies) and sizes (diffusivities) of the microscopic diffusion tensors, rather than by their relative orientations.

In this study, we describe a new framework that simplifies the measurement and analysis of diffusion heterogeneity in microscopic water pools within gray matter using a non-parametric DTD. Specifically, if the voxel size is small enough compared to the radius of curvature of the cortex, we can constrain all the microscopic (subvoxel) diffusion tensors to share the same principal reference frame determined, for instance, by the dyadic of the principal diffusion eigenvectors, ϵ1, ϵ2, ϵ3, measured with DTI. With this constraint, the DTD is completely characterized by the voxel reference frame ϵ1ϵ1 T,ϵ2ϵ2 T,ϵ3ϵ3 T, and by the 3D joint distribution of corresponding subvoxel principal diffusivities, λ1, λ2, λ3, which are random variables. This joint probability distribution can be estimated with a 3D Inverse Laplace Transform analysis using only single diffusion encoded (SDE) MR measurements. This practical, non-parametric framework for mapping DTDs, called COnstrained Reference frame diffusion TEnsor Spectroscopic (CORTECS) MRI, has the potential to quantify a wide range of cortical diffusion heterogeneity in healthy or diseased brains.

2. Methods 2.1. Higher spatial resolution reduces the intravoxel orientational dispersion of dMRI signals

The net diffusion signal in an imaging voxel containing complex tissue microstructure can be described generally using an ensemble of subvoxel (i.e., microscopic) diffusion tensors with different sizes, shapes, and orientations assumed to be in slow exchange, i.e., the diffusion tensor distribution (DTD). Ordinarily, we can quantify DTDs by analyzing diffusion-weighted images (DWIs) acquired with multidimensional diffusion encoding (MDE) (Westin et al., 2016; Topgaard, 2017; Magdoom et al., 2021). The net dMRI voxel signal, S, is a function of the tensor-valued encoding variable called the b-tensor, b, computed by integrating the time-dependent diffusion gradient waveforms amplitudes, and is related to the underlying DTD, p(D):

S(b)=∫M+p(D)e-b·DdD,    (2)

where the integral runs over the space or domain of all positive definite matrices, M+. Since the random variable D has 6 degrees of freedom, p(D) is essentially a 6-dimensional joint probability distribution (or correlation spectrum) of the diffusion tensor elements. The high dimensionality and the inherent challenge of defining the subspace of positive-definite random tensor-valued variables, D, make solving this problem infeasible in practice, as no closed-form solution exists. Measuring p(D) requires a prohibitively large number of measurements with a very high signal-to-noise ratio (SNR) and MDE. Previously, approximations to p(D) have been proposed either by assuming parametric models and/or by using statistical reconstruction algorithms (Jian et al., 2007; Westin et al., 2016; Topgaard, 2017; Magdoom et al., 2021).

In cortical GM the orthogonal coordinate axes along which diffusive fluxes align at the microscopic scale of cellular and subcellular structures (i.e., diffusion length scale) are propagated at larger mesoscopic scales guiding the assembly of these structures into orthogonal tissue architectural patterns of cortical laminae and columns (Nieuwenhuys, 2013; Rubenstein and Rakic, 2020). If the voxel size of dMRI data is significantly smaller than the minimum radius of the curvature of the underlying anatomy (i.e., cortical folding) the orientational variance of subvoxel (microscopic) diffusion processes can be neglected (Figure 1). Microscopic diffusion processes are coincident with the axes of the local microstructural reference frame determined by the cortical cyto- and myeloarchitecture. For a continuously varying cortical anatomy with a minimum radius of curvature, R, the range of orientational misalignment between the microscopic diffusion tensors and the voxel reference frame, ±θmax, in a cubic voxel of side length, x, is:

θmax=tan-1(x32R)    (3)

Figure 1B shows that θmax decreases rapidly at low spatial resolutions, Rx, but changes slowly at higher values of Rx (Figure 1B). At a spatial resolution of a few hundred micrometers the voxel size is much smaller than the cortical radius of curvature (R = 5 mm) leading to very small values of θmax. Under these circumstances, it is reasonable and practical to constrain all diffusion tensor processes in microscopic water pools throughout the voxel (i.e., the DTD) to be described using the same local orthogonal reference frame (Figures 2A,B).

Figure 1. (A) As we decrease the voxel size, x, relative to the radius of curvature of the tissue (e.g., due to cortical folding), R, the intravoxel orientational variance of the continuously varying microstructural reference frame also decreases. For a voxel with an arbitrary orientation relative to the underlying microstructure, the range of intravoxel orientational variation due to tissue curvature is ±θmax. (B) The value of θmax decreases rapidly at low spatial resolutions, R/x, but changes very slowly at higher spatial resolutions, R/x. (C) A quantitative comparison of θmax at different voxel sizes assuming a cortical radius of curvature R = 5mm shows the significant reduction in intravoxel orientational variance due to the effects of anatomical curvature at high spatial resolutions.

Figure 2. At a mesoscopic length scale cortical cyto- and myeloarchitecture is organized preferentially along the axes of an orthogonal frame of reference (A). If the dMRI spatial resolution is sufficiently small (Figure 1) we can measure DTDs efficiently using the constraints of the CORTECS MRI framework (B). If we constrain all microscopic diffusion tensors to have the same principal axes of diffusion (C) we can quantify the DTD as the 3D correlation spectrum of the corresponding principal diffusivities (D). If the microarchitecture varies along a single radial orientation, we can further constrain the DTD to contain only axisymmetric tensors (F) and quantify the 2D correlation spectrum of the corresponding radial and tangential diffusivities (G). We can also quantify the shape-size (i.e., microscopic FA-MD) correlation spectra of microscopic tensors from the 3D (E) or 2D (H) constrained reference frame DTDs (cDTDs).

2.2. COnstrained Reference frame diffusion TEnsor Correlation Spectroscopic (CORTECS) MRI

Fixing the local reference frame for all subvoxel tensors has several surprising advantages. First, it significantly reduces the dimensionality of p(D) and decouples the statistical random variables needed to describe p(D). Specifically, the 6D vector/tensor random variable, D, corresponding to the 6 components (or degrees of freedom) needed to describe the general DTD is reduced to a 3D random variable comprising the three principal diffusivities, λ1, λ2, λ3 along the axes of the fixed voxel frame of reference, ϵ1ϵ1 T,ϵ2ϵ2 T,ϵ3ϵ3T, respectively, which are sufficient to describe the constrained DTDs (cDTDs) within the CORTECS MRI framework (Figures 2C–E). Using the eigenvalue decomposition of the diffusion tensor (Equation 1) we can re-write Equation (2) as a more tractable 3D Inverse Laplace transform (ILT) problem:

S(b)=∫0∞∫0∞∫0∞ p(λ1,λ2,λ3)e −λ1ϵ1 T bϵ1−λ2ϵ2 Tbϵ2−λ3ϵ3 Tbϵ3dλ1dλ2dλ3,    (4)

where ϵiTbϵi is a non-negative scalar weighting that represents the reciprocal Laplace variable corresponding to λi, and b is the measurement b-tensor (Westin et al., 2016). We can estimate non-parametric DTDs in the reduced dimensional space of the principal diffusivities by applying the ILT to dMRI data acquired with SDE and/or MDE. Besides the drastic reduction in the computational complexity due to the dimensionality reduction, the CORTECS framework inherently enforces positive definiteness of diffusion tensors by requiring positivity for λi.

Another very important advantage of constraining the reference frames of the DTD tensor random variable is that we can measure p(λ1, λ2, λ3) using only DWIs acquired with single pulse-field gradient (sPFG), or SDE, a.k.a. linear tensor encoding with rank-1 b-tensors. For a conventional SDE DWI with an arbitrary b-value, b, and diffusion gradient direction given by the unit vector g=[gx,gy,gz]T, the encoding b-tensor has rank-1, b = bggT. We can rewrite the signal equation above with respect to the components of g expressed in the voxel frame of reference, g′=[g′1,g′2,g′3]T=[ϵ1ϵ2ϵ3]g:

S(b)=∫0∞∫0∞∫0∞ p(λ1,λ2,λ3)e-λ1bg′12e-λ2bg′22e-λ3bg′32dλ1dλ2dλ3    (5)

The factors bg′i2 are the non-negative weighting parameters of the principal diffusivities, λi, in the Laplace Transform representation of the signal. We can generate a wide range of joint weighting parameters bg′i2 by varying the b-value and diffusion gradient orientations in conventional SDE preparations. Subsequently, from multiple SDE DWIs we can estimate, in each voxel, the correlation spectrum of principal diffusivities, p(λ1, λ2, λ3) which quantifies the properties of all microscopic diffusion tensors. Compared to MDE-DWIs, the conventional SDE-DWI can be acquired efficiently using product single pulsed-field gradient (sPFG) spin-echo (SE) diffusion MR sequences (Stejskal and Tanner, 1965) available on all microimaging and clinical MRI scanners. In general, SDE-DWIs can achieve higher b-values, shorter echo times (TEs), higher spatial resolution, and/or better SNR than MDE-DWIs using double or triple diffusion encoding (Sjölund et al., 2015). Moreover, the spectral reconstruction of p(λ1, λ2, λ3), henceforth referred to as 3D cDTD, does not require statistical methods to enforce positive definiteness but can still benefit from various techniques that may be used to solve ILT-like problems, such as L2- or L1-norm regularization, compressed sensing (Bai et al., 2015), or constrained optimization (Benjamini et al., 2016), etc.

The microarchitecture of certain tissues can be described more economically and effectively using a single dominant axis. For example, in cortical gray matter, microscopic diffusion processes align with cell and tissue structures (Budde and Annese, 2013; Seehaus et al., 2013, 2015) that are preferentially oriented radially or tangentially with respect to the cortical surface (Kleinnijenhuis et al., 2013; McNab et al., 2013; Leuze et al., 2014). Under these circumstances, we can further simplify the problem and assume that the DTD comprises only tensors with cylindrical symmetry (Figures 2F–H). Thus, the voxel reference frame is determined by a single orientation, ϵ1ϵ1T, i.e., normal to the cortical surface, which implicitly determines the orthogonal, tangential component described by the rank-2 tensor ϵ2ϵ2T+ϵ3ϵ3T=I3-ϵ1ϵ1T, where I3 is the 3x3 identity matrix. We can relate the signal in a voxel with fixed principal axis ϵ1ϵ1T to a two-dimensional correlation spectrum of principal diffusivities along radial (cortical columns) and tangential (cortical layers) orientations with respect to the cortical surface, p(λr, λt):

S(b)=∫0∞∫0∞p(λr,λt)e-λrbcos2ϕge-λtbsin2ϕgdλrdλt    (6)

The parameter ϕg=arccos(ϵ1Tg) represents the angle between the applied gradient direction, g, and the radial direction of the underlying reference frame, ϵ1ϵ1T, while p(λr, λt) completely determines the corresponding distribution of cylindrically symmetric diffusion tensors, henceforth referred to as the 2D cDTD.

Lastly, in a final simplifying step, if all subvoxel diffusion processes are isotropic, the correlation spectrum of diffusion tensor eigenvalues reduces to a distribution of a single scalar diffusivity random variable, λ0, which can be viewed as a 1D cDTD:

S(b)=∫0∞p(λ0)e-λ0bdλ0    (7)

As an aside, we should point out an important connection between 1D cDTD MRI and our previously proposed methods for one- and multidimensional MD spectroscopic MRI using isotropic diffusion encoding (IDE) (Avram et al., 2019, 2021). Mapping non-parametric spectra of MD values in microscopic tissue water pools using multiple IDE measurements does not require that diffusion in these pools is isotropic. Meanwhile, the 1D cDTD MRI spectral reconstruction using Equation (7) correctly quantifies the spectra of water mobilities only if all diffusion processes within the voxel are isotropic, in which case the two methods will provide congruent results.

2.3. Mapping distributions and correlation spectra of microscopic fractional anisotropy and mean diffusivity

From the measured cDTD within each voxel, we can compute non-parametric distributions and correlation spectra of any DTI parameters derived from the microscopic diffusion tensors, such as fractional anisotropy (FA) or mean diffusivity (MD). Specifically, we can define a new random variable, α, that quantifies the FA of each microscopic diffusion tensor in the cDTD:

α=12(λ1-λ2)2+(λ2-λ3)2+(λ3-λ1)2λ12+λ22+λ32    (8)

From p(λ1, λ2, λ3) we can then derive the probability density function (one-dimensional spectrum) of the microscopic tensor FAs, pFA(α), which quantifies the cDTD shape heterogeneity non-parametrically. The statistical moments of pFA(α) provide important microstructural parameters, such as the microscopic anisotropy, μFA, computed as the mean of pFA(α) (Lasic et al., 2014; Westin et al., 2016; Magdoom et al., 2021). Similarly, we can define another cDTD-derived random variable that quantifies the mean diffusivity of each microscopic tensor, μ = (λ1+λ2+λ3)/3, and compute its probability density function pMD(μ) to describe the spectrum the microscopic tissue water mobilities non-parametrically.

Finally, from p(λ1, λ2, λ3) we can also compute non-parametric multidimensional correlation spectra of two or more microscopic DTI metrics. For example, we can quantify non-parametrically the correlations between the shapes and sizes of the diffusion ellipsoids corresponding to the microscopic diffusion tensors by computing the joint probability density function of the two random variables α and μ, pFA−MD(α, μ). This practical and efficient decomposition of tissue heterogeneity based on diffusion anisotropy and mean diffusivity correlations in microscopic water pools may reveal specific microstructural motifs or patterns potentially relevant to many clinical applications.

2.4. A generalization of several multi-tensor diffusion signal models

The CORTECS framework can describe a wide range of heterogeneous diffusion processes in healthy and diseased tissues and subsumes several diffusion tensor signal models. For example, if we constrain p(λ1,λ2,λ3)=δ(λ1-λ′1,λ2-λ′2,λ3-λ′3), 3D cDTD simplifies to conventional DTI with the three mean eigenvalues λ′1,λ′2,λ′3. In this way, 3D cDTD can be viewed as a generalization of high-resolution DTI that quantifies intravoxel diffusion heterogeneity as a non-parametric correlation spectrum of the principal diffusivities in microscopic water pools. To describe multi-exponential or multi-tensor signal decays in heterogeneous tissues (Stanisz et al., 1997; Mulkern et al., 1999; Avram et al., 2020) we can assume that p(λ1, λ2, λ3) can be represented as a sum of delta functions (point masses) (Avram et al., 2020). Moreover, the spectroscopic decomposition of the net voxel signal in cDTD makes it easy to disentangle partial volume contributions, such as those from cerebrospinal fluid (CSF), or free water in tissues caused by edema or other processes (Pasternak et al., 2009).

2.5. Monte Carlo simulations

We conducted Monte Carlo (MC) simulations to evaluate the numerical stability and accuracy of the voxel-wise estimation of 3D and 2D cDTDs from noisy data. Specifically, starting from ground truth DTDs constrained with fixed voxel reference frames (2D and 3D cDTDs), defined analytically using multidimensional log-normal distributions, respectively, we computed the dMRI signals expected from an experiment using conventional single-diffusion encoded (SDE) DWI measurements with the same gradient orientations and b-values as in our fixed-brain experiment described below. Next, from these ground truth signals, we generated 500 instances of noisy measurements by adding Rician noise to simulate real measurements with different SNR levels. From each set of noisy measurements, we computed the corresponding normalized 3D correlation spectra of principal diffusivities, or normalized 2D correlation spectra of radial and tangential diffusivities and compared the statistics of these spectra (mean and standard deviation) to the corresponding ground truth 3D and 2D DTDs, respectively.

2.6. Ultra high-resolution dMRI of a fixed macaque monkey brain

The brain of a healthy young adult rhesus macaque monkey (Macaca mulatta) weighing 13.55 kg was prepared using a well-controlled perfusion fixation process, as described in Saleem et al. (2021). In brief, the animal was deeply anesthetized with sodium pentobarbital and perfused transcardially with heparinized saline, followed by 4% paraformaldehyde in 0.1 M phosphate buffer (pH 7.4). After perfusion, the brain was removed from the cranium and post-fixed for 8 h in the same buffered paraformaldehyde solution. Following the post-fixation, the brain was transferred into 0.1 M phosphate-buffered saline with sodium azide before the MRI data acquisition. All procedures were carried out under a protocol approved by the Institutional Animal Care and Use Committee of the National Institute of Mental Health (NIMH) and the National Institute of Health (NIH) and adhered to the Guide for the Care and Use of Laboratory Animals (National Research Council).

Based on a preliminary structural MRI scan of the specimen, we fabricated a three-dimensional (3D) brain mold in