Generalized Richardson-Lucy (GRL) for analyzing multi-shell diffusion MRI data

Spherical deconvolution is a widely used approach to quantify fiber orientation distribution from diffusion MRI data. The damped Richardson-Lucy (dRL) is developed to perform robust spherical deconvolution on single shell diffusion MRI data. While the dRL algorithm could in theory be directly applied to multi-shell data, it is not optimised to model the signal from multiple tissue types. In this work, we introduce a new framework based on dRL - dubbed Generalized Richardson Lucy (GRL) - that uses multi-shell data in combination with user-chosen tissue models to disentangle partial volume effects and increase the accuracy in FOD estimation. The optimal weighting of multi-shell data in the fit and the robustness to noise and partial volume effects of GRL was studied with synthetic data. Subsequently, we investigated the performances of GRL in comparison to dRL on a high-resolution diffusion MRI dataset from the Human Connectome Project and on an MRI dataset acquired at 3T on a clinical scanner. The feasibility of including intra-voxel incoherent motion (IVIM) effects in the modelling was studied on a third dataset. Results of simulations show that GRL can robustly disentangle different tissue types at SNR above 20 and improves the angular accuracy of the FOD estimation. On real data, GRL provides signal fraction maps that are physiologically plausible and consistent between datasets. When considering IVIM effects, high blood pseudo-diffusion fraction is observed in the medial temporal lobe and in the sagittal sinus. In comparison to dRL, GRL provides sharper FODs and less spurious peaks in presence of partial volume effects and results in a better tract termination at the grey/white matter interface or at the outer cortical surface. In conclusion, GRL offers a new modular and flexible framework to perform spherical deconvolution of multi-shell data.


Introduction
Diffusion MRI (dMRI) is an established technique to study the microstructure of brain's white matter (WM). With fiber tractography (Mori et al., 1999), it is possible to reconstruct the main WM bundles given knowledge on the direction of the underlying diffusion process. When the process has a single diffusion orientation, the diffusion tensor imaging (DTI) model (Basser et al., 1994;Pierpaoli and Basser, 1996) can be employed to efficiently derive the main diffusion direction from six or more diffusion weighted measurements. However, the in-vivo brain is characterized by the presence of multiple fiber crossings (Behrens et al., 2007;Jeurissen et al., 2013, Dell'Acqua et al. 2013, which cannot be disentangled with DTI. To reconstruct the directions of multiple crossing fibers, more general approaches such as spherical deconvolution (Anderson, 2005;Dell'Acqua et al., 2010;Dell'Acqua and Tournier, 2019;Tournier et al., 2004) or q-space inversion methods (Tuch et al., 2002;Wedeen et al., 2005) are commonly employed.
Spherical deconvolution approaches estimate the fiber orientation distribution (FOD) of WM using data acquired diffusion gradients along multiple spatial orientations distributed over a sphere. When the acquisition is performed at a single diffusion weighting, the data is commonly referred as singleshell data, whereas multi-shell data implies that the sampling has been performed at two or more diffusion weightings. Once the FOD has been derived, it can be used directly to perform fiber tractography as well as for generalized voxel-wise fiber specific statistics (Dell'Acqua et al., 2013a;Raffelt et al., 2010). In general, a suggested data acquisition scheme for single shell spherical deconvolution approaches consists of at least 60 gradient directions with a diffusion weighting of around b-value = 3000 s/mm 2 . Data at such strong diffusion weighting are generally characterized by degraded SNR (Descoteaux et al., 2011;Jones et al., 2013), but allows to resolve crossing fibers with higher angular resolution as compared to lower diffusion weightings. Several spherical deconvolution approaches (Dell'Acqua et al., 2007;Tournier et al., 2007) have been proposed to estimate the FOD from single shell diffusion weighted signals. Among others, the Richardson-Lucy spherical deconvolution method (Dell'Acqua et al., 2007) is a well-established method to estimate the FOD from single-shell dMRI data. In a later formulation (Dell'Acqua et al., 2010), the approach was extended as damped Richardson-Lucy (dRL) to increase the robustness of the FOD estimation to partial volume effects with isotropic tissues and to prevent overfitting to noise.
Recent improvements in hardware have strongly increased the acquisition efficiency, making it feasible to collect dMRI data at multiple diffusion weighting (multi-shell) in a single acquisition within clinical achievable time. Multi-shell dMRI has potential to further improve the performance of the FOD estimation, as previously shown by Jeurissen et al.  in the context of constrained spherical deconvolution (Tournier et al., 2007). In that study, the authors modelled the three main tissue types of the brain (WM, GM and CSF) to improve the accuracy of WM FODs in presence of partial volume effects in comparison to single shell CSD (Tournier et al., 2007), while also computing the signal fractions of each compartment. Different tissue types, as for instance WM, cerebrospinal-fluid (CSF) and grey matter (GM) are characterized by distinct attenuation profiles over increasing diffusion weighting, which allows to model their contributions, a problem otherwise ill-posed with single shell data. Once signal fraction maps of different tissue classes have been derived, they could potentially be used to guide the fiber tractography procedure, in analogy to what is done with T1 derived segmentations . This unexplored application would knot only avoid the extra-acquisition time of structural imaging, but also prevent residual misalignments between T1w data and dMRI. Further, the reconstruction of WM FODs could potentially benefit of multi-shell data by leveraging the high angular resolution of strong diffusion weighting in conjunction with the higher SNR of lower diffusion weighting signals.
In this study, we introduce the Generalized Richardson Lucy (GRL) framework, a modular spherical deconvolution that aims at improving the reconstruction accuracy of WM FODs by accounting for an arbitrary -used defined -number of signal models representing different tissue types. We show through simulations that GRL can disentangle up to 4 different tissue types, and accommodate multiple choices to represent the WM signal, as the DTI model with a single set of eigenvalues, the DTI model with shell-dependent eigenvalues, the diffusion kurtosis imaging (DKI) model, and the neurite orientation dispersion and density imaging (NODDI) model. GRL allows to perform fiber tractography with higher accuracy than dRL without the need to define FOD thresholds for the fiber tractography procedure but rather uses the derived signal fraction maps to terminate the tractography procedures, delivering high quality termination at the GM/WM interface or at the outer GM surface.

