Diffusion tensor distribution imaging.

Conventional diffusion MRI yields voxel-averaged parameters that suffer from ambiguities for heterogeneous anisotropic materials such as brain tissue. Using principles from solid-state NMR spectroscopy, we have previously introduced the shape of the diffusion encoding tensor as a separate acquisition dimension that disentangles isotropic and anisotropic contributions to the observed diffusivities, thereby allowing for unconstrained data inversion into diffusion tensor distributions with "size," "shape," and orientation dimensions. Here we combine our recent non-parametric data inversion algorithm and data acquisition protocol with an imaging pulse sequence to demonstrate spatial mapping of diffusion tensor distributions using a previously developed composite phantom with multiple isotropic and anisotropic components. We propose a compact format for visualizing two-dimensional arrays of the distributions, new scalar parameters quantifying intra-voxel heterogeneity, and a binning procedure giving maps of all relevant parameters for each of the components resolved in the multidimensional distribution space.

diffusion tensor distribution

orientation distribution function

INTRODUCTION

Diffusion MRI1, 2 has enabled non‐invasive studies of the microscopic structure of the central nervous system3, 4, 5, 6, 7, 8, 9 and other parts of the human body.10 While conventional diffusion MRI observables, such as the mean diffusivity and fractional anisotropy,11 are exquisitely sensitive to changes in tissue structure, they are difficult to relate to specific tissue properties,12 in particular for heterogeneous voxels containing multiple nerve fiber bundles with different orientations,13 varying amounts of free water,14 or mixtures of gray and white matter.15 Methods capable of dealing with tissue heterogeneity are of great interest for investigations of brain plasticity,16 cancer treatment,17 and neurodegeneration.18

Within the diffusion MRI field, it is customary to approximate the signal S(b) from heterogeneous tissues as the sum of contributions from a few discrete components,14, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43 or more generally as an integral transform according to44, 45, 46 where S 0 is the signal intensity at vanishing diffusion weighting, D is the diffusion tensor, P(D) is a continuous diffusion tensor distribution (DTD), b is the diffusion encoding tensor,47, 48, 49, 50 and the colon denotes a generalized scalar product. When considering only scalar diffusivities, Equation (1) reduces to the more familiar one‐dimensional multi‐exponential form.51, 52, 53

The concept of tissue heterogeneity should be interpreted in light of the time‐ and length‐scales defined by the measurement process.40, 54 Diffusion MRI measurements on clinical scanners are usually performed at a fixed diffusion time of about 0.1 s, during which the observed water molecules undergo translational diffusion on a length‐scale of 10 μm. Tissue environments exchanging water molecules on time‐scales much shorter than the observational one cannot be resolved and yield apparent diffusivities that include averaged effects of all structures and exchange processes that are too small and fast to be detected on the length‐ and time‐scales defined by the measurement. The weights and diffusivities of the DTD components may depend on the choice of experimental parameters on account of the influence of the observation time on the effects from exchange55, 56 and restricted diffusion,57, 58 which even on its own gives rise to non‐exponential signal decay for both closed compartments59 and interconnected porous structures.60, 61

Limiting our scope to the popular multi‐tensor paradigm14, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46 where the signal is given by Equation (1), an accurately estimated DTD can be seen as a crude but experimentally accessible characterization of intra‐voxel tissue heterogeneity. Unfortunately, Equation (1) belongs to a class of integral equations that are notoriously difficult to invert,62, 63, 64, 65 giving rise to non‐unique results where multiple distributions are consistent with the same experimental data. The inversion stability can be improved by imposing constraints on the obtained distributions, e.g. by fixing the diffusion tensor eigenvalues and solving for the orientation distribution function (ODF)44, 66 or estimating the weights of a few components with pre‐defined or heavily constrained diffusivities.37, 67 Although the constraints facilitate data inversion, they often rely on unverified assumptions, which may give rise to erroneous conclusions whenever the true microstructure is different from the preconceptions.39, 67

