A robust deconvolution method to disentangle multiple water pools in diffusion MRI.

The diffusion-weighted magnetic resonance imaging (dMRI) signal measured in vivo arises from multiple diffusion domains, including hindered and restricted water pools, free water and blood pseudo-diffusion. Not accounting for the correct number of components can bias metrics obtained from model fitting because of partial volume effects that are present in, for instance, diffusion tensor imaging (DTI) and diffusion kurtosis imaging (DKI). Approaches that aim to overcome this shortcoming generally make assumptions about the number of considered components, which are not likely to hold for all voxels. The spectral analysis of the dMRI signal has been proposed to relax assumptions on the number of components. However, it currently requires a clinically challenging signal-to-noise ratio (SNR) and accounts only for two diffusion processes defined by hard thresholds. In this work, we developed a method to automatically identify the number of components in the spectral analysis, and enforced its robustness to noise, including outlier rejection and a data-driven regularization term. Furthermore, we showed how this method can be used to take into account partial volume effects in DTI and DKI fitting. The proof of concept and performance of the method were evaluated through numerical simulations and in vivo MRI data acquired at 3 T. With simulations our method reliably decomposed three diffusion components from SNR = 30. Biases in metrics derived from DTI and DKI were considerably reduced when components beyond hindered diffusion were taken into account. With the in vivo data our method determined three macro-compartments, which were consistent with hindered diffusion, free water and pseudo-diffusion. Taking free water and pseudo-diffusion into account in DKI resulted in lower mean diffusivity and higher fractional anisotropy values in both gray and white matter. In conclusion, the proposed method allows one to determine co-existing diffusion compartments without prior assumptions on their number, and to account for undesired signal contaminations within clinically achievable SNR levels.

apparent diffusion coefficient

cerebrospinal fluid

diffusion kurtosis imaging

diffusion‐weighted magnetic resonance imaging

diffusion tensor imaging

fractional anisotropy

free water

gray matter

maximum gradient output

Gaussian mixture model

intra‐voxel incoherent motion

non‐negative least squares with L2 regularization

mean diffusivity

Marchenko-Pastur Principal component analysis

non‐linear least squares

non‐negative least squares

non‐negative least squares with data‐driven L2 regularization

sensitivity encoding

signal‐to‐noise ratio

T 1‐weighted/T 2‐weighted

echo time

repetition time

white matter

INTRODUCTION

Diffusion‐weighted magnetic resonance imaging (dMRI) is an imaging technique capable of estimating the diffusion displacement of water molecules in vivo. Although several models have been proposed to quantify parameters from the dMRI signal, the apparent diffusion coefficient (ADC)1 model and its generalization, the diffusion tensor imaging (DTI) model2, are still among the most popular choices, especially in clinical applications. These approaches assume that the diffusion signal can be described by a single Gaussian process. However, at the imaging resolution typically achieved in a dMRI experiment – on the order of millimeters – this assumption is often violated.3, 4, 5

The vast majority of studies on multiple diffusion components are based on multi‐compartment models, which require prior knowledge of the number and structure of the compartments. Originally, Le Bihan et al.6 proposed a bi‐exponential formulation of the isotropic diffusion signal defining the intra‐voxel incoherent motion (IVIM) model, which quantifies the effect of blood pseudo‐diffusion on the dMRI signal. Similarly, Pierpaoli et al.7 and Pasternak et al.8 extended the DTI model with a fixed isotropic compartment to account for the intra‐voxel contamination of free water (FW). Likewise, other models9, 10, 11, 12, 13, 14 have been proposed to quantify the intrinsic diffusion parameters of the tissues, whilst accounting for partial volume effects. However, such models generally rely on restrictive physiological assumptions and do not include the pseudo‐diffusion and FW components.

The pseudo‐diffusion component has been investigated previously in conjunction with ADC,15 DTI16 and diffusion kurtosis imaging (DKI).17 The extension of classic models with pseudo‐diffusion and FW compartments has provided valuable biomarkers in several applications.18, 19, 20 Although such partial volume effects are an important bias in DTI and DKI analysis, taking them into account results in improved signal characterization, reduced estimation biases in DTI and DKI21, 22 metrics, and additional specificity, as it potentially allows intra‐ and extracellular‐related changes to be disentangled. Nonetheless, the design of models to account for the correct number of components is challenging, as these are voxel and acquisition dependent. For this reason, the number of components needed to appropriately describe the dMRI signal remains a topic of active research.14, 23, 24

Recent work22 has shown that both the pseudo‐diffusion and FW compartments should be considered when including the non‐weighted image in the model fitting. This suggests the need for at least three compartments for an accurate disentanglement of partial volume effects in the brain, which is in line with other studies on brain25 and non‐brain26, 27, 28 tissue. However, the appropriate fit of multiple exponentials with non‐linear least‐squares (NLLS) approaches is a difficult and time‐consuming numerical problem.8 A number of studies22, 25, 29, 30 have proposed the segmentation of the fit, i.e. to determine part of the parameters on a subset of the data and to fix their values in the global fit. Although effective in stabilizing the fit, this assumes that the contribution of some compartments is negligible in subsets of the data, which might not hold in experimental conditions.

Other studies have proposed methods that do not explicitly impose the number of compartments, but rather derive them from the data.31, 32, 33, 34 Although prone to overfitting,35 such methods have been shown to provide good results in problems in which the number of components is difficult to establish a priori, such as crossing fibers36, 37, 38 or the joint diffusion–T 2 relaxometry quantification.39

Recently, Keil et al.40 have proposed to fit the data acquired at multiple diffusion weightings in three orthogonal directions with a deconvolution method to separate pseudo‐diffusion from hindered diffusion in tissues. Different from other approaches, their solution decomposes the dMRI signal into an arbitrary number of compartments, which are then grouped into macro‐compartments. This method results in more physiologically plausible estimates of pseudo‐diffusion compared to a bi‐exponential model fit using an NLLS estimator. However, in their spectral analysis method, only hindered diffusion and blood pseudo‐diffusion were considered in combination with predefined hard thresholds. Moreover, their deconvolution method requires high signal‐to‐noise ratio (SNR) despite regularization.41

In this work, we propose a spectral analysis framework that automatically disentangles diffusion domains beyond pseudo‐diffusion, whilst enforcing robustness to noise and to spurious measurements. In addition to our methodological innovations, we show how this framework can be used to robustly estimate FW and pseudo‐diffusion partial volume effects. Further, we show how accurate estimation allows the mitigation of their effects in DTI and DKI42 analysis.

The reliability of our proposed method is investigated using simulations and in vivo MRI data. First, simulations are used to compare the methods with regard to their ability to correctly estimate signal fractions and diffusion coefficients at different SNR levels. Second, a numerical phantom simulating a human brain is employed to evaluate whether the proposed methods can mitigate biases in DTI and DKI analysis as a result of the presence of co‐existing diffusion processes. Finally, our proposed method is evaluated with in vivo multi‐shell dMRI data and is compared with other state‐of‐the‐art approaches.

METHODS

In addition to the background and theory presented in Keil et al.,40 the following sections present the main contributions of this work, which include a modified non‐negative least‐squares (NNLS) regularization term based on a data‐driven prior, an iterative NNLS estimator with outlier rejection and a Gaussian mixture model (GMM) to automatically subdivide the NNLS estimation into multiple diffusion components without prior assumptions. Figure 1 shows a schematic representation of the method, which is detailed in the following paragraphs.

Figure 1

Flow chart of the proposed method (non‐negative least squares with data‐driven L2 regularization, P‐L2‐NNLS) from data acquisition to parametric maps. The diffusion‐weighted magnetic resonance imaging (dMRI) data acquired at multiple b values (a) are deconvolved with a signal dictionary (b) after geometric averaging (c) at two different precision levels (ε). The low‐precision solution (d) is voxel‐wise averaged and fitted with a Gaussian mixed model (GMM) to detect the number of macro‐compartments (f). The high‐precision solution (e) is voxel‐wise averaged to define the prior p used for the P‐L2‐NNLS fit (h). The operation is performed voxel‐wise (i) to compute the signal fraction (f) and diffusion coefficient (D) maps

Dictionary‐based representation of the signal

Assuming an arbitrary number of Gaussian components (n) with diffusion tensor , the diffusion signal S measured at diffusion weighting b along gradient direction g can be expressed as: where is the non‐weighted signal and is the signal fraction of each component.43

When only the trace of ( is of interest, Equation (1) can be simplified considering the geometric average of the data at each diffusion weighting b (Figure 1a–c). This expression can be restructured using a two‐dimensional dictionary (W), in which each column represents an isotropic mono‐exponential decaying signal with different diffusion coefficients D 1,…,, as shown in Figure 1b: where b 1,…, are the unique applied diffusion weightings (i.e. the b values). Similar to Equation (1), the signal for this model can be expressed in a linear matrix formulation: In this notation, vector p represents the contribution of each diffusion component (columns of W) to the measured signal. In addition to acting as a basis function for the linear decomposition of S, W can also be seen as a discrete pseudo‐Laplace operator, that transforms the diffusion‐weighted signal from the space of b values to the space of diffusion coefficients (Figure 1d,e). Therefore, following Keil et al.,40 we refer to p as the “diffusion spectrum”.40, 44

Robust NNLS fit

As p has only positive values, a non‐negative least‐squares estimator can be employed to efficiently solve the deconvolution operation (or the inverse pseudo‐Laplace transform) in its standard (NNLS) or L2‐regularized (L2‐NNLS) formulation45: The regularization term in Equation (5) enforces the robustness to noise, but penalizes sparse solutions, therefore smoothing the solution p with increasing values of .

When prior knowledge is available on the distribution of p, Equation (5) can be modified to penalize the divergence from such a prior, leading to prior enforced L2‐NNLS (P‐L2‐NNLS): where is the diffusion spectrum prior. In this article, we used a data‐driven prior p 0, which was derived independently for each dataset by averaging the voxel‐wise diffusion spectrum obtained with Equation (4), and thus represents the observed probability of each diffusion component (Figure 1 g) within each dataset. The L2 regularization term in Equation (6) also depends on the parameter , whose effect on the solution p is further explained in Supporting Information.

Deconvolution methods are known to be sensitive to outliers; therefore, we implemented an iterative deconvolution approach to discard erroneous data points from the fit. For this we used a framework that alternates between an heteroscedastic and an homoscedastic fit of the dMRI signal, which was previously introduced in REKINDLE.46 The weights used for the homoscedastic fit were iteratively adapted based on the model estimates as proposed in Veraart et al.47. Once outlier points had been identified and removed, the final fit was performed with Equations (4)–(6) (Figure 1e).

In the following analysis, the regularization parameter was set to for P‐L2‐NNLS and for L2‐NNLS to avoid excessive smoothing, but to maintain the robustness to noise.

Diffusion spectrum subdivision (GMM)

For each component i of the diffusion spectrum, p (Figure 1e) is the contribution of a spectral component with diffusion coefficient D to the measured signal. The direct interpretation of p can therefore be difficult when the spectrum is non‐sparse, as a result of noise, regularization or physiological complexity.

For this reason, Keil et al.40 proposed to divide the spectrum into two macro‐compartments (Cx), one describing hindered diffusion for diffusion coefficients between 0.5 × 10−3 and 3 × 10−3 mm2/s, and a second describing blood pseudo‐diffusion for diffusion coefficients greater than 3 × 10−3 mm2/s. To avoid prior assumptions, we propose an automated method to define the number of macro‐compartments and their thresholds. Assuming that, in a noisy experiment (both acquisition noise and biological noise), the diffusivities of each compartment are distributed around a common value, we can use a GMM to identify the number of biological “macro”‐compartments from the voxel‐wise spectral components.

However, in noisy experiments, the variance of each spectrum may be too large to sufficiently separate the macro‐compartments. For this reason, a fast fit of the data was performed with Equation (4), limiting the tolerance of the NNLS solver (lsqnonneg, MATLAB, The Mathworks, Natick, Massachussets, USA) to 10−2 (Figure 1d), which in practice resulted in a less accurate and noise‐sensitive solution. This fit was then voxel‐wise averaged over the entire brain and fitted with a GMM (Figure 1f) to determine the number and location of compartments, where the maximum number was limited to 300. Individual Gaussian components with area overlap above 15% were merged, after which the voxel‐wise diffusion spectrum was subdivided into the resulting macro‐components (Figure 1h). The number of resulting macro‐compartments is dependent on the overlap setting. An additional explanation on the area overlap setting and how it can affect the results can be found in Supporting Information.

Combining the spectral components of each voxel in macro‐compartments allowed us to sufficiently approximate the non‐Gaussian distributions of each component. Finally, for each macro‐compartment Cx, the signal fraction f(Cx) was defined as the sum of the spectral components p in the specified range. The diffusion coefficient D was defined as the average of the diffusion coefficients D weighted by their amplitude p . According to this formulation, the value of D represents a summary statistic of the spectral components included in each macro‐compartment.

Simulations

With simulations, the performance of P‐L2‐NNLS was compared with that of NNLS and the regularized L2‐NNLS approach. The first simulation shows the robustness to noise of the three methods. The second simulation investigates P‐L2‐NNLS and L2‐NNLS using a realistic brain phantom. Further details are provided in the following sections.

Simulation I

The purpose of this simulation was to show that L2 regularization, in conjunction with a prior obtained from NLLS fitting, can further improve the robustness of the method proposed by Keil et al.40 Data for simulation I were based on an acquisition scheme used for the in vivo MRI data (see Section 2.5). The signal was generated as the sum of three isotropic exponential decaying signals with respective fractions of 0.7, 0.2 and 0.1 and corresponding diffusion coefficients of 0.7 × 10−3, 3.0 × 10−3 and 200.0 × 10−3 mm2/s. These diffusion coefficients were chosen to resemble a mixture of hindered tissue diffusion, FW and blood pseudo‐diffusion, respectively, with typical literature values. A thousand realizations of Rician noise were added to the synthesized signal to simulate SNR levels of 10, 20, 30, 50, 70 and 90, with SNR defined as the average over standard deviation of the signal value corresponding to b = 0 s/mm2.

For each noise realization, the diffusion spectrum was computed with NNLS and L2‐NNLS (Equations (4) and (5)). Next, the average spectrum from NNLS was used as a prior for the P‐L2‐NNLS fit (Equation (6)). For this simulation, we set the regularization parameter to , which is a lower value than that used in Keil et al.40 The effect of on P‐L2‐NNLS estimates is further explained in Supporting Information. The number of macro‐compartments was automatically derived from the spectrum with the GMM fit, and then the 25th, 50th and 75th percentiles of the signal fractions and diffusion coefficient associated with each macro‐compartment were computed. The relative distance of the median from the predefined values was derived for each SNR level as a quantitative measure of error.

Simulation II

An in vivo MRI acquisition (see Section 2.5) was used as a template for this simulation. Its T 1‐weighted (T1W) data and the corresponding probabilistic segmentation into cerebrospinal fluid (CSF), gray matter (GM) and white matter (WM), derived with FSL FAST,48 were non‐linearly aligned to the corresponding dMRI space.49, 50 The derived tissue segmentations were used to simulate voxel‐wise partial volume effects plausible for the three classes. For each voxel, the diffusion tensor model was fitted to the data at b = 5 and 1000 s/mm2,46 from which the main eigenvector was derived. The dMRI signal of each tissue class (S ) was generated as22: where b is the b value, g is the corresponding gradient orientation, and and refer to rotation matrices around the y and z axes to align the tensor with V 1. In this model, only the hindered diffusion component was assumed to be anisotropic.

The diffusivity coefficients associated to FW (D FW) and to pseudo‐diffusion (D*) were set to × 10−3 mm2/s and × 10−3 mm2/s, respectively, with indicating a normal distribution with mean and standard deviation . The signal fractions , as well as the diffusion tensor ( refers to its average value) and the mean kurtosis , were predefined differently in GM and WM to take into account the physiological variability normally observed in the brain. In particular, values were set as follows: , , .

CSF: , , ;

GM: and , , ;

WM: and , where I is a 3 × 3 identity matrix;

An additional phantom with K = 0 for all components was generated. Finally, voxel‐wise realizations of Rician noise were added to simulate an SNR of 30 and 100 for the signal at b = 0 s/mm2. Rician noise was preferred over Gaussian noise to account for signals at high b values, given that, in voxels characterized by D ≥ 1.1 × 10−3 mm2/s, the signal at b = 2500 s/mm2 would be within two standard deviations of the noise floor.

The dictionary W was built using 300 mono‐exponential isotropic Gaussian signals with diffusion coefficients log‐spaced in the diffusion range DR = [0.1, 1000] × 10−3 mm2/s. The simulated data were voxel‐wise geometric averaged for each b value and fitted with the L2‐NNLS and P‐L2‐NNLS methods. In line with simulation I, the voxel‐wise diffusion spectrum computed with NNLS was averaged and set as the prior p 0 used for the P‐L2‐NNLS fit.

After the GMM fit had been performed, the signal fraction and the diffusion coefficients of each macro‐component were derived for L2‐NNLS and P‐L2‐NNLS in each voxel. The agreement between predefined and recovered signal fractions was evaluated with scatterplots including a line of equality.

To evaluate the performances of the model at mitigating partial volume biases, data were fitted47 with the DTI (b = 0, 500, 750, 1000 s/mm2) and DKI (with isotropic kurtosis K 21, b = 0, 500, 750, 1000, 1750, 2500 s/mm2) models, taking into account macro‐compartments larger than hindered diffusion. Subsequently, the mean diffusivity (MD), fractional anisotropy (FA, only in WM) and K (DKI only), before and after correction, were compared with the ground‐truth values.

MRI data

The feasibility of P‐L2‐NNLS on in vivo data and its performances were assessed on a clinically compatible MRI dataset. The data were acquired with a Philips Achieva 3 T (dStream, Philips Medical System NV, Best, The Netherlands) at IRCCS Eugenio Medea (Bosisio Parini, Italy) using a Stejskal–Tanner spin‐echo sequence and a 32‐channel head coil. The acquisition parameters were as follows: echo time (TE) = 80 ms; repetition time (TR) = 6.9 s; sensitivity encoding (SENSE) = 2.5; voxel size, 2.5 × 2.5 × 2.5 mm3; G max = 66 mT/m; raw matrix, 96 × 96.

The acquired gradient scheme, also used in simulations, featured a total of 124 volumes, with six volumes at b values of 0, 1, 5, 10, 20, 50, 80, 110, 150, 200, 250, 500 and 750 s/mm2, 15 volumes at b = 1000 s/mm2, 15 volumes at b = 1750 s/mm2 and 15 volumes at b = 2500 s/mm2, for a total acquisition time of around 15 minutes.

In addition to dMRI, a T1W image (TE/TR = 3.7/8.3 ms; SENSE = 2; voxel size, 1 × 1 × 1 mm3) and a T2W image (TE/TR = 80 ms/4.2 s; SENSE = 1.9; resolution, 1.4 × 1.4 × 1.4 mm3) were acquired. Preprocessing of dMRI data included eddy currents and motion correction51 with ExploreDTI (http://exploredti.com).52 The T1W volume was preprocessed according to the standard FSL anatomic pipeline fsl_anat 53 to derive probabilistic tissue‐type segmentation into CSF, GM and WM.48 The first non‐weighted volume of each dMRI dataset was non‐linearly registered to the T1W image54 using the T2W image as midpoint,50 and then the two transformations were combined to co‐register GM and WM segmentation to dMRI space. The standard deviation of noise at b = 0, 1000, 1750 and 2500 s/mm2 was determined for each dataset with the Marchenko‐Pastur Principal component analysis (MP‐PCA) method55; then, the SNR of each diffusion‐weighted volume was computed as the average signal of the brain tissue divided by the standard deviation of noise.

Data were fitted as in simulation II. The average whole‐volume diffusion spectra of each subject obtained with L2‐NNLS and P‐L2‐NNLS were computed. In addition, statistics (average and standard deviation) on the number of non‐zero spectral components in the voxel‐wise P‐L2‐NNLS spectra were computed for each subject. The evaluation of the effect of the regularization parameter and of the overlap threshold of the spectral components and resulting macro‐compartments for the in vivo data can be found in Supporting Information.

To show the robustness of the proposed method, we also performed the analysis on a subset of the data corresponding to 5‐minute acquisition time, which is shown in Supporting Information.

P‐L2‐NNLS was compared to a three‐exponential compartment model, which represents the most popular alternative for separating hindered diffusion from FW and pseudo‐diffusion.22 The three‐exponential model was fitted with NLLS and with a segmented NLLS fit, as suggested previously.25, 56 Then, the average signal fraction and diffusion coefficients were computed in GM and WM for comparison with P‐L2‐NNLS.

Fractional maps derived from P‐L2‐NNLS were compared with those derived from L2‐NNLS. In addition, average values and standard deviations of signal fractions and diffusion coefficients were evaluated in GM and WM. The correction for perfusion and FW contamination was evaluated with a DKI fit (with isotropic K 21) of the data at b = 0 s/mm2 and b ≥ 500 s/mm2 with a weighted least‐squares approach taking into account the spectral components that are faster than hindered diffusion. The average value ± standard deviation of MD, FA and K before and after correction were evaluated.

RESULTS

Simulation I

Figure 2 reports the results of simulation I, which compared the 25th, 50th and 75th percentiles of the D and f values estimated with L2‐NNLS and P‐L2‐NNLS at SNR levels between 10 and 90. The GMM fit automatically subdivided the average spectrum into three macro‐components, which were approximately centered on the ground‐truth values.

Figure 2

Performances of L2‐NNLS (non‐negative least squares with L2 regularization, red) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization, blue) in disentangling signal mixtures at different signal‐to‐noise ratio (SNR) levels (simulation I). Columns 1–3 correspond to each of the three simulated diffusion components. Error bars represent the 25th, 50th and 75th percentiles of the estimated values of the diffusion coefficients (first row) and signal fractions (second row) from 1000 realizations, and the black broken line represents the predefined value. P‐L2‐NNLS resulted in the most consistent estimates of compartments 1 and 2 from SNR = 30. All three methods had similar performances on component 3

According to the simulation results, all methods asymptotically converged to the ground truth. However, P‐L2‐NNLS resulted in less biased and dispersed estimates from SNR = 30, with relative errors equal to 0.13% and 0.80% for f(C1) and D(C1), respectively, compared with 1.93% and 7.77% for L2‐NNLS. P‐L2‐NNLS outperformed L2‐NNLS in terms of the estimation of f(C2) and f(C3), with respective errors at SNR = 30 equal to 0.18% and 1.51%, compared with 6.99% and 2.50% for L2‐NNLS. The simulation also showed that compartments 2 and 3 were characterized by high uncertainty with L2‐NNLS for SNR below 70. D(C2) and D(C3) were characterized by higher uncertainty than D(C1), as visible in Figure 1.

Simulation II

The first row of Figure 3(a,b,e,f) shows the average diffusion spectrum obtained after a voxel‐wise L2‐NNLS fit of the phantoms without and with normally distributed kurtosis at SNR = 30 and 100. At SNR = 100, peaks were observed in concordance with the ground‐truth values D = 0.7 × 10−3, 3.0 × 10−3 and 100 × 10−3 mm2/s, with an additional peak at D = 0.8 × 10−3 mm2/s when kurtosis was simulated. At SNR = 30, the same peaks were observed, but part of the component at D = 0.7 × 10−3 mm2/s was decomposed into two peaks at D = 0.1 × 10−3 and 0.9 × 10−3 mm2/s. The second row of Figure 3(a,b,e,f) shows the same spectrum with P‐L2‐NNLS, which was sparser. The volume‐averaged spectrum was subdivided by GMM into three compartments with division points of 1.75 × 10−3 and 11.14 × 10−3 mm2/s (red broken lines in Figure 3(a,b,e,f)). Figure 3(c,d,g,h) shows the ground‐truth signal fractions and those computed with the two tested methods for an example slice of the same phantoms. L2‐NNLS resulted in visibly higher values of f(C1) and lower values of f(C2) compared with those from P‐L2‐NNLS, which were conversely closer to the ground‐truth values.

Figure 3

(a,b,e,f) The average deconvolution spectrum (full black line) from the voxel‐wise non‐negative least squares (NNLS) fit of the phantom without kurtosis (a,e) and with normally distributed kurtosis (b,f) at signal‐to‐noise ratio (SNR) = 30 (a,b) and SNR = 100 (e,f) with L2‐NNLS (non‐negative least squares with L2 regularization, first row) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization, second row). The broken red lines represent the separation points of the three macro‐compartments as determined by the Gaussian mixture fit. (c,d,g,h) Simulated fraction maps for macro‐compartments 1–3 (first column) and maps estimated with L2‐NNLS and P‐L2‐NNLS on the phantom without kurtosis (c,g) and with normally distributed kurtosis (d,h) at SNR = 30 (c,d) and SNR = 100 (g,h). P‐L2‐NNLS resulted in better estimation of C1 and C2 than L2‐NNLS, whereas C3 was comparable

Figure 4 shows the scatter plot analysis performed on the phantoms at SNR = 30 and 100 with and without kurtosis. With K = 0 (left), the relation between predefined and estimated values was overall linear (red broken line) for C1 and C2, with accumulation points corresponding to very small fractions. Values from P‐L2‐NNLS were closer to the ideal line for these two compartments. No visible differences were noticed for C3 among the two methods. The presence of normally distributed kurtosis (right) resulted in slightly higher dispersion for both methods.

Figure 4

Scatter plots of simulated versus estimated values (blue dots) with L2‐NNLS (non‐negative least squares with L2 regularization) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization) on the synthetic phantoms without (left) and with (right) kurtosis at signal‐to‐noise ratio (SNR) levels of 30 and 100. The red broken line shows the ideal agreement. Signal fractions of C1 and C2 estimated with P‐L2‐NNLS were closer to the ideal line for both SNR levels, with and without kurtosis. Values of f(C3) were similarly estimated by the two methods

The DTI fit of the phantom with K = 0 without any correction resulted in MD overestimation of 26.4% in GM and 13.9% in WM, and an FA underestimation of around 8% in WM, independent of SNR. When C2 and C3 were taken into account with L2‐NNLS, biases in MD at SNR = 30 and 100 were reduced to 11.94% and 7.14% in GM, respectively, and 8.48% and 5.67% in WM, respectively. However, P‐L2‐NNLS resulted in the smallest errors, which, at SNR = 30, were 8.03% and 6.50% in GM and WM, respectively, for MD and 4.08% in WM for FA. Similar results for FA and MD were observed for the DKI fit with normally distributed K, as reported in Table 1. In addition, K overestimation in GM was greatly reduced by taking into account C2 and C3 with P‐L2‐NNLS before DKI fitting, but only small improvements were noticed in WM.

Table 1

Percentage errors of diffusion tensor imaging (DTI) and diffusion kurtosis imaging (DKI) estimates computed in gray matter (GM) and white matter (WM) at two signal‐to‐noise ratio (SNR) levels (simulation II). Errors were computed when no correction was applied and after removal of macro‐components C2 and C3 with L2‐NNLS (non‐negative least squares with L2 regularization) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization). Signal filtering with P‐L2‐NNLS resulted in the smallest errors

	Correction	GM		WM
		SNR = 30	SNR = 100	SNR = 30	SNR = 100
Simulated K = 0	DTI MD (simulated value 0.77 × 10−3 mm2/s)
	Not corrected	26.41	26.44	13.90	13.87
	C2 + C3 (L2‐NNLS)	11.94	7.14	8.48	5.67
	C2 + C3 (P‐L2‐NNLS)	8.03	1.01	6.50	1.08
	DTI FA (simulated value 0 in GM, 0.80 in WM)
	Not corrected	NA	NA	8.02	8.13
	C2 + C3 (L2‐NNLS)	NA	NA	5.10	3.63
	C2 + C3 (P‐L2‐NNLS)	NA	NA	4.08	0.93
Normally distributed K	DKI MD (simulated value 0.77 × 10−3 mm2/s)
	Not corrected	40.62	40.62	20.08	20.14
	C2 + C3 (L2‐NNLS)	17.37	13.83	13.44	10.72
	C2 + C3 (P‐L2‐NNLS)	15.88	4.51	9.63	8.22
	DKI FA (simulated value 0 in GM, 0.80 in WM)
	Not corrected	NA	NA	10.96	10.89
	C2 + C3 (L2‐NNLS)	NA	NA	7.95	6.25
	C2 + C3 (P‐L2‐NNLS)	NA	NA	5.65	5.03
	DKI K (simulated value 0.2 in GM, 0.8 in WM)
	Not corrected	216.10	217.70	14.54	13.86
	C2 + C3 (L2‐NNLS)	133.26	102.54	19.09	9.82
	C2 + C3 (P‐L2‐NNLS)	122.86	55.35	11.83	11.24

FA, fractional anisotropy; MD, mean diffusivity; NA, not applicable.

Average SNR ± standard deviation values of the acquired data were 49 ± 7 at b = 0 s/mm2, 16 ± 2 at b = 1000 s/mm2, 8 ± 3 at b = 1750 s/mm2 and 5 ± 1 at b = 2500 s/mm2. MRI data from the two time points of each subject were fitted with Equations (4)–(6). Figure 5 shows the average voxel‐wise diffusion spectrum computed with NNLS and P‐L2‐NNLS for each time point. The NNLS spectra were similar between the two time points for each subject, and comparable across subjects, with peaks in the NNLS spectrum at diffusion coefficients of (0.92 ± 0.06) × 10−3, (2.47 ± 0.19) × 10−3 and (57.64 ± 3.55) × 10−3 mm2/s. The broken red lines shown in Figure 5 represent the automatically computed subdivisions of the diffusion spectrum into three compartments at (1.72 ± 0.03) × 10−3 and (13.03 ± 1.10) × 10−3 mm2/s, which were consistent across acquisitions. The spectra computed with P‐L2‐NNLS were sparser and sharper compared to those from NNLS. The average number of non‐zero components in the diffusion spectrum of all subjects was 18, with intra‐subject standard deviation equal to 9.39 and inter‐subject standard deviation equal to 0.69.

Figure 5

Average whole‐brain diffusion spectra (full black line) obtained with L2‐NNLS (non‐negative least squares with L2 regularization, top) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization, bottom) of subjects 1–3 (rows S1–S3) at time points 1 (TP1, left) and 2 (TP2, right). The red broken lines represent the separation thresholds between the three macro‐compartments automatically detected by the Gaussian mixed model (GMM) fit, which were consistent across subjects. The second main peak of all subjects was very close to the diffusion value of free water at 37°C (3 × 10−3 mm2/s, black asterisk). P‐L2‐NNLS resulted in sharper peaks

The results of the NLLS and segmented NLLS fit of one subject are shown in Figure 6. Voxel‐wise maps were very similar among the two NLLS methods, but different from P‐L2‐NNLS. Considering the average values in GM and WM, which are reported in Table 2, f(C1) estimated with the NLLS methods was around 70% lower than with P‐L2‐NNLS, whereas values of f(C2) were noticeably higher (around 20%) for the two NLLS methods. Compared with P‐L2‐NNLS, values of f(C3) obtained with NLLS were also higher (close to 10%).

Figure 6

Signal fractions (f, columns 1, 3 and 5) and diffusion coefficients (D, columns 2, 4 and 6) obtained using a three‐compartment model fit with NLLS (non‐linear least squares, columns 1 and 2), a segmented NLLS fit (columns 3 and 4) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization, columns 5 and 6). The two NLLS fits resulted in higher estimates of f(C2) (free water, FW) compared with P‐L2‐NNLS. Differences between the method proposed in this work and NLLS were also visible in the pseudo‐diffusion range (C3), with higher diffusion coefficient for the latter, and in C1, which was characterized by a lower signal fraction

Table 2

Average values ± standard deviation of the signal fractions (f) estimated with a tri‐exponential segmented non‐linear least‐squares (NLLS) fit (SegTriExp), a tri‐exponential NLLS fit (TriExp), L2‐NNLS (non‐negative least squares with L2 regularization) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization)

Tri‐exponential NNLS
	GM		WM
	SegTriExp	TriExp	SegTriExp	TriExp
f(C1)	0.69 ± 0.01	0.70 ± 0.01	0.69 ± 0.02	0.71 ± 0.02
f(C2)	0.21 ± 0.02	0.20 ± 0.02	0.26 ± 0.01	0.22 ± 0.01
f(C3)	0.10 ± 0.02	0.10 ± 0.02	0.06 ± 0.01	0.08 ± 0.01
D(C1) (mm2/s × 10−3)	0.59 ± 0.01	0.61 ± 0.02	0.45 ± 0.10	0.50 ± 0.01
D(C2) (mm2/s × 10−3)	NA	NA	NA	NA
D(C3) (mm2/s × 10−3)	549.40 ± 66.30	496.50 ± 57.30	776.70 ± 50.80	714.30 ± 43.30

GM, gray matter; NA, not applicable; WM, white matter.

Example images of the signal fractions and diffusion coefficients computed for each class with L2‐NNLS and P‐L2‐NNLS are shown in Figure 7. Values of f (C1) were closer to unity when estimated with L2‐NNLS relative to P‐L2‐NNLS, whereas values of f(C2) were almost zero in brain tissue voxels for the former, in line with simulation II. Maps of macro‐compartments 1–2 from P‐L2‐NNLS were less noisy and smoother than with L2‐NNLS, whereas f(C3) and D(C3) did not visibly differ between the two methods. The average values and standard deviations observed for the two methods are reported in Table 2.

Figure 7

Voxel‐wise estimates of signal fractions (f, rows 1–3) and diffusion coefficients (D, rows 4–6) provided by L2‐NNLS (non‐negative least squares with L2 regularization, left) and P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization, right) fitting of the dataset. Maps derived from P‐L2‐NNLS appeared less noisy and more spatially homogeneous. f(C2) from L2‐NNLS had lower values than those from P‐L2‐NNLS, consistent with simulation II. D(C1) from P‐L2‐NNLS had more plausible and less spurious values than those from L2‐NNLS

The results of the data DKI fit with and without correction for C2 and C3 (P‐L2‐NNLS only), and their voxel‐wise differences, are shown in Figure 8 for an example slice. After correction, MD values were reduced by −15.50% and −8.56% in GM and WM, respectively, whereas FA estimated values were increased by 15.80% and 7.38%, respectively. Lower mean kurtosis values were observed after correction: −11.47% in GM and −5.47% in WM. The average values and standard deviations before and after correction for each time point are reported in Table 3.

Figure 8

Example axial images of mean diffusivity (MD), fractional anisotropy (FA) and mean kurtosis (K) maps obtained with a diffusion kurtosis imaging (DKI) fit of the data without any correction (NC) and after removal of macro‐compartments 2 and 3 (C2 + C3) with P‐L2‐NNLS (non‐negative least squares with data‐driven L2 regularization). The difference between corrected and uncorrected metrics is indicated as Δ. The fit of the corrected data resulted in lower MD and K values, especially in gray matter

Table 3

Average values ± standard deviation of mean diffusivity (MD), fractional anisotropy (FA) and mean kurtosis (K) in gray matter (GM) and white matter (WM) for the two time points (TP1 and TP2) in each subject

	S1		S2		S3
	TP1	TP2	TP1	TP2	TP1	TP2
GM
MD NC (mm2/s) × 103	0.96 ± 0.27	0.98 ± 0.26	0.91 ± 0.20	0.92 ± 0.19	0.94 ± 0.26	0.93 ± 0.24
MD C2 + C3 (mm2/s) × 103	0.81 ± 0.16	0.82 ± 0.12	0.78 ± 0.13	0.78 ± 0.13	0.82 ± 5.72	0.78 ± 0.13
FA NC	0.16 ± 0.08	0.15 ± 0.08	0.16 ± 0.08	0.15 ± 0.07	0.17 ± 0.08	0.16 ± 0.08
FA C2 + C3	0.19 ± 0.10	0.18 ± 0.10	0.19 ± 0.10	0.18 ± 0.09	0.19 ± 0.09	0.19 ± 0.10
K NC	0.59 ± 0.16	0.61 ± 0.17	0.61 ± 0.18	0.62 ± 0.18	0.54 ± 0.16	0.64 ± 0.16
K C2 + C3	0.51 ± 0.18	0.55 ± 0.18	0.53 ± 0.19	0.55 ± 0.19	0.48 ± 0.18	0.57 ± 0.19
WM
MD NC (mm2/s) × 103	0.83 ± 0.16	0.85 ± 0.14	0.82 ± 0.14	0.83 ± 0.14	0.81 ± 0.13	0.81 ± 0.13
MD C2 + C3 (mm2/s) × 103	0.76 ± 0.14	0.79 ± 0.12	0.73 ± 0.13	0.75 ± 0.13	0.75 ± 0.12	0.74 ± 0.12
FA NC	0.40 ± 0.15	0.36 ± 0.15	0.39 ± 0.16	0.36 ± 0.15	0.39 ± 0.16	0.37 ± 0.16
FA C2 + C3	0.43 ± 0.17	0.39 ± 0.16	0.43 ± 0.18	0.39 ± 0.17	0.41 ± 0.17	0.39 ± 0.17
K NC	0.87 ± 0.20	0.85 ± 0.20	0.88 ± 0.21	0.87 ± 0.21	0.82 ± 0.22	0.91 ± 0.21
K C2 + C3	0.81 ± 0.24	0.83 ± 0.21	0.81 ± 0.23	0.83 ± 0.24	0.77 ± 0.24	0.87 ± 0.24

NC, no correction.

DISCUSSION

Previous work40 introduced the spectral analysis as an investigative tool for the dMRI signal acquired at multiple b values. However, the proposed approach could only account for two predefined diffusion processes under the demand of high SNR acquisitions. The method presented in this work addresses these shortcomings, introducing outlier rejection in conjunction with data‐driven regularization in the diffusion spectrum computation, and proposing a method to automatically determine the number of macro‐compartments from the spectrum itself. Our findings show that: (1) the prior‐based regularization term introduced here improves the diffusion spectrum computation compared with the original approach; (2) it is possible to automatically and consistently divide the spectrum into multiple macro‐compartments without prior assumptions on their number; and (3) this approach can be used to reduce biases in DTI and DKI estimates.

The method proposed by Keil et al.40 showed that SNR levels of at least 100 are needed to obtain reliable estimates of the diffusion spectrum, which, although lower than previous reports of a minimum required SNR equal to 150,57 is still difficult to achieve in clinical acquisitions. In this work, we have shown that the addition of a modified regularization term included in P‐L2‐NNLS allows the estimation of the diffusion spectrum from clinically achievable SNR levels of 30. The usage of priors in conjunction with deconvolution methods has been proposed previously, for example following Bayesian approaches36 or enforcing sparsity constraints.38 Here, we rely on an L2 regularization term to account for the distance between the observed spectrum and its observed probability. Besides being computationally very efficient, our work shows that, among the three tested methods, it performs best with low‐SNR data, with remarkable improvements over L2‐NNLS.

For in vivo data, a known prior is not available, and therefore we have proposed the use of the average voxel‐wise spectrum computed with NNLS as diffusion spectrum prior. As it is a whole‐volume average, the diffusion spectrum prior is inherently not voxel‐wise correct. Nonetheless, the method showed superior performances to L2‐NNLS, in particular in the estimation of f(C2).

In simulations, the automated subdivision of the diffusion spectrum with the GMM fit was consistently capable of correctly identifying the number and location of the predefined diffusion components. For the in vivo data, 18 components were used on average to model the voxel‐wise signals, which were then consistently grouped into three automatically identified macro‐compartments in all datasets. This is in line with the recent dMRI literature investigating pseudo‐diffusion and FW contamination, which suggests the need for at least three distinct components to properly model the signal in the brain as well as in other tissues, e.g. the liver and the prostate.26, 27, 58 Furthermore, the diffusion coefficients of the three observed components were in agreement with previously reported literature values of brain tissue, FW and blood pseudo‐diffusion.3, 8, 25 Recent studies have suggested the existence of multiple components in the pseudo‐diffusion part itself (i.e. micro‐ and macrovascular network).3 Disentangling these components is possible with the proposed method if data with sufficient SNR are provided. However, this may have minimal impact when only the diffusion properties of the tissue are of interest.22

Compared with the more common tri‐exponential isotropic NLLS fit, the proposed method showed remarkably lower values in the estimation of the FW signal fraction (C2) and of the pseudo‐diffusion fraction (C3). For the tri‐exponential isotropic NLLS fit, the average values were 21% and 8.75%, respectively. In contrast, for P‐L2‐NNLS, the average values of C2 and C3 were 12% and 2.3%, respectively, which is similar to previously reported values.20, 22, 25 In addition, pseudo‐diffusion maps from the NLLS fit were noticeably noisier than those from P‐L2‐NNLS. These differences might be caused, for example, by convergence of the NLLS fitting procedure to local minima.29, 40, 59, 60, 61

The disentanglement of the FW compartment from hindered water pools is known to be a difficult (and ill‐posed) problem.8 Nonetheless, disentanglement of these two domains was achieved with P‐L2‐NNLS as shown by the simulations and in vivo data. In addition, the framework introduced here can be easily extended to accommodate multi‐contrast data, e.g. multiple diffusion weightings in combination with multiple TE acquisitions,62 which has been shown to further enhance the decoupling performances exploiting the different distributions of T 2 values between FW and tissues.

Considering the results obtained with in vivo data, the fractional maps were consistent and similar across time points and subjects, and showed good spatial correspondence to known physiology for both L2‐NNLS and P‐L2‐NNLS. Nonetheless, the latter resulted in more plausible D(C1) values and higher f(C2) estimates, in line with those observed in the simulations. Maps of D(C3) were the noisiest amongst the diffusion coefficient maps, in line with previous literature pointing towards the high variability of the pseudo‐diffusion coefficient.63 Consistent estimates were also obtained on a highly subsampled dataset (Supporting Information) corresponding to a 5‐minute‐long acquisition, which suggests the clinical applicability of the method.

The feasibility of our proposed method to correct for undesired partial volume biases was demonstrated using simulations. In both simulations, P‐L2‐NNLS outperformed L2‐NNLS, also at SNR = 100, where the standard L2 norm is likely to cause a minimum imposed blurring of the solution. With correct estimation and consideration of the various diffusion components, estimated values of MD, FA and K were much closer to their ground‐truth values. For example, the MD estimation error decreased from 26.41% in GM and 13.90% in WM to 8.03% and 6.50%, respectively. The correction for partial volume effects was also investigated with in vivo data, where it resulted in lower values and more homogeneous MD maps. The changes as a result of bias correction were in line with our simulations and previous findings of lower MD and K values after pseudo‐diffusion and FW correction.17, 21

Although our proposed method was successful in recovering the predefined values, it remains ill‐posed and ill‐conditioned. Therefore, it requires great care in data acquisition, fitting and interpretation. The basis functions used in this framework assumed Gaussian diffusion, which might be incorrect for the description of data over a wide range of b values.35, 64 Jbabdi et al.35 suggested the use of gamma distributions to limit data over‐fitting. However, although such a basis can be of interest in the disentanglement of fiber crossings, it is not trivial to derive quantitative metrics, such as diffusion coefficients, from such signal representation. In our work, we observed a rather large number of voxel‐wise components, 18, because of the effect of L2 regularization. Although other regularization techniques, e.g. the SADD approach,38 might help to reduce possible over‐fitting, the grouping into macro‐compartments offered a viable and effective solution. The diffusion coefficient of each macro‐compartment was determined as the weighted average of the corresponding spectral components, which, with low SNR, might not converge to the real underlying value. Although such a value represents mainly a summary statistic of the diffusion processes within each macro‐compartment, in our data it resulted in plausible diffusion maps (Figure 7).

The results of simulation I suggest that the method should be applied to data with an SNR of at least 30 (with respect to the b = 0 s/mm2 image). Further analysis, reported in Supporting Information, showed that P‐L2‐NNLS could be beneficial also at lower SNR, if a good prior distribution can be estimated. In this work, we used voxel‐wise NNLS fits of the whole brain to initialize the diffusion spectrum prior. This prior penalizes the less likely values for the diffusion spectrum, but still allows the fitting of local heterogeneity. However, it should be emphasized that this prior is only valid for SNR values above 20. Deriving the prior from other estimation techniques, such as L2‐NNLS with variable regularization settings, could potentially further improve the method performance. We observed small improvements in P‐L2‐NNLS over L2‐NNLS in the estimation of f(C3) in simulation I. This is again dependent on the NNLS prior, which was characterized by small and broadly distributed values in the pseudo‐diffusion range. Under this condition, the P‐L2‐NNLS regularization term reduced to a normal L2‐NNLS term, thus not further improving the estimation results. A further difficulty in the pseudo‐diffusion estimation might arise from the crusher gradients included in the diffusion sequence, which result in a minimum applied gradient around 2 s/mm2. Here, we do not explicitly account for their effect. However, we do not expect them to affect signal fraction estimations, but eventually the associated pseudo‐diffusion coefficient.

Furthermore, in the proposed method, we used the geometric mean, which does not allow us to capture the intrinsic anisotropy of tissues, FW and pseudo‐diffusion. Nonetheless, this is generally well accepted for FW,6 has been shown to be adequate on brain data at multiple diffusion values65 and is not expected to be detrimental on pseudo‐diffusion estimates.66, 67 In addition, as shown in this work, if part of the data is acquired with higher angular resolution, anisotropic models can still be applied after estimation of the isotropic components.

The representation of the diffusion signal using the diffusion spectrum allows us to easily assess the distribution of the main water compartments in tissue. As there are no tissue‐specific assumptions about the number of compartments, this method is widely applicable. Furthermore, the parameters obtained might be valuable to study pathology,18, 19, 68, 69 where changes in either diffusion coefficients and/or partial volume effects are expected. Further dedicated studies exploiting the potential of this technique in such a context, in conjunction with the tuning of parameters, such as the threshold overlap (see Supporting Information), are nevertheless required.

CONCLUSIONS

This work presents a robust method to disentangle multiple components from dMRI data acquired at multiple diffusion weightings, and to reduce the estimation biases caused by partial volume effects in common techniques, such as DTI and DKI. Unlike classic multi‐compartment modeling techniques, this approach does not require a knowledge of the number of components a priori, with the potential to be a tissue‐independent framework.

Figure S1 – Performances of P‐L2‐NNLS in disentangling signal mixtures at different SNR levels (simulations I) with different regularization parameter settings (γ). Columns 1–3 correspond to each of the three simulated diffusion components. Error bars represent the 25th, 50th and 75th percentile of the estimated values of diffusion coefficients (first row) and signal fractions (second row) from 1000 realizations, while black dashed line represent the predefined value.

Figure S2 – Performances of P‐L2‐NNLS compared to L2‐NNLS in disentangling signal mixtures at different SNR levels (simulations I) when employing the simulated diffusion spectrum as prior. Columns 1–3 correspond to each of the three simulated diffusion components. Error bars represent the 25th, 50th and 75th percentile of the estimated values of diffusion coefficients (first row) and signal fractions (second row) from 1000 realizations, while black dashed line represent the predefined value.

Figure S3 – Voxel‐wise estimates of signal fractions (f, columns 1–3) and diffusion coefficients (D, columns 4–6) provided by P‐L2‐NNLS fit of the first time‐point of subject 1 with different regularization parameter settings (γ).

Figure S4 – voxel‐wise estimates of signal fractions (f) provided by P‐L2‐NNLS fit of the first time‐point of subject 1. Each row corresponds to a different threshold value used in the automatic GMM fit. A higher number of classes were detected with higher threshold values, up to a maximum of 5.

Figure S5 – voxel‐wise estimates of signal fractions (f, columns 1–3) and diffusion coefficients (D, columns 4–6) provided by P‐L2‐NNLS fit of the first time‐point of subject 1 considering the full acquisition presented in the main text (first row), and the same data subsampled by a factor 3 (second row).

Click here for additional data file.