Diffusion-tensor-based method for robust and practical estimation of axial and radial diffusional kurtosis.

OBJECTIVES: A new method that can estimate diffusional kurtosis image (DKI), estimated DKI (eDKI), parallel and perpendicular to neuronal fibres from greatly limited image data was designed to enable quick and practical assessment of DKI in clinics. The purpose of this study was to discuss the potential of this method for clinical use.
METHODS: Fourteen healthy volunteers were examined with a 3-Tesla MRI. The diffusion-weighting parameters included five different b-values (0, 500, 1,500, 2,000 and 2,500 s/mm(2)) with 64 different encoding directions for each of the b-values. K values were calculated by both conventional DKI (convDKI) and eDKI from these complete data, and also from the data that the encoding directions were abstracted to 32, 21, 15, 12 and 6. Error-pixel ratio and the root mean square error (RMSE) compared with the standard were compared between the methods (Wilcoxon signed-rank test: P < 0.05 was considered significant).
RESULTS: Error-pixel ratio was smaller in eDKI than in convDKI and the difference was significant. In addition, RMSE was significantly smaller in eDKI than in convDKI, or otherwise the differences were not significant when they were obtained from the same data set.
CONCLUSION: eDKI might be useful for assessing DKI in clinical settings. KEY POINTS: • A method to practically estimate axial/radial DKI from limited data was developed. • The high robustness of the proposed method can greatly improve map images. • The accuracy of the proposed method was high. • Axial/radial K maps can be calculated from limited diffusion-encoding directions. • The proposed method might be useful for assessing DKI in clinical settings.

Introduction

The diffusional kurtosis imaging (DKI) technique is being gradually accepted in clinics. The advantage of DKI over classical apparent diffusion coefficient (ADC) is that the method provides information about non-Gaussian water diffusion in in vivo tissues [1]. Generally speaking, in vivo diffusion does not follow a Gaussian distribution because of the interactions between water molecules and cell structures (especially cell membranes). Thus, important information concerning the tissue microstructure would be lost if the data were fitted to a Gaussian model as it would be in ADC [1-3]. In this regard, many diffusion models have been evaluated, and DKI has been accepted as a relatively simple model that can be utilized in clinics. Its significance as a diagnostic tool has been mostly discussed in terms of the central nervous system (CNS) [4-7], but it is gradually being adapted for other lesions such as tumours [8, 9].

However, the long acquisition time of DKI resulting from the requirements for multiple b-values and the large number of diffusion encoding directions (NEDs) of the raw diffusion-weighted images (DWIs) remains a problem [3, 10]. A typical conventional protocol can consist of 30 uniformly distributed encoding directions in three to six different b-values between 0 and 3,000 s/mm2 [3]. This is quite large compared to the widely accepted diffusion-tensor imaging (DTI) that usually consists of six to 20 encoding directions and two different b-values.

Another problem is that the parameters assessed in DKI were mostly limited to mean diffusional kurtosis (mean K) in previous clinical reports. Axial and radial kurtosis (axial K, radial K), which in the CNS are diffusional kurtosis values of the directions parallel and perpendicular to neuron fibres, have not been much evaluated for use in clinics. This might due to the severe noise and errors derived from the quite tedious post-processing operation to calculate these parameters. Mean K was commonly calculated by averaging the apparent K value over all the acquired encoding directions [3]. This method is robust because it does not require calculation of the kurtosis tensor, which is generally the most time-consuming and noise-generating part. This method can also smooth the images and decrease noise and errors via the averaging step. On the other hand, obtaining axial and radial K was not very robust because the step of calculating the kurtosis tensor could not be skipped. Some computing algorithms have been proposed to overcome this problem [11-13], but the required image data are still too large for daily clinical use.

However, considering the importance of the interactions of water molecules and the cell membrane, which was the basic principal of non-Gaussian diffusion, diffusional anisotropy will represent important information, especially in the CNS. The significance of DKI in clinics might become more promising if axial and radial DKI can be easily obtained with high quality from limited image data.

In this study, we designed a new method, eDKI (estimated DKI) based on diffusion-tensor imaging (DTI) and quadratic function fitting, which can robustly estimate axial and radial K from even smaller data than the general DKI protocol. The purpose of this study was to assess the potential capability of eDKI compared to the standard conventional DKI method regarding robustness and accuracy of the estimation.

Materials and methods

This volunteer study was approved by our local ethics committee. Written informed consent was obtained from all subjects.

Subjects and imaging procedure

We recruited 15 healthy female volunteers for the study. None had a past history of CNS disease. Brain magnetic resonance imaging (MRI) scans were obtained for each volunteer with a 3-Tesla clinical MRI (Skyra 3T; Siemens, Germany) with a 20-channel head coil. First, T1-weighed imaging (three-dimensional magnetization-prepared rapid gradient-echo imaging (MPRAGE)) was performed to gather and assess anatomical information. Second, diffusion-weighted imaging (DWI) with five different b-values (0, 500, 1,000, 1,500 and 2,500 s/mm2) and 64 different motion-probing gradient (MPG) encoding directions, which were uniformly distributed for each b-value, was obtained. The readout-segmented echo-planar imaging multi-shot (RESOLVE®) technique was used for DWI acquisition. The duration time of MPG (δ) and the interval time between MPGs (Δ) were kept identical at 37.8 and 47.3 ms, respectively. The other major parameters of this DWI sequence were as follows: repetition time/echo time, 4,700/100 ms; field of view: 220 × 220 mm; matrix: 126 × 126; section thickness: 3 mm; 16 slices with 0.9-mm gap. The images were acquired in a horizontal slice direction perpendicular to the line connecting anterior and posterior commissures. The slices were located between the frontal cranial base and the top of the brain parenchyma. Images were acquired three times and averaged for b = 0 s/mm2. The other images were acquired once. Scan time for this DWI sequence was approximately 1 h.

No volunteers had abnormal findings in their brain structures or abnormal signal intensities in MPRAGE and b = 0 images. One volunteer (23 years old) was excluded from this study because of excessive motion during the scan. Finally, 14 female volunteers (age range 20–39 years, average 24.4 years) were selected for this study.

Image processing

The pixels outside the area of brain parenchyma were excluded semi-automatically from the images as follows. First, regions of interest (ROIs) were roughly designed around the brain parenchyma on b = 0 images manually. Second, to exclude the cerebrospinal fluid (CSF) pixels, the pixels with a mean ADC [2] calculated from data at signal intensities of b = 0 and 1,000 s/mm2 of over 1.5 × 10−3 mm2/s were excluded.

The images were also smoothed with a median filter (eight connected neighbourhoods).

The raw data consisting of 64 diffusion-encoding directions were abstracted in five different patterns to create new partial data sets with 32, 21, 15, 12 and six encoding directions, respectively. For example, to make a data set with 21 encoding directions, data from 43 encoding directions were excluded from the original full data with 64 encoding directions. The excluded directions were identical for all b-values, and they were selected evenly so that the remaining directions were not biased. From each of these data sets, axial and radial DKI images were calculated by the conventional method and also estimated by the proposed method (eDKI) as described in the following sections. All the image processions in this study were made by our in-house software working on a commercial analysis package (MATLAB® version 2014a, MathWorks Inc., Natick, MA, USA) (See Appendix for a stand-alone version).

Calculation of axial and radial diffusional kurtosis imaging (DKI) by conventional method

Conventional DKIs (convDKI) were calculated in axial and radial directions of each subject in four different sets of encoding directions (64, 32, 21 and 15 encoding directions, respectively; abbreviated in this study as convDKI-64, convDKI-32, etc.). They were calculated pixel-by-pixel according to previous reports [13-15]. The other data sets consisting of 12 and six encoding directions could not be applied to this procedure, because the conventional method requires at least 15 different encoding directions. The brief calculation procedure for each pixel was as follows. First, quadratic-function curve fitting was performed for the b-value-related signal changes in each encoding direction, separately (a closed-form method was applied) [10, 14].

The proper apparent diffusional kurtosis (K) and corrected diffusion coefficient (D) of each encoding direction were obtained by this step. Second, the kurtosis tensor, which has 15 independent variables, was calculated by the least-squares method from these results and the parameters of the applied encoding directions. Third, the diffusion tensor was calculated at the b-value pair of b = 0 and 1,500 s/mm2 to define the proper axial and radial direction of the pixel. Finally, axial and radial convDKIs were generated by projecting the kurtosis tensor to the axial and radial directions defined above.

Axial and radial K calculated by this conventional method with 64 encoding directions (axial and radial convDKI-64) were regarded as standards in this study.

Estimation of axial and radial DKI by the proposed method (eDKI)

The estimated images of axial and radial DKI were also calculated for each subject in six different encoding numbers (e.g. 64, 32, 21, 15, 12 and 6; abbreviated in this study as eDKI-64, eDKI-32, etc.) by the proposed method. The essence of the method is to simulate the virtual b-value-dependent signal changes regarding the proper axial and radial directions of the pixel on the basis of the diffusion tensor imaging method (estimated DWI [16]), and then apply the quadratic-function fit to it to calculate apparent diffusion kurtosis for each diffusion encoding direction). For this virtual imaging, the following steps were carried out.

The diffusion tensor and its eigenvalues (λ1b, λ2b and λ3b) were calculated pixel-by-pixel for every b-value accompanied by b = 0 (e.g. the pair of b = 0 and 500, the pair of b = 0 and 1,000… and the pair of b = 0 and 2,500) (Fig. 1A). As the first eigenvalue and the average of the second and third eigenvalues of each tensor indicate the supposed axial and radial diffusion coefficients of the pixel at the applied b-value, they were substituted in Db of the equation,

Fig. 1

Methods for obtaining axial/radial estimated diffusional kurtosis imaging (eDKI). (A) From the raw data with multiple b-values and diffusion encoding directions, diffusion tensor was calculated for each b-value. The largest eigenvalue of this tensor (λ1) was regarded as the axial diffusion coefficient, and the average of the other two eigenvalues (λ2, λ3) was regarded as the radial diffusion coefficient. (B) For each axial and radial direction, the b-value-dependent signal change was calculated and plotted at each b-value using the diffusion coefficients obtained above. Axial and radial diffusional kurtosis were calculated by fitting this signal change curve to the quadratic function [14]. Linear correction was also performed to bring the contrast of eDKI closer to conventional DKI (convDKI) (see the section Estimation of axial and radial DKI by the proposed method)

to obtain the virtual signal intensity (S ) in each pixel at the b-value reflecting the axial and radial diffusions (Fig. 1B). Here, S0 is the signal intensity at b = 0 and b indicates the applied b-value. Finally, estimated axial and radial K were obtained by applying the quadratic-function curve fit, similar to Eq. 1 in the previous section, to these virtual b-value-dependent signals of estimated DWI (in Eq. 3: K and D indicate the estimated axial and radial K and D that was obtained simultaneously in this fitting).

In a previous phantom study and a simple computer simulation (not shown), the values of eDKI (K in Eq. 3) were shifted linearly in comparison to those of convDKI. Therefore, the eDKI values were corrected by a linear function to bring the values and the contrast closer to those of convDKI.

ConvDKI-64 and eDKI-64 were used as follows to define the correcting parameters p and q for axial and radial directions, respectively. First, the correcting parameters p and q were changed from 0.6 to 1.4 and from −0.5 to 0.5 step-by-step (0.02 for each step), respectively. For each pair of p and q, the root mean square error (RMSE) between convDKI-64 and corrected eDKI-64s for that pair (RMSEcor(p, q)) were calculated.

N corresponds to the number of pixels included in the brain parenchymal area, and the variable i indicates each of the pixels included there. Second, the pair of p and q achieving the smallest RMSEcor(p, q) was defined as the temporary-correcting parameter pair of that subject. Third, for the final correction of each subject, the temporary-correcting parameter pairs of p and q were averaged among the other 13 subjects and applied, instead of directly adopting the pair defined in each subject. This was done to prevent biased results for the later statistical assessments (leave-one-out cross-validation). The same p and q were also applied in each subject to correct the eDKIs obtained from other data sets of less encoding directions (eDKI-32, eDKI-21, eDKI-15, eDKI-12 and eDKI-6).

Comparison of conventional and proposed methods

Comparisons between convDKI and eDKI were made to assess robustness and accuracy. First, for robustness, the pixels of which the calculated value of DKI was not in a certain range were regarded as error-pixels. For this purpose, the lower limit was set to zero and the upper limit to 1.5 and 3 for axial and radial directions, respectively, according to previous reports [17, 18]. The ratio of the error-pixels was calculated for both convDKIs and eDKIs, respectively, as follows.

The ratio of error-pixels was statistically compared between convDKI and eDKI.

Second, to assess accuracy, RMSE values between the conventional standard DKI (convDKI-64) and the other DKI images (convDKI-32 to 15 and eDKI-32 to 6) were calculated and analysed between convDKIs and eDKIs. In this analysis, the pixels defined above as error-pixels in convDKI and/or eDKI were excluded. The differences in RMSEs between convDKI and eDKI were evaluated statistically. The Wilcoxon signed-rank test was used for the statistical comparisons. A P-value <0.05 was considered significant.

Results

The motions of the volunteers were assessed by using automated registration of SPM8 (Statistical Parametric Mapping, www.fil.ion.ucl.ac.uk/spm) [19]. In the one excluded subject, the translation of head-feet direction among the scans was more than 10 mm while those of the other volunteers were less than 3 mm.

Average and standard deviations of the parameters used for correcting the contrast of eDKI were: p = 0.92 ± 0.07 and q = 0.14 ± 0.06 for axial eDKI, and p = 0.90 ± 0.04 and q = 0.07 ± 0.03 for radial eDKI.

The noise and errors of the map images were well reduced in eDKI compared to convDKI (Fig. 2).

Fig. 2

Axial and radial conventional diffusional kurtosis imaging (convDKI) and estimated DKI (eDKI). Images of the axial and radial convDKI and eDKI maps obtained from a representative volunteer. The eDKI maps are much less noisy and more consistent among the different numbers of encoding directions than those of convDKI. The left-right asymmetries seen in convDKI are thought to be errors, of which less are found in eDKI. ConvDKI cannot be calculated from less than 15 encoding directions (for theoretical reasons) while eDKI could be obtained even from six encoding directions

The ratios of error-pixels became elevated in both convDKI and eDKI as the NED decreased (Fig. 3). The elevation was smaller in eDKI. The median error-pixel ratio of eDKI was constantly smaller than that of convDKI in any NED. The differences were statistically significant except between convDKI-64 and eDKI-6 in a radial direction (Fig. 3). In addition, the difference between the methods in the same NED became larger as the NED decreased.

Fig. 3

Analysis of error-pixel ratio for assessment of robustness of conventional diffusional kurtosis imaging (convDKI) and estimated DKI (eDKI). The median error-pixel ratio of eDKI was constantly smaller than that of convDKI in any numbers of encoding directions. The differences were statistically significant (Wilcoxon signed-rank test, P < 0.05) except for * pairs between convDKI-64 and eDKI-6 in radial direction. The number after convDKI and eDKI indicates the number of encoding directions

In the RMSE analysis, the values of convDKI and eDKI both elevated as the NED decreased, as well as the ratio of error-pixels. The elevation was also dramatically smaller in eDKI (Table 1, Fig. 4). In the same NED, the median RMSE of eDKI was smaller than that of convDKI with a significant difference except in the radial DKI of 32 encoding directions (Table 1, Fig. 4). In addition, the difference became larger as the NED decreased. The median RMSE of eDKIs compared with that of convDKI-32 was smaller by a significant difference, or the difference was not significant when the NED for eDKI was 15 or more in axial direction, and 21 or more in radial direction (Fig. 4). The median RMSE values of eDKIs were smaller with a significant difference than those of convDKI-21 and convDKI-15 in any NED (Fig. 4).

Table 1

Root mean square error (RMSE) calculated with standard images (convDKI-64)

Number of encoding directions		64	32	21	15	12	6
Required scan time (min : sec)		61 : 08	30 : 55	20 : 32	14 : 52	12 : 02	6 : 22
Axial	ConvDKI	-	0.185	0.249	0.461	-	-
	(min, max)		(0.173, 0.199)	(0.240, 0.270)	(0.451, 0.471)
	eDKI	0.174	0.178	0.185	0.190	0.193	0.198
	(min, max)	(0.161, 0.186)	(0.165, 0.189)	(0.173, 0.195)	(0.176, 0.200)	(0.182, 0.206)	(0.182, 0.278)
	P-value a	-	<0.01	<0.01	<0.01	-	-
Radial	ConvDKI	-	0.277	0.396	0.986	-	-
	(min, max)		(0.253, 0.296)	(0.372, 0.418)	(0.965, 1.004)
	eDKI	0.269	0.279	0.288	0.296	0.309	0.312
	(min, max)	(0.232, 0.287)	(0.239, 0.292)	(0.251, 0.302)	(0.259, 0.310)	(0.270, 0.323)	(0.276, 0.424)
	P-value a	-	0.27	<0.01	<0.01	-	-

aWilcoxon signed-rank test. P-value under 0.05 was considered significant in this study

Numbers denote medians and ranges

DKI diffusional kurtosis imaging, convDKI conventional DKI, convDKI-64 convDKI obtained from 64 encoding directions, eDKI estimated DKI

Fig. 4

Analysis of root mean square error (RMSE) for assessment of accuracy of conventional diffusional kurtosis imaging (convDKI) and estimated DKI (eDKI). RMSEs were obtained between the standard image (convDKI-64) and the other images. The median RMSE values were smaller in eDKI than in convDKI with a significant difference except for * and ** pairs between convDKI-32 and some eDKIs.*difference was not significant, **eDKI was larger than convDKI with a significant difference; Wilcoxon signed-rank test, P < 0.05. The number after convDKI and eDKI indicates the number of encoding directions

Discussion

Various techniques to analyse water diffusion in vivo had been developed and improved recently. Some of those techniques such as Mean Apparent Propagator (MAP) MRI [20] may enable us to attain more knowledge about the molecular mechanisms involved in this topic. However, most of these methods will potentially require highly developed scanning that includes, for example, higher magnetic field strength, higher gradient strength, a long scan time and so on. On the other hand, some other techniques may not provide such advanced information but could be used in clinical MRI settings to provide conventional information more easily and/or in higher quality than previously. The eDKI method proposed in the present study is one of the techniques that belong to the latter group. It was designed to estimate axial and radial K robustly and with practical accuracy from limited image data. Its purpose was to facilitate the evaluation of axial and radial K in clinics, a process providing insufficient data by the conventional method. The map images of eDKI (Fig. 2) as well as the results of the assessment of the error-pixel ratio (Fig. 3) and RMSE (Fig. 4) may suggest that eDKI is superior to the conventional method in this regard: first because eDKI can improve the quality of the map image with even higher accuracy, and second because eDKI maps could be generated from very limited encoding directions, such as 12, which requires only approximately 12 min of scanning with our settings. Note that eDKI subsumes DTI, so a further additional scan time for DTI is not needed. Potentially, scan time can be further shortened by reducing the number of b-values (e.g., 3 min in total when the imaging sequence consists of three b-values and six encoding directions), but this might be better discussed in further studies including pathological imaging.

The basic idea for gaining robustness is to skip the steps of mathematical fitting as much as possible. The convDKI requires a quadratic function fitting of each encoding direction (at least 15 times), followed by a fitting calculation to obtain rank-4 kurtosis tensor that has 15 independent variables [13, 14, 18]. This complicated fitting process generates a lot of noise and errors. On the other hand, eDKI requires only a fitting calculation to obtain diffusion tensor in each b-value, and needs only one quadratic function fitting for each axial and radial direction. As diffusion tensor has only six independent variables, it could be robustly calculated from a relatively limited data size. The concept of estimating the virtual b-value-dependent signal change in a certain direction before curve fitting was previously discussed in terms of combining diffusion-tensor imaging and the two-compartment model of water diffusion in vivo to assess the normal-appearing white matter of multiple sclerosis [21]. This method may be useful for combining other multiple b-value analyses and diffusion-tensor imaging as well.

The robustness of acquiring axial and radial K was greatly improved by the eDKI method. The image quality of the map images was improved (Fig. 2), and this was also statistically confirmed by comparison of the error-pixel ratios (Fig. 3).

The accuracy of the estimation of eDKI might be high enough for practical use, based on the result of the analysis of RMSE, compared with the standard convDKI (Table 1, Fig. 4). RMSE was significantly smaller in eDKI than in convDKI in both axial and radial directions in most of the cases with the same NED. The only exception was that between convDKI-32 and eDKI-32 in radial direction, but the difference was not significant (Table 1). From this result, eDKI may generally be substituted for convDKI when the NED is 32 or less.

The eDKI method should be strongly recommended for data with 21 encoding directions (which might represent modern clinical standard with relatively poor data) because both the error-pixel ratio and RMSE of eDKIs were always significantly smaller than those of convDKI-21, even when there were only six encoding directions for eDKI. The advantage of eDKI grows as the NED decreases. Viewed from another perspective, axial and radial K could reasonably be discussed at a contemporary level by using eDKI if there are data with six encoding directions.

The robustness of eDKI might also be reflected in its higher consistency among different NEDs. The elevation of the error-pixel ratio and RMSE as the NED decreased were both dramatically smaller in eDKI than in convDKI (Figs. 3 and 4). This consistency is important for the accumulation of evidence, because otherwise the K values cannot be compared with each other when the NED is different. From another point of view, with the consistency of eDKI, scan time can be flexibly modified for each clinic by changing the NED.

Another merit of eDKI is that the restriction in the numbers and combinations of encoding directions among different b-values is small, as diffusion tensor is calculated for each b-value separately: six or more NEDs for each b-value is the only requirement. In convDKI, the set of encoding directions must be identical for all b-values, because the calculation of apparent K for each direction is required first. Diffusion imaging and analysis methods that apply different sets of encoding directions for each b-value are increasingly being evaluated currently [22, 23]. It may be important that eDKI can also be used for the data acquired by those newly emerging imaging sequences.

The linear correction of eDKI in each case was performed by applying the temporary-correcting parameters averaged among the other subjects to avoid the statistical results being biased. For practical use in clinics, the correction can be done by applying the overall averaged temporary-correcting parameters that are indicated in the results section.

This study has several limitations. First, the significance of eDKI was discussed only in healthy volunteers. Its significance for detecting pathological changes should be discussed in further studies. Second, the sample size was not large. However, the statistical significances were so stable and strong that the size seemed sufficient to support our conclusions.

In conclusion, the proposed eDKI method could calculate axial and radial K robustly with practical accuracy from limited data. The method might be useful for assessing diffusional kurtosis in clinical settings.