The problem of estimating DTDs is further exacerbated by the fact that most data acquisition protocols are based on the Stejskal‐Tanner pulse sequence,68 which is only capable of generating linear b‐tensors and sampling a limited region of the available acquisition space where the information about the isotropic, anisotropic, and orientational properties of the diffusion tensors is inextricably entangled.50 As a remedy to this problem, we have recently developed a set of multidimensional diffusion MRI methods69 using principles from multidimensional solid‐state NMR spectroscopy.70 Our methods rely on time‐varying amplitudes and orientations of the magnetic field gradient to produce diffusion encoding as a function of not only the conventional b‐value and b‐vector, but also the shape of the b‐tensor.48, 49 Assuming that the component diffusion tensors are axisymmetric, such detailed data can be converted to DTDs without imposing further constraints on the individual eigenvalues or number of components.69 In a recent experimental demonstration of the approach,71 we were able to resolve two isotropic components and one anisotropic component in a composite phantom mimicking the Stanisz model of nerve tissue.72 This study was limited to spectroscopic signal read‐out and did not provide directional or spatially resolved information.

Here, we merge our previously developed spectroscopic DTD acquisition protocol,71 imaging pulse sequence with single‐shot signal read‐out,73 and data inversion algorithm that retains the directional information.69 Proof‐of‐principle measurements are made on microimaging hardware using a phantom with multiple isotropic and anisotropic components.71, 73 In order to visualize the rich information of the spatially resolved multidimensional DTDs, we propose a compact 2D graph format and new scalar metrics for quantifying intra‐voxel heterogeneity. The pulse sequence has been used on whole‐body scanners in several recent publications,67, 73, 74, 75, 76, 77, 78 and we consider the present contribution as a crucial step towards implementing DTD imaging in vivo.

THEORY

A comprehensive overview of the theory for multidimensional diffusion MRI is given in a recent book chapter50 and review article.69 Below follows a succinct summary of the relevant theory for the present paper.

Diffusion encoding with axisymmetric b‐tensors

An axisymmetric diffusion tensor D can be parameterized with the axial and radial diffusivities, D || and D ⊥, as well as the polar and azimuthal angles, θ and ϕ, giving the orientation of its symmetry axis within the laboratory frame of reference. The isotropic diffusivity D iso and normalized diffusion anisotropy D Δ are given by49 respectively. Correspondingly, the diffusion encoding tensor b can be parameterized with the axial and radial eigenvalues, b || and b ⊥, and its orientation (Θ, Φ), which yield the conventional b‐value and normalized anisotropy b Δ through where b Δ is the “shape” of the diffusion encoding with special values b Δ = −1/2, 0, and + 1 for planar, spherical, and linear diffusion encoding.49 Using Equations (2) and (3), the tensor scalar product in Equation (1) can be expressed as69 where and P 2(x) = (3x 2–1)/2 is the second Legendre polynomial. Equation (4) clearly shows that the signal obtained with conventional linear diffusion encoding (b Δ = +1) depends on D iso, D Δ, and (θ, ϕ), while spherical diffusion encoding (b Δ = 0) isolates the contribution from D iso. Signal acquired in the 4D space [b, b Δ, Θ, Φ] enables estimation of DTDs expressed in the corresponding 4D space [D iso, D Δ, θ, ϕ].69

Data inversion

For the purpose of data analysis, Equation (1) is rewritten in discrete form as where S is the signal for the mth measurement with b‐tensor b and w is the weight of the nth microscopic diffusion tensor D . Equation (6) can be expressed in matrix form as where s is a M × 1 (rows × columns) vector with signal amplitudes for M diffusion encoding tensors, w is a N × 1 vector with weights for N discrete microscopic diffusion tensors, and K is an M × N matrix with elements given by the exponential in Equation (6).

