Evidence for microscopic kurtosis in neural tissue revealed by correlation tensor MRI.

PURPOSE: The impact of microscopic diffusional kurtosis (µK), arising from restricted diffusion and/or structural disorder, remains a controversial issue in contemporary diffusion MRI (dMRI). Recently, correlation tensor imaging (CTI) was introduced to disentangle the sources contributing to diffusional kurtosis, without relying on a-priori multi-gaussian component (MGC) or other microstructural assumptions. Here, we investigated µK in in vivo rat brains and assessed its impact on state-of-the-art methods ignoring µK. THEORY AND METHODS: CTI harnesses double diffusion encoding (DDE) experiments, which were here improved for speed and minimal bias using four different sets of acquisition parameters. The robustness of the improved CTI protocol was assessed via simulations. In vivo CTI acquisitions were performed in healthy rat brains using a 9.4T pre-clinical scanner equipped with a cryogenic coil, and targeted the estimation of µK, anisotropic kurtosis, and isotropic kurtosis.
RESULTS: The improved CTI acquisition scheme substantially reduces scan time and importantly, also minimizes higher-order-term biases, thus enabling robust µK estimation, alongside Kaniso and Kiso metrics. Our CTI experiments revealed positive µK both in white and gray matter of the rat brain in vivo; µK is the dominant kurtosis source in healthy gray matter tissue. The non-negligible µK substantially were found to bias prior MGC analyses of Kiso and Kaniso .
CONCLUSIONS: Correlation Tensor MRI offers a more accurate and robust characterization of kurtosis sources than its predecessors. µK is non-negligible in vivo in healthy white and gray matter tissues and could be an important biomarker for future studies. Our findings thus have both theoretical and practical implications for future dMRI research.

INTRODUCTION

Diffusion MRI (dMRI) has become one of the most important methods for non‐invasively probing microstructural features in health and in disease. , , Single diffusion encoding (SDE) experiments , (Figure 1A), probing diffusion along a single axis, have been widely used to measure the directional apparent diffusion coefficient, , , the diffusion propagator from q‐space experiments, , , , the diffusion tensor , , and its time‐dependence, , , and the diffusional kurtosis (via diffusional kurtosis imaging, DKI ). Together with microstructural and biophysical models, , , , , SDE has provided important insight into microstructural tissue.

FIGURE 1

Experiments for kurtosis source estimation. A, Parameters of a SDE pulse sequence, where Δ is the diffusion time, is the gradient pulse duration, and g is the gradient pulse amplitude (alternatively, it is sometimes convenient to describe the diffusion encoding by a gradient direction n and b‐value ). B, Parameters of a DDE pulse sequence, where Δ1 and Δ2 are the diffusion encoding blocks’ diffusion times, and are their gradient pulse durations, g 1 and g 2 are their gradient amplitudes and is the mixing time (ie, the time between the two diffusion encoding bocks). At the long mixing time regime, DDE experiments for powder‐averaged systems can be fully described by two b‐values ( and ) and the angle describing the relative orientation between the directions of the two diffusion encoding modules n 1 (blue) and n 2 (red). C, Acquisition strategy for resolving the different kurtosis sources via CTI. The diffusion encoding profiles for all dastasets are shown in panels C1.1, C2.1, C3.1, and C4.1. The gradient directions used for signal powder‐average calculation are shown in panels C1.2, C2.2, C3.2, and C4.2. The equivalent b‐tensor shapes for each experiment are shown in panels C1.3, C2.3, C3.3, C4.3. Abbreviation: CTI, correlation tensor imaging

DKI in particular provides important information on the degree of non‐Gaussian diffusion by quantifying the diffusional excess‐kurtosis. DKI has been shown to be very sensitive to, inter‐alia, age‐related microstructural changes, , , , , ischemia, , tumors, , , traumatic brain injury, , and Parkinson and Alzheimer diseases. , However, the biological interpretation of DKI remains limited because non‐Gaussian diffusion may arise from different sources. , , Several compartment models have been proposed to relate non‐Gaussian diffusion with more specific biological underpinnings , , , , ; however, the specificity of these models can be severely compromised by their strong assumptions and constraints. , , , ,

Unique information about the non‐Gaussian nature of diffusion can be resolved using multidimensional diffusion encoding (MDE) strategies. , , In contrast to SDE methods, MDE probes diffusion correlations across different dimensions by either including additional pairs of pulse gradients , , , , , , , , or using continuous gradient waveforms with 3D trajectories. , , , , , , Under the strict multiple Gaussian components (MGC) assumption (no time‐dependence and no kurtosis arising from restricted diffusion or structural disorder), Westin et al. showed that diffusion encoding in MDE can be generally described by tensor‐valued information. Using different b‐tensor shapes, MDE can resolve anisotropic and isotropic kurtosis sources ( and ), which represent the shape and size variances the diffusion tensors comprising the entire system. , , , Such MGC analyses are, however, doubly prone to bias arising from restricted diffusion (eg, upon interaction with microscopic boundaries): First, the MGC analyses of MDE data biases and when continuous gradient waveforms are long compared to the time it takes to probe the boundaries. , , Second, as will be shown in this work, even for fixed diffusion times and tuned MDE sequences, MGC analyses can be biased by microscopic non‐Gaussian effects (namely, microscopic kurtosis [µK]) that arise from restricted diffusion or systems characterised by complex microstructural features with diffusion path lengths in the order of the probed scales, , , , , eg, intra‐cellular cross‐sectional size variance, extra‐cellular tortuosity, etc.

Double diffusion encoding (DDE , , , , , , , , , ) is an MDE variant probing diffusion via two pairs of diffusion pulsed gradients (Figure 1B). Although previously used without MGC analyses to resolve microscopic anisotropy, , , , , DDE’s inherent capability to provide both linear and planar encodings was recently harnessed to resolve and based on MGC analyses. Moreover, DDE at the long mixing time regime can be used to minimize diffusion time‐dependence effects. , , , , , Although prior DDE studies have attempted to measure µK, orientation dispersion was inherently conflated in these studies, , making the importance of microscopic kurtosis in tissues unclear. Going beyond the MGC framework, the correlation tensor imaging (CTI) approach was recently introduced for µK measurements. The CTI framework allows the simultaneous decoupling of , from µK effects without resorting to multi‐Gaussian component assumptions and without conflation with other mesoscopic effects (such as orientation dispersion). However, the initial CTI approach can suffer from higher‐order effects, and was quite time‐consuming to acquire, thereby limiting its in vivo applicability.

Here, we aimed to investigate the existence of µK in in vivo neural tissues and its impact on the increasingly popular MGC approaches. We first develop a highly improved and accelerated CTI acquisition scheme, which is more robust towards µK estimation, with minimized high‐order‐term biases. Our data reveal that µK cannot be ignored in in vivo neural tissues, and that the highly popular MGC approaches provide biased information due to µK. Our results provide a new window for quantifying microstructure in health and disease, and show that µK must be considered in future dMRI studies.

THEORY

Total kurtosis estimates from single diffusion encoding

The SDE signal attenuation () can be expressed (with Einstein summation convention) as the following 2nd order cumulant expansion : where b is the b‐value defined by , is the diffusion gradient direction, D and W are the diffusion and excess‐kurtosis tensors, and is the mean diffusivity.

To quantify non‐Gaussian diffusion decoupled from confounding effects of tissue dispersion, it is also useful to consider the cumulant expansion of powder‐averaged SDE signal decays (ie, signals averaged across multiple gradient directions) , , : where is the powder‐averaged SDE signal decay, D and K are the isotropic diffusivity () and isotropic excess‐kurtosis of powder‐averaged signals . In the absence of exchange, the total kurtosis K can be described by the sum of three different sources : where K is related to system’s miscrocospic anisotropy µA (), , , , and K is related to the variance of components’ apparent mean diffusivities D (, with representing the variance across the mean diffusivities of components). , The µK, which was previously referred to as intra‐compartmental kurtosis, is a weighted sum of different microscopic sources of non‐Gaussian diffusion µK with representing the average over tissue components. Here µK can be related to non‐Gaussian diffusion arising from restricted diffusion , , or tissue disorder due to the presence of microscopic hindrances to water molecules, eg, membranes, organelles, axonal caliber variations, etc. , , , , Although the total kurtosis K can be estimated by fitting Equation (2) to data acquired with at least two non‐zero b‐values, it is important to note that the kurtosis sources in Equation (3) cannot be decoupled from SDE experiments in a model free manner.

Correlation tensor imaging kurtosis source estimation

Recently, the CTI methodology was proposed to resolve different kurtosis sources from DDE signals. Figure 1B shows an illustration of the DDE sequence which probes diffusion using two pairs of pulsed gradients with magnitudes g 1 and g 2, widths and , separations time Δ1 and Δ2, and mixing time (Figure 1B). Note that the DDE pairs can also be applied along different directions, 1 and 2. To probe kurtosis for fixed timing parameters, CTI uses and Δ1 = Δ2 = Δ. Moreover, to avoid diffusion time‐dependent biases, CTI is applied to DDE data acquired at long mixing time. In this regime and up to 2nd order in b, the DDE signal attenuation () can be expressed as , , , : where and are the b‐values associated with the two DDE gradient wavevectors, and Z is a tensor that approaches the covariance tensor () at long mixing times. We previously showed that K, K, and µK can in theory be extracted from the tensors of Equation (5). However, our preliminary validation showed that high‐order‐terms can introduce biases on the different kurtosis estimates that depend on dispersion levels. To suppress this dependence, powder‐averaged DDE signals () are used: where is defined as the angle between the gradient directions 1 and 2. Note that several different pairs of 1 and 2 with constant evenly sampling a 3D unit sphere are required for powder‐averaging , , —specific pairs of gradient directions used on this study are shown below. From the parameters of Equation (6), µK can be estimated by (c.f. Equation 3).

Accelerating CTI and increasing its robustness towards higher order effects

To accelerate the CTI acquisition, we find that only the four different sets of DDE experiments illustrated in Figure 1C (in addition to acquisitions without diffusion sensitization, ie, ), are required to extract CTI’s metrics. The four sets are as follows:

Powder‐averaged signals with and . Note that these experiments are equivalent to SDE experiments (Figure 1C1);

Powder‐averaged symmetric DDE with diffusion weighting and parallel gradient directions (, Figure 1C2).

Powder‐averaged symmetric DDE with diffusion weighting and perpendicular gradient directions (, Figure 1C3);

Powder‐averaged symmetric and parallel DDE as 3) but with a different total b‐value (Figure 1C4). Note that all previous sets (1‐3) have the same total b‐value equal to

To ensure homoscedastic signals, all four experiment sets should be acquired with an equal number of gradient direction pairs for powder‐averaging.

Specificity of the improved protocol to different kurtosis sources

The analysis to resolve different kurtosis sources proceeds as follows:

µK can be extracted from the log difference of powder‐averaged signals from the experiments’ set 1 and 2:

as pointed in previous studies (eg, Refs. [ , , ]), K can be extracted from the log difference of powder‐averaged signal from sets 2 and 3:

to decouple D̅, K, and K, powder‐averaged signals also require at least two non‐zero total b‐values b. Therefore, DDE experiments with symmetric intensities and parallel directions for a lower total b‐value b = b are acquired (set 4).

Diffusion tensor variance approach

Under the MGC assumption, previous studies showed that K and K can be estimated from signals measured using any MDE sequence that probes different b‐tensor magnitudes b and shapes. For the sake of simplicity, here we only consider axial tensor‐valued experiments, , , where the b‐tensor shape is characterized by a single parameter . Thus, the powder‐averaged signal is

Note that and can be estimated from Equation (8), using the same dMRI experiments as for CTI, since these correspond to data acquired with at least two b‐tensor shapes ( = 1 for sets 1, 2, and 4, and = ‐1/2 for set 3) and two non‐zero b‐tensor magnitudes ( for sets 1‐3, for set 4).

METHODS

Simulations

The robustness of the different kurtosis source estimation strategies (CTI and MGC) was first assessed via simulations. The full details of the simulations can be found in the Supporting Information, Section A, which is available online. Briefly, synthetic signals were generated for different models with known ground truth kurtosis sources, which included:

Sum of multiple isotropic Gaussian diffusion components (Figure 2A).

FIGURE 2

Simulations for synthetic diffusion‐weighted signals from three different systems containing only a single microenvironment “classes”. A, Isotropic Gaussian diffusion components with different mean diffusivities sampled from a Gaussian distribution with mean 0.65 µm2/ms and standard deviation 0.21 µm2/ms. B, Isotropically oriented anisotropic Gaussian diffusion components with axial and radial diffusivities of 1.45 µm2/ms and 0.25 µm2/ms. C, Non‐Gaussian diffusion due to microscopic disorder with ground truth and D = 0.65 µm2/ms. From left to right, each panel shows: a schematic representation of the models (A1, B1, C1, D1); the signal decays for three different DDE experiment types in which signals of the improved CTI protocol are marked by the blue triangles (A2, B2, C2); the kurtosis ground truth values (A3, B3, C3); the kurtosis estimates obtained from MGC (A4, B4, C4); the kurtosis estimates obtained from the CTI using its previous “old” protocol (A5, B5, C5); the kurtosis estimates obtained from the CTI using its improved “new” protocol (A6, B6, C6). Note that when µK exists (panel C), the MGC approaches conflate µK with non‐existing anisotropic and isotropic sources, while CTI clearly identifies the µK in the system. Also note that the new CTI protocol better resolves anisotropic and isotropic sources compared to its old counterpart (panels A and B)

Sum of multiple uniformly oriented anisotropic Gaussian diffusion components with identical axial and radial diffusivities (Figure 2B).

Single compartment with non‐vanishing µK (Figure 2C). Although in Figure 2C1 a single component with positive µK is sketched as the extra‐compartmental medium that encompasses randomly oriented anisotropic compartments, ground truth positive µK in both intra‐ and extra‐“cellular” components can arise according to effective medium theory due to, eg, cross sectional size variance and packing degree. , , , The exact µK value for a medium represented in Figure 2C will depend on the volume fraction, anisotropy, size, and packing of the anisotropic compartments as well as on the acquisition parameters. To simplify, the DDE signal decay for the single isotropic compartment with positive µK is here numerically computed using the signal representation with and µK ground truths set to an arbitrary value of 0.65 µm2/ms and 1, respectively. Note that µK can also arise due to restricted diffusion. In section B of Supporting Information, simulations for restricted diffusion inside spheres are also produced using the MISST package, , which could represent, for example, neural soma and other quasi‐spherical objects such as boutons.

A system comprising different components and with non‐zero contributions for all different kurtosis sources (Figure 3A). For this system, we consider a sum of all compartment types used in the previous simulations with equal weights. As the mean diffusivities of the simulations 1, 2 and 3 are equal, this ensemble model can assess the robustness of estimates for different kurtosis sources individually by varying concrete model parameters. We thus varied ground truth , ground truth , and ground truth µK values. Please see the Supporting Information, Section A, for full details on these simulations.

FIGURE 3

Simulations of synthetic diffusion‐weighted signals generated in a complex system containing a sum of different microenvironment “classes” (isotropic Gaussian + anisotropic Gaussian + isotropic non‐Gaussian components). A, The microstructural scenario (isotropic component mean diffusivities sampled from a Gaussian distribution with mean 0.65 µm2/ms and standard deviation 0.21 µm2/ms, anisotropic Gaussian diffusion components with axial and radial diffusivities of 0.45 µm2/ms and 0.25 µm2/ms; and non‐Gaussian diffusion with K = 1 and D = 0.65 µm2/ms) (A1), signal decays (A2), and results from the model fit (A3‐A6 for the various sources), showing that CTI describes the ground truth much more accurately than its MGC counterpart. B, CTI and MGC kurtosis source estimates and their sensitivity to varying ground truth K (by changing the mean diffusivity variation of isotropic Gaussian components). C, CTI and MGC kurtosis estimates sensitivity to varying ground truth K by changing the axial and radial diffusivities of anisotropic Gaussian compartments. D, sensitivity of CTI and MGC estimates to varying ground truth µK by changing the µK levels in the non‐Gaussian component. These simulations suggest a very good stability for CTI (especially the new protocol), such that a change in one direction of any of the metrics will induce a systematically correct change in the estimated values

For all models, powder‐averaged signals were generated for the four different sets of DDE acquisition parameters (c.f. Figure 1) for total b‐values b = 2.5 ms/µm2 (sets 1, 2, 3) and b = 1 ms/µm2 (set 4) ( and for all experiments)—note the maximum b‐value of 2.5 ms/µm2 was selected since this was showed to provide an optimal trend between signal contrast to measure diffusional kurtosis and minimization of biases from high order‐terms. For all four sets, the 45 directions of a 3D spherical eight‐design were used for the single encoding of set 1 (Figure 1, panel C1.2), for the double diffusion encodings of sets 2 and 4 (Figure 1, panels C2.2 and C4.2), and for the first diffusion encoding of set 3 (Figure 1, panels C3.2). The directions for the second diffusion encoding of set 3 were repeated for three equidistant orthogonal directions relative to each direction of the 3D spherical eight‐design (Figure 1, panels C3.2), yielding a total of 135 pairs of directions. To ensure homoscedasticity of powder‐averaged signals, the acquisition of the 45 directions for sets 1, 2, and 4 were repeated three times (45 × 3 = 135 pairs of directions).

For reference, signals were also produced for an adapted version of the previous (old) DDE protocol suggested for CTI. The , , and values, as well as gradient directions schemes used for the improved (new) and old CTI protocols (CTInew and CTIold) are summarized in Table 1 (for more information on the old protocol, c.f. Supporting Information, Section C). In addition to the diffusion‐weighted signals for the different CTI sets, 135 signal replicas for were incorporated in both protocols; these data are treated as an independent set. To assess the robustness of estimates towards noise, all synthetic signals were corrupted by Rician noise with a nominal signal‐to‐noise ratio (SNR) of 40 before powder‐averaging. Consideration about the precision of CTI for other SNR levels are discussed in Supporting Information, section D (c.f. Supporting Information Figure S3).

TABLE 1

Summary of the DDE parameter combination used for the “new” CTI protocol (CTInew) and the reference “old” CTI protocol (CTIold). Parameters b1, b2, bt are expressed in ms/µm2

Set	b1	b2	bt	θ	bΔ	Direction scheme
New CTI protocol—CTInew
1	2.5	0	2.5	0o	1	45 directions of the 8‐design (x3 repetitions)
2	1.25	1.25	2.5	0o	1	45 directions of the 8‐design (x3 repetitions) for both diffusion encodings
3	1.25	1.25	2.5	90o	‐1/2	45 directions of the 8‐design for the 1st encoding, repeated for 3 orthogonal directions for the 2nd encoding
4	0.5	0.5	1	0o	1	45 directions of spherical 8‐design (x3 repetitions) for both diffusion encodings
Old CTI protocol—CTIold
1	2.5	0	2.5	0o	1	45 directions of the 8‐design (x3 repetitions)
2	2.5	2.5	5	0o	1	45 directions of the 8‐design (x3 repetitions) for both diffusion encodings
3	2.5	2.5	5	90o	‐1/2	45 directions of the 8‐design for the 1st encoding, repeated for 3 orthogonal directions for the 2nd encoding
4	1	0	1	0o	1	45 directions of the 8‐design (x3 repetitions)
5	1	1	2	0o	1	45 directions of the 8‐design (x3 repetitions) for both diffusion encodings
6	1	1	2	90o	‐1/2	45 directions of the 8‐design for the 1st encoding, repeated for 3 orthogonal directions for the 2nd encoding

Data processing

CTI metrics were obtained by fitting Equation (6) to the log of the powder‐averaged signals of both new and old protocols using an ordinary linear‐least‐squares (OLLS) procedure. , and were also estimated from the MGC approach by fitting Equation (9) to the powder‐average of CTI’s new protocol using an OLLS procedure. Mean and standard deviation of each kurtosis estimates were computed by repeating simulations for 1000 iterations.

MRI experiments

All animal experiments were preapproved by the institutional and national authorities and carried out according to European Directive 2010/63. Data were acquired from N = 3 female Long Evans rats (ages = 22/19/22 weeks old, weights = 354.6/260.4/334.5g in a 12 h/12 h light/dark cycle with ad libitum access to food and water) under anesthesia (~2.5% isoflurane in 28% oxygen) using a 9.4T Bruker Biospec scanner equipped with an 86 mm quadrature transmission coil and four‐element array reception cryocoil (Bruker, Fallanden, Switzerland).

Auxiliary sagittal T2‐weighted images were acquired using a RARE sequence (c.f. Supporting Information Section E for experimental parameters). Diffusion MRI datasets were then acquired using a DDE‐EPI pulse sequence (Δ = = 12 ms, δ = 3.5 ms) for three evenly spaced coronal slices placed using the sagittal T2‐weighted images as a reference (c.f. Supporting Information Figure S4). Per‐slice respiratory gating was applied. These acquisitions followed the first CTI protocol reported in Table 1 along with 24 acquisitions per DDE set. Other acquisition parameters included: repetition time/echo time (TR/TE) = 3000/50.9 ms, in‐plane resolution = 200 × 200 μm2, slice thickness = 0.9 mm acquisition bandwidth = 400 kHz, number of shots = 1, partial Fourier factor = 7/8 (partial Fourier acceleration = 1.25). Acquisition time for the diffusion‐weighted data was around 40 mins. To empirically check if our measurements satisfied the long mixing time regime, additional DDE data were acquired in rat 1 with parallel and antiparallel directions and an intermediate total b‐value of 2 ms/µm2 (b 1 = b 2 = 1 ms/µm2) (Supporting Information, Section F).

C.f. Supporting Information Section G for full details. Briefly, data were denoised and Gibbs unrung per channel, , and the four channels were combined using sum‐of‐squares. Data were registered using a sub‐pixel algorithm.

Both CTI and MGC approaches were used to estimate kurtosis metrics by voxel‐wise fitting of Equations (6) and (9) to the set powder‐averaged data (masked to avoid regions distorted due to b0 inhomogeneities (Figure 4A) using an OLLS procedure. Moreover, for a quantitative analysis of the kurtosis sources, 10 regions of interest (ROIs) were manually defined on the images of rat 1, including the left and right cortical gray matter (GM), the right and left corpus callosum genu (CCg), the right and left corpus callosum body (CCb), the right and left corpus callosum splenium (CCs), and the right and left internal capsule (IC). To ensure consistency across animals, the ROIs for rats 2 and 3 were automatically generated for the other animals using non‐parametric registration of fractional anisotropy maps. ,

FIGURE 4

Robustness of diffusion‐weighted MRI data. A, Representative b = 0 image for all animals (coronal slices) with delineated ROIs: brain mask (MASK); right and left cortical GM; right and left CCg; the right and left CCb; right and left CCs; and right and left IC. B, Representative nominal SNR map that was computed from all b = 0 acquisitions of rat 1 (slice 2). C, Representative diffusion‐weighted images (rat 1, slice 2) for a single diffusion gradient direction (marked in blue) for the maximum total b 1 used in this study (ie, b 1 = µm2/ms and b 2 = 0) before and after PCA denoising (right and left respectively). Note the high signal to noise, and the lack of imaging artifacts. Abbreviation: den, denoised

RESULTS

Figure 2 shows the simulation results for systems containing single component types (models 1‐3). For isotropic Gaussian diffusion components with different mean diffusivities (Figure 2A1), all DDE signal sets reveal identical log‐signal dependencies with (Figure 2A2). The non‐linearity of the log‐signal decays was, thus, correctly identified as non‐zero by all strategies (Figure 2A3‐6). For uniformly distributed anisotropic Gaussian components (Figure 2B1), perpendicular DDE signals evidenced stronger diffusion‐weighted attenuations, as expected (Figure 2B2). The kurtosis estimated from MGC and the new CTI approaches was again correctly attributed to (Figure 2B4,B6). Identically biased magnitudes were observed for both strategies for likely due to higher‐order‐term biases which are expected to be larger for systems with high tissue dispersion. As expected, the old CTI protocol reports non‐zero µK in this case (Figure 2B5), while the new CTI approach overcomes this bias and correctly attributes a zero µK for such systems (Figure 2B6).

For a system with microscopic disorder (Figure 2C1, positive µK), asymmetric DDE signals (ie, ) differ from their symmetric DDE counterparts (ie, and , respectively, Figure 2C2). The finite µK strongly biases both and estimates from the MGC analyses (Figure 2C4). On the other hand, CTI—both old and new protocols—correctly estimate the finite µK. These results are consistent with the non‐vanishing (positive/negative) µK values of restricted diffusion inside spheres—see Supporting Information Figure S1.

In a realistic voxel, all kurtosis sources could exist simultaneously; therefore it is instructive to assess whether the CTI framework can disentangle the kurtosis sources with specificity when all sources are present (Figure 3). For this system, log‐signal DDE decays are different for the three conditions (Figure 3A2). The new CTI protocol successfully estimates all kurtosis sources (Figure 3A3‐6), while the old protocol overestimates µK.

We then investigated how changes in ground truth kurtosis sources would impact the different kurtosis source estimates (see the Methods section). Particularly, when changing the ground truth (Figure 3B) or (Figure 3C), the new CTI and MGC approaches correctly report larger changes in and respectively. The old CTI protocol has limited specificity for µK as varying clearly affects µK. and from the MGC approach also show limited specificity when the ground truth µK is varied (Figure 3D).

Interestingly, we find that CTI new protocol correctly tracked the specific ground truth sources, ie, major changes in , , and µK estimates from the new CTI protocol are only observed when , , and µK ground truths are varied, respectively (Figure 3B3,C2,D4), even if some offsets exist. Note that in Figure 3, due to the biases in MGC estimates introduced by µK, these only match the improved (new) CTI estimates when ground truth µK is zero (c.f. Figure 3D).

Coronal images for all three rats are shown in Figure 4A. Nominal SNR was ~30 for the raw data (Figure 4B). Figure 4C shows a representative diffusion‐weighted image for a dMRI experiment acquired with the maximum b‐value used in this study, before and after PCA denoising. An SNR gain of ~1.3 was noted upon denoising (the nominal SNR of the denoised data was ~40).

Figure 5A shows the powder‐averaged signal decays for the four sets of the improved CTI protocol. The log difference between the powder‐averaged data from set 1 and set 2 (Figure 5B1) shows the sensitivity to µK (c.f. Equation 4), while the log difference between the powder‐averaged data from set 2 and set 3 (Figure 5B4) shows CTI’s sensitivity to (c.f. Equation 5). Positive log differences between set 1 and set 2, reveal a non‐vanishing, positive µK for both GM and WM (c.f. Figure 5B2,B3).

FIGURE 5

Powder‐averaged experimental data for all slices of a representative animal (rat 1). A, Powder‐averaged images for all four DDE sets required for the CTI methodology: set 1 (A1); set 2 (A2); set 3 (A3); and set 4 (A4). B, Raw sensitivity to µK and anisotropic kurtosis (K): raw µK sensitivity maps are quantified by the log difference between powder‐averaged signals of set 1 and 2 (B1); histogram of the log difference between powder‐averaged signals of set 1 and 2 for all concatenated GM and WM regions of interest (B2‐3); raw K sensitivity maps are quantified by the log difference between powder‐averaged signals of set 2 and 3 (B4); histogram of the log difference between powder‐averaged signals of set 2 and 3 for all concatenated GM and WM regions of interest (B5‐6). Note that the data have sufficient power to resolve µK and K even from the raw data itself

The CTI‐driven kurtosis estimates for all slices of rat 1 and the second slice of rats 2 and 3 are shown in Figure 6. Both and were higher in WM regions (Figure 6A‐B). and µK maps show noisier spatial profiles than and maps (Figure 6C‐D). Nevertheless, µK is a prevalent source of kurtosis, with higher values found in GM and lower values found in WM (Figure 6D, red arrows). The ROI analysis supported the trends observed in the kurtosis maps (Figure 7A); in general, µK shows to explain 64 ± 6% and 30 ± 14% of the total kurtosis in GM and WM regions, respectively (Figure 7, panel A4). The higher error bars for and µK is consistent with the lower precision observed in their parametric maps (the errors in our experimental µK estimates are in line to the estimation error predicted in Supporting Information Figure S3, which also shows that if µK was due to noise, the mean would be centered around 0). Figure 7B shows the histograms of the different CTI‐driven kurtosis estimates for combined GM and combined WM ROIs. The mean values of GM and WM are significantly different and above zero (two‐sample t‐test with unequal variances, P < .001 for kurtosis estimates and all three animals).

FIGURE 6

Kurtosis source maps from CTI. A, Total kurtosis (K) map. B, Anisotropic kurtosis (K) map. C, Isotropic kurtosis (K) map. D, µK map. The red arrow in panel D points to lower µK estimates observed particularly in WM brain regions. In each panel, maps are presented for all slices of rat 1 to show the spatial distribution of metrics, and the 2nd slice of rats 2 and 3 to show the reproducibility of the findings in the different rats. Note how the total kurtosis maps, which have the largest values, are decomposed into their underlying sources. Anisotropic kurtosis is highest in WM, while isotropic kurtosis is highest in areas with larger partial volume effects between different diffusivities (namely areas with cerebrospinal fluid). Although the µK maps are slightly noisier, they do display clear GM/WM contrast

FIGURE 7

ROI analysis of CTI metrics in various brain regions. A, Mean and standard deviation (SD) of K (A1), K (A2), K (A3), and µK (A4) estimated from the CTI approach for different regions of interest. The three black intervals in each bar represents the mean and SDs for the three individual rats, respectively. The colored bars show the mean values across all three animals. Regions of interest for these panels include: the right and left cortical GM (GMr and GMl); the right and left CCs (CCsr and CCsl); the right and left CCb (CCbr and CCbl); the right and left CCg (CCgr and CCgl); and the right and left IC (ICr and ICl). Note the excellent right‐left symmetry, and the excellent reproducibility over the different animals. B, Histograms of the K (B1), K (B2), K (B3), and µK (B4) estimates for GM (histograms in green) and WM (histograms in red) ROIs; in each panel, the differences between the mean values from GM and WM ROIs are statistically tested using a two‐sampled t‐test with unequal variances (* for P < .05, ** for P < .01, and *** for P < .001 ***). The histograms reveal that µK is in fact the dominant contributor to the total kurtosis in the GM.

We then turned to assess the impact of the finite µK on MGC analysis, where only diffusion tensor variance is assumed. Figure 8A‐C shows the kurtosis maps obtained from the MGC approach, while Figure 8D shows the histograms of MGC kurtosis values from WM and GM ROIs. In comparison to their CTI counterparts, MGC‐derived was lower in both GM and WM regions (Figure 8A), while MGC and values were higher. As for CTI, MGC mean and appear higher in WM; however, non‐significant differences between the WM and GM voxels were observed for the MGC estimates (Figure 8D3).

FIGURE 8

Analysis of kurtosis bias when assuming MGC. A, MGC total kurtosis () maps. B, MGC anisotropic kurtosis maps. C, MGC isotropic kurtosis maps. D, Histograms of the K (D1), K (D2), and K (D3) MGC estimates for all GM (histograms in green) and WM (histograms in red) regions of interest; in each panel, the differences between the mean values of GM and WM ROIs are statistically tested using a two‐sampled t‐test with unequal variances (* for P < .05, ** for P < .01, and *** for P < .001 ***). When forcing an MGC analysis on the data, the isotropic and anisotropic kurtosis metrics “absorb” the ignored µK, leading to biased maps (compare the MGC‐driven maps in this Figure with those shown in Figure 6A‐C, respectively, and the histograms in this figure with those shown in Figure 7B)

To further investigate the correlation between these metrics, scatter plots of MGC and CTI kurtosis estimates are shown in Figure 9. Points in the scatter plots are color‐coded according to CTI’s µK estimates, showing that higher differences between CTI and MGC estimates are associated with higher degrees of µK.

FIGURE 9

Scatter plots between MGC and CTI estimates. A, Scatter plots between MGC total kurtosis () and CTI total kurtosis () estimates. B, Scatter plots between MGC anisotropic kurtosis () and CTI anisotropic kurtosis () estimates. C, Scatter plots between MGC isotropic kurtosis () and CTI isotropic kurtosis () estimates. Points in the scatter plots are colour coded according to CTI’s µK values. Note that as µK increases, the bias in the MGC estimates becomes more dramatic, especially for isotropic and anisotropic kurtosis sources

DISCUSSION

The conflation of underlying kurtosis sources in SDE was a major motivation in developing multidimensional diffusion encoding approaches. Even under the MGC assumption, and contrasts have shown great promise for, eg, distinguishing different tumor types and grades, , depicting healthy and pathological age‐related microstructural alterations, mapping multiple sclerosis lesions, , and characterizing body organs. , However, the MGC assumption implicitly ignores diffusion‐time dependence , and µK effects, thereby risking the conflation of µK effects into the metrics, which in turn can complicate the interpretation of the metrics. We therefore sought in this study to characterize the insofar ignored µK in vivo and assess its impact on the more conventional MGC approaches.

The CTI approach was recently introduced for µK mapping. CTI goes beyond the tensor‐valued framework and simultaneously estimates , and µK, albeit at the expense of a larger number of acquisitions. The original CTI framework was lengthy and included DDE measurements that could bias µK estimation in some scenarios (eg, Figure 2B5) due to higher order terms, a common issue for techniques based on the truncation of the signal cumulant expansion. , , To alleviate these drawbacks and accelerate the acquisitions, a new and improved CTI strategy was here developed. We found a sparser set of DDE acquisitions for robustly resolving kurtosis sources and found that they can much more accurately estimate µK compared with the older CTI protocol. By balancing the total b‐values used for different DDE sets (c.f. Section C of the Supporting Information), the higher order term effects are greatly diminished (eg, Figures 2B6 and 3). The four CTI quantities (, ) can be fully resolved from only four different combinations of DDE parameters (, , ) (c.f. DDE the sets in Figure 1). In addition, we managed to accelerate CTI from 2 h to under 40 min, thereby making it applicable for in vivo preclinical and even clinical mapping. It is interesting to note that µK alone can be estimated from the log signal differences of two different DDE experiments (c.f. Equation 7), specifically using (1) parallel symmetric DDE gradient waveforms and (2) experiments analogous to SDE with the same total b‐value, much like how microscopic anisotropy is estimated from parallel and perpendicular waveforms , , , , ). Therefore, experiments aiming to resolve only µK could be even further accelerated.

In Ref. [40], we found significant positive µK estimates both in GM and WM brain regions in the ex vivo and in vivo rodent brain, suggesting it cannot be ignored in MGC approaches. Our new results with the improved CTI protocol confirmed the overall non‐vanishing positive µK effects in both WM and GM (Figures 6 and 7) with higher confidence (due to the attenuation of the potentially remaining higher order effects in the previous study performed with the older CTI protocol) and further highlighted significant µK differences between GM and WM brain regions (Figure 7). Positive µK is consistent with non‐Gaussian diffusion effects due to intra‐cellular cross‐section size variance , , and/or the presence of obstacles in tortuous extra‐cellular environments. , Therefore, the µK differences between GM and WM could perhaps be explained by differences in both intra‐ and extra‐cellular microstructural configurations (eg, different compartmental cross‐sectional variance, different degree of cellular packing, etc.), or by the presence of a more negative µK contributions from restricted diffusion in WM. Although µK will depend on a weighted sum of all above‐mentioned effects, positive µK contributions are expected to prevail over negative µK contributions from completely restricted diffusion, as the latter are typically associated with low apparent diffusivities which strongly affects the signal contributed from these compartments (note that the total µK measured by CTI is a weighted average of all its contributions, where the weights depend on the squared apparent diffusivities of each contribution .

Another significant result in this study, is that this non‐vanishing µK can have a dramatic effect on kurtosis sources computed from tensor‐valued and MGC framework (Figures 8 and 9). In general, both and derived from the MGC approach were biased towards higher values compared with their more accurate CTI counterparts. The color‐coded scatter plots between MGC and CTI estimates revealed that differences can be fully explained by µK biases on MGC estimates (Figure 9). In is important to note that the influence of µK on and can obscure microstructural differences, such as the differences between GM and WM brain regions (Figure 9), and can thus bias the interpretation of metrics arising from MGC approaches in the presence of finite µK.

Limitations and future work

Although it does not rely on the Gaussian diffusion assumption, CTI is still a cumulant expansion of DDE signals, which induces some implicit assumptions. Namely, disregarding higher‐order cumulant terms; assuming the long mixing time regime (which can be [and was in this study] empirically evaluated); and ignores exchange. Higher‐order‐term biases were here minimized by the new CTI protocol. The long mixing time regime effectively suppresses unwanted time‐dependent diffusion correlations from the Q and S‐tensors , and guarantees that the Z‐tensor in Equation (5) approaches the covariance tensor. , , , This regime was empirically verified by measuring parallel and antiparallel DDE signals and showing that they produce identical decays , , , , (c.f., Supporting Information Figure S5). Last, since exchange is not yet explicitly modelled by CTI, it may affect the µK metrics. , , , Future studies will investigate the biological underpinning of µK and how exchange between biological components can affect its magnitude.

In this study, simple models were used to investigate the origin of CTI‐driven kurtosis sources, and to illustrate the impact of finite µK on previous MGC approaches. These simple models are, however, are likely not sufficiently complex to fully represent biological tissues. Future studies should expand such in‐silico experiments towards more complex scenarios allowing the assessment of the relationships between different kurtosis sources and concrete (sub)cellular features (eg, cellular cross‐sectional variance, cellular packing, exchange, etc.). , , , , , Another limitation of the current study is the reduced number of animals used. Although the power of (measures have a precision of ~0.2), future studies using a larger number of animals could be performed to explore for detailed regional differences, test‐retest, and inter‐animal variance of µK.

Lundell et al. (2019) have recently shown that continuous diffusion gradient waveforms experiments probing identical b‐tensors but different power spectra can provide different information about diffusion time‐dependence and µK. One could, therefore, argue that MGC‐driven biases could be reduced by adjusting the diffusion parameters of our CTI acquisition protocol, or by entirely excluding some acquisitions. In section H of the Supporting Information, we report MGC kurtosis estimates obtained by fitting Equation (9) to only DDE sets 2, 3, and 4, which correspond to an acquisition with identical waveform profiles, similar to Lundell's idea. Under this condition, MGC and CTI are identical (c.f., Supporting Information Figure S6). However, the MGC estimates from this modified protocol are still a combination of isotropic and µK effects. The comparison between CTI and MGC approaches could assist to define the regimes in which multidimensional MGC estimates are accurate.

Here, we managed to accelerate the CTI scan times to about 40 mins. Although this is still not (yet) sufficiently rapid for clinical translation, we note that the objective was only to identify the minimal acquisition set requirements for the extraction of all CTI quantities. Indeed, we acquired a large and likely redundant number of directions (135 per experiment set) to enhance the precision of our kurtosis estimates. In future studies, further acceleration could be obtained by reducing the number of directions acquired for powder‐averaging. , We note that µK can also be estimated directly from only two sets of the new CTI protocol, albeit at the expense of not resolving and , as was done for the raw µK sensitivity analysis in Figure 5B. Future studies aiming to measure µK should also consider the desired estimation precision/accuracy when designing their experiments; some considerations on the relationship between µK precision and acquisition parameters are described in section D of the Supporting Information.

CONCLUSIONS

The accelerated CTI approach developed here facilitated more accurate in vivo acquisitions, and revealed that the commonly neglected µK is a significant source of total kurtosis in the brain, both in GM and in WM. In fact, µK is the dominant kurtosis source in GM. Ignoring µK leads to significant bias in MGC approaches, underscoring the importance of accounting for µK in multidimensional diffusion encoding approaches. Our findings suggest promising new biomarkers in health and disease.

FIGURE S1 CTI estimates of a single spherical compartment. A) Non‐Gaussian diffusion inside restricted spheres with radius of 3 µm and intrinsic diffusivity of 2 µm2/ms: A1) model schematic representation; A2) the signal decays for three different DDE experiment types in which signals of the improved CTI protocol are marked by the blue triangles; A3) the kurtosis ground truth values; A4) the kurtosis estimates obtained from MGC; A5) the kurtosis estimates obtained from the CTI using its previous “old” protocol; A6) the kurtosis estimates obtained from the CTI using its improved “new” protocol. B) Results from noise free synthetic signals as a function of r: B1) ground truth mean diffusivity; B2) CTI measured mean diffusivity; B3) ground truth microscopic kurtosis; B4) CTI measured microscopic kurtosis. C) Results from noise corrupted signals (SNR=40, mean and std values from 1000 iterations) as a function of r: C1) CTI total kurtosis ; C2) CTI anisotropic kurtosis ; C3) CTI microscopic kurtosis µK; C4) CTI measured isotropic kurtosis K

FIGURE S2 Kurtosis estimates of non‐optimized CTI protocol with parameters b = 1.25 ms/µm2 and b = 0.5 ms/µm2 for two simulations (model 2 and 4 described in main manuscript); A) Evenly oriented anisotropic Gaussian diffusion components; B) Sum of different compartment types (isotropic Gaussian + anisotropic Gaussian + isotropic non‐Gaussian components). From left to right, each panel shows: a schematic representation of the toy models (A1, B1); the kurtosis ground truth values (A2, B2); and the kurtosis estimates obtained from the CTI using its “old” non‐optimized protocol (A3, B3)

FIGURE S3 Microscopic kurtosis estimation uncertainty. A) Microscopic kurtosis uncertainty prediction for different ground truth mean diffusivities D and microscopic kurtosis values and given the acquisition parameters of this study (N=135 and b = 2.5 ms/µm2) and a single acquisition nominal SNR of 40 (ie, ); B) prediction for different ground truth D and values and given the acquisition parameters of this study (N=135 and b = 2.5 ms/µm2) and a single acquisition nominal SNR of 20 (ie, ); C) Required SNR to obtain a =0.05 and given the acquisition parameters of this study (N=135 and b = 2.5 ms/µm2). In all panels, the mean D and values observed on the white matter and grey matter ROIs of the rat brain data is marked in by the red and green points, respectively. D) histogram of estimates obtained from signals computed using numerical simulations (signals from model 3) and corrupted with Rician noise at a nominal SNR of 40 for three ground truth set of parameters representing: 1) D = 0.8 µm2/ms, and (histogram marked by the grey line); 2) D = 0.76 µm2/ms, and (histogram marked by the green line); 3) D = 0.82 µm2/ms, and (histograms marked by the red line). These latter histograms were produced for 1000 simulation iterations and the ground truth values are marked by the dashed lines

FIGURE S4 Midsagittal T2‐weighted images of all three animals. These images were used as a reference for placing the 3 coronal slices for diffusion MRI acquisition – positions of the coronal slices are marked in red

FIGURE S5 Powder‐averaged data for the extra experiments performed on Rat #1. A) Powder‐averaged signal decays for DDE parallel acquisitions (repetition #1); B) Powder‐averaged signal decays for DDE parallel acquisitions (repetition #2); C) Powder‐averaged signal decays for DDE antiparallel acquisitions. D) Map of the log differences values between powder‐averaged signals of set #5 repetitions #1 and #2 (brain mask region delineated in blue, and B0 inhomogeneities artifacts are pointed by the red arrow). E) Map of the log differences values between powder‐averaged signals of set #5 and #6 brain mask region delineated in blue, and B0 inhomogeneities artifacts are pointed by the red arrow). F) Histograms of the values inside the brain mask for the log differences values between powder‐averaged signals of set #5 repetitions #1 and #2 (grey) and the log differences between powder‐averaged signals of set #5 and #6 (blue). Note parameters was set to 1 ms/µm2 (i.e total b‐value = 2 ms/µm2)

FIGURE S6 MGC kurtosis estimates using only data from sets #2, #3, #4 and its comparison to CTI estimates using all data. A) MGC total kurtosis () maps. B) MGC anisotropic kurtosis maps. C) MGC isotropic kurtosis maps. D) Scatter plots between MGC anisotropic kurtosis () and CTI anisotropic kurtosis () estimates; E) Scatter plots between MGC isotropic kurtosis () and CTI isotropic kurtosis () estimates

TABLE S1 Summary of the complete DDE parameter combination used for the CTI estimates in Henriques et al. (2020). The 12 parallel directions of the 5‐design and the 60 perpendicular directions of Jespersen’s scheme are reported by (Jespersen et al., 2013)

TABLE S2 Summary of the reference non‐optimized CTI protocol used on the current study

Click here for additional data file.