Toward more robust and reproducible diffusion kurtosis imaging.

PURPOSE: The general utility of diffusion kurtosis imaging (DKI) is challenged by its poor robustness to imaging artifacts and thermal noise that often lead to implausible kurtosis values. THEORY AND METHODS: A robust scalar kurtosis index can be estimated from powder-averaged diffusion-weighted data. We introduce a novel DKI estimator that uses this scalar kurtosis index as a proxy for the mean kurtosis to regularize the fit.
RESULTS: The regularized DKI estimator improves the robustness and reproducibility of the kurtosis metrics and results in parameter maps with enhanced quality and contrast.
CONCLUSION: Our novel DKI estimator promotes the wider use of DKI in clinical research and potentially diagnostics by improving the reproducibility and precision of DKI fitting and, as such, enabling enhanced visual, quantitative, and statistical analyses of DKI parameters.

INTRODUCTION

Despite a growing interest in biophysical models of diffusion in white matter to develop specific biomarkers of microstructural changes, signal representations, for example, diffusion tensor imaging (DTI) or diffusion kurtosis imaging (DKI), retain the potential to become invaluable tools in diagnostic and clinical research settings. Such signal representations provide metrics that are highly sensitive to microstructural changes associated with development, , aging , , disease and disorder without adopting (often contested , , , ) model assumptions. An extensive literature has demonstrated the sensitivity of DKI parameters to white and gray matter alterations in, for example, aging, , , , stroke, , , , , traumatic brain injury, , , multiple sclerosis, schizophrenia, , autism spectrum disorder, , epilepsy, , migraine, Parkinson’s Disease, and Alzheimer’s disease. , Despite a focus on neuroapplications, DKI also has emerging applications in body MRI. , However, sensitivity is not enough; diagnostics and clinical research may only adopt biomarkers that are also accurate, precise, and robust.

Technically, DKI forms a straightforward extension of DTI and provides, aside from the diffusion tensor, an estimate of the diffusion kurtosis tensor which quantifies the degree of directional non‐Gaussian diffusion. , , DKI parameters, for example, the mean kurtosis (), radial kurtosis (), and axial kurtosis (), , , have been shown to yield clinically relevant information that is not captured by a more conventional DTI representation. Moreover, DTI parameters (eg, mean diffusivity, radial diffusivity, axial diffusivity and fractional anisotropy ) are themselves estimated more accurately when evaluated within the DKI framework (ie, estimating the diffusion tensor and kurtosis tensor simultaneously). Compared to DTI, DKI requires a slightly more extensive scan protocol in the sense that there must be at least 3 distinct b‐values (typically one of these is set to ), which differ only in their gradient magnitude. , Given the widespread availability of accelerated image acquisition techniques such as simultaneous multiband imaging, , data compatible with whole‐brain DKI‐analysis can be acquired in a few of minutes, thereby facilitating clinical utility. Moreover, many recent large cohort studies, for example, the Human Connectome Project, , UK Biobank, Adolescent Brain Cognitive Developement, and the Cambridge Centre for Ageing Neuroscience (CamCAN), , provide now large‐scale multishell data that are well‐suited to DKI analyses.

Unfortunately, DKI has been challenged ever since its introduction by a poor robustness to imaging artifacts and thermal noise that often leads to non‐physical kurtosis values during the fitting, especially in voxels with very low radial diffusivity. , , , Indeed, implausible negative values are ubiquitous in many kurtosis maps, hampering visual, quantitative, and statistical analysis of the data. Advances in artifact correction, noise removal, and constrained parameter estimation have reduced, but not eliminated the problem. , , , , , Alternatively to the full tensor estimations, a scalar kurtosis quantity can be obtained from diffusion‐weighted signals averaged across different isotropically distributed gradient directions (and for each individual b‐value): that is, the powder kurtosis , , , distinct from mean kurtosis. Due to the higher signal‐to‐noise ratio (SNR) of the powder signals and decreased number of estimated model parameters, can be computed with higher precision than other kurtosis tensor derived metrics. We will demonstrate in this work that can be used to provide a robust prediction of the mean kurtosis without the need to estimate the kurtosis tensor. However, the powder‐averaged signal does not allow for the estimation of directional kurtosis values, which have been shown to provide unique information in well‐aligned structures. , , Moreover, a robust estimation of the full kurtosis tensor is required for biophysical modeling , , and tractography. , ,