Methods for solving ill‐posed linear equation systems resembling Equation (7) for the unknown vector w have been investigated in detail in the literature.62, 63, 64, 65, 79, 80, 81, 82 Here we use an iterative approach based on the inversion algorithm introduced in Reference 71 and described in detail in Reference 83. A set of D is generated by randomly selecting sampling points in the [log10(D ||), log10(D ⊥), cos(θ), ϕ]‐space, the only constraint being that the values of D || and D ⊥ are within a range that can be probed with the used b‐values. Weights w for the members of the set are obtained by inverting Equation (7) using the non‐negative least squares algorithm.79 All points with non‐zero weight survive to additional rounds of inversions, where they compete with newly generated random points. After a pre‐defined number of rounds, the surviving points are split and subjected to minor random mutations. The offspring points then compete for an additional pre‐defined number of rounds, and the values of w and D of the final survivors are taken as one possible solution to the problem. The inversion uncertainty is estimated by a bootstrapping procedure generating an ensemble of possible solutions.71, 83, 84, 85

DTD visualization

The symmetric 3 × 3 matrices of the diffusion tensors are often visualized as ellipsoid86 or superquadric tensor glyphs,87 where the lengths and orientations of the glyph semi‐axes correspond to the eigenvalues and eigenvectors of the matrices. Within such a geometric interpretation, it is intuitive to report a distribution of axisymmetric diffusion tensors in a 4D space of “sizes” D iso, “shapes” D Δ, and orientations (θ, ϕ), as explained in Figure 1.50 The shape parameter D Δ derives from the conventions of solid‐state NMR spectroscopy88 and leads to compact mathematical expressions as in Equation (4). Diffusivities in biological materials often vary over several orders of magnitude and the corresponding distributions are more clearly visualized on a logarithmic scale. Consequently, we choose to express the tensor shape as log10(D ||/D ⊥) in graphs while reserving the D Δ notation for equations. Mapping between the two representations is straightforward.50 The directional color‐coding is given by which returns black for spherical tensors where the orientation is undefined. The scatter plot in Figure 1B is sufficiently compact to be plotted for each voxel of a 2D image with moderate matrix size.

Figure 1

Illustration of a DTD. A, schematic voxel with an ensemble of superquadric tensor glyphs representing axisymmetric microscopic diffusion tensors D with different “sizes” D iso, “shapes” D Δ, and orientations (θ, ϕ). The parameters D iso and D Δ are defined from the axial and radial diffusivities, D || and D ⊥, in equation (2). B, C, bubble charts reporting the voxel composition as 4D size‐shape‐orientation distributions. The horizontal and vertical axes show the size and shape dimensions while the circle areas correspond to the statistical weight w of each component n. The orientational color‐coding is given by equation (8). B and C show different representations of the same information: The [log10(D iso/m2 s−1), log10(D ||/D ⊥)] style in B is used for plotting DTDs in Figure 3 and Figure 4, while the [D iso, D Δ 2] version in C is more directly related to the expectation value E[x], variance Var[x], and covariance Cov[x, y] metrics defined in Equations (9)−(13) and shown as parameter maps and histograms in Figure 5 and Figure 6

Figure 3

Signal data and DTDs for three representative voxels of the composite phantom comprising a 5 mm glass tube with a randomly oriented liquid crystal inserted into a 10 mm tube containing an aqueous polymer solution. A, acquisition protocol with magnitude b, normalized anisotropy b Δ, and orientation (Θ, Φ) of the b‐tensor as a function of data point index. B, S 0‐map with circles indicating three voxels with liquid crystal (green), liquid crystal + polymer solution (red), and polymer solution (blue). C, experimental (colored circles) and fitted (points) normalized signals S/S 0 versus index. All graphs in panels A and C are drawn with the same horizontal axis. D, DTDs obtained by unconstrained Monte Carlo inversion of equation (7) for the signal data in C. The DTDs are visualized as bubble charts according to the convention in Figure 1B

Figure 4

DTD imaging of the liquid crystal/polymer solution phantom. A, Array of DTD bubble charts for every 0.6 × 0.6 × 1.0 mm3 voxel of the 16 × 16 × 1 image. The boundaries between voxels (squares) coincide with the axes of the bubble charts, which have the same limits as in Figure 3D. B, global DTD obtained by superimposing the spatially resolved DTD bubble charts in panel A. The numbered rectangles (gray) mark the bins for calculating component‐resolved parameters (1, “sticks,” water in liquid crystal; 2, “small spheres,” polymer; 3, “large spheres,” water in polymer solution; 4, “planes,” artifact from data inversion)

Figure 5

Color‐coded parameter maps derived from the spatially resolved DTDs. Columns: Initial signal intensity S 0, expectation value E[x], variance Var[x], and covariance Cov[x, y] of the isotropic diffusivity D iso and the square normalized diffusion anisotropy D Δ 2. Rows: Parameter calculations from the total and four bins of the DTD space as indicated in Figure 4B. The numerical values of the linear color scales are identical to the horizontal axes of the histograms in Figure 6

Figure 6

Uncertainty estimation of the DTD‐derived parameters for the liquid crystal (green), liquid crystal + polymer solution (red), and polymer solution (blue) voxels indicated in Figure 3B. The histograms represent the set of parameters calculated from the 96 individual realizations of the DTDs in the bootstrapping procedure. The matrix of histograms is ordered in the same way as the parameter maps in Figure 5, and the horizontal axes of the histograms are the same as the linear color scale of the maps

Statistical measures

In order to display the rich information contained within the DTDs for larger image matrices, it is useful to extract statistical measures such as expectation values E[x], variances Var[x], and covariances Cov[x, y] over all relevant dimensions of the distribution. Here we use the [D iso, D Δ 2] space shown in Figure 1C to define and where is the initial signal. The parameter E[D iso] is identical to the conventional mean diffusivity and E[D Δ 2] carries similar information to other anisotropy measures from the literature, namely the microscopic anisotropy,89 microscopic fractional anisotropy,73, 90, 91 and fractional eccentricity.92 The parameters Var[D iso], Var[D Δ 2], and Cov[D iso, D Δ 2] report on various aspects of intra‐voxel heterogeneity. Of these three measures, only the first one has previously appeared in the literature under the symbols μ 2 iso,73, 90 V I,74, 76 V MD,75 and V(D iso),50, 69 and has been shown to be correlated with cell density heterogeneity in brain tumors.76

Bin‐resolved parameters

The measures defined in Equations (9)−(14) can be separated into contributions from different components by defining bins in the DTD space. The [D iso, D ||/D ⊥] representation of the size‐shape dimensions of the DTD in Figure 1B makes it natural to express the boundaries of the bins as

This subdivision of the distribution space is conceptually analogous to the binning of 1D T 2‐distributions for estimating the fraction of bound fluid in porous media84 or myelin water in brain tissue.93

MATERIALS AND METHODS

A phantom with separate polymer solution and liquid crystal compartments was assembled from 5 and 10 mm NMR tubes as described by Lasič et al.73 The aqueous polymer solution was prepared by adding H2O (Milli‐Q quality) and the polymer poly(ethylene glycol), with a molecular weight of 400 g/mol, to D2O (99.8 mol% 2H, Armar Chemicals, Döttingen, Switzerland), giving a final concentration of 10 wt% for both H2O and polymer. The liquid crystal was of the reverse hexagonal type,94, 95 comprising water channels with diameters on the nanometer scale in a continuous matrix of detergent and hydrocarbon. The liquid crystal was prepared as described in Reference 78 using the chemical composition 43 wt% H2O, 5 wt% D2O, 38 wt% sodium 1,4‐bis(2‐ethylhexoxy)‐1,4‐dioxobutane‐2‐sulfonate, and 14 wt% 2,2,4‐trimethylpentane, the two latter chemicals being purchased from Sigma‐Aldrich, Stockholm, Sweden. The nuclear relaxation properties of the components of the liquid crystal are such that only water contributes to the detected signal.78 The orientations of the liquid crystalline domains were randomized by melting the sample at 50°C and subsequently reforming the liquid crystal at 15°C, giving microscopic crystallites with dimensions larger than the approximately 10 μm displacements during the diffusion encoding gradients and smaller than the approximately 1 mm voxel size.96 A liquid crystal with domains aligned in a single direction was prepared by exposing the phantom to an 11.7 T magnetic field for seven days.

MRI measurements were made at 22°C on a Bruker Avance‐II 500 MHz spectrometer equipped with an 11.7 T standard bore magnet and a MIC‐5 triple‐axis gradient probe with a 10 mm RF coil insert (Bruker, Karlsruhe, Germany). Images were acquired with the pulse sequence of Lasič et al.73 (see Figure 2), consisting of spin echo signal preparation and single‐shot RARE image read‐out,97 with 35 ms echo train duration, 16 × 16 matrix size, 0.6 × 0.6 mm2 in‐plane spatial resolution, 1 mm slice thickness, and 10 s repetition time. Tensor‐valued diffusion encoding48, 49, 69 was implemented with a pair of τ = 10 ms modulated gradient waveforms bracketing the 180° pulse in the TE = 35.3 ms spin echo. The waveforms were calculated as described in detail in Reference 96 using 3 T m−1 maximum gradient amplitude. The acquisition protocol comprised an array of eight b‐values (b = 0.022, 0.054, 0.132, 0.319, 0.776, 1.886, 4.583, and 11.136·109 s m−2) and four b‐tensor anisotropies (b Δ = −0.5, 0.0, 0.5, 1.0). For each combination of (b, b Δ), images were acquired at 15 gradient waveform rotations given by the electrostatic repulsion scheme,98 resulting in a total number of 480 images and 80 min scan time. A graphical representation of the acquisition protocol is displayed in Figure 3A.

Figure 2

Pulse sequence for DTD imaging. RARE image read‐out97 is preceded by a spin echo block (90° and 180° RF pulses) of length TE. The 180° pulse is bracketed by a pair of modulated gradient waveforms (red, green, and blue), each having duration τ. Waveforms calculated with the equations in reference 96 yield diffusion weighting as a function of the magnitude b, normalized anisotropy b Δ, and orientation (Θ, Φ) of the diffusion encoding tensor b. (adapted from Topgaard96 with permission)

The acquired k‐space data was converted to images using 0.9 × 0.9 mm2 in‐plane Gaussian smoothing. Per‐voxel DTDs were estimated with the algorithm described in Section 2.2 using 200 randomly selected points in the space [−11 < log10(D ||/m2 s−1) < −8.3, −11 < log10(D ⊥/m2 s−1) < −8.3, −1 < cos(θ) < 1, 0 < ϕ < 2π], 20 rounds of generating new random points, and 20 rounds of mutations. Bootstrapping with replacement was used to generate 96 different DTDs consistent with the data. For each DTD, scalar parameters were calculated with Equations (9)−(14) for the total distribution as well as four bins with the following limits: (i) “sticks,” −11 < log10(D iso/m2 s−1) < −8.3 and 0.5 < log10(D ||/D ⊥) < 3.3; (ii) “small spheres,” −11 < log10(D iso/m2 s−1) < −9.5 and − 0.5 < log10(D ||/D ⊥) < 0.5; (iii) “large spheres,” −9.5 < log10(D iso/m2 s−1) < −8.3 and − 0.5 < log10(D ||/D ⊥) < 0.5; and (iv) “planes,” −11 < log10(D iso/m2 s−1) < −8.3 and − 3.3 < log10(D ||/D ⊥) < −0.5. These bins are expected to contain the following: (i) water in the liquid crystal; (ii) polymer; (iii) water in the polymer solution; and (iv) spurious planar components, which sometimes appear in the unconstrained data inversion.83 The directional information in Bin 1 was converted to an ODF by mapping the discrete components onto a 3994‐node spherical mesh using a 10° Gaussian smoothing kernel.

Bruker TopSpin pulse sequence and gradient waveforms, MATLAB code for image reconstruction and data inversion, and experimental data are available at https://github.com/daniel‐topgaard/.

RESULTS AND DISCUSSION

Experimental results from three representative voxels are shown in Figure 3. The selected voxels contain liquid crystal in the inner 5 mm tube and aqueous polymer solution in the outer 10 mm tube, as well as a “mixed” voxel with contributions from both the liquid crystal and the polymer solution. Visual inspection of the signal data reveals distinct differences between the voxels: the liquid crystal gives rise to a pattern sensitive to mainly the acquisition variables b and b Δ, and only some minor dependence on Θ and Φ, consistent with anisotropic diffusion in orientationally disordered microdomains; the polymer solution data is insensitive to b Δ, Θ, and Φ, indicating isotropic diffusion; and the mixed voxel appears to be a superposition of the two other voxels with a major contribution from the polymer solution. These observations of the raw signal data are consistent with the DTDs obtained by unconstrained inversion: the liquid crystal gives rise to a DTD having components with log10(D iso/m2 s−1) ≈ −9.4, log10(D ||/D ⊥) ≈ 2.0, and multiple orientations (θ, ϕ); the polymer solution yields isotropic components at log10(D iso/m2 s−1) ≈ −9.8 and −9.0, corresponding to the polymer and the water, respectively; and the mixed voxel gives a superposition of the liquid crystal and polymer solution results. The liquid crystal, polymer, and water components are expected to have different values of the relaxation rate constants R 1 and R 2, which are not quantified or compensated for in the current acquisition protocol and inversion procedure. Consequently, the relative weights of the DTD components cannot be directly related to the known chemical compositions. As discussed in Reference 78, the signal from the liquid crystal originates only from water since the detergent and hydrocarbon components have too high values of R 2 to be detected at the echo times of about 50 ms used herein.

Figure 4 compiles results for all voxels as both a square grid and a superposition of the per‐voxel DTDs, showing that most voxels are similar to the representative ones in Figure 3. The liquid crystal water component is located in voxels corresponding to the interior of the 5 mm tube, while the water and polymer components of the aqueous polymer solution can be found in the area between the walls of the 5 and 10 mm tubes. The chemical compositions are homogeneous for both the liquid crystal and the polymer solution, and the overlay of all spatially resolved DTDs in Figure 4B would ideally comprise delta‐functions at the three discrete components. The global DTD features three dense clusters of points at the expected components as well as a sparse cloud of spurious low‐weight points that we attribute to experimental noise and the stochastic nature of the unconstrained data inversion procedure. Since the chemical compositions and microstructures are the same throughout each compartment of the phantom, the spread of points within each cluster reflects the measurement and inversion uncertainty rather than spatial variation of the component diffusivities. Any real spread of diffusivities within each component, which sometime occurs for polymers with a distribution of molecular weights,99 would be at least partially masked by the measurement uncertainty. Taking a closer look at the boundary between the liquid crystal and polymer solution in Figure 4A, it is clear that the expected three‐component DTDs are noisier than the one‐ and two‐component ones in the neighboring voxels.

The rich information in the DTDs can be visualized in a more compact and noise‐resistant way using the parameters defined in Equations (9)−(14). While the S 0 and E[D iso] parameters are identical to ones obtained with conventional diffusion MRI, the E[D Δ 2], Var[D iso], Var[D Δ 2], and Cov[D iso, D Δ 2] parameters can only be determined with tensor‐valued diffusion encoding and give easily interpretable information about intra‐voxel microstructure and heterogeneity. First focusing on the parameters calculated from the total DTD space (see Row 1 in Figure 5), we note that E[D Δ 2] and Var[D iso] clearly distinguish between locally anisotropic and multicomponent isotropic diffusion as first demonstrated by Lasič et al.73 The maps of E[D Δ 2] and Var[D iso] correspond to the locations of the liquid crystal and polymer solution, respectively. Heterogeneous voxels having overlapping contributions from isotropic and anisotropic components are captured with the novel Var[D Δ 2] parameter which gives high values in a circular region at the boundary between the liquid crystal and polymer solution. The Cov[D iso, D Δ 2] parameter carries unique information about the correlation between the isotropic and anisotropic dimensions. Even though the individual DTDs appear noisy for the boundary region in Figure 4A, the values of Cov[D iso, D Δ 2] are finite and negative, showing that there is a correlation between high values of D iso and low anisotropy.

Moving on to Rows 2–5 in Figure 5 with parameter maps calculated from the bins in the DTD space, we note that the heterogeneity measures Var[D iso], Var[D Δ 2], and Cov[D iso, D Δ 2] are nearly zero, indicating that the chosen bins are indeed capable of separating the three expected components. Bin 1 contains the “stick” component from the liquid crystal in the inner tube, while Bins 2 and 3 correspond to the slow and fast isotropic components from the polymer and water in the outer tube. The spurious “plane” component in Bin 4 is only visible as a faint band at the border between the tubes. Finite values in the Var[D iso], Var[D Δ 2], and Cov[D iso, D Δ 2] maps occur only for voxels and bins with low S 0 and are therefore attributed to noise rather than to any true heterogeneity of the material.

The DTD uncertainty was estimated with a bootstrapping procedure generating multiple solutions consistent with the data. For every voxel, the range of DTDs obtained with bootstrapping is similar to the spread of solutions across the different voxels with identical chemical composition in Figure 4A. The uncertainty of the scalar parameters is reported as histograms in Figure 6 for the three representative voxels labeled in Figure 3B. From the widths of the peaks in the histogram, we note that there is a general trend of increasing uncertainty in the order S 0, E[D iso], E[D Δ 2], Var[D iso], Var[D Δ 2], and Cov[D iso, D Δ 2], as well as a correlation between high values of S 0 and low uncertainty. The widths of the clusters of points in Figure 4B depend on both the “true” widths Var[x] and the uncertainties in E[x] and Var[x] for each of the bins. The narrow histogram peaks for in particular E[D iso] in Figure 6 illustrate that this parameter can be determined with fairly high precision even though the DTD scatter plots in Figure 4 appear noisy. The mixed voxel (red) with an expected three‐component DTD yields higher uncertainty than the simpler DTDs from the liquid crystal (green) and polymer solution (blue). The very high uncertainty of the Var[D iso] and Var[D Δ 2] parameters for the Bin 1 mixed voxel confirms that the finite values in the corresponding maps in Figure 5 should be interpreted as resulting from noise rather than from heterogeneity.

The directional information of the Bin 1 “stick” component is visualized as arrays of ODF glyphs in Figure 7, showing that the liquid crystal with domain orientations randomized by a temperature cycle indeed has a wide range of orientations in every voxel. Conversely, the liquid crystal exposed to an 11.7 T magnetic field for 7 d features domain alignment in the direction of the field. With the exception of the ODFs, the randomized and aligned samples give nearly identical results in the size and shape aspects of the DTDs, as illustrated in the Supplementary Material with figures corresponding to Figures 3, 4, 5, 6 for the aligned sample.

Figure 7

Voxel‐resolved ODFs for the “stick” bin 1 in Figure 4B. The ODFs are displayed as directionally color‐coded surfaces (RGB = xyz) with radius scaled according to the values of the ODF and S 0 for bin 1. A, randomized domain orientations. B, domains aligned by 7 d immersion in an 11.7 T magnetic field in the z‐direction. The black squares in the upper panels outline the 8.6 × 8.6 mm2 field of view in the xy‐plane. The lower panels show 4× magnifications of the upper ones

The present acquisition protocol for microimaging applications has recently been implemented and validated with phantom measurements78 on the whole‐body Connectome scanner100 at the Athinoula A. Martinos Center at Massachusetts General Hospital. The two protocols include the same number of acquired images and cover similar ranges of b‐values,78 indicating that there is potential for non‐parametric DTD imaging also for human in vivo applications. The current data inversion does not include the effects of time‐dependent diffusion57 which have been observed for human in vivo using oscillating gradients101, 102 and stimulated echo pulse sequences103 allowing for the extended ranges of diffusion times necessary to observe the often subtle effects. Future attempts to apply DTD imaging in vivo will thus require dedicated measurements to test if effects of time‐dependent diffusion are observed within the narrow window of diffusion times, approximately 100–150 ms, afforded by our sequence with spin echo signal preparation and single‐shot EPI read‐out.73 Such tests could be performed by investigating the signal modulation as a function of the spectral content104 or the direction of gradient waveforms giving spherical b‐tensors.105 If the effects can be detected on a voxel‐by‐voxel basis within the limitations of our sequence, we expect that the time dependence could be incorporated into our analysis as an additional measurement dimension carrying information about compartment sizes and local diffusivities,106 for instance by using the confinement tensor approach.107

CONCLUSIONS

We have demonstrated that imaging of 4D size‐shape‐orientation DTDs can be achieved by unconstrained inversion of signal data acquired in the corresponding size‐shape‐orientation space of the b‐tensor. The data inversion algorithm recovered all expected DTD components and gave reproducible results across voxels having the same chemical composition and microstructure. Conversion of the rich multidimensional information of the DTDs into a compact 2D format and novel scalar parameters allowed simultaneous visualization of the results for an entire imaging slice. The recent validation of our acquisition protocol on a whole‐body scanner78 is encouraging for future attempts to apply DTD imaging in vivo.

Supporting Figure 1. Signal data and DTDs for three representative voxels of the composite phantom with a 5 mm glass tube with a magnetically aligned liquid crystal inserted into a 10 mm tube with an aqueous polymer solution. (a) Acquisition protocol with magnitude b, normalized anisotropy , and orientation (Θ, Φ) of the b‐tensor as a function of data point index. (b) S 0‐map with circles indicating three voxels with liquid crystal (green), liquid crystal + polymer solution (red), and polymer solution (blue). (c) Experimental (colored circles) and fitted (points) normalized signals S/S 0 vs. acquisition index. All graphs in panels (a) and (c) are drawn with the same horizontal axis (d) DTDs obtained by unconstrained Monte Carlo inversion of main article Equation 8 for the signal data in panel (c). The DTDs are visualized as color‐coded scatter plots according to the convention in main article Figure 1.

Supporting Figure 2. DTD imaging of the magnetically aligned liquid crystal/polymer solution phantom. (a) Array of DTD scatter plots for every 0.6 × 0.6 × 1.0 mm3 voxel of the 16 × 16 × 1 image. The boundaries between voxels (squares) coincide with the axes of the scatter plots which have the same limits as in Figure 1(d). (b) Global DTD obtained by superimposing the spatially resolved DTD scatter plots in panel (a). The numbered rectangles (gray) mark the bins for calculating component‐resolved parameters (1 “sticks”: water in liquid crystal, 2 “small spheres”: polymer, 3 “big spheres”: water in polymer solution, 4 “planes”: artifact from data inversion).

Supporting Figure 3. Color‐coded parameter maps derived from the spatially resolved DTDs. Columns: initial signal intensity S 0, expectation value E[x], variance Var[x], and covariance Cov[x,y] of the isotropic diffusivity D iso and the square normalized diffusion anisotropy . Rows: parameter calculation from the total and four bins of the DTD space as indicated in Figure 1(b). The numerical values of the linear color scales are identical to the horizontal axes of the histograms in Figure 4.

Supporting Figure 4. Uncertainty estimation of the DTD‐derived parameters for the magnetically aligned liquid crystal (green), liquid crystal + polymer solution (red), and polymer solution (blue) voxels indicated in Figure 1(b). The histograms represent the set of parameters calculated from the 96 individual realizations of the DTDs in the bootstrapping procedure. The matrix of histograms is ordered in the same way as the parameter maps in Figure 3, and the horizontal axes of the histograms are the same as the linear color scale of the maps.

Click here for additional data file.