Theory
In the following sections we introduce the essential notions of the original damped Richardson Lucy implementation (section 2.1), and show how GRL accommodates multi-shell diffusion MRI data to perform multi-tissue decomposition (section 2.2). Further, in sections 2.3 and 2.4, we show how different WM models can be used to generate the deconvolution kernel, as well as how to take into account intra-voxel incoherent motion (IVIM) effects to quantify signal fractions of blood pseudodiffusion.

The original dRL spherical deconvolution approach
The Richardson-Lucy approach (Dell'Acqua et al., 2007) derives the FODs directly in the signal domain, which avoids artefacts that can arise from alternative representations as spherical harmonics, including truncation and Gibbs ringing. The Richardson-Lucy framework requires to define a deconvolution kernel H, or H-matrix, that projects the diffusion MRI signals measured along each gradient direction onto the unit sphere. Given the signal s acquired at a single diffusion weighting (bvalue) applied along multiple gradient directions, the FOD of each voxel can be estimated with an iterative expectation maximization approach, as reported in Eq. [1].
In the equation above, FOD (k) is the FOD estimated at the k th iteration along the i-th direction sampled on the unit sphere, H is the m by n deconvolution kernel that maps the m (e.g., 60) measured signals onto n (e.g., 300) arbitrarily discretized directions uniformly sampled on the unit sphere, and H T is the transpose of H. Eq.
[1] is semi-convergent, thus its solution practically requires a stopping criteria.
In Eq.
[2], the regularization term u (k) is updated at each iteration k based on the amplitude of the current FOD estimation and the standard deviation of the diffusion weighted signals std(s), as shown in Eq.
In this work, the constant was set to =8 (Dell'Acqua et al., 2010), whereas was set to two times the maximum amplitude of the FOD describing an isotropic diffusion process with diffusion coefficient 0.7x10 -3 mm 2 /s. Such value of is chosen to minimize the contribution of GM to the FOD of WM.

The GRL approach for multi-shell diffusion MRI data
In section 2.2.1 we describe the GRL framework to describe the diffusion signal from multiple tissue types by using multi-shell dMRI. [5] Here p is the total number of considered diffusion weighting (i.e., shells). In practice, the H-matrix introduced with GRL consists of the stacking of multiple H-matrices corresponding to each of the acquired shells. To account for mixtures of multiple tissue types (partial volume effect), the problem can be formulated by adding further columns to the H-matrix to describe the tissues of interest. In name of simplicity, here we will consider a formulation similar to a previous study introducing multishell CSD . Therefore, we will consider three tissue types when creating the Hmatrix, but other choices, e.g. considering intra-vs extra-cellular signals (Alexander et al., 2002), could be also taken into account. [6] In the following sections, we refer to GRL as the modeling of three tissue types representing WM, GM and CSF, where not otherwise specified. The columns of H-matrix corresponding to isotropic GM and CSF contributions can be generated with the apparent diffusion coefficient (ADC) model and using typical values of the two compartments, e.g. 0.7 x 10 -3 mm 2 /s and 3.0 x 10 -3 mm 2 /s. Regarding the modeling of WM, the following sections discuss some possible model choices to generate the H-matrix, such as the diffusion tensor model with a single fixed tensor (GRL-FT) (HDTIs), DTI with a shell-dependent variable tensor (GRL-VT), the diffusion kurtosis imaging model (GRL-DKI), and a simplified white matter model from the NODDI framework (GRL-NODDI), as explained in the following sections.

DTI-based H-matrix with fixed tensor (GRL-FT)
The original formulation of the single-shell dRL method (Dell'Acqua et al., 2010) used the DTI model and a predefined set of eigenvalues to initialize the H-matrix. Similarly, we generated the columns of the H-matrix of GRL-FT using eigenvalues equal to [1.7, 0.2, 0.2] x 10 -3 mm 2 /s. This choice of the eigenvalues implicitly assumes non-Gaussian effects to be negligible.

DTI-based H-matrix with variable eigenvalues per shell (GRL-VT)
When fitting the DTI model to individual shells, different tensor values are determined as a result of the heterogeneous organization of in-vivo tissues and of non-Gaussian diffusion effects. For this reason, estimates obtained with data at strong diffusion weighting (e.g., b = 3000 s/mm 2 ) generally differ from those obtained with data at lower diffusion weighting (e.g., b = 1000 s/mm 2 ). The columns of the Hmatrix of GRL-VT were generated using the DTI model but different tensor values for each shell independently, which implicitly allows for non-Gaussian effects in the diffusion process. The values of the tensors were determined from the data itself by fitting each shell and the non-weighted images with the DTI model independently and averaging the values in a mask obtained by thresholding the fractional anisotropy mask above 0.7.

DKI-based H-matrix (GRL-DKI)
The diffusion signal is known to be non-Gaussian in WM at strong diffusion weighting (Jensen et al., 2005). Previous studies (Jensen et al., 2005;Jensen and Helpern, 2010) introduced the DKI model, an extension of DTI to account for the restriction of water molecules observed at strong diffusion weighting, as described in Eq. [7].
In the above equation, SDWI and S0 are the diffusion weighted and the non-diffusion weighted signals, respectively, D is the diffusion coefficient, b the diffusion weighting, and K the value of mean kurtosis. The columns of the H-matrix of GRL-DKI were generated with Eq.
[7] using tensor values and mean kurtosis K values estimated from the data. All the data was fit with a DTI model extended to account for isotropic kurtosis (De Luca et al., 2018b) and using a weighted linear least squares estimator (Veraart et al., 2013). The values of the tensor and of K were averaged within a mask defined by fractional anisotropy values above 0.7.

Microstructure-based H-matrix (GRL-NODDI)
The columns of the H-matrix describing WM can be modelled not only with statistical descriptions of the diffusion signal, such as DTI and DKI, but also with microstructural models. As a proof of concept, here we consider the neurite model defined in the neurite orientation dispersion and density imaging (NODDI) (Zhang et al., 2012). We generated the columns of the H-matrix of GRL-NODDI by setting the intracellular volume fraction to 1, the intrinsic free diffusivity of the medium to 1.7 x 10 -3 mm 2 /s, the concentration parameter of the Watson distribution to 3.5, and the isotropic volume fraction to 0. These values are meant to represent the typical an average NODDI parameter set observed in WM, with the exception of CSF for which we account explicitly in GRL.

A two-step iterative fit approach
The solution of the deconvolution problem with the above-mentioned H-matrices would result in unstable solutions if using directly the iterative dRL scheme. For this reason, we decoupled the problem into two steps, namely the estimation of the tissue fractions and the FOD calculation. The new GRL approach can be written as an iterative process alternating the dRL algorithm and a non-negative leastsquares estimation, described in Eq.
In the above equation, fWM, fGM and fCSF represent the signal fractions associated to each compartment, s is the measured signal (normalized by the non-weighted image), H is the H-matrix describing the fiber response, and Y is an m by 2 matrix containing the signal representation of the isotropic contributions (GM, CSF). The * operator symbolizes a spherical convolution operation to compute the signal given the estimated FOD. In eq. [9] we use an intrinsic regularization function for the FOD function to further penalize spurious peaks, which is the median thresholding of the FOD, denoted as FOĎ in Eq. [9]. The values of fractions fGM and fCSF were set to zero in the first iteration.
The modularity of the GRL framework allows to easily take into account more isotropic compartments by simply expanding the deconvolution matrix along the columns dimension. As a proof of concept, here we show that in addition to WM, GM and CSF, it is also possible to account for blood pseudodiffusion using the intra-voxel incoherent motion (IVIM) model (Le Bihan et al., 1988Bihan et al., , 1986. The diffusion coefficient of the IVIM component was set to an arbitrary value of 50.0 x 10 -3 mm 2 /s to generate an additional column of the isotropic deconvolution matrix Y. The adjusted fit procedure to account for WM, GM, CSF and IVIM is reported in Eq.
To numerically solve Eq.
[10] and [11], the number of acquired shells must be larger than the number of modelled compartments, and the columns of H not rank-deficient. In the case of Eq.
[10] and [11], this means that at least 4 unique diffusion weightings are needed, of which at least one where the IVIM contribution is not null.

Methods
Single-shell dRL and GRL were fit to both simulated and in-vivo data. The H-matrix used for dRL was generated in agreement with previous work (Dell'Acqua et al. 2010), e.g. with a diffusion tensor described by the eigenvalues [1.7, 0.2, 0.2] x 10 -3 mm 2 /s.

Weighting of multi-shell data
Data acquired at strong diffusion weighting has high angular resolution but generally poor SNR, whereas data at low diffusion weighting has high SNR but poor high angular resolution. In the GRL formulation, the contribution of each shell to the fit can be weighted ) by multiplying both the signal and the rows of the H-matrix corresponding to a specific shell by a scalar . To investigate the effect of low diffusion weighting data on the angular resolution, we simulated a single-voxel crossing fiber configuration with 47 degrees crossing angle and 1000 realizations of Rician noise at SNR 50. We performed the fit with GRL-FT by varying the weight of all shells but the outer between 0.1 and 1.

Simulations
The robustness to noise of GRL was tested with numerical simulations. Single fiber populations were generated with a three compartments model accounting for WM, GM and CSF. The DTI model was employed for all three compartments, setting the eigenvalues of the WM tensor equal to [1.7, 0.2, 0.2] x 10 -3 mm 2 /s, while GM and CSF were considered isotropic with mean diffusivity values of 0.7 x 10 -3 mm 2 /s and 3.0 x 10 -3 mm 2 /s, respectively. Three different volume fractions (fWM, fGM, fCSF) mixtures were generated under the constraint of f WM + f GM + f CSF = 1, including: 1. fGM = 0, fWM assuming values between 0 and 1 with step 0.1; 2. fCSF = 0, fWM assuming values between 0 and 1 with step 0.1; 3. fWM assuming values between 0 and 1 with step 0.1, fGM = fCSF = 0.5 x (1 -fWM).
Realizations of Rician noise were added to the signals to simulate the SNR levels at b = 0 s/mm 2 equal to [20,30,50, 1000] with 1000 repetitions each. As a measure of performance, the derived signal fractions and the FOD peak orientations were compared to the simulated values.

Data acquisition
Three in-vivo datasets with different acquisition protocols were used in this work to demonstrate the performance of GRL.

Data Analysis
Data1 and Data2 were fit with GRL using the above-described H-matrices (section 2.2.2) to evaluate the effect of voxel resolution on the fractional maps. The WM FOD obtained on the two datasets were compared to those derived with a dRL fit of the largest shell of each dataset. The number of fibers (NuFO) (Dell'Acqua et al., 2013b) of each FOD derived with dRL and mDRL was derived and compared between methods. In particular, we counted only peaks with magnitude above 20% the average amplitude of the first peak in pure WM, which was defined as voxels with a value of fractional anisotropy above 0.7. Finally deterministic fiber tractography was performed in ExploreDTI (Leemans et al., 2009) with the derived FODs for both dRL and GRL. The tracking procedure was performed on dRL using angle-threshold 45 degrees, step size equal to half of the voxel-size. Tractography of the FODs derived with dRL was performed first with an FOD threshold equal to 0.1, which corresponds to the FOD amplitude of a GM voxel and is meant to stop the tractography at the WM/GM interface. Subsequently, the tractography was repeated with an FOD threshold equal to 0.001 to allow the continuation of the tracts into cortical GM. For the tracking of FODs derived with GRL, equal parameters but no FOD threshold were used. The signal fraction maps associated to WM and GM, which are derived with GRL but not with dRL, were used as termination criteria to stop tracking at the GM/WM interface, and at the outer GM surface, respectively.

Results
The effect of including shells at multiple diffusion weightings in the FOD estimation is shown in Fig. 1.
Increasing the weight of the shells at b = 1000, 2000 s/mm 2 in the fit reduces the angular resolution of data at b = 3000 s/mm 2 , inducing angular deviations on the first peak up to 7.3. Including data at multiple diffusion weightings is, however, essential to account for partial volume effects, which can be optimally estimated at = 0.2 and results in a partial volume estimation error of 1012% with a corresponding angular deviation of 3.8. This weighting value was therefore used throughout this manuscript.
The signal fraction maps estimated with GRL-FT in Simulation I-II-III and their standard deviations are shown in Fig. 2. At SNR 20, the precision of the estimation of small CSF-like partial volume effects drops (fWM > 0.9, fCSF < 0.1). For SNR levels above 20, GRL can effectively disentangle partial volume effects with high accuracy. For linearly decreasing partial volume effects, linearly increasing values of fWM are estimated, although a small but consistent overestimation of fWM -underestimation of fGM -is revealed.
The average angular error of the main FOD orientation estimated with dRL and GRL-FT in Simulations I-II-III are shown in Fig. 3. While dRL can effectively suppress the effect of partial volume effects on the FOD for fWM values above 0.5, angle errors up to 50 degrees are observed for lower fWM values, in particular at SNR 30. Conversely, FODs estimated with GRL have less angular error in presence of strong partial volume effects for SNR above 20, with a maximum angular error of 9 degrees for fWM = 0.2. At lower SNR levels (SNR 20), the application of GRL-FT is still beneficial in comparison to dRL, but the angular error of the FOD remains large enough to discourage its use. For high SNR levels (SNR  50), dRL results in minimal angular error for fWM  0.2. GRL can be used to take into account more than three tissue models, including the IVIM effect observed at low diffusion weightings. This is shown in Fig. 4 for an axial slice of the MASSIVE dataset.
In addition to signal fraction maps associated to WM, GM and CSF remarkably similar to those observed with Data1 and Data2 (Fig. 3), GRL estimated an fIVIM map showing high perfusion values in the middle temporal lobe, and in locations compatible with the middle cerebral artery. The estimated values of fIVIM in tissues is 116%.
The FODs estimated with dRL and GRL with Data1 and Data2 in an example coronal slice are shown in Fig. 6 and Fig. 7, respectively. Generally, GRL-FT and GRL-DKI perform similarly to dRL within WM, providing excellent separation of up to three crossing fiber configurations (red box). In proximity of the ventricles (yellow box), GRL-DKI visibly reduces the estimation uncertainty of the FODs, and further suppresses FODs clearly belonging to the CSF area. GRL-NODDI also reduces the effect of partial volume effects on the WM FOD as compared to dRL but is less effective at separating crossing fiber configurations. The use of GRL-VT results in numerically unstable FOD estimation, possibly in combination with the chosen weighting of the multi-shell data. Fig. 8 and Fig. 9 show an example axial slice of the FODs estimated with dRL and GRL with Data1 and Data2 in a region where the WM enters the cortical folding (red box). Within WM, the FODs estimated with dRL, GRL-FT and GRL-DKI show remarkable similarity. On Data1, both GRL-FT and GRL-DKI estimate sharper FODs with a reduced number of peaks in comparison to dRL, particularly in proximity of GM/CSF (yellow box). As previously observed, the FODs estimated with GRL-NODDI contain a lower number of crossing fibers also in this region. On Data 2, which has lower spatial resolution than Data1 and includes a diffusion weighting below b = 1000 s/mm 2 , the FODs estimated with GRL in proximity of the cortex are strongly attenuated as compared to those estimated with dRL.
To investigate to which extent GRL resulted in a different number of spatial orientations as compared to dRL, we derived the FOD peaks of both methods throughout the brain which amplitude was above 20% that of a typical WM peak. Results, which are reported in Fig. 10, show that GRL-FT performs very similarly to dRL with the exception of a slight reduced number of voxels with 3 or more crossing fiber configurations, which visually are located in proximity of the cortex. GRL-DKI results in a reduced number of three fiber crossings throughout the WM, but this result might be driven by the peak thresholding as such fiber configurations can be appreciated in the FODs shown in Fig. 6 and Fig. 7. In GM, GRL-DKI and GRL-NODDI essentially reduced the number of observed peaks to 1.
The whole brain fiber tractography of Data1 and Data2 are shown in Fig. 11 and Fig. 12. The figure shows the fiber tractography results obtained with dRL for a threshold equal to the average FOD amplitude in GM (t=0.1) -which should therefore terminate the tracking at the WM/GM interfaceand without using constraints on the FOD amplitude (t=0.001), allowing the tracking to continue to the cortical surface (and beyond). When tracking with GRL, we did not use any FOD-based termination criteria but rather stopped the tracking when fWM=0 and fGM=0, respectively. On both datasets, we observe that GRL provides tractography results that are practically identical to those provided by dRL within WM, as exemplified by the red tracts shown in Fig. 10 and Fig. 11. On Data2, however, the WM tractography provided by dRL showed some spurious -likely artefactual -fiber reconstructions in the frontal lobe, which were not observed with GRL. When lowering the FOD threshold to allow tracts, for instance, to enter in cortical grey matter, tractography results provided by dRL (t=0.001) are swamped by erroneous reconstructions. Conversely, the GM-constrained fiber tractography with GRL show much cleaner results with a physiologically plausible delineation of the cortical ribbon, and a generally plausible continuation of WM bundles into the cortex (yellow box).

Discussion
In this study, we generalize the Richardson Lucy framework (GRL) to leverage the information content of multi-shell diffusion MRI data and model multiple tissue types. GRL can be used in combination with dMRI models of choice -including DTI, DKI and NODDI -to disentangle partial volume effects of WM with GM, CSF and other intra-voxel incoherent motion effects (IVIM) to attenuate their contamination on WM FODs.
Previously, Dell'Acqua et al. (Dell'Acqua et al., 2010) introduced the dRL framework to reduce the number of spurious peaks in the FOD due to partial volume effects, showing the importance of accounting for isotropic contaminations of cerebrospinal fluid in spherical deconvolution. More recently, Jeurissen et al.  have shown that by using multi-shell diffusion MRI data and modeling explicitly the contributions of CSF and GM in the spherical deconvolution problem, it is possible to effectively improve the FOD estimation accuracy as compared to deconvolution methods using only single-shell diffusion MRI data (Roine et al., 2014). In analogy with the work of Jeurissen et al. in the context of constrained spherical deconvolution, here we revisit the dRL approach to improve the accuracy of the FOD estimation by taking into account multiple tissue types and the information content of multi-shell diffusion MRI data. Differently from other approaches, and in particular from the work of Jeurissen et al., GRL performs the spherical deconvolution operation in the original signal space rather than in the spherical harmonics basis (Cheng et al., 2014;Christiaens et al., 2017;Tournier et al., 2004b), and uses a model-based response function rather than estimating it from the data (Tax et al., 2014;Tournier et al., 2007). While a more extensive comparison between the two methods can be found elsewhere (Dell'Acqua and Tournier, 2018;Parker et al., 2013), the choice of using a modelbased approach in GRL potentially allows to apply the method to tissues which response function is difficult to objectively isolate, such as the kidneys (De Luca et al., 2018a), and to data acquired with non-spherically distributed diffusion gradients without performing any q-space interpolation (Wedeen et al., 2005).
The SNR of dMRI data typically degrades with the strength of the applied diffusion gradient, which conversely enhances the angular resolution of the data . Performing multi-shell deconvolution might allow to combine the benefits of both, but this has not been previously well investigated. In the MSCSD approach  for instance, only the larger shell is used to perform the spherical deconvolution operation, whereas only the average value of the lower diffusion weightings is used to identify partial volume effects. In Fig. 1, we investigated a new way to include data at lower diffusion-weighting in the spherical deconvolution without worsening the angular resolution. We established that weighting the intermediate shells at about 20% represents an optimal choice in our study, trading off angular resolution in exchange of the ability to model tissue types over different diffusion weightings, a required step to allow for signal formulations as NODDI (Zhang et al., 2012). It should be noted, however, that this choice of the weighting represents an optimal scenario only for Data1, and that the weighting optimization should be repeated in sight of the effectively considered gradient scheme.
With simulations, we investigated the performances of GRL at disentangling partial volume effects of WM-like, GM-like and CSF-like at different SNR levels. The results shown in Fig. 1 indicate that GRL can estimate three class mixtures with high accuracy when the SNR is equal or above 30. A small but consistent bias is observed that leads to overestimation of fWM and underestimation of fGM. Nevertheless, the linear relation between simulated and estimated fWM values is not affected by such effect, which we dim negligible to the final application of GRL. At SNR 20, the uncertainty on the estimation of small fCSF becomes relatively large, suggesting that partial volume effects below 5% might be undetected if a sufficient SNR level is not met. This is unsurprising, as the accurate and reliable estimation of mixture of exponentials with little signal fraction is a well-known problem in diffusion MRI (De Luca et al., 2018b;Pasternak et al., 2009;Pierpaoli and Jones, 2004). Partial volume effects are known to introduce spurious peaks in the FODs and are therefore expected to introduce angular errors in the main FOD peaks. Fig. 2 shows that GRL can substantially increase the robustness of the FOD estimation to partial volume effects in comparison to dRL. For very low WM-like signal fractions, e.g. fWM  0.2, the angular deviation of the FOD is reduced of a factor about 2 as compared to dRL, which might allow future applications outside pure white matter. At SNR 20, the application of GRL still improves the angular accuracy of the FOD in comparison to GRL, but the angular error might remain too large in correspondence of small WM-like signal fractions.
Results on in-vivo MRI data suggest that GRL can be applied to datasets characterized by different diffusion weightings and imaging resolution, computing high quality FODs on data from the HCP project (Data1), a dataset acquired with a 3T clinical scanner at a lower spatial resolution (Data2), and the MASSIVE dataset (Data3). The signal fractions estimated with GRL on Data1 and Data2, shown in Fig. 4, are anatomically plausible and show high correspondence to the tissue types ideally included in the modelling, and similarly to what was previously shown with multi-shell CSD . The choice of the H-matrix has a critical impact on the estimated signal fractions. In this work we have shown the applicability of GRL in combination with common models, such as the diffusion tensor imaging (Basser et al., 1994), the diffusion kurtosis imaging (Jensen et al., 2005) and the neurite orientation dispersion and density imaging models (Zhang et al., 2012). The purpose of considering multiple models in our analysis is to enforce the generalizability of the GRL framework, showing its application with both statistical representation of the diffusion MRI signal as well as with microstructural models. With recent work showing the inhomogeneity of the WM microstructure (Schilling et al., 2019), GRL might serve as a supporting framework to perform spherical deconvolution while accounting for microstructure specific models, eventually informed by acquisitions decoding additional details on the tissue microstructure, such as neurite density (Nilsson et al., 2018;Westin et al., 2014).
In our results, the CSF signal fraction maps estimated using the NODDI based H-matrix (GRL-NODDI) have higher values throughout the brain as compared to other H-matrix choices, which suggests that the parameters used to initialize the NODDI models might need further optimization. The number of isotropic components to be considered in brain tissue modelling remains debated. In the context of spherical deconvolution, it has become common to consider WM, GM and CSF separately to improve the performance of FOD estimation. Outside the context of spherical deconvolution, however, the quantification of signal fraction maps associated to different water compartments is of high interest. The CSF fraction map quantified by GRL (and other previous approaches) can be compared to the free water mapping technique (Pasternak et al., 2009), which has been shown to provide additional insights as compared to DTI in neurodegenerative diseases as Parkinson's disease. With diffusion MRI data acquired at low diffusion weighting, it is also possible to quantify pseudo-diffusion signal fractions, which provide valuable information in diseases as brain tumors (Gong et al., 2018). Fig. 5 shows that, if data at sufficient diffusion weightings is available, GRL can suppress partial volume effects on the FOD while quantifying signal fraction maps of WM, GM, CSF and IVIM (Fig. 5). These signal fractions, which are in line with previous reports (Hare et al., 2017), pose a first step to bridge the gap between quantitative multi-compartment models and multi-shell spherical deconvolution. The values of fIVIM estimated in this work, 116% are relatively high as compared to previous literature on pseudodiffusion, and probably reflect an under calibration of the pseudo-diffusion coefficient used in the modelling. Fig. 6 and Fig. 7 show that GRL estimates high quality FODs in WM, suggesting that the state-of-theart angular performance of dRL is preserved in this new formulation. When partial volume effects become prominent, however, as in proximity of the ventricles or in cortical regions (Fig. 8, Fig. 9), GRL improves the FOD estimation by reducing the number of spurious peaks as compared to dRL, which results in sharper -more plausible -fiber orientations. As expected, the choice of the model (H-matrix) used to represent WM has a major impact on the estimated FOD. We observe that GRL-FT and GRL-DKI preserve the ability of dRL to efficiently resolve WM crossing fibers, although the results shown in Fig. 10 suggest that GRL-DKI might tradeoff the ability to disentangle a third peak for robustness. However, it should be noted that in Fig. 10 we used an arbitrary threshold on the WM FOD equal to 20% the average amplitude of the first FOD peak in WM. Lowering this threshold might reveal a large number of 3+ crossing fibers in GRL-DKI. In contrast to the above-mentioned, the use of both GRL-VT and GRL-NODDI worsened the spherical deconvolution performances in comparison to dRL at resolving crossing fibers. Regarding GRL-VT, a large number of spurious peaks was also observed, which hints toward instability of this signal representation in the context of spherical deconvolution. In first instance, the diffusion tensor model itself might not be suited to fit data at strong diffusion weighting (e.g. b = 3000 s/mm 2 ) due to the effect of membrane restrictions in specific spatial directions (Jensen and Helpern, 2010). Furthermore, the fit might become unstable as a result of the multi-shell weighting strategy explained in Fig. 1, which is however needed to preserve the angular resolution of the outer diffusion shell. Regarding GRL-NODDI, the NODDI model takes into account the effect of axonal dispersion (Zhang et al., 2012), which might the angular resolution in the deconvolution context. Optimizing the initialization parameters of GRL-NODDI is therefore expected to improve the FOD estimation. Besides the angular resolution, identifying the correct number of peaks is critical when computing measures such as the apparent fiber density (Raffelt et al., 2012) or the Number of Fiber Orientations (Dell'Acqua et al., 2013b), which are becoming used in clinical studies (Genc et al., 2018). Among the four tested H-matrices, GRL-DKI seems to provide the best choice (Fig. 10), providing good performances in the solution of the fiber crossing problem while improving the fit quality by accounting for the non-Gaussian behavior of the signal at stronger diffusion weighting (Jensen et al., 2005).
Fiber tractography is one of the most common applications of spherical deconvolution methods. Once the FOD has been estimated, fiber tractography generally requires defining some termination criteria, including a threshold on the amplitude of the FODs. If the study interest is to track purely the white matter, defining an FOD threshold amplitude equal to that of GM (t=0.1 in this study) can provide an effective solution, as shown for the tractography results of dRL FODs in Fig. 11 and Fig. 12. Nevertheless, it should be noted that such approach might also lead to the removal of genuine WM peaks, and is prone to be affected by inhomogeneity in the image, due to either B1 or receive field performance. The problem worsens when attempting to study tracts entering in GM, as commonly done in connectivity studies. To allow tracts entering the cortex, the FOD threshold should be lowered, but no explicit criteria can be easily established onto which threshold value to use. Removing the FOD threshold termination criteria is generally also not an option, as it would generate a large number of false positives, as shown for dRL in Fig. 11 and Fig. 12 (condition t=0.001). The GRL framework offers an elegant solution, as it allows to disregard FOD thresholding, and use the information of the signal fraction maps to constrain the tracking procedure. As shown in Fig. 11 and Fig. 12, following this approach tractography can be stopped (approximately) at the GM/WM interface by using the WM signal fraction map as termination criteria, or allowed to enter the cortex by using the GM signal fraction map. Reducing the number of spurious tracts can potentially contribute to reduce the false positive problem, which has been shown being critical in current fiber tractography (Maier-Hein et al., 2017). Furthermore, by improving the accuracy of WM fiber tracking in proximity of the cortex, GRL can potentially provide better results than dRL when computing the structural connectome (Tau and Peterson, 2009). Indeed, in connectomics, graph representations of the brain are constructed by assigning each tract to the closest cortical end-point by dilation of the cortical region, and inaccurate tracking in the GM/WM interface is likely to affect the process (Wei et al., 2017). Beyond the connectome, fiber tractography is commonly applied in other fields such as neurology (Reijmer et al., 2017(Reijmer et al., , 2013 or surgical planning (Essayed et al., 2017). Brain tumors, among other conditions, can cause presence of edema in the surrounding healthy tissues, which can worsen the performances of fiber tracking unless properly considered (Gong et al., 2018). Similar considerations apply also to conditions such as Parkinson (Ofori et al., 2017) or psychosis (Lyall et al., 2019), where increases of free-water fractions in WM have been reported. We here speculate that another isotropic contamination effect, e.g. the intra-voxel incoherent motion of blood in the micro-vascular network (Le Bihan et al., 1988), might act as an additional confounder in spherical deconvolution along with free water. While the extent of such effect -if any -has not yet been proved, in this work we have shown that it is possible to take it into consideration in the GRL framework (Fig. 5).
Some limitations of our work should be acknowledged. All the test dataset used in this study are brain data. While the brain is the most common target application of spherical deconvolution techniques, the acquisition of multi-shell diffusion MRI protocols outside the brain is feasible, e.g. to resolve crossing fibers in the human tongue muscles (Voskuilen et al., 2019) or in the kidneys (De Luca et al., 2018b). GRL is in principle applicable to such districts by considering an appropriate tissue model, but its feasibility remains to be proven. Regarding the application to the brain, recent works suggest to account for three tissue types (Dhollander et al., 2016;Jeurissen et al., 2014) as also considered in this work. However, differently from the mentioned studies we did not estimate the signal representation of each tissue type by using tissue masks, but rather used models with plausible values to represent the considered tissues. Gray matter and cerebrospinal fluid were modeled using the ADC equation and typical representative values taken from literature (Chiapponi et al., 2013;Pasternak et al., 2009). The use of literature-based values has the advantage not to require any additional data and processing beyond the diffusion dataset, standardizing the analysis, but can be inaccurate if any change in the diffusivity of the two tissues would occur. While this was not taken into account in this study, it is straightforward to estimate the diffusivity values in tissue specific masks prior to applying GRL, in analogy with what is done in other methods . With regards to the modeling of WM, the parameters of the DTI (GRL-VT) and DKI based (GRL-DKI) H-matrices were estimated from the data itself using a WM mask with coherent WM (fractional anisotropy > 0.7), which might be suboptimal to select only single population WM bundles. Similarly, the parameters used for the singleshell DTI (GRL-FT) and the NODDI (GRL-NODDI) based H-matrices were assumed from previous studies (Dell 'Acqua et al., 2007;Zhang et al., 2012), and might be sub-optimal when applied to the considered datasets. Nevertheless, we have here shown that GRL allows to perform spherical deconvolution while allowing large freedom to the user in terms of model choice, which parameters can be further optimized in specific applications.
In summary, GRL improves the performance of the Richardson-Lucy spherical deconvolution scheme while also estimating signal fraction maps of multiple water compartments. The derived fractional maps can be used as surrogate to existing techniques for multi-compartment quantification, as free water elimination (Berlot et al., 2014;Henriques et al., 2017;Pierpaoli and Jones, 2004) , which have been shown promising in previous clinical studies (Coelho et al., 2018;Hoy et al., 2017). Further, the estimated signal fraction maps of WM and GM can be used as termination criteria for fiber tractography, relieving from the definition of arbitrary thresholds on the fiber orientation distributions.

Conclusions
In this work we presented the GRL framework, a method to perform spherical deconvolution of multishell diffusion MRI data while accounting for partial volume effects to provide better estimation of fiber orientation distributions in white matter. We have shown that GRL can disentangle signal fractions from multiple tissue components, such as white matter, grey matter, cerebrospinal fluid and, if sufficient data is acquired, pseudo-diffusion. The presented framework allows flexibility in terms of acquisition protocols and can accommodate different choices of diffusion models to perform the spherical deconvolution, making it potentially applicable to any tissue for which a signal model is available.

Figure 1.
The effect of multi-shell weighting on a simulated fiber configuration including GM-like partial volume (20%) and 2 WM-like fibers crossing with an angle of 47 degrees (red lines). GRL-FT was applied 1000 times for different noise realizations at SNR 50, then the main peaks were extracted (yellow lines). The angular error increased with increasing weighting of the lower diffusion-weighted shells, but weighting values below = 0.2 did not allow to effectively separate partial volume effects.

Figure 2.
The mean (solid lines) and the standard deviation (the shadowed area) of the signal fractions estimated with GRL-FT in simulations. The WM-like signal fraction was increased from 0.1 to 1 and mixed with CSF-like signal (Simulation I), GM-like signal (Simulation II), and both CSF and GM-like signals (Simulation III). GRL-FT recovered the simulated signal mixtures with minimal bias for all SNR levels, a small but consistent overestimation of the WM signal fraction (fWM) and underestimation of the GM signal fraction was observed (fGM).

Figure 3.
The mean (solid line) and the standard deviation (shaded area) of the angular deviation of the FODs estimated with dRL (first row) and GRL-FT (second row) in correspondence of multiple SNR levels. WMlike signal was mixed with decreasing GM-like partial volume (Simulation I, first column), decreasing CSF-like partial volume (Simulation II, middle column), and both effects simultaneously (Simulation III, last column). If sufficient SNR is provided (e.g. SNR above 20), GRL-FT can remarkably reduce the angular error as compared to dRL in presence of strong partial volume effects.  An example axial slice of the signal fractions estimated with GRL-FT on Data3 when also including pseudo-diffusion effects into the H-matrix formulation to quantify IVIM effects in addition to WM, GM and CSF. The estimated IVIM contributions are non-zero throughout the brain, and are in-line with known anatomy, with signal fraction values around 5% in tissues.   An example axial view of the FODs estimated with dRL and GRL on Data2 with focus in a frontal region where the WM enters the GM. The FODs are colored encoded according to the conventional diffusion directional color scheme and overlaid on the first non-weighted image. GRL generally suppressed the effect of partial volume effects onto the FOD with higher efficacy than dRL. Figure 10.
The number of peaks detected with the FODs computed with dRL and GRL (in correspondence of different H-matrix choices) on Data1. A) an example coronal slice showing the number of detected peaks and B-C) the peak frequency in WM and GM, respectively. GRL-FT resulted in slightly reduced number of 3+ fiber crossing as compared to dRL, in particular in peri-cortical regions. The presented results are computed with a minimum peak value threshold equal to 20% the average amplitude in WM of an FOD derived with dRL. Figure 11.
Rows 1-2) the results of whole brain fiber tractography obtained with dRL and GRL on Data1. Fiber tractography with dRL is terminated using FOD thresholding, whereas tractography with GRL is free of thresholds and uses the WM-like (first row) or the GM-like (second row) signal fraction maps as termination criteria. The third row shows a coronal slab of the tracking, which for GRL was colored according to the traversed signal fraction. Tracking with GRL-FT was comparable to dRL in WM. Conversely, with dRL it is not possible to selectively allow tracts to enter the cortex without producing a large number of spurious tracts (second row). The choice of the H-matrix has a major impact on the fiber tractography.

Figure 12.
Rows 1-2) the results of whole brain fiber tractography obtained with dRL and GRL on Data2. Fiber tractography with dRL is terminated using FOD thresholding, whereas tractography with GRL is free of thresholds and uses the WM-like (first row) or the GM-like (second row) signal fraction maps as termination criteria. The third row shows a coronal slab of the tracking, which for GRL was colored according to the traversed signal fraction. Tracking with GRL-FT was comparable to dRL in WM. The lower resolution of Data2 is likely to cause more severe partial volume effects, which when taken into account when GRL-FT allow a finer transition into the cortical layer. Even considering tractography only in WM, the first row shows that GRL removes spurious fibers in the frontal lobe otherwise reconstructed with dRL.