Here we propose and evaluate a novel technique for more robust and precise estimation of the full kurtosis tensor and derived metrics. We introduce a regularized DKI estimator in which the estimated mean kurtosis is evaluated against a robust prediction of the mean kurtosis, which in turn is derived from the powder kurtosis. We will describe the technical details and demonstrate how they improve the reproducibility, robustness and precision of DKI parameters.

METHODS

Diffusion kurtosis imaging

DKI provides a representation of the diffusion‐weighted signal S in terms of the second and fourth cumulant as a function of the diffusion‐weighting strength, or b‐value, b. In 1D or in the case of isotropic diffusion,with and the apparent diffusion and apparent kurtosis coefficients, respectively. Since proton diffusion in biological tissue, e.g. brain white matter, is typically anisotropic, 3D diffusion and kurtosis tensors are needed to describe the orientational dependence of the diffusion‐weighted signal adequately , :Here, is the unit vector along which the diffusion‐weighting gradient is applied, and represents the element of the fully symmetric second‐order diffusion tensor D for which a third of the trace equals the mean diffusivity . In addition, denotes the element of the fully symmetric fourth order diffusion kurtosis tensor W. Because both tensors are fully symmetric, D and W have and degrees of freedom, respectively.

From D and W, various diffusion and kurtosis parameters can be derived. , , , , Fractional anisotropy (FA), mean diffusivity (), radial diffusivity (), axial diffusivity (), mean kurtosis (), radial kurtosis (), and axial kurtosis () are amongst the most widely adopted DKI parameters. , , , Without loss of generality we here limit ourselves to the original definition of mean kurtosis , computed by evaluating the apparent kurtosis from the along a large number of directions followed by averaging.

Powder kurtosis

Analogous to Equation (1), the decay of the powder‐averaged signals as a function of b can be approximated using the cumulant expansion , :where and are the scalar diffusivity and excess‐kurtosis of powder signals, respectively. Since powder‐averaged signals are independent of the orientation distribution of microscopic components, , , , the scalar quantity extracted from these signals are decoupled from mesoscopic properties such as tissue dispersion or fiber architecture configurations.

Robust prediction of mean kurtosis

In this section, we will present 3 strategies to predict the mean kurtosis without relying on the full tensor. First, from the powder kurtosis, can be computed analytically from the signal in the limit :withThe correction term mainly depends on the anisotropy captured by the elements of the diffusion tensor . Importantly, the derivation of this expression relies on the mean kurtosis tensor definition proposed by Hansen et al, , which was shown to present nearly identical contrast to the original mean kurtosis definition .

Second, for an isotropic diffusion tensor D, the correction term is zero and yields an accurate prediction of the mean kurtosis:Third, for finite b‐values, a mean kurtosis prediction can be obtained from and the diffusion tensor D using a polynomial regression model in which the thousands of non‐problematic voxels, that is, positive apparent kurtosis in each direction, serve as training data. Note that due to the strong dependence of DKI parameters on scan settings, b‐values, and subject‐specific alteration of the underlying tissue microstructure, the training voxels must be sampled from the same subject or from datasets with the same acquisition parameters. In this strategy, coined voxel quality transfer, the polynomial coefficients are estimated using multivariate polynomial regression in which the parameters , , and are included up to the order. Hence,with a multivariate polynomial function of order N. In this work .

To summarize, the predicted mean kurtosis can be estimated robustly from the diffusion tensor D and the powder kurtosis using 3 strategies: (a) powder kurtosis with analytical correction, (b) without analytical correction, or (c) voxel quality transfer using polynomial regression. Where deemed relevant, the strategy will be specified using the superscript index. If not specified, as motivated by our results (see below).

Regularized kurtosis fitting

After establishing a robust prediction of the mean kurtosis (using 1 of the 3 strategies described above), we will use this metric to regularize the DKI fitting. Here we propose the following regularized nonlinear least squares (NLS) estimator:with the regularization weight and h the operator that computes the mean kurtosis from the tensor coefficients . , Note that is an alternative formulation of the DKI model in Equation (2) in which B represent the extended b‐matrix. For , the estimator reduces to the ordinary NLS estimator, a widely adopted estimator for DTI and DKI. , ,

We hypothesize that the prediction of the mean kurtosis from the powder‐averaged data is a much more robust and reproducible metric than estimated from the fitted kurtosis tensor and, as such, the norm of the difference between and regularizes and stabilizes the DKI fit.

The nonlinear fitting is initiated by a starting point obtained by fitting the DKI model using the ordinary NLS estimator. In the few cases where the estimator does not converge to a plausible solution, the estimation is repeated with a starting point that is the result of a constrained DKI fitting in which positivity of is imposed.

Alternative fitting strategies

We compare performance of the regularized NLS estimator to commonly used fitting strategies: (a) ordinary NLS, (b) ordinary NLS after smoothing the data with a 2D isotropic Gaussian filter with a kernel and full‐width‐half‐maximum of 1.25 times the voxel size, and (c) constrained NLS. The constrained NLS estimator solves the same objective function as the ordinary NLS estimator, but imposes one or more constraints on the estimated parameters. Here, we only adopted the constraint that the apparent kurtosis coefficient must be positive in each gradient direction .

Data

Simulated data

Monte Carlo simulations with 2500 trials were performed to evaluate the performance of the regularized DKI estimator. For each trial, diffusion‐weighted signals were generated by evaluating the DKI signal for the diffusion encoding scheme of our study‐specific data (vide infra). The ground truth diffusion tensor and diffusion kurtosis tensor coefficients were randomly sampled from the plausible DKI estimates of the HCP data. Gaussian noise was added to the simulated data with a corresponding SNR of 30 for the non‐diffusion‐weighted signal. We opted for Gaussian noise instead of complex Gaussian noise to avoid the Rician noise bias as a confound in our interpretation. The effect of Rician bias on the estimation of DKI parameters has previously been studied and documented in detail. , ,

Study‐specific MRI data

Data were collected under the approval of the Cardiff University School of Psychology Ethics Committee. Five healthy volunteers were recruited and data were collected on 2 different scanning sessions with exactly the same imaging protocol on a Siemens Connectom 3T MR scanner using a 32‐channel receiver coil. For each volunteer, the 2 test‐retest scanning sessions were performed on the same day interleaved by a short break. In both sessions, subjects were re‐positioned by the same operator. The repeated images were analyzed individually, without data averaging.

The diffusion gradients were characterized by and maximal gradient amplitude of ; note that in this work the full power of the Connectom gradient system was not exploited to ensure translation of the findings to clinical scanners (where gradient amplitudes of 80 mT/m are prevalent). Diffusion weighting was applied along 30 isotropically‐distributed gradient directions for 0.5, 1, and 2.5 , with , matrix: , and 54 slices with a spatial resolution of . Data acquisition was accelerated using simultaneous multiband (SMS = 2) and GRAPPA (R = 2), but partial Fourier encoding was turned off. In addition, 4 non‐diffusion‐weighted images were acquired with the same and reversed phase encoding to enable susceptibility‐induced geometrical correction. The total acquisition time was 6 m 15 seconds per session.

Additional MRI data

Additional experiments are based on arbitrarily selected data sets from various public neuroimaging repositories (eg, Openneuro) that include multi‐shell diffusion‐weighted MRI data. We selected data from a wide variety of (clinical) research projects: MGH Adult Diffusion Data from the Human Connectome Project (HCP, , , ), the MASiVar project, the Mexican substance use disorder database neuroimaging dataset (SUDMEX), , , and Cambridge Centre for Ageing Neuroscience (CamCAN) dataset inventory. , Although all datasets meet the minimal requirements for DKI analysis, various scan and subject parameters vary widely. A comprehensive summary of the protocols and subjects is presented in Table 1—full protocol details can be found on each project’s resource pages.

TABLE 1

Overview of the scan parameters

	HCP	MASiVar	SUDMEX	CamCAN	Study‐specific data
Scan protocol
T_E (ms)	57	79	127	104	66
b‐values (ms/μm^2)	0, 1, 3	0, 1, 2	0, 1, 3	0, 1, 2	0, 0.5, 1, 2.5
# directions	5, 64, 64	16, 40, 56	8, 32, 96	3, 30, 30	5, 30, 30, 30
Voxel size (mm^3)	1.25×1.25×1.25	2.1×2.1×2.2	2×2×2	2×2×2	2.5×2.5×2.5
Subject information
Age (years)	40‐44	23	43	66	26
Gender	F	M	M	M	F
Preprocessing
	HCP 77	PreQual 87	Designer 52	Designer 52	see Section 2.6.2

Image processing

The study‐specific data were corrected for Gibbs ringing, eddy current and susceptibility‐induced geometrical distortions, , and gradient nonlinearities prior to DKI analyses. As a part of the eddy current correction, signal outlier detection to identify, for example, motion‐corrupted slices was performed. Seven regions of interest (ROIs) were automatically segmented in the native spaces of the images to minimize data interpolation. The ROIs included major white matter tracts, that is, genu and splenium of the corpus callosum (GCC and SCC), corticospinal tract (CST), arcuate fasciculus (AF), inferior fronto‐occipital fasciculus (IFO), superior longitudinal fasciculus (SLF), and optic radiation (OR). The automated segmentation of all white matter bundles was performed using the TractSeg algorithm.

If (minimally) preprocessed data were available for the cohort data, then we used those data sets without applying additional processing steps for DKI analyses. Otherwise, the DESIGNER pipeline was applied to correct for thermal noise, Gibbs ringing, and eddy current distortions. This information is also listed in Table 1.

Statistics

The test‐retest reliability of the estimated ROI‐averaged was evaluated for each fitting strategy using the test‐retest variability (TRV). Note that the test and retest data were not aligned to each other, since reproducibility was only evaluated for ROIs that were segmented in the native spaces of the images to minimize data interpolation. The TRV was computed across N subjects by averaging the ratio of the absolute difference and the average of the test‐retest estimates over the subjects. A scaling factor of was applied to derive the standard deviation from the mean absolute difference for Gaussian distributed variables.

RESULTS

In Figure 1, we show the scatter plot of predicted mean kurtosis against the corresponding in simulations and experimental data, using the 3 prediction strategies. In the simulations, is the ground truth value, whereas for the experimental data is estimated using an ordinary NLS estimator. For the experimental data, we only include all 441,271 voxels of the 5 study‐specific retest data with positive . The similarity between and (or ) is high. The voxelwise percentage error is and (mean ± standard deviation) for the simulations and all experimental data, respectively. The correction using the analytical correction term () reduces the accuracy significantly, leading to an increase of the mean percentage errors with an order of magnitude; and for the simulations and all experimental data, respectively. The third strategy, that is, correction using polynomial dictionary learning (), provides more accurate predictions of mean kurtosis than the other two. The mean percentage error is only for the simulated data and for the experimental data. Of importance, the polynomial regression model was trained on all voxels of the 5 study‐specific test data with positive . The Pearson correlation coefficients between and , , and was 0.70, 0.95, and 0.95, respectively, for the simulated data, and 0.82, 0.99, and 0.99, respectively, for the experimental data.

FIGURE 1

Simulations: a scatter plot shows the relationship between the actual mean and predicted mean kurtosis for the simulated and study‐specific data

In Figure 2, it is apparent that the predicted mean kurtosis is more robust than . Indeed, the prevalence of the ubiquitous “black voxels” in major white matter structures is reduced significantly, even almost nullified. The percentage of voxels in the whole brain with negative varies from 0.95 to 3.27%—with an average of 1.84%—in the ten study‐specific data sets. In contrast, for , this percentage dropped on average to 0.07%. Therefore, the quality of the maps is increased, allowing for a less confounded visual and statistical analysis. The robustness is similar for all 3 strategies to predict the mean kurtosis.

FIGURE 2

An representative map of the mean kurtosis and the predicted mean kurtosis ‐as estimated using the third strategy: voxel quality transfer using polynomial regression. The grayscale intensities are scaled between 0 and 2

Based upon the above results, we opt for the third strategy to compute the predicted mean kurtosis . Hence, for the remainder of this work, .

Accuracy of regularized tensor fitting

Figure 3 summarizes the simulation results by showing the ground truth reference of DKI parameters versus their estimates obtained using the ordinary (red) and regularized NLS (blue) estimator. The average percentage difference between the estimated and ground truth mean kurtosis is and for the ordinary and regularized NLS estimators. For the radial kurtosis, the average percentage difference decreases from to when regularizing the tensor fitting. For axial kurtosis, the percentage difference is minimally altered: and .

FIGURE 3

Simulations: Scatter plots show the relationship between the ground truth (“reference") and estimates of mean kurtosis (), radial kurtosis (), and axial kurtosis () when fitting the DKI model using the ordinary (red) and regularized (blue) nonlinear least squares estimator

For the study‐specific data, the ground truth kurtosis values are missing. However, for the majority of the gray and white matter voxels, the ordinary NLS estimator yields plausible and non‐problematic estimates with positive kurtosis values. In Figure 4A, we show the correspondence between estimated kurtosis parameters using the ordinary and regularized NLS estimator for those voxels. The average percentage difference between the ordinary and regularized NLS estimators for , , and is , , and , respectively. The Pearson correlation coefficients between the ordinary and regularized NLS are 0.995 0.993, and 0.995 for , , and , respectively. The effect of the regularized NLS estimator on the shape of the kurtosis tensor from voxels near the mid‐sagittal plane of the genus of the corpus callosum is highlighted in Figure 4B. The regularized NLS estimator shows to resolve the higher kurtosis values perpendicularly to the diffusion tensor main direction for both non‐problematic and problematic voxels.

FIGURE 4

Study‐specific MRI data: A, Scatter plots show the relationship between the ordinary and regularized NLS estimates of mean kurtosis (), radial kurtosis (), and axial kurtosis () for all non‐problematic voxels in the study‐specific MRI data. Note that a “problematic voxel” has at least one negative when fitted with the ordinary NLS. B, The 3D kurtosis geometry, ie, apparent kurtosis coefficient evaluated along different directions, as estimated with the ordinary and regularized NLS for a non‐problematic and problematic voxel. These 3D kurtosis geometry were extracted from a voxels near the mid‐sagittal plane of the genus of the corpus callosum (direction marked by the black line corresponds to the principal direction of the diffusion tensor and are truncated

Reproducibility

Figure 5 shows the mean within each ROI as computed in 5 subjects using various fitting strategies for the test and retest data. In addition, we include the median within each ROI for the ordinary NLS estimator.

FIGURE 5

The tract‐averaged using various fitting strategies and (bottom row) for the test (filled marker) and retest data (open marker) each subject (labeled by marker shape). The graphs on the right column show the same data, but windowed differently for enhanced contrast. Seven major white matter tracts were evaluated: genu and splenium of the corpus callosum (GCC and SCC), corticospinal tract (CST), arcuate fasciculus (AF), inferior fronto‐occipital fasciculus (IFO), superior longitudinal fasciculus (SLF), and optic radiation (OR)

In various tracts, that is, the SCC, GCC, and CST, the test‐retest reproducibility is poor when using ordinary NLS estimators, with and without spatial smoothing, because in at least one of the scans the negative kurtosis outliers dominate the ROI‐averaged . The median operator is robust to such outliers and can be used as a reference target for reproducibility analysis.

The test‐retest reproducibility is greatly improved by the use of the regularized fitting algorithm. The TRV of the ROI‐averaged metrics, computed over the 5 subjects, is tabulated per ROI in Table 2.

TABLE 2

The test‐retest variability in the estimation of for various tracts and fitting strategies

	GCC	SCC	CST	AF	IFO	SLF	OR
Ordinary NLS	‐	‐	‐	1.66	2.21	1.62	2.38
Median	4.20	2.19	2.17	1.93	2.60	1.68	2.03
Constrained NLS	5.30	2.89	2.16	1.75	1.75	1.59	1.99
Smoothing + ordindary NLS	‐	‐	‐	1.65	7.51	1.63	2.37
Regularized NLS	4.76	2.08	2.27	1.65	2.10	1.49	1.77
K^	4.73	2.08	2.28	1.65	2.10	1.49	1.76

Notes: If the test‐retest variability is dominated by outliers, that is, , the value is not listed. To illustrate the reproducibility of the predicted mean kurtosis , we also show those results (bottom row).

The percentage of negative estimates per tract strongly impacts the TRV. The regularized NLS estimator decimated, if not nullified, this number in comparison to the ordinary NLS estimator, with and without smoothing. Importantly, 80‐90% of the voxels with negative using the ordinary NLS estimator converged to plausible, positive, values with regularized NLS estimator if the optimization was initiated with the outcome of the ordinary NLS estimator. For the remaining 10‐20%, the regularized fit needed to be repeated with a starting point that was obtained from the constrained NLS estimator.

General applicability

In Figures 6 and 7, we show the maps of , and for a single, but representative slice of various data sets. The data quality varied highly from 1 dataset to another due to variations in SNR and number of diffusion‐weighted images. However, in each dataset, traditional kurtosis estimation suffered from low robustness and many negative in structures such as the CC. Overall, the quality of the maps improved drastically for all datasets, with a few remaining black voxels (negative kurtosis values), mostly around some GM/WM boundaries, where the kurtosis and their reference values are not well understood yet. In Figure 6, we also show, for comparison, the map of the predicted mean kurtosis .

FIGURE 6

The maps for the various dataset are shown for the ordinary and regularized NLS in the top and middle row, respectively. Moreover, we show the map of the predicted mean kurtosis (bottom row) to illustrate the similarity in contrast

FIGURE 7

The and maps for the various dataset are shown for the ordinary and regularized NLS in the top and bottom row of each panel

DISCUSSION

Improving the robustness of DKI fitting has been an active topic of research since the introduction of the technique. Black voxels have been intrinsic to DKI and challenge the visual and statistical analysis of potentially clinically relevant biomarkers of tissue integrity. Despite increasing evidence of the potential value of DKI during the past decade, the lack of a robust parameter estimation prevents widespread clinical/diagnostic adoption of the kurtosis biomarkers.

It might be impossible to attribute a single cause of the vulnerability of DKI to outliers. It has been shown that thermal noise alone can result in extremely negative kurtosis values. However, imaging artifacts, such as signal voids, Gibbs ringing, CSF pulsation, and/or misalignment of the diffusion‐weighted images might also contribute to the problem, even when the artifacts are largely corrected using state‐of‐the‐art image processing tools. , , The wide variety of signal fluctuations that might cause the DKI estimator to fail motivated the use of a strong image smoothing prior to parameter fitting and/or constraining the parameter estimation by enforcing positive kurtosis estimates. , , ,

The biophysically implausible kurtosis estimates have previously been associated with artifactually low non‐diffusion weighted signals. Motivated by this, Zhang et al recently proposed an approach in which the non‐diffusion‐weighted signal is altered in a selective data‐driven range. Although the approach yields visually appealing kurtosis maps, the alteration is strongly dependent on a user‐dependent tuning variable, while a selection criterion is missing.

However, for voxels with negative mean kurtosis in our reproducibility data, only of voxels have at least one diffusion‐weighted signal exceeding the non‐diffusion weighted signal, coined as “physically implausible signal (PIS)” , ; see Figure 8C. Moreover, only of all such PIS voxels resulted in a negative mean kurtosis. Hence, techniques focusing on the correction of the non‐diffusion‐weighted signal might be limited to a subset of the voxels only or bias the outcomes by forcing a correction on the non‐diffusion‐weighted signal instead of considering spurious signal fluctuations across the whole of the diffusion‐weighted data.

FIGURE 8

The map of the GCC of a single subject of the reproducibility data is shown for the various fitting strategies (A, B, E, F). The grayscale intensities are scaled between 0 and 2. In (C), the voxels in which at least one diffusion‐weighted signal was larger than the non‐diffusion‐weighted signal, “physically implausible signal (PIS)” are highlighted in red. In (D), the voxels that do not satisfy the Pearson’s goodness‐of‐fit criterion are indicated in red. Many problematic voxels in terms of plausibility of pass the PIS or goodness‐of‐fit criterion. G, A scatter plot shows the relationship between the as estimated using the constrained and regularized NLS, for the voxels that were characterized with a negative using the ordinary NLS estimator. The probability density estimate is color‐encoded

Moreover, when comparing the fit residuals of the ordinary and regularized NLS in the voxels of interest, that is, the voxels with negative mean kurtosis, we observed that the signal prediction of the regularized NLS estimated was significantly different for each b‐value with both positive and negative signal offsets. This observation (data not shown) also suggests that only correcting non‐diffusion‐weighted signals is likely to be insufficient for a robust and unbiased estimator.

Up to of all voxels with negative mean kurtosis within a single subject of our reproducibility data satisfied the Pearson’s goodness‐of‐fit criterion (Figure 8D shows in red voxels in which signals did not satisfied the Pearson’s goodness‐of‐fit for a zoomed brain region of the subject reproducibility data). Hence, despite the implausible kurtosis estimate, the fit residuals were normally distributed and no signal outliers were detected using a residual analysis. Robust estimators, such as RESTORE and related techniques, , would not make any difference to the outcome. Indeed, those methods are designed to improve the robustness of the fit to signal outliers, but are not necessarily tackling the problem described here.

Overall, smoothing has proven to be a very effective way to suppress the signal fluctuations that might lead to negative kurtosis. , , , However, such a brute force strategy introduces image blur and partial volume artefacts. Analogous to other artifacts, for example, Gibbs ringing, , , there is a recent trend toward the development of image processing tools that are more specifically targeted to correct a particular artifact and to avoid smoothing. The development of a robust kurtosis estimator that does not require image smoothing prior to fitting is more in line with such a strategy and, ultimately, will result in sharper images with higher effective spatial resolution.

The constrained NLS estimator only yields positive kurtosis values by design. Hence, the constrained NLS typically appears to be more accurate and robust in the estimation of DKI parameters. However, imposing constraints is likely to bias the estimator in any single voxel. In of all voxels with negative using the ordinary NLS estimator, the constrained NLS estimator provides a solution that lies on the bounds of the search landscape. Although such estimates technically satisfy the imposed constraints, they are not necessarily closer to the most biophysically plausible solution, see Figure 8.

All of these considerations are visualized and summarized in one representative figure that shows the Genu of the CC of one the study‐specific MRI data sets, Figure 8.

In this work, we demonstrated that the mean kurtosis can be very well predicted without the need to estimate the kurtosis tensor. Indeed, the powder kurtosis is already a good approximation of the mean kurtosis, despite a significant approximation error that depends on the anisotropy of the underlying process. We showed that this approximation error can be reduced significantly by using a polynomial regression model in which the mapping between the powder kurtosis, the mean diffusivity, an anisotropy index, and the mean kurtosis can be learned from the hundreds of thousands of non‐problematic voxels in the same or similar data sets. By doing so, we transfer the excellent quality of the DKI estimation in the majority of the gray and white matter voxels to the few, but persistently problematic voxels.

The accurate and robust prediction of the mean kurtosis might be relevant and sufficient for various studies. However, the regularized fitting is a necessary additional step for all studies that have interest in directional kurtosis metrics (eg, radial or axial kurtosis), DKI‐derived biophysical modeling (ie, WMTI ), or tractography. ,

Similar to the negative kurtosis, various models developed to analyze diffusion MRI data suffer from an abundance of biophysically implausible outcomes or, more subtly, a multitude of biophysically plausible solutions for which the goodness‐of‐fit is not significantly different. Although beyond the scope of this work, we hypothesize that using the regularization term presented in this work might increase the stability and precision of parameter estimators of such models. Similarly, the regularization term can be added to various variations of the loss function, including the log‐linearized or maximum likelihood function.

The implementation of this regularized estimator is publicly available on https://github.com/jelleveraart/RobustDKIFitting/.