Estimation of Tensors and Tensor-Derived Measures in Diffusional Kurtosis Imaging

doi:10.1002/mrm.22655

Magn Reson Med. Author manuscript; available in PMC 2012 March 20.

Published in final edited form as:

Magn Reson Med. 2011 March; 65(3): 823–836.

Published online 2010 October 28. doi: 10.1002/mrm.22655

PMCID: PMC3042509

NIHMSID: NIHMS247734

Estimation of Tensors and Tensor-Derived Measures in Diffusional Kurtosis Imaging¹

Ali Tabesh,^a² Jens H. Jensen,^a^b Babak A. Ardekani,^c^d and Joseph A. Helpern^a^b^c^d

^aCenter for Biomedical Imaging, Department of Radiology, New York University School of Medicine, 660 First Avenue New York, NY 10016 USA

^bDepartment of Physiology and Neuroscience, New York University School of Medicine, 550 First Ave, New York, NY 10016 USA

^cCenter for Advanced Brain Imaging, The Nathan S. Kline Institute for Psychiatric Research, 140 Old Orangeburg Rd, Orangeburg, NY 10962 USA

^dDepartment of Psychiatry, New York University School of Medicine, 550 First Ave, New York, NY 10016 USA

² Corresponding author. Phone: +1-212-263-3344; fax: +1-212-263-4830; Email: ali.tabesh/at/nyumc.org

The publisher's final edited version of this article is available at Magn Reson Med

See other articles in PMC that cite the published article.

Abstract

This paper presents two related advancements to the diffusional kurtosis imaging (DKI) estimation framework to increase its robustness to noise, motion, and imaging artifacts. The first advancement substantially improves the estimation of diffusion and kurtosis tensors parameterizing the DKI model. Rather than utilizing conventional unconstrained least squares (LS) methods, the tensor estimation problem is formulated as linearly constrained linear LS, where the constraints ensure physically and/or biologically plausible tensor estimates. The exact solution to the constrained problem is found via convex quadratic programming methods or, alternatively, an approximate solution is determined through a fast heuristic algorithm. The computationally more demanding quadratic programming-based method is more flexible, allowing for an arbitrary number of diffusion weightings and different gradient sets for each diffusion weighting. The heuristic algorithm is suitable for real-time settings such as on clinical scanners, where run time is crucial. The advantage offered by the proposed constrained algorithms is demonstrated using in vivo human brain images. The proposed constrained methods allow for shorter scan times and/or higher spatial resolution for a given fidelity of the DKI parametric maps. The second advancement increases the efficiency and accuracy of the estimation of mean and radial kurtoses by applying exact closed-form formulae.

Keywords: MRI, diffusion, kurtosis, non-Gaussian, brain

Introduction

Diffusion of water molecules in biological tissues is conventionally quantified via diffusion tensor imaging (DTI) (1). DTI is a valuable tool for noninvasive characterization of tissue microstructural properties (2–3). The DTI model uses a Gaussian approximation to the probability distribution governing the random displacement of water molecules. In many biological tissues, however, the displacement probability distribution can deviate considerably from a Gaussian form. Several techniques have been proposed to characterize non-Gaussian diffusion, with diffusion spectrum imaging (DSI) (4) being arguably the most comprehensive. Although it is being used in human research, the application of DSI as a routine clinical protocol has been limited by its extended acquisition time and hardware requirements.

Diffusional kurtosis imaging (DKI) is a clinically feasible extension of DTI which enables the characterization of non-Gaussian diffusion by estimating the kurtosis of the displacement distribution (5–8), in addition to the estimation of the standard DTI-derived parameters. DKI has shown promising results in studies of human brain aging (9), tumor characterization in gliomas (10) and head and neck cancers (11), and rodent brain maturation (12). Moreover, the additional information provided by DKI has been exploited to resolve intravoxel fiber crossings (13), which could help to improve upon DTI-based fiber tracking methods.

The DKI model is parameterized by the diffusion tensor (DT) and kurtosis tensor (KT) from which several rotationally-invariant scalar measures are extracted. The most common DT-derived measures are mean, axial, and radial diffusivity (2,14), as well as fractional anisotropy (FA) (2); and the KT-derived measures are axial, radial, and mean kurtoses (6–7,15). The interpretability of these metrics is influenced by the estimation accuracy of the tensors. Noise, motion, and imaging artifacts can introduce errors into the estimated tensors. Sufficiently large errors can cause the tensor estimates to be physically and/or biologically implausible. For instance, the directional diffusivities may become negative, that is, the DT may become non-positive definite (NPD). A well-known consequence of NPD DT estimates is that FA values, which in theory should range between 0 and 1, may exceed 1, particularly in high FA regions of the brain such as the corpus callosum (16). Inaccuracies in the estimated tensors may also drive the directional kurtoses outside of an acceptable range. Empirical evidence in the brain as well as idealized multi-compartment diffusion models suggest that directional kurtoses should typically be positive and should not exceed a certain level depending on tissue complexity (6). The maximum allowable kurtosis is also influenced by the maximum b-value used in image acquisition (8,13).

In our previous work, the DT and KT were estimated using unconstrained nonlinear least squares (UNLS) (7). In a related work, higher-order DTs were estimated using unconstrained linear least squares (ULLS) (17). These unconstrained schemes do not guarantee acceptable tensor estimates. To address this drawback in the context of DTI, the Cholesky decomposition has been utilized to impose the non-negative definiteness constraint on the DT (16,18) using either the ULLS or UNLS algorithms. Moreover, a parameterization has been proposed in a UNLS framework to guarantee a positive diffusivity function in a fourth-order tensor-only model of the diffusion signal (19). While these methods outperform other algorithms for imposing the positive diffusivity function in the context of a second order-only and a fourth order-only diffusion signal model, their generalization to impose the constraints on the DKI signal model is not straightforward.

We propose that the tensor estimation problem be cast as linear least squares (LS) subject to linear constraints. The constraints ensure that the directional diffusivities and kurtoses along the imaged gradient directions remain within a physically and biologically plausible range. The proposed constrained linear LS (CLLS) formulation yields a convex objective function that permits efficient solutions via convex quadratic programming or a fast heuristic algorithm. The two algorithms proposed for solving the CLLS problem strike different tradeoffs between the speed and exactness of the solution as well as algorithm flexibility. The quadratic programming-based (CLLS-QP) algorithm exactly satisfies the constraints and can also handle an arbitrary number of diffusion weightings and different gradient sets for each diffusion weighting, but it does so at a moderately high computational cost. On the contrary, the heuristic (CLLS-H) algorithm produces an approximation to the optimal solution, but it accomplishes this at almost no computational overhead compared to the ULLS. Both constrained algorithms substantially improve the estimation fidelity of the tensors.

Once the tensors are estimated, they are utilized to determine the scalar diffusion and kurtosis measures. In our previous work, mean kurtosis (MK) was being approximated as the average of directional kurtoses (5–7). Here, we present an exact analytical formula for the MK, expressing it as an explicit function of the DT and KT. We also present an exact analytical expression for the radial kurtosis, defining it as the diffusional kurtosis averaged over all directions perpendicular to the DT eigenvector with the largest eigenvalue. Our definition is consistent with the definition of radial diffusivity in DTI and the definitions given in (8,20), but is different from that of (15). The analytic formulae for the MK and radial kurtosis were previously quoted in (8) but here for the first time we present their derivations and discuss important practical details of their implementation. The application of these formulae improves both the accuracy and efficiency of DKI parameter estimation.

Theory

Constrained LS Formulation

For direction vector n and diffusion weighting b, the DKI approximation to diffusion signal intensity S(n,b) is given by (6)

[1]

where S₀ is the signal intensity for b = 0, D_ij and W_ijkl are the elements of the second-order DT D and fourth-order kurtosis tensor W, respectively, and = (1/3)tr(D) is mean diffusivity, where tr(.) denotes matrix trace. For the classic Stejskal-Tanner sequence (21), the b-value is given by b = (γδg)²(Δ − δ/3), where Δ is the diffusion time, δ is the pulse duration, g is the diffusion-sensitizing gradient amplitude, and γ is the nuclear spin gyromagnetic ratio. Diffusivity D(n) and kurtosis K(n) along direction n are given by

[2]

and

[3]

Both D and W are symmetric, that is, they are invariant to permutations of their indices. Consequently, D has 6 and W has 15 independent parameters.

To ensure that the diffusivity and kurtosis parameters are physically and biologically plausible, D(n) and K(n) must satisfy

[4]

and

[5]

and

[6]

where n

N, with N denoting the set of gradient directions used in acquisition. The constraints in Eq. [4] ensure that D satisfies the necessary conditions for non-negative definiteness. The constraints in Eq. [5] arise from the fact that biologically relevant tissue geometries impose restrictions on physically acceptable directional kurtoses (6). While the theoretically feasible minimum kurtosis of a probability distribution is −2, multi-compartment diffusion models and empirical evidence in the brain suggest a super-Gaussian displacement distribution; that is, K_min = 0 (6). Although scenarios are conceivable under which S(n,b) can be an increasing function of b (22), this has never been observed in biological tissues. Thus, it is reasonable to require that the estimated S(n,b) be strictly decreasing functions of b for all n [set membership]

N. The maximum K(n) in Eq. [6] with C = 3 ensures that this condition is satisfied in the range of b-values used for data acquisition (13). Smaller C values may be used to further constrain the directional kurtoses.

The objective in the estimation problem is to find D and W such that the right-hand side of Eq. [1] closely matches its left-hand side, while the constraints in Eqs. [4] and [5] are satisfied. The estimation problem may therefore be formulated as

[7]

where

[8]

is a 21×1 vector of unknowns, and

[9]

The change of variables in Eq. [9] is made to permit a linear LS formulation. Note that V is reminiscent of the fourth-order tensor in the formulation of (17).

Matrix A is given by

[10]

where M is the number of nonzero b-values, and row k of equation M11

and

are given by

[11]

where

. Note that each diffusion weighting b_m may correspond to a distinct gradient set equation M15

. The number of rows of A is thus given by

Vector B is given by

[12]

where

is the diffusion signal intensity for gradient direction equation M19

and diffusion weighting b_m, and

Matrix C representing the linear constraints is given by

[13]

where

[14]

and

[15]

Vector d is given by

[16]

where the zero vectors are 1 × N -dimensional.

Two comments regarding the above formulation are in order:

At least two nonzero b-values and 15 gradient directions per b-value are required for the solution of Eq. [7]. The need for two b-values can also be seen from Eq. [1], where to find the two unknowns D(n) and K(n), at least two signal measurements (obtained with two different b-values) are needed.
For the typical choice K_min = 0, Eq. [7] provides an exact formulation of the estimation problem. If K_min ≠ 0, the terms in d need to be approximated. This may be accomplished via estimating D once and plugging that estimate into Eq. [2]. The estimate of D may be obtained from the ULLS estimate of X, denoted as , given by
[17]
where A⁺ denotes the pseudoinverse of A. The terms can be further refined using the CLLS-based estimate of D.

Quadratic Programming (CLLS-QP) Algorithm

The problem defined by Eq. [7] is a convex quadratic programming problem whose solution can be determined via standard algorithms. Two common classes of algorithms for solving this problem are active set and interior-point methods (23). Here we used the active set method which works as follows. First, the algorithm is initialized with a feasible point, if one exists. A feasible point, defined as a point that satisfies the constraints in Eq. [7], may be found via linear programming. In our implementation, the ULLS solution of Eq. [17] is used first as a potential feasible point. If this point is feasible, it is also the solution to the CLLS problem. If the point is not feasible, then the algorithm uses the feasible point found via linear programming.

In the second phase, the feasible point is iteratively refined until it converges to the optimal solution. At each iteration, an estimate of the active constraints in Eq. [7] is maintained. Active constraints are a subset of the constraints that are actually enforced, that is, they are the ones that become equality constraints. As the algorithm proceeds, the set of active constraints is updated by sequentially adding the active and removing the inactive constraints from the set.

Heuristic (CLLS-H) Algorithm

Alternatively, the heuristic procedure described below may be used to estimate D and W such that the constraints of Eqs. [4] and [5] are approximately satisfied. This algorithm assumes that exactly two nonzero b-values are used, and that the gradient directions are the same for both b-values; that is, N⁽¹⁾ = N⁽²⁾ = {n₁,…,n_N}. The central idea in this algorithm is to first estimate D(n_i) and K(n_i) along individual gradient directions, and then to threshold the ones that violate the constraints of Eqs. [4] and [5] to within their prescribed limits (steps D1, K1, and K2 below). More specifically, D(n_i) and K(n_i) are first determined from S(n_i,b₁) and S(n_i,b₂) (Eqs. [18] and [19]). Then, D(n_i) are thresholded such that they satisfy the constraints on both D(n_i) and K(n_i) (steps D1-D4). Next, the constrained D(n_i) are used to estimate D (Eq. [21]), which in turn is used to re-estimate D(n_i) (Eq. [22]). The re-estimated diffusivities D^(R) (n_i) are then used to obtain the re-estimated K^(R) (n_i) (Eq. [23]). Finally, K^(R) (n_i) are thresholded (steps K1 and K2), and are then used along with D^(R) (n_i) to compute W.

Along direction n_i, D_i = D(n_i) and K_i = K(n_i) can be estimated using S(n_i,b₁) and S(n_i,b₂) as

[18]

and

[19]

where

[20]

Note that Eqs. [18] and [19] result when exactly two nonzero b-values are used.

Next, a set of rules is applied to constrain D_i:

D1.
If D_i ≤ 0, set D_i = 0.
D2.
Otherwise, if , set D_i = 0.
D3.
Otherwise, if D_i > 0 and K_i < K_min,
- If K_min = 0, set .
- Otherwise, set
D4.
Otherwise, if D_i > 0 and K_i > K _{max_i}, where K_{max_i} = K_max (n_i), set
where b_max = max(b₁,b₂). Note that the above equation results by plugging the expression for K_{max_i} = K_max(n_i) given by Eq. [6] into Eq. [1], setting b = b₁, and solving for D_i. The diffusion signal for b = b₁ is utilized here, as it is less sensitive to noise and is a stronger function of D_i.

From the constrained D_i set, D may be estimated as

[21]

where _D is the ULLS estimate of X_D = [D₁₁ D₂₂ D₃₃ D₁₂ D₁₃ D₂₃]^T, and b_D = [D₁ … D_N]^T.

The re-estimated equation M36

is then obtained as

[22]

The re-estimation step has a noise removal effect on b_D by forcing it to be consistent with _D. Next, equation M38

are used to obtain equation M39

[23]

Next, a set of rules is applied to constrain equation M41

K1.
If , set .
K2.
If , set .

Finally, the kurtosis tensor is estimated as

[24]

where _K is the ULLS estimate of X_K = [W₁₁₁₁ … W₁₁₁₂ …W₁₁₂₂ … W₁₂₃₃]^T, and equation M47

Tensor-Derived Kurtosis Measures

The MK is defined as the diffusional kurtosis averaged over all gradient directions. In a reference frame that diagonalizes D, it may be calculated as (see Appendix for the derivation)

[25]

where λ₁ > λ₂ > λ₃ denote the eigenvalues of D, and _ijkl are the elements of the kurtosis tensor in the rotated frame and are given by

[26]

with R_ij being the j-th component of the eigenvector of D corresponding to λ_i. Functions F₁ and F₂ are given by

[27]

and

[28]

where R_F and R_D are Carlson’s elliptic integrals (24) defined as

[29]

and

[30]

As there are highly efficient numerical algorithms for calculating Carlson’s elliptic integrals (24–25), Eqs. [27] and [28] may be used for fast computation of the MK. Note that when two or more of the eigenvalues of D coincide, F₁ and F₂ have removable singularities whose treatment is discussed in the Appendix.

The axial kurtosis K_║ is the diffusional kurtosis in the direction of the highest diffusion. In a reference frame that diagonalizes D, it is simply given by

[31]

The radial kurtosis K is defined as the mean diffusional kurtosis perpendicular to the direction of highest diffusion, consistent with the definition of radial diffusivity. In a reference frame that diagonalizes D, it is given by (see the Appendix for the derivation)

[32]

where

[33]

and

[34]

The treatment of the singularities of G₁ (λ₁,λ₂,λ₃) and G₂ (λ₁,λ₂,λ₃) is discussed in the Appendix.

Method

Magnetic Resonance Imaging Procedure

DKI scans were performed on three healthy volunteers using a 3 T Siemens Tim Trio system with a 12-channel head coil. This study was approved by the New York University School of Medicine Institutional Review Board and written informed consent was obtained from all subjects. Diffusion-weighted images (DWI’s) were acquired along 30 gradient directions with a twice-refocused spin-echo sequence. Scan parameters were TR = 5900 ms, TE = 96 ms, matrix = 82 × 82, FOV = 222 × 222 mm², 45 slices, slice thickness = 2.7 mm with no gap, partial Fourier encoding = 3/4, NEX = 44 for b = 0 (NEX = 45 for Dataset 2), and NEX = 4 for b = 500, 1000, 1500, and 2000 s/mm² (NEX = 5 for Dataset 2). The total acquisition time was 53:20 minutes (65:32 minutes for Dataset 2). We emphasize that standard DKI scans (discussed below) require 7:29 minutes.

Processing Pipeline

In the first stage of processing, the DWI’s are optionally smoothed using a Gaussian kernel to reduce the impact of noise and misregistration. Next, D and W are estimated using the CLLS-QP algorithm, and finally, the scalar measures are obtained from the tensors. The pipeline was implemented in-house in the MATLAB environment (http://www.mathworks.com). The program, designated as Diffusional Kurtosis Estimator (DKE), is available upon request from the corresponding author.

The parameters used with DKE were as follows: the full width at half maximum of the Gaussian kernel was 3.375 mm; K_min = 0; and C = 3. The tensors were also estimated using the UNLS, ULLS and CLLS-H algorithms following smoothing of DWI’s. Three processing protocols were considered where different subsets of the acquired DWI’s were used to estimate the parametric maps. In the reference protocol, the maps were obtained using all of the available DWI’s. In the standard protocol, the maps were estimated using NEX = 11 for b = 0 images, and NEX = 1 for b = 1000, 2000 s/mm² DWI’s. Finally, in the fast protocol, the maps were obtained with NEX = 11 for b = 0 images, and NEX = 1 with 15 of b = 1000 s/mm² and all of b = 2000 s/mm² DWI’s. We used the maps obtained using the reference protocol and ULLS as baseline for evaluating the maps estimated with the other protocols. The standard protocol is currently used for clinical research acquisitions, and the fast protocol was conjured up as a potential candidate for future clinical research acquisitions. The estimated acquisition time for the fast protocol is 5:54 minutes.

Results

Table 1 shows the processing times for generating the scalar tensor-derived measures using different tensor estimation methods, and summarizes the features of the algorithms. Figure 1, Figure 2, and Figure 3 show the MK, MD, and FA maps obtained using the standard protocol. The regions with the most noticeable change in the estimated MK across all datasets are in the splenium and the genu of the corpus callosum, where the voxels with the incorrectly estimated MK values are in black. The CLLS-H algorithm corrects the MK values for some of these voxels, whereas the CLLS-QP method corrects all of the voxels. Although constrained methods also improve the estimates of the MD and FA, this improvement is visually less discernible. Figure 4 shows the voxels where tensor estimates obtained via ULLS and the standard and reference protocols violate the constraints in Eqs. [4] or [5]. With the standard protocol, very few voxels violate the Eq. [4] constraint in general. Such voxels are typically located in the corpus callosum. Minimum kurtosis constraint violations also typically occur in the corpus callosum. These voxels are less likely to violate the maximum kurtosis constraints, although this is not uncommon. The cerebrospinal fluid (CSF) voxels show a prominent pattern of violating the maximum kurtosis constraints.

Table 1

Feature comparison between the UNLS, ULLS, CLLS-H, and CLLS-QP algorithms. Run times are for all slices in Dataset 2 obtained with the standard protocol and include the input/output time.

Figure 1

MK maps obtained using the standard protocol and (a) UNLS; (b) ULLS; (c) CLLS-H; (d) CLLS-QP; and (e) the reference protocol and ULLS.

Figure 2

MD maps obtained using the standard protocol and (a) UNLS; (b) ULLS; (c) CLLS-H; (d) CLLS-QP; and (e) the reference protocol and ULLS.

Figure 3

FA maps obtained using the standard protocol and (a) UNLS; (b) ULLS; (c) CLLS-H; (d) CLLS-QP; and (e) the reference protocol and ULLS.

Figure 4

Voxels (indicated in red) for which the diffusion and kurtosis tensors estimated using standard (rows 1, 3, and 5) and reference (rows 2, 4, and 6) protocols and ULLS violate the constraints on (a) directional diffusivities; (b) minimum directional kurtoses; (more ...)

Table 2 shows the root mean-square error (RMSE) between the maps obtained using the standard protocol and the maps estimated using the reference protocol. RMSE values are shown for the whole brain as well as for voxels with constraint violations. The whole-brain RMSE values averaged over all subjects indicate that the constrained methods outperform the unconstrained methods in estimating all three scalar quantities. The RMSE improvement offered by the constrained methods for the voxels with constraint violations is more dramatic, and is more than 35% for the MK, 7% for the MD, and 8% for the FA. Note that a lower threshold of −2, the theoretical minimum kurtosis, was applied to MK values estimated with the unconstrained methods to prevent their RMSE values from becoming excessively large. The difference between the RMSE’s of the two unconstrained methods is fairly small, and so is the difference between the constrained methods.

Table 2

RMSE between the parametric maps estimated using the standard protocol and different estimation methods, and the maps obtained using the reference protocol and ULLS; (a) MK; (b) MD; and (c) FA. In each cell, two RMSE values are reported: the first RMSE (more ...)

Figure 5, Figure 6, and Figure 7 show the MK, MD, and FA maps estimated using the fast protocol, while Figure 8 shows the voxels for which the constraints are violated. In this case, the voxels with constraint violations are distributed more widely in the brain parenchyma, although their most consistent location is the corpus callosum. Table 3 shows the RMSE values for the maps obtained using the fast protocol. As with the standard protocol, the RMSE values of the unconstrained methods are comparable. However, with the fast protocol, the whole-brain RMSE gain due to the CLLS-QP method is larger than with the standard protocol. For the voxels with constraint violations, the RMSE improvement due to the constrained methods exceeds 40% for the MK, 10% for the MD, and 19% for the FA.

Figure 5

MK maps obtained using the fast protocol and (a) UNLS; (b) ULLS; (c) CLLS-QP; and (d) the reference protocol and ULLS. Note that due to the different numbers of gradient directions, the CLLS-H algorithm is not applicable.

Figure 6

MD maps obtained using the fast protocol and (a) UNLS; (b) ULLS; (c) CLLS-QP; and (d) the reference protocol and ULLS.

Figure 7

FA maps obtained using the fast protocol and (a) UNLS; (b) ULLS; (c) CLLS-QP; and (d) the reference protocol and ULLS.

Figure 8

Voxels (indicated in red) for which the diffusion and kurtosis tensors estimated using the fast protocol and ULLS violate the constraints on (a) directional diffusivities; (b) minimum directional kurtoses; and (c) maximum directional kurtoses. The shade (more ...)

Table 3

RMSE between the parametric maps estimated using the fast protocol and different estimation methods, and the maps obtained using the reference protocol and ULLS; (a) MK; (b) MD; and (c) FA. In each cell, two RMSE values are reported: the first RMSE value (more ...)

Figure 9 shows the directional diffusivities and kurtoses for two Dataset 1 voxels with constraint violations. In the corpus callosum voxel, the UNLS solution violates 5 out of 90 constraints, and the ULLS solution violates 11 constraints. The CLLS-H solution reduces the total number of violations from 11 to 5, and the estimate based on CLLS-QP satisfies all constraints as expected. In the prefrontal white matter voxel, the unconstrained algorithms both violate one constraint, whereas both constrained algorithms satisfy all constraints.

Figure 9

Voxels in the (a) corpus callosum and (b) prefrontal white matter from Dataset 1 and their corresponding directional diffusivities (top right), kurtoses (bottom left), and the product b_maxD(n)K(n) (bottom right), estimated using different unconstrained (more ...)

Finally, Figure 10 shows the maps estimated using the new expressions for the MK and K, along with the map for K_║.

Figure 10

Kurtosis maps for Dataset 2 estimated using the standard protocol, utilizing the CLLS-QP algorithm and the proposed exact expressions; (a) MK; (b) K_║; and (c) K. Note that the maps are displayed with different scales.

Discussion

The voxels in Figure 1 and Figure 5 where CLLS-H and CLLS-QP produce improved MK estimates are those for which ULLS yields an NPD or a nearly NPD estimate of D and/or the constraints on the minimum and maximum kurtoses are violated. In these voxels, which most consistently occur in highly anisotropic regions such as the corpus callosum, λ₁ is large and λ₃ is small. This situation increases the sensitivity of unconstrained methods to noise. One effect of noise is to create the appearance of non-decreasing diffusion signal intensity with increasing b. This typically occurs along gradient directions with high diffusivity, but can also happen along directions with very low diffusivity. At high b-values, a small disturbance due to noise can artificially increase the signal intensity with increasing b-value, yielding falsely large estimates of directional kurtoses. A prominent instance of this effect can be observed in the CSF voxels in Figure 4 and Figure 8.

The other effect of noise is to create the appearance of the diffusion signal being a concave function of b along low-diffusivity directions, yielding negative estimates of directional kurtoses. In this case, the measured diffusion-weighted signal intensity may even exceed the signal intensity at b = 0, in which case a negative estimate of diffusivity, that is, an NPD DT, arises.

It is well-known that imaging noise causes λ₁ to be overestimated and λ₃ to be underestimated (26). The constraints in Eqs. [4] and [5] force a smaller λ₁ and larger λ₃, yielding more reasonable estimates of D and W. The top right plot of Figure 9(a) shows an instance of this, where the CLLS-QP estimates of directional diffusivities along directions with high diffusivity are slightly lower than those obtained with ULLS, whereas the opposite is the case along low-diffusivity directions.

Table 2 and Table 3 provide further confirmation of the improvement offered by the constrained algorithms. Note that the whole-brain RMSE improvement is actually due to only 10–20% of voxels for which there are constraint violations, hence the RMSE improvement within the set of voxels with constraint violations is quite larger than the improvement over the whole brain. The constrained algorithms achieve comparable RMSE, and both outperform the unconstrained algorithms consistently across all subjects, measures and protocols. However, the RMSE improvement is more evident with the fast protocol. Note also that the average whole-brain RMSE in estimating the MK and MD with the fast protocol and the constrained methods is smaller than or is comparable to the whole-brain RMSE with the standard protocol and unconstrained methods. This means that the constrained methods can achieve a savings of more than 20% in acquisition time. (Note that a total of 56 b = 0 images and DWI’s were utilized in the fast protocol, whereas in the standard protocol 71 images were used.) The RMSE improvement offered by the constrained methods is less pronounced in estimating the FA. A possible explanation for this may be that the FA is more heavily influenced by the Rician noise corrupting the DWI’s than the MD, consistent with previous reports in the literature (27). Therefore, to obtain more reliable FA estimates, constrained methods may need to be combined with schemes to compensate for the Rician noise bias.

The CLLS-H algorithm strikes a balance between complexity and accuracy by correcting the MK in most of the voxels that violate the constraints, doing so at a negligible 2.5% overhead compared to the ULLS (cf. Table 1). Thus, the CLLS-H method can be beneficial in real-time situations where run time is crucial. On the other hand, the CLLS-QP approach removes all such voxels, but it has a much longer run time. It should be noted that while the improvement offered by CLLS-QP comes at the expense of increased algorithm complexity and consequently higher run time, it is straightforward to speed up the computations via parallelization.

A limitation of the formulation in Eq. [7] is that it ensures that constraints in Eqs. [4] and [5] are satisfied along all imaged gradient directions, rather than all possible directions. Imposing the latter constraints is more challenging and requires a deeper understanding of the kurtosis tensor. A practically viable alternative is to estimate the tensors using either of the constrained algorithms and then check if the solution satisfies the constraints along other directions. The directions of the eigenvectors of D are of the most interest as Eq. [25] and Eq. [32] are explicit functions of the kurtosis tensor elements along these directions. If the solution does not satisfy the constraints along these directions, then these constraints can be added to the set of constraints in Eqs. [4] and [5] and the problem re-solved. This procedure may be repeated until the solution satisfies the constraints.

Finally, as the expressions for K and MK used for estimating the maps in Figure 10 represent the exact values of these parameters, they more accurately quantify the diffusional kurtosis characteristics of the tissue. The advantage of K as defined here and in (8,20) over that of (15) is that ours gives all radial directions equal weight, rather than selecting out two particular radial directions.

Conclusions

We presented two related advancements to the DKI estimation framework. The first advancement significantly improves the fidelity of the tensors parameterizing the DKI model by casting the estimation problem as CLLS. The CLLS formulation ensures that the estimated diffusivities and kurtoses along the imaged gradient directions remain within a physically and biologically plausible range. Two algorithms were provided for minimizing the CLLS cost function: a convex quadratic programming (CLLS-QP) method and a heuristic (CLLS-H) algorithm. The CLLS-H algorithm imposes a negligible computational overhead compared to the ULLS algorithm, making it suitable for real-time applications with two b-values and equal numbers of gradients per b-value, whereas the computationally more demanding CLLS-QP algorithm permits an arbitrary number of b-values and different gradient sets per each b-value. In vivo image data demonstrated the significant improvement offered by the proposed constrained tensor estimation methods over unconstrained algorithms. The advantage of the constrained methods is particularly evident when fewer DWI’s are available. The RMSE improvement due to constrained methods may be translated into shorter scan times and/or increased spatial resolution for a given fidelity of the parametric maps; our results indicate a potential time savings of more than 20%. The second advancement allows for more accurate and efficient calculation of MK and K by utilizing the exact closed-form formulae derived in this paper. The proposed tensor estimation methods along with the exact expressions for the kurtosis measures improve the robustness of the DKI metrics, thus facilitating the clinical application of DKI.

Acknowledgement

The authors would like to thank Caixia Hu, Els Fieremans, and Vitria Adisetiyo for their technical assistance.

Appendix

The MK may be calculated from the formula

[A1]

where

[A2]

and n_i is the i-th element of unit direction vector n. MK may be obtained by numerically evaluating the integral in Eq. [A2]. However, this can be computationally inefficient if a high degree of accuracy is desired. Below we give an analytic formula for A_ijkl that enables fast computation of the MK.

In a reference frame that diagonalizes D (see Eq. [26]), we can write

[A3]

where A_ijkl are given by

[A4]

and

[A5]

It can be seen from Eq. [A4] that A_ijkl is invariant with respect to permutations of its indices. Thus, we need consider only the 15 components A₁₁₁₁, A₁₁₁₂, A₁₁₁₃, A₁₁₂₂, A₁₁₂₃, A₁₁₃₃, A₁₂₂₂, A₁₂₂₃, A₁₂₃₃, A₁₃₃₃, A₂₂₂₂, A₂₂₂₃, A₂₂₃₃, A₂₃₃₃, and A₃₃₃₃. From symmetry of the integrand in Eq. [A4], one may show

[A6]

This leaves only 6 nonzero components to consider. We can further reduce this number by noting the symmetries

[A7]

and

[A8]

Hence we only need to find analytic expressions for A₁₁₁₁ and A₂₂₃₃. From these, and by applying conventional integration techniques and also utilizing = (λ₁ + λ₂ +λ₃)/3, the MK can be calculated as prescribed by Eq. [25], where F₁(λ₁, λ₂, λ₃) = A₁₁₁₁ (λ₁, λ₂, λ₃) and F₂(λ₁, λ₂, λ₃) = 6A₂₂₃₃(λ₁, λ₂, λ₃).

When λ₁ = λ₂ or λ₁ = λ₃, Eq. [27] has a removable singularity, which can be resolved with

[A9]

When λ₂ = λ₃, Eq. [28] has a removable singularity. This can be resolved with

[A10]

where

[A11]

Finally, for the special case when λ₁ = λ₂ = λ₃, we simply have

[A12]

The radial kurtosis K is the mean diffusional kurtosis perpendicular to the direction of the highest diffusion. In a reference frame that diagonalizes D, it is given by

[A13]

where

[A14]

and

[A15]

It can be seen from Eq. [A14] that C_ijkl is invariant with respect to permutations of its indices. Thus, we need consider only the 15 components C₁₁₁₁, C₁₁₁₂, C₁₁₁₃, C₁₁₂₂, C₁₁₂₃, C₁₁₃₃, C₁₂₂₂, C₁₂₂₃, C₁₂₃₃, C₁₃₃₃, C₂₂₂₂, C₂₂₂₃, C₂₂₃₃, C₂₃₃₃, and C₃₃₃₃. From symmetry of the integrand in Eq. [A14] and the fact that n₁ (ϕ) = 0, one may show

[A16]

This leaves only 3 nonzero components to consider. We can further reduce this number by noting the symmetry

[A17]

Hence we only need to find analytic expressions for C₂₂₂₂ and C₂₂₃₃. From these, K can be obtained as prescribed by Eq. [32], where the integrals

[A18]

and

[A19]

are elementary forms which lead to Eqs. [33] and [34], and G₁(λ₁,λ₂,λ₃) = C₂₂₂₂(λ₁,λ₂,λ₃) and G₂(λ₁,λ₂,λ₃) = 6C₂₂₃₃(λ₁,λ₂,λ₃).

For the special case of λ₂ = λ₃, Eq. [33] and Eq. [34] reduce to

[A20]

and

[A21]

Footnotes

¹Grant support: NIH 1R01AG027852, 1R01EB007656, and Litwin Foundation for Alzheimer’s Research.

References

1. Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophysical Journal. 1994;66:259–267. [PubMed]

2. Basser PJ. Inferring microstructural features and the physiological state of tissues from diffusion-weighted images. NMR in Biomedicine. 1995;8:333–344. [PubMed]

3. Basser PJ, Pierpaoli C. Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. Journal of Magnetic Resonance B. 1996;111:209–219. [PubMed]

4. Wedeen VJ, Hagmann P, Tseng W-YI, Reese TG, Weisskoff RM. Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magnetic Resonance in Medicine. 2005;54(6):1377–1386. [PubMed]

5. Jensen JH, Helpern JA. Quantifying non-Gaussian water diffusion by means of pulsed-field-gradient MRI; Toronto, Canada. Proceedings of the International Society for Magnetic Resonance in Medicine Annual Meeting; 2003. p. 2154.

6. Jensen JH, Helpern JA, Ramani A, Lu H, Kaczynski K. Diffusional kurtosis imaging: The quantification of non-Gaussian water diffusion by means of magnetic resonance imaging. Magnetic Resonance in Medicine. 2005;53:1432–1440. [PubMed]

7. Lu H, Jensen JH, Ramani A, Helpern JA. Three-dimensional characterization of non-gaussian water diffusion in humans using diffusion kurtosis imaging. NMR in Biomedicine. 2006;19:236–247. [PubMed]

8. Jensen JH, Helpern JA. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine. 2010 Published online: May 19, 2010. DOI:10.1002/nbm.1518.

9. Falangola MF, Jensen JH, Babb JS, Hu C, Castellanos FX, Di Martino A, Ferris SH, Helpern JA. Age-related non-Gaussian diffusion patterns in the prefrontal brain. Journal of Magnetic Resonance Imaging. 2008;28:1345–1350. [PMC free article] [PubMed]

10. Raab P, Hattingen E, Franz K, Zanella FE, Lanfermann H. Cerebral gliomas: Diffusional kurtosis imaging analysis of microstructural differences. Radiology. 2010;254:876–881. [PubMed]

11. Jansen JFA, Stambuk HE, Koutcher JA, Shukla-Dave A. Non-Gaussian analysis of diffusion-weighted MR imaging in head and neck squamous cell carcinoma: A feasibility study. AJNR Am J Neuroradiol. 2009 ajnr.A1919.

12. Cheung MM, Hui ES, Chan KC, Helpern JA, Qi L, Wu EX. Does diffusion kurtosis imaging lead to better neural tissue characterization? A rodent brain maturation study. NeuroImage. 2009;45:386–392. [PubMed]

13. Lazar M, Jensen JH, Xuan L, Helpern JA. Estimation of the orientation distribution function from diffusional kurtosis imaging. Magnetic Resonance in Medicine. 2008;60:774–781. [PMC free article] [PubMed]

14. Song S-K, Sun S-W, Ramsbottom MJ, Chang C, Russell J, Cross AH. Dysmyelination revealed through MRI as increased radial (but unchanged axial) diffusion of water. NeuroImage. 2002;17(3):1429–1436. [PubMed]

15. Hui ES, Cheung MM, Qi L, Wu EX. Towards better MR characterization of neural tissues using directional diffusion kurtosis analysis. NeuroImage. 2008;42:122–134. [PubMed]

16. Koay CG, Carew JD, Alexander AL, Basser PJ, Meyerand ME. Investigation of anomalous estimates of tensor-derived quantities in diffusion tensor imaging. Magnetic Resonance in Medicine. 2006;55:930–936. [PubMed]

17. Liu C, Bammer R, Acar B, Moseley ME. Characterizing non-Gaussian diffusion by using generalized diffusion tensors. Magnetic Resonance in Medicine. 2004;51:924–937. [PubMed]

18. Wang Z, Vemuri BC, Chen Y, Mareci TH. A Constrained variational principle for direct estimation and smoothing of the diffusion tensor field from complex DWI. IEEE Transactions on Medical Imaging. 2004;23:930–939. [PubMed]

19. Barmpoutis A, Hwang MS, Howland D, Forder JR, Vemuri BC. Regularized positive-definite fourth order tensor field estimation from DW-MRI. NeuroImage. 2009;45:S153–S162. [PMC free article] [PubMed]

20. Poot DHJ, den Dekker AJ, Achten E, Verhoye M, Sijbers J. Optimal Experimental Design for Diffusion Kurtosis Imaging. IEEE Transactions on Medical Imaging. 2010;29:819–829. [PubMed]

21. Stejskal EO, Tanner JE. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. Journal of Chemical Physics. 1965;42:288–292.

22. Callaghan PT, Coy A, MacGowan D, Packer KJ, Zelaya FO. Diffraction-like effects in NMR diffusion studies of fluids in porous solids. Nature. 1991;351:467–469.

23. Nocedal J, Wright SJ. Numerical Optimization. New York: Springer; 1999.

24. Carlson BC. Computing elliptic integrals by duplication. Numerische Mathematik. 1979;33:1–16.

25. Carlson BC, Notis EM. Algorithm 577: Algorithms for incomplete elliptic integrals [S21] ACM Transactions on Mathematical Software. 1981;7(3):398–403.

26. Anderson AW. Theoretical analysis of the effects of noise on diffusion tensor imaging. Magnetic Resonance in Medicine. 2001;46(6):1174–1188. [PubMed]

27. Basu S, Fletcher T, Whitaker R. Lecture Notes in Computer Science. Volume 4190. Berlin: Springer; 2006. Rician noise removal in diffusion tensor MRI; pp. 117–125.