Microscopic anisotropy misestimation in spherical-mean single diffusion encoding MRI.

PURPOSE: Microscopic fractional anisotropy (µFA) can disentangle microstructural information from orientation dispersion. While double diffusion encoding (DDE) MRI methods are widely used to extract accurate µFA, it has only recently been proposed that powder-averaged single diffusion encoding (SDE) signals, when coupled with the diffusion standard model (SM) and a set of constraints, could be used for µFA estimation. This study aims to evaluate µFA as derived from the spherical mean technique (SMT) set of constraints, as well as more generally for powder-averaged SM signals.
METHODS: SDE experiments were performed at 16.4 T on an ex vivo mouse brain (Δ/δ = 12/1.5 ms). The µFA maps obtained from powder-averaged SDE signals were then compared to maps obtained from DDE-MRI experiments (Δ/τ/δ = 12/12/1.5 ms), which allow a model-free estimation of µFA. Theory and simulations that consider different types of heterogeneity are presented for corroborating the experimental findings.
RESULTS: µFA, as well as other estimates derived from powder-averaged SDE signals produced large deviations from the ground truth in both gray and white matter. Simulations revealed that these misestimations are likely a consequence of factors not considered by the underlying microstructural models (such as intercomponent and intracompartmental kurtosis).
CONCLUSION: Powder-averaged SMT and (2-component) SM are unable to accurately report µFA and other microstructural parameters in ex vivo tissues. Improper model assumptions and constraints can significantly compromise parameter specificity. Further developments and validations are required prior to implementation of these models in clinical or preclinical research.

INTRODUCTION

Diffusion MRI is perhaps one of the most valuable reporters for dimensions much smaller than the MRI voxel size, and thus it has been extensively utilized to characterize microstructural changes both in health and in disease.1, 2, 3, 4, 5 Numerous methods, mainly relying on the Stejskal‐Tanner technique6 – later dubbed more generally single diffusion encoding (SDE)7 – were developed for characterizing diffusion anisotropy,8, 9 leading to robust white matter (WM) orientation mapping in vivo.5 Other methods deriving diffusion anisotropy include, among others, quadratic‐form analysis10 and analysis of displacement profiles in q‐space imaging (QSI).11, 12, 13

Diffusion anisotropy is highly dependent on the mesoscopic tissue organization.14, 15 An emerging diffusion MRI field of research strives to decouple microscopic diffusion properties from the mesoscopic tissue organization. Methods based on double diffusion encoding (DDE)16 have been proposed for decoupling microscopic diffusion anisotropy (or its normalized form, microscopic fractional anisotropy (µFA)) from orientation dispersion using displacement correlations.17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 Isotropic diffusion encoding can also estimate μFA if multiple Gaussian components are assumed.30, 31, 32, 33, 34, 35, 36

An alternative approach for decoupling microscopic diffusion properties from the mesoscopic tissue organization using the more commonly available SDE methodology involves microstructural modeling.37, 38, 39, 40 This was pioneered by Stanisz et al. in the characterization of water diffusion in bovine optical nerve41 and used early on to quantify directional uncertainty in tractography42, 43, 44 and to estimate fiber caliber.45, 46, 47 Typically, biophysical models assume that biological tissues can be represented by a sum of nonexchanging Gaussian diffusion components,48, 49, 50, 51, 52, 53 and a specific instance of this class of models is the so‐called standard model (SM).40, 54 In this physical picture, the models can be expressed as a convolution between the mesoscopic orientation distribution function (ODF) and a kernel containing scalar microstructural parameters (e.g., neurite density, neurite diffusivity).49, 50, 51, 52, 53, 54, 55, 56

The SM fitting landscapes are rather flat,56, 57 and therefore the number of scalar parameters that can be estimated can be limited, depending on the b value regimes. For a small number of relatively low b‐values, strong constraints must be imposed. For example, the Neurite Orientation Dispersion and Density Imaging (NODDI) model constrains the diffusion coefficients to constant values and assumes a predetermined function (Watson52 or Bingham53 distributions) for the ODF. If higher b‐values are reached, softer constraints are usually imposed: the Neurite Density Model (NDM) truncates spherical harmonics at order 4 to represent the ODF,49 while LEMONADE estimates all SM parameters using signal moments up to the sixth order.56

It was realized early on that in powder‐averaged systems,58, 59, 60, 61, 62, 63, 64 the ODF can be treated as a constant, and thus scalar microstructural parameters can be estimated independently to specific ODF constraints. Water diffusing in endosperm tissue,59 3He gas diffusing in lungs,60 and metabolites diffusing (predominantly) in randomly oriented61, 62, 63, 64 neuronal dendrites and axons, were all assumed to constitute powder‐averaged systems. More recently, it was noticed that, rather than using powder‐averaged systems, one could use powder‐averaged signals to remove the orientational complexity for estimating scalar signal parameters, such as µFA.28, 29, 30, 31, 32, 33, 34, 35, 36 Recently, the signal powder‐averaging approach was extended to SDE signals and coupled with a constraint of a single Gaussian component for estimating µFA in neural tissues.65 The framework was later expanded to accommodate two (Gaussian) components representing intracellular and extracellular domains, such that the “axonal” volume fraction could be measured.66 However, further constraints become necessary in such expansions (due to the relatively low and few b‐values used): specifically, it was assumed that (1) extracellular diffusion follows the tortuosity model67; (2) axial diffusivities of intracellular and extracellular components are equal; and (3) intracellular radial diffusivity is zero (Table 1). The coupling between powder‐averaged SDE signals and imposition of these specific constraints was termed spherical mean technique (SMT),65, 66 and we refer to them hereafter as SMT1 and SMT2, respectively, where 1 or 2 represents the number of components assumed by each model. Some of these constraints can be released when higher b‐values are reached.56

Table 1

Summary of the four SM set of constraints used in this studya

Name	Free model parameters	n	Constraints
SMT1	λ_‖ and λ_⊥	2	1) 1 component
SMT2	λ and f	2	1) 2 components
			2) λ⊥e=(1-f)λ
			3) λ‖i=λ‖e=λ
			4) λ⊥i=0
SM3	λ⊥e,λ, and f	3	1) 2 components
			2) λ‖i=λ‖e=λ
			3) λ⊥i=0
SM4	λ‖i,λ‖e, λ⊥e and f	4	1) 2 components
			2) λ⊥i=0

Apart from the constraints mentioned, it is assumed that diffusion components have only 2 unique eigenvalues, and .

Imposing constraints on microstructural models may be a practical necessity, but it may also corrupt the parameter specificity.40 Indeed, the assumption that biological tissues can be well described by a small number of Gaussian diffusion components might be overly simplistic, given the large heterogeneity of microscopic structures within a typical voxel containing axons, neurites, neuronal bodies, glia cell bodies, glial processes, myelin, etc.1 Moreover, even in relatively homogeneous voxels (e.g. voxels predominantly containing WM), diffusion heterogeneity might exist as a result of size polydispersity (e.g., distribution of axon diameter)68, 69, 70, 71 or nonzero cell body fractions.72 In principle, such effects will lead to signal kurtosis, with two main physical origins: (1) for Gaussian diffusion, the diffusivity variance from multiple components can induce nonexponential signal decays73, 74, 75, 76, 77 (intercomponent kurtosis); (2) restricted, time‐dependent diffusion,78, 79, 80, 81, 82, 83, 84, 85, 86 which has been observed in biological tissues, can contradict the Gaussian assumption within individual microstructural compartments (intracompartmental kurtosis). Diffusion kurtosis can also arise from exchange,75 and any combination of the kurtosis factors mentioned can be envisioned in a given voxel.

Here, we aimed to investigate the validity and specificity of µFA derived from powder‐averaged SDE signals. We present a few simple theoretical considerations of current SMT and powder‐averaged SM assumptions and the associated parameter constraints vis‐à‐vis kurtosis. Experiments in ex vivo mouse brains show large discrepancies between DDE‐derived µFA and SM‐driven SDE counterparts. We explore in simulation how violation of the constraints mentioned can affect µFA estimation. Our findings suggest that powder‐averaged SM estimates do not capture the ground truth, and that intravoxel heterogeneity may be a culprit.

THEORY

Definition of microscopic fractional anisotropy

Diffusion in biological tissues is often represented by an ensemble of microscopic components (or compartments),30, 40, 77 with the ith component characterized by its own individual diffusion tensor (Figure 1). The ensemble‐averaged diffusion tensor then equals the average of the individual diffusion tensors of all components contributing signal to the ensemble:

Figure 1

Schematic representation of five different environments and respective individual diffusion tensors and ensemble diffusion tensor . All 5 microenvironments consist of cylindrical structures where water can diffuse in the interior and/or exterior spaces (random trajectories of water molecules diffusing in the interior and exterior of cylinders are illustrated by the lines marked by 1 and 2, respectively). A‐C, Considering only water inside the structures, diffusion might be represented by homogeneous individual tensor with constant microscopic fractional anisotropy (μFAi). While the fractional anisotropy (FA) of tensor depends on both μFAi (differences between A and B) and dispersion orientation of the individual tensors (differences between A and C), the ensemble μFA depends only on μFAi. D and E, Considering water inside and outside cylindrical microstructures, environments are better characterized by 2 types of tensors . For these latter cases, μFA is still independent of orientation dispersion and depends on the ensemble anisotropy of the diffusion tensors

where brackets represent an ensemble average. The fractional anisotropy of (equivalent to DTI’s FA) can be expressed as31

where is the trace of , and is the variance of the eigenvalues of . The µFAi is the FA of each individual diffusion tensor , given by31, 34

For a better understanding of the quantities defined by Equations (1), (2), (3), individual diffusion tensors and their respective ensemble diffusion tensor for 5 different environments are illustrated in Figure 1. The environment represented in Figures 1A and 1B consists of well‐aligned structures with intracylindrical signal only, but with .80 or .65, respectively. Figure 1C shows that FA differs substantially from for finite orientation dispersion.

For the sake of simplicity, all microscopic diffusion tensors in Figure 1A‐C are assumed to have constant . A different and often‐used representation of WM consists of two diffusion tensors representing signals from intracylindrical and extracylindrical spaces, respectively (Figure 1D‐E). In this case, it might be impractical to measure the microscopic fractional anisotropy of individual microscopic components. Therefore, it is useful to define the effective μFA as31, 34

where is the average eigenvalue variance of microscopic diffusion tensors and is the average trace of microscopic diffusion tensors It is important to note that the effective µFA is not the simple average of (i.e. ). Moreover, according to Equation 4, the effective μFA is expected to be invariant to tissue organization (as illustrated in Figure 1D, E). The preceding definitions are general for any environment containing an indefinite number of different types of diffusion components.

µFA estimates from powder‐averaged SDE signals

In the original SMT work, the powder‐averaged signal was constrained to one Gaussian component65; i.e. it is assumed that any voxel can be fully described by one axial and one radial diffusivity (Figure 1A‐C). Under this constraint (SMT1), the powder‐averaged signal is given by65

where is the powder‐averaged signal decay, b is the b‐value, and erf() is the error function. and can then be obtained by fitting Equation 5 to measured with at least 2 nonzero b‐values. From the estimated and and under the same single‐component assumption, μFA from SMT1 can be computed as65

Although Equation 6 is equivalent to Equations 2, 3 and 4 for systems containing exactly one component with constant and values (Figure 1A‐C), this is not the case for heterogeneous environments (e.g. Figure 1D‐E). Equation 6 can be easily expanded to the standard model comprising two components (e.g. intracellular and extracellular domains). The powder‐averaged signal for the 2‐component standard model (SM) is

where and are the axial and radial diffusivities of the intracellular component, and are the axial and radial diffusivities of the extracellular components, and is volume fraction of the intracellular component. Introducing these model parameters into Equation 4, μFA for the 2‐component SM can be calculated as

As previously demonstrated,56, 57 Equation 7 presents a flat fitting landscape. To increase the fit robustness, Kaden et al.66 proposed to impose the following constraints: (1) tortuosity constraint,67 ; (2) equal axial diffusivities constraint, ; and (3) stick constraint, = 0. We term this set of constraints SMT2 (Table 1). Equations 7 and 8 can then be rewritten as

and

where .

Equation 10 shows that, under the constraints imposed by Kaden et al.,66 μFA only depends on the “neurite” fraction . This theoretically reveals that Equation 10 cannot properly quantify µFA for arbitrary tissue configurations. For instance, for loosely packed randomly oriented fibers, the low extracellular diffusion anisotropy should naturally decrease the effective µFA. However, since only depends on f, Equation 10 will fail to capture such effects.

Cumulant expansion of the standard model powder‐averaged signals

To assess how kurtosis contributes to SM, it is useful to express first other easily interpretable quantities under the assumption that SMT assumptions are valid. For instance, the second‐order cumulant expansion74, 87 can be used to derive the effective diffusion and effective kurtosis of the powder‐averaged signals:

For the 1‐component SM (SMT1), the theoretical diffusion and kurtosis of the powder‐averaged signals can be calculated by expanding Equation 5 to second‐order in b:

and

Equation 13 shows that SMT1 only considers kurtosis from the variance across the eigenvalues of a single diffusion component.

Analogous to Equations 12 and 13, the diffusion and kurtosis for the 2‐component powder‐averaged SM can be computed by expanding Equation 7 to second order in b:

and

where the diffusivities and correspond to the mean diffusivities of the intracellular and extracellular components (i.e. and . Note that the right‐side terms of Equation 15 are related to the eigenvalue variances of both intracellular and extracellular components and the isotropic variance between and .

Under SMT2 constraints, Equations 14 and 15 can be simplified to

and

Equation 17 shows that depends solely on . Other kurtosis effects, e.g. arising from eigenvalue variance, are no longer accounted for.

µFA estimates from powder‐averaged DDE signals

In contrast to the SDE methodology, estimates from DDE ideally do not rely on specific microstructural model assumptions or parameter constraints. Jespersen et al. showed that the average eigenvalue variance can be computed from signal decays acquired for 2 parallel diffusion gradients and signal decays acquired for 2 perpendicular diffusion gradients 28:

where and are the powder‐averaged signals of and , using the 5‐design orientation scheme.28 The average eigenvalue variance extracted from Equation 18 can then be applied to Equation 4 to obtain a model‐free estimate of In this study, of Equation 4 is set to the effective diffusion computed from Equation 11.

METHODS

Specimen preparation and MRI experiments

All animal experiments were preapproved by the local ethics committee operating under local and EU law. A mouse brain was perfused intracardially from a healthy adult animal (N = 1) and was then immersed in 4% Paraformaldehyde (PFA) solution for 24 h, followed by immersion in a Phosphate‐buffered saline (PBS) solution for at least 48 h. The specimen was placed in a 10‐mm NMR tube filled with Fluorinert (Sigma Aldrich, Lisbon, PT) and scanned at 37 oC using a 16.4 T Bruker scanner (Karlsruhe, Germany) equipped with a Micro5 probe and a gradient system capable of producing up to 3000 mT/m in all directions.

For the generation of the gold‐standard μFA map, diffusion‐weighted signals were acquired using DDE‐MRI for 19 evenly sampled b‐values from 0 to 4.5 ms/μm2 (Δ = 12, δ = 1.5 ms, and mixing time = 12 ms) and 72 pairs of directions per b‐value following Jespersen et al.’s 5‐design.28 To achieve a high SNR, acquisitions were repeated 14 times and averaged. Moreover, data for each b‐value were further denoised using a Marchenko‐Pastur‐PCA denoising procedure (with an 8 × 8 sliding window),88 while Gibbs ringing artifacts were suppressed using a subvoxel shift algorithm.89 The μFA maps are then computed from the powder averages of and and using Equation 18 (higher‐order corrections90 were applied). It is important to note that the DDE‐MRI sequence was recently validated and shown to be robust against gradient artifacts (such as concomitant fields) under the experimental parameters used in this study.90

The SDE experiments were then acquired for 5 evenly spaced b‐values from 0.5 to 2.5 ms/μm2 (2 averages), 5 b‐values from 3 to 5 ms/μm2 (4 averages), 5 b‐values from 5.5 to 7.5 ms/μm2 (10 averages), and 3 b‐values from 8 to 9 ms/μm2 (14 averages). The increasing number of averages were designed to compensate for the inherently lower SNR of high b‐value data. For each b‐value, SDE experiments were performed along 72 directions. Other MRI parameters common to both DDE and SDE experiments were the following: navigated 2‐shot EPI acqusitions with double‐sampling and a bandwidth of 384615 Hz; 3 coronal slices with thickness of 0.7 mm; FOV = 11×8.64 mm2, matrix size = 92 × 72 leading to an in‐plane resolution of 0.12 × 0.12 mm2; partial Fourier factor = 1.25; TR/TE = 3000/50.4 ms.

The SDE data underwent the same preprocessing as in DDE (MPPCA and Gibbs unringing). The standard model was then used to extract μFA estimates from powder‐averaged signals. For this, the 4 sets of constraints (SMT1, SMT2, SM3, and SM4) summarized in Table 1 were tested. According to Novikov et al.56 SM4 has 2 plausible branch solutions ( and ) for low b values, and we therefore fit data up to rather high b values of 9 ms/µm2. All models were fitted using grid search sampling parameters over 30 evenly spaced values within their plausible ranges (diffusivities between 0 and 3 ms/μm2 and between 0 and 1). The set of parameters with lowest sum of squared residuals were then used as the initial guess of a nonlinear least‐squares fitting procedure, which was implemented in Matlab® (The Mathworks, Natick, MA, USA).

Simulations

The effects of the model assumptions and parameter constraints were evaluated using 5 experiments based on synthetic data, where ground‐truth μFA values are known a priori. To assess the robustness of μFA estimates toward violation of model assumptions separately from the violation of parameter constraints, Simulations 1‐3 were performed only for SMT1, which estimates both and of a single component without imposing further constraints. The robustness of the 2‐component SM‐based μFA estimates (SMT2, SM3, and SM4) was then assessed in Simulations 4 and 5. All simulations were generated noise‐free. Synthetic signals were produced along the same gradient directions and for the same b‐values as the ex vivo SDE experiments. Each simulation is described in detail in the following.

1. Multiple randomly oriented diffusion tensors (intercomponent kurtosis > 0): Intercomponent kurtosis was assessed by simulating an ensemble of 10,000 randomly oriented microscopic components characterized by different tensors (no restricted diffusion, dispersion in the magnitude of ). The diffusion tensor of each microscopic component was generated by different combinations of both and giving constant μFAi values (μFAgt = μFAi) (gt, ground truth), but varying MDi. The MDi values of the 10,000 components were randomly drawn from a Gaussian distribution with mean MDi of 0.8 μm2/ms and a given MDi standard deviation (negative MDi values were zeroed). This experiment was also repeated for MDi values sampled from a bimodal Gaussian distribution with modes set to 0.8 and 1.7 μm2/ms.

2. Multiple diffusion tensors to represent a distribution of axons with nonzero radius (intercomponent kurtosis > 0, no restricted diffusion): One way of treating WM diffusion is to assume that diffusion transverse to axons of finite diameters is adequately described by the Gaussian approximation, hence producing a mapping between axon sizes and apparent transverse diffusivities. To assess kurtosis in this “intercomponent” case, diffusion‐weighted signals of cylinders with different radii were generated using the MISST package.91 Apparent transverse diffusivities were then set to the radial diffusivity obtained by fitting the standard DTI to individual MISST signals. Radii from a log‐normal distribution with mean (m) and standard deviation (std) values were adjusted according to known diameter distributions in rat corpus callosum genu (m = 1 µm, std = 0.5 µm) and body (m = 1.5 µm, std = 0.7 µm), as well as in rat spinal cord WM (m = 3.0 µm, std = 0.6 µm).68, 71 The intrinsic diffusivity of these simulations was set to 2.5 µm2/ms. Simulations were repeated for different cylinder lengths between 5 and 40 µm and repeated for an added extracellular Gaussian component with volume fraction 0.3, 1.7 µm2/ms, and = 0.4 µm2/ms. Ground‐truth values were computed from a simulated DTI experiment on the same perfectly aligned cylinders.

3. Multiple restricted cylindrical compartments to represent a distribution of axons with nonzero radius (intercompartmental and intracompartmental kurtosis > 0): Next, we simulated restricted diffusion in realistic compartment size distributions of WM, containing both intercompartmental and intracompartmental kurtosis. For this, the signal of log‐normally distributed cylinders was generated by directly averaging the signals obtained from the MISST package.91 The dimensions of the cylinders and their log‐normal distributions were identical to those described previously.68, 71 This experiment was also repeated with an added extracellular Gaussian component.

4. Testing the robustness of the SMT2 constraints: To assess the robustness of SMT2 to the violation of its model constraints, synthetic signals based on the 2‐component model, which represents extraneurite and intraneurite domains as two different diffusion tensors (no intracomponent kurtosis nor time‐dependent diffusivities), are generated. The deviation of the tortuosity model was first tested by performing simulations for extraneurite values, sampled between 0 and 2.0 μm2/ms independently of f, and while and were set to 2 μm2/ms. This experiment was then repeated for = 2.3 μm2/ms and = 1.7 μm2/ms (testing the violation of equal axial diffusivities).

5. Multiple Gaussian components and restricted compartments to test the robustness of 2‐component SM constraints: The 2‐component simulations were then applied to all 3 different sets of the 2‐component SM contraints (SMT2, SM3, SM4). To assess intercomponent kurtosis effects not considered by SM, synthetic signals were also generated for multiple components (as described in experiment 2) representing a distribution of axons with different radii together and an extracellular component with parameters set to = 1.7 µm2/ms, = 0.4 µm2/ms, and 0.3. For the sake of simplicity, these simulations were performed only for log‐normal radii distribution with m = 1.5 µm and std = 0.7 µm. To assess the added effect of intracompartmental kurtosis, the latter simulation was then repeated by replacing the sum of Gaussian components with restricted diffusion cylindrical compartments (as described in experiment 3). Moreover, since previous studies suggested that brain tissues might also contain an isotropic restricted diffusion compartment (e.g. representing cell bodies),41, 47, 92, 93 the SM fits were further tested on a simulation where a fraction of isotropic compartments was added. For the sake of simplicity, the diffusivity of these isotropic compartment was set to 0 (also known as a “dot” compartment).93 Ground‐truth values were computed based on the apparent diffusivities of individual compartments computed from individual simulated DKI experiments90.

To assess the robustness of the SM nonlinear least‐squares fitting procedure to its initial guess, the simulations of experiment 5 were processed using a random initializer. This procedure was repeated for 100 sets of randomly sampled initial guesses and compared to the grid search sampling fitting procedure.

RESULTS

Experimental validation

Figure 2A shows the maps obtained for the 3 DDE data slices (hereafter referred to as maps). Under the assumption that DDE at the long mixing time regime provides a relatively model‐free estimate of µFA, these maps can be considered as a “ground truth” for this tissue. As expected, WM areas such as corpus callosum or internal capsule exhibit very high , while gray matter (GM) areas exhibit lower (typically lower than 0.5). The µFA maps calculated from SDE signals using the SMT1 constraints () for the same slices are shown in Figure 2B. The values in the entire brain are much higher than . Notice that is still higher in WM compared with GM and approaches unity. Figure 2C similarly shows the µFA maps for SMT2 constraints (). In this case, is clearly much lower than in GM, while seems to be overestimated in WM. Figures 2D and 2E plot µFA from powder‐averaged SDE fit using SM3 and SM4 constraints ( and , respectively). Both and appear to be overestimated compared to in GM and in WM. For reference, plots assessing the quality of the raw SDE and DDE data and standard DTI FA maps are shown in Supporting Information.

Figure 2

Results for the ex‐vivo mouse brain experiments: A, μFA gold‐standard maps of the 3 DDE data slices; B, μFA maps of the SDE data obtained from SMT1; C, μFA maps of the 3 SDE data slices obtained from SMT2; D, μFA maps of the 3 SDE data slices obtained from SM3; E, μFA maps of the 3 SDE data slices obtained from SM4

To quantify these aspects better, Figure 3 displays µFA estimates for different sets of SM constraints plotted against the “ground‐truth” from all voxels. For reference, identity and regression lines are marked by red dashed and black solid lines, respectively. For a better assessment of differences between WM and GM estimates, vertical lines that roughly segment these brain tissues estimates are shown. Compared to , appears to be overestimated always (Figure 3A), while appears to be mainly underestimated, particularly for GM voxels (Figure 3B). is also overestimated in higher µFA regions (mainly in WM). and are closer to (Figure 3C and D, respectively); however, they are still higher than the gold standard.

Figure 3

Correlations between μFA estimates from SDE and DDE experiments. A, μFA estimates from SMT1 plotted as a function of the μFA estimates from DDE; B, μFA estimates from SMT2 plotted as a function of the μFA estimates from DDE; C, μFA estimates from SM3 plotted as a function of the μFA estimates from DDE; D, μFA estimates from SM4 plotted as a function of the μFA estimates from DDE. For reference, the identity and regression lines are marked by the red dashed and black solid lines, while the DDE value that roughly segments WM and GM is marked by the gray vertical line

Parametric maps for all SMT2, SM3, and SM4 model parameters are shown in Figure 4. For all sets of constraints, estimates are close to 1 in WM and close to 0 in GM regions (Figure 4A1, B1, and C1 for SMT2, SM3, and SM4, respectively). maps derived from SM3 exhibit low contrast between gray and white matter regions (Figure 4B2), while maps calculated from SM4 (Figure 4C2) show lower values than maps (Figure 4C3). estimated from SM3 (Figure 4B3) is similar to extracted from SM4 (Figure 4C4). The relationships between different model parameters are assessed in the lower panels of Figure 4. In general, extracted from both SM3 and SM4 is lower than the values predicted by the tortuosity model (Figure 4D1 and D2). No dependence was observed between and obtained from SM4 (Figure 4D3).

Figure 4

SMT2, SM3, and SM4 parameters for the ex vivo mouse brain experiments. A, Parameter maps from SMT2: A1, axonal volume fraction f; A2, intrinsic diffusivity ; B, parameter maps from SM3: B1, axonal volume fraction f, B2, axial diffusivity , B3, extracellular radial diffusivity ; C, parameter maps from SM4: C1, axonal volume fraction f; C2, extracellular axial diffusivity , extracellular axial diffusivity ; C4, extracellular radial diffusivity . D, relationship across parameters: D1, vs. for SM3; D2, vs. for SM4; D3, vs. for SM4. In these latter plots, values corresponding to DDE higher than 0.5 are plotted in blue (mainly WM voxels), while voxels with DDE lower than 0.5 are plotted in red (mainly GM voxels)

1. Multiple randomly oriented diffusion tensors (intercomponent kurtosis > 0): Simulations for multiple randomly oriented diffusion tensor distributions are shown in Figure 5. Figure 5A1 shows unimodal distributions, as would be encountered for example in the intraneurite space alone. The corresponding estimates of , , and are plotted as a function of ground‐truth values. Even for modest intercomponent kurtosis, and , are overestimated and underestimated, respectively, leading to a large overestimation of . To represent another distinct component, as would for example emerge from extracellular space, Figure 5B1 plots bimodal tensor distributions. Figure B2‐4 shows that the trends mentioned for overestimates in and underestimates in are much exacerbated, leading to a very large overestimation of .

Figure 5

Results from simulations based on multiple Gaussian components: A1, Unimodal Gaussian distribution of the MDi values sampled for different standard deviations; A2, SMT1 estimates for the unimodal distributions plotted as function of ground‐truth μFA values; A3, SMT1 estimates for the unimodal distributions plotted as function of ground‐truth μFA values; A4, SMT1 μFA estimates for the unimodal distributions plotted as function of its ground‐truth values; B1, bimodal Gaussian distribution of the MDi values sampled for different standard deviations; B2, SMT1 estimates for the bimodal distributions from SMT1 plotted as function of ground‐truth μFA values; B3, SMT1 estimates for the bimodal distributions from SMT1 plotted as function of μFA ground‐truth values; B4, SMT1 μFA estimates for the bimodal distributions from SMT1 plotted as function of its ground‐truth values

2. Multiple diffusion tensors to represent a distribution of axons with nonzero radius (intercomponent kurtosis > 0, no restricted diffusion): Figure 6 shows simulations for multiple Gaussian components, distributed under the assumption that axon sizes could be translated to Gaussian diffusivities. Simulated ground‐truth distributions from selected tissues are shown in Figure 6A. Figure 6B‐C plots , , and as a function of the ground‐truth . When a hindered extracellular component is not considered (upper panel), , , and (solid lines) are quite close to their ground‐truth values (dashed line, Figure 6B). However, as soon as another component is added, , , and are again subject to very large errors (Figure 6C).

Figure 6

Results from simulations based on multiple Gaussian components representing log‐normal distributed axons. A, The 3 log‐normal distributions of ground‐truth radial diffusivities adjusted according to previous histological diameter measures of corpus callosum genu (CCg), corpus callosum body (CCb), and spinal cord (SC); B1, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diffusion distributions plotted as a function of ground‐truth values; B2, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diffusion distributions and plotted as a function of ground‐truth values; B3, μFA estimates from SMT1 (solid lines) and μFA ground‐truth values (dashed lines) for the different diffusion distributions and plotted as a function of ground‐truth values; C1, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diffusion distributions involved by a single extracellular diffusion component; C2, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diffusion distributions involved by a single extracellular diffusion component; C3, μFA estimates from SMT1 (solid lines) and μFA ground‐truth values (dashed lines) for the different diffusion distributions involved by a single extracellular diffusion component

3. Multiple cylindrical compartments to represent a distribution of axons with nonzero radius (intercompartmental and intracompartmental kurtosis > 0): To investigate whether restricted diffusion – i.e. intracompartmental kurtosis – effects may be significant, Figure 7 presents simulated signals from realistic axon size distributions under the assumption of fully restricted diffusion within the axons. Figure 7A plots the distributions, which are identical to those used in Figure 6A. Figure 7 plots , , and (solid lines) and their ground‐truth values (dashed lines), as a function of (finite) cylinder length. When only intracellular signals are considered, does not match their ground‐truth values (Figure 7B2), leading to underestimated for cylinders with lengths smaller than 10 µm and overestimated for larger lengths (Figure 7B3). When an extracellular component is added, and are overestimated and underestimated, respectively (Figure 7C1 and 7C2). Consequently, is overestimated independent of compartment length (Figure 7C3).

Figure 7

Results from simulations generated based on multiple cylindrical compartments representing log‐normal distributed axons: A, The 3 log‐normal distributions of ground‐truth diameters adjusted according to previous histological diameter measures of rat corpus callosum genu (CCg), rat corpus callosum body (CCb), and rat’s spinal cord (SC); B1, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diameter distributions plotted as a function of ground‐truth compartment length; B2, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diameter distributions and plotted as a function of ground‐truth compartment length; B3, μFA estimates from SMT1 (solid lines) and μFA ground‐truth values (dashed lines) for the different diameter distributions and plotted as a function of ground‐truth compartment length; C1, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diameter distributions involved by a single extracellular diffusion component; C2, estimates from SMT1 (solid lines) and ground‐truth values (dashed lines) for the different diameter distributions involved by a single extracellular diffusion component; C3, μFA estimates from SMT1 (solid lines) and μFA ground‐truth values (dashed lines) for the different diameter distributions involved by a single extracellular diffusion component

4. Robustness of the SMT2 parameter constraints: Figure 8A shows a perfect match between (solid lines) and ground‐truth values (dashed lines) for a 2‐component system produced exactly according to SMT2 constraints. However, as soon as the values of and deviates from the tortuosity constraint (), significant errors in , , and arise (Figure 8B). The red line signifies that SMT2’s tortuosity constraint is satisfied, and = = 2 μm2/ms. Figure 8C shows the errors in , , and for arbitrary and ground‐truth values and for 2.3 μm2/ms and μm2/ms. Notice the strong negative biases on , , and even when the tortuosity assumption is met (dashed red line).

Figure 8

Robustness of the SMT2 parameter constraints. A, μFA estimates from SMT2 when all SMT2 assumptions are met; B1, errors of the estimates from SMT2 when the tortuosity constraint is not met; B2, errors of the λ estimates from SMT2 when the tortuosity constraint is not met; B3, errors of the μFA estimates from SMT2 when the tortuosity constraint is not met. C1, Errors of the estimates from SMT2 when the tortuosity and axial diffusivities constraints are not met; C2, errors of the λ estimates from SMT2 when the tortuosity and axial diffusivities constraints are not met; C3, errors of the μFA estimates from SMT2 when the tortuosity and axial diffusivities constraints are not met. In panels B1‐3 and C1‐3, the ground‐truth and f pair of values that match SMT2’s tortuosity constraints are marked by the red dashed line

5. Multiple Gaussian components and restricted compartments to test the robustness of 2‐component SM constraints: The scatterplots of Figure 9 show pairs of µFA and estimates obtained for different SM constraint sets and for different tissue scenarios. For each set of constraints (SMT2, SM3, SM4), estimates obtained from the grid search initialization approach are plotted by cross‐marked points, while estimates obtained from the random initialization are plotted by the dot‐marked points. For all tissue scenarios, estimates from random initializations overlapped their respective estimate from the grid search approach. All estimates are, however, lower than the ground‐truth value (black cross‐marked point in Figure 9A‐D). The wide spread of and estimates obtained from the randomly initialized fitting approaches shows that SM3 and SM4 are associated to flat fitting landscapes. Nevertheless, and estimates obtained from the grid search approach are close to the µFA ground‐truth value when intracellular and extracellular domains are represented by single diffusion tensors (Figure 9A). When both intercomponent kurtosis (Figure 9B) and intracompartmental kurtosis (Figure 9C) (e.g., arising from intracellular domains) are considered, biases are present in µFA/ estimates of SM3 and SM4. These biases are dramatically higher when a fraction of zero‐diffusion isotropic compartments is added (Figure 9D).

Figure 9

Robustness of the standard model given different set of constraints (SMT1, SMT2, SM3, and SM4): A, Results for two Gaussian components representing intracellular ( 2.3 μm2/ms, μm2/ms, 0.7) and extracellular ( 1.7 μm2/ms, μm2/ms, 0.3) diffusivities; B, Results for Gaussian components representing WM axons ( 0.7 and with a log‐normal radius distribution characterized by m = 1.5 µm and std = 0.7µm) and a single component representing the extracellular ( 1.7 μm2/ms, μm2/ms, 0.3) diffusivity; C, Results for cylindrical compartments representing WMaxons ( 0.7 and with a log‐normal radius distribution characterized by m = 1.5 µm and std = 0.7 µm) and a single component representing the extracellular ( 1.7 μm2/ms, μm2/ms, 0.3) diffusivity. D, Results for cylindrical compartments representing WM axons ( 0.56 and with a log‐normal radius distribution characterized by m = 1.5 µm and std = 0.7µm), a single component representing the extracellular ( 1.7 μm2/ms, μm2/ms, 0.24) diffusivity, and a fraction of zero‐diffusion isotropic compartments ( 0.2). For each panel, the pairs of µFA and values obtained from the grid search sampling fit procedure are plotted by cross‐markers “+”; the pairs of µFA and values obtained from parameters initialized randomly are plotted by the dot markers]

DISCUSSION

Quantifying diffusion properties independent of mesoscopic tissue organization is of high interest in many applications since it provides a more informative representation of underlying microstructural features.28, 29, 30, 31, 32, 33, 34, 35, 36, 49, 50, 51, 52, 53, 54, 55, 56, 57 The µFA, typically obtainable from DDE acquisitions, quantifies diffusion anisotropy in the object’s eigenframe, which is (ideally) not affected by potentially deleterious orientation dispersion effects.16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 Recently, spherically averaged SM signals with different constraints have been proposed for quantifying µFA using the more commonly available SDE‐MRI without having to estimate the ODF.65, 66 This study’s purpose was focused on assessing whether µFA could be accurately estimated using the recently proposed SMT frameworks65, 66 in an ex vivo mouse brain. Additionally, we explored more generally whether 2‐component SM powder‐averaged fits could be appropriate for mapping µFA accurately.

In the long mixing time regime, DDE is widely established to report on µFA accurately without requiring prior assumptions.28, 29 Potentially confounding higher‐order term effects90 were accounted for in this study, making estimation more reliable. However, it should be noted that the long mixing time regime approximation needs to be fulfilled, and that exchange between different components is not accounted for. Still, can be (cautiously) considered as a “ground truth” for the experimental study.

Our experiments clearly showed that, in a realistic tissue such as a mouse brain, is crucially overestimated compared with (Figures 2 and 3). The source of this overestimation was revealed from the theory: it arises from any source of diffusion kurtosis (Equation 13), which “masquerades” as the nonmonoexponential decay underlying SMT1’s parameter estimation. Our simulations fully corroborated these experimental findings, and theoretical predictions showed that kurtosis leads to overestimation. The exception is diffusion in very elongated structures with a large difference between and . Whether represented by Gaussian or non‐Gaussian diffusion, is rather accurate so long as “intracellular” diffusion is the only source of signal (cf. Figures 6B and 7B); in this case, the deviations in estimates are not sufficiently large to skew . However, as soon as another component is introduced with a somewhat different diffusivity (e.g. “extracellular” diffusion), becomes unreliable (Figures 6C and 7C). Hence, SMT1 can play a role in characterizing, e.g., diffusion of intracellular metabolites, as already done61, 63 but not diffusion of ubiquitous water.

We then tested whether the more advanced framework incorporating (exactly) 2 components,66 termed SMT2, could provide better estimates of the underlying microstructure. Our experimental findings clearly showed that was significantly underestimated compared with (Figures 2 and 3). From a theoretical perspective, this has been shown to originate from the direct relationship imposed by SMT2’s tortuosity constraint: exclusively depends on the stick volume fraction (Equation 10), unlike the (Equation 8), which, as would be expected, depends also on the eigenvalues. If diffusivity is larger/smaller than that predicted by the tortuosity model, can be either overestimated or underestimated, respectively (Figure 8B). Figure 4 suggests that in our experimental data (and in agreement with human in vivo data),35, 56, 94 the diffusivity is lower than that predicted by the tortuosity model. If other SMT2 assumptions are violated, such as fixing the axial diffusivities for both intra‐ “neurite” and extra‐“neurite” (stick) domains to a single intrinsic diffusion coefficient , more than 50% error can be incurred (Figure 8C), in line with our experimental findings (Figures 2 and 3).

Another interesting question is whether powder‐averaged signals subject to the commonly used 2‐component SM fits could better correspond to the microstructural information reported by the model‐free DDE. In the vast majority of our voxels, both SM3 and SM4 (acquired with b values up to 9 ms/µm2, which aids in avoiding branch selection and stabilizes the fitting landscape) overestimated µFA. Before drawing any conclusions, it is worth pointing out that we have used powder‐averaged SM fits to keep consistent with SMT1 and SMT2 frameworks; however, such powder‐averaged signals have an a priori relatively flat fitting landscape, arising from the truncation of spherical harmonics at L = 0,56 which may make the fit more difficult.56, 57 In the future, directional SM fits will be studied and reported, although it is clear that for powder‐averaged systems,59, 60, 61 this approach will not work since the higher‐order terms are inherently zero. Nevertheless, assuming that the grid search of our powder‐averaged SM fits was relatively effective (cf. Figure 9A‐C), our results likely suggest that the powder‐averaged SM description – at least in ex vivo tissues scanned under the experimental parameters in this study – may be incomplete for neural tissues. One reasonable explanation is that 2 components simply do not capture the tissue heterogeneity sufficiently well (Figure 9D). For example, cell bodies exist in both white and gray matter72 and have nonzero volume fractions,92 and ex vivo samples might present fully restricted diffusion (a non‐zero fraction of “dot” compartments).47, 93 These features could perhaps suggest that a third component may be required for models of diffusion in tissues. Another potential confound of SM is the stick assumption, which ignores axon size distributions in WM or dendrites and astrocytic branches in GM, whose radii can be far from zero.

One perhaps alarming property of fitting powder‐averaged SDE signals is that the maps they produce may seem very agreeable (Figure 2). Without independent validation, it would be difficult to infer that there is anything wrong with such maps. This underscores the importance of method‐based vs. model‐based comparisons. In fact, much can be learned from the exercise of comparing such signals, including, perhaps, deciding on which biophysical model is correct for a given voxel; alternatively, DDE (or other method‐driven) signals can be used to constrain biophysical model fits.

CONCLUSIONS

Spherically averaged SDE signals – whether strongly constrained by parameter fixing, such as SMT1 and SMT2, or more weakly constrained such as SM3 or SM4 – produce discrepant µFA compared to their DDE counterparts in neural tissue. The origin of these discrepancies is traced to the imposed constraints and in SMT1 and SMT2, mainly to signal kurtosis arising from any source of heterogeneity, such as distributions of tensors or restricted diffusion. The 2‐component, minimally constrained powder‐averaged SM framework also produces different estimates compared to , suggesting either that the powder‐averaged approach is not robust enough, or that it is insufficient for describing the underlying microstructure sufficiently well. Metrics extracted from the powder‐averaged standard models should be taken with great caution and validated on the basis of data from histology or by further comparisons with model‐free measures obtained from nonconventional diffusion acquisition schemes.

FIGURE S1 The ex vivo mouse brain raw data and standard DTI FA maps: A.1, b‐value = 0 map displaying the WM and GM regions used to plot the signal decays in panels A.2 and A.3; A.2, raw signal decays for different SDE gradient directions and for WM and GM regions of interest; A.3, raw powder‐averaged signal decays for parallel and perpendicular DDE gradient directions and for WM and GM regions of interest; B, DTI FA maps of the three SDE data slices [Color figure can be viewed at wileyonlinelibrary.com]

Click here for additional data file.