NODDISH: a computational efficient NODDI extension for fODF estimation in diffusion MRI
Abstract
Diffusion Magnetic Resonance Imaging (DMRI) is the only noninvasive imaging technique which is able to detect the principal directions of water diffusion as well as neurites density in the human brain. Exploiting the ability of Spherical Harmonics (SH) to model spherical functions, we propose a new reconstruction model for DMRI data which is able to estimate both the fiber Orientation Distribution Function (fODF) and the relative volume fractions of the neurites in each voxel, which is robust to multiple fiber crossings. We consider a Neurite Orientation Dispersion and Density Imaging (NODDI) inspired single fiber diffusion signal to be derived from three compartments: intracellular, extracellular, and cerebrospinal fluid. The model, called NODDISH, is derived by convolving the single fiber response with the fODF in each voxel. NODDISH embeds the calculation of the fODF and the neurite density in a unified mathematical model providing efficient, robust and accurate results. Results were validated on simulated data and tested on invivo data of human brain, and compared to and Constrained Spherical Deconvolution (CSD) for benchmarking. Results revealed competitive performance in all respects and inherent adaptivity to local microstructure, while sensibly reducing the computational cost. We also investigated NODDISH performance when only a limited number of samples are available for the fitting, demonstrating that 60 samples are enough to obtain reliable results. The fast computational time and the low number of signal samples required, make NODDISH feasible for clinical application.
1 Introduction
The diffusion process of water molecules in brain tissues reflects the cytoarchitecture of the observed tissue. Neurites density and axon fibers topology are two of the main factors which influence the diffusion process in the white matter. With Diffusion Magnetic Resonance Imaging (DMRI) it is possible to indirectly observe the Ensemble Average Propagator (EAP) of the water molecules movements and reconstruct the white matter fiber Orientation Distribution Function (fODF) [1] as well as numerical indexes reflecting the neurites density in each voxel [2]. The proposed approach provides a novel method for the joint estimation of the fiber orientation and neurite densities also in presence of complex topologies (crossing, fanning) while keeping low computational complexity which ensures fast processing. In addition, robustness to the reduction in the number of samples enables clinical exploitability providing a powerful computational tool for microstructure characterization.
So far, different models have been proposed to solve the two problems of fODF and volume fraction estimation individually. The most clinically used model is the Diffusion Tensor Imaging (DTI) [3], [4] which considered the EAP as a single multivariate Gaussian function. DTI enables the identification of only one single principal direction of diffusion and thus is not suitable to model complex WM fiber configurations such as fiber crossings, kissing, and fanning [5].
Spherical Deconvolution (SD) methods [6], [7], [8] consider the diffusion signal as the result of the convolution of the fODF times a single fiber response. By deconvolving the diffusion signal with a given response function (usually obtained from single fiber areas of the brain) it is possible to estimate the fODF. However, the deconvolution process is an illconditioned problem which can be affected by noise or error in the estimation of the response function. Constrained Spherical Deconvolution (CSD) [1] added nonnegative constraints to the fODF estimation, increasing the resolution of the fODF which is extremely important for tractography application [9].
The EAP represents only an indirect measure of the different components of the tissue. Several models have been proposed to explicitly quantify the contribution of the different compartments present in each voxel in order to provide a map of the different tissues. A noteworthy example of multicompartment model is the Composite Hindered and Restricted Model (CHARMED) [10], which considers two different types of substrate in which the water molecules can diffuse: intracellular and extracellular. The intracellular compartment is the one in which the water diffusion is restricted by the cell walls, e.g. axons and cell bodies, and it is modeled as an axially symmetric cylinder. The extracellular compartment captures the contributions of all the water molecules trapped inbetween different neurites and in which the diffusion is only hindered, and it is modeled as a multivariate Gaussian. Two important modifications were added to the CHARMED by [11] and [12]. [11] extended CHARMED with the addition of a third compartment in order to fit the cerebrospinal fluid (CSF), while [12] replaced the coherently oriented cylinders of CHARMED with a Watson distribution of cylinders in order to model fiber dispersion. Such distribution is limited by the fact that it presents only one principal direction of diffusion and cannot be used to model crossing fibers. Histological evidence shows that the human axonal diameter is normally less than two micrometers [13], [14]. Cylinders presenting such diameters exhibit almost no signal decay using gradient strength of mT/m in Pulse Gradient Spin Echo (PGSE) acquisitions with a single diffusion time [15], [16]. Therefore, in the Neurite Orientation Dispersion and Density Imaging (NODDI) [2] model the cylinders of [12] were replaced with a Watson distribution of “sticks” (cylinder of zero radius, equivalent to a Gaussian function with zero perpendicular diffusion). NODDI was able to provide good results on simulated data for retrieving the volume fractions of the different compartments. Improvements in computational efficiency were reached by the reformulation within a convex optimization framework proposed by [17], called Accelerated Microstructure Imaging via Convex Optimization (AMICO). Moving from the original nonlinear optimization to a linear optimization greatly reduced the computational time of both techniques, without downgrading the performance of the models. The Watson distribution suffers from the fact that it can only model isotropically dispersed distributions, and not anisotropic fiber fanning and crossings. In order to overcome this limitation, [18] replaced the Watson distributions with a Bingham distribution, allowing to model fiber fanning, but not crossings. A mixture of Bingham distributions was previously adopted also in [19] and [20], but the problem of finding the number of Bingham distributions to use in each voxel, and the hardness of the fitting make these techniques less exploitable for practical cases. A general framework for calculating any kind of probability density function of fibers distribution on the sphere, given the model of a single fiber response, was proposed in Fiber ORientation Estimated using Continuous Axially Symmetric Tensors (FORECAST), [21]. FORECAST exploits SH to model any distribution of axially symmetric tensors on the sphere. To our knowledge, FORECAST is also one of the first models exploiting the mean value of the signal to estimate rotation invariant fiber specific microstructural parameters (in this case parallel and perpendicular diffusivities). Jespersen et al. (2007) [22] advanced the FORECAST model considering instead of a single axially symmetric tensor, a twocompartmental model composed of an axially symmetric tensor and an isotropic tensor for modeling the single fiber response. Neurite density estimation obtained using this model was compared with light microscopy and electron microscopy images in [23]. FORECAST parameters estimation which uses the mean of the signal on each bvalue for estimating microstructural parameters has been recently called Spherical Mean Technique (SMT), by [24] and [25]. In particular in [25] employ the same strategy as FORECAST in order to estimate the parallel and perpendicular diffusivity in each voxel. [26] extended the SMT to other rotational invariant features, combining them with machine learning to estimate several microstructural parameters from the diffusion signal.
In this work, we expanded the FORECAST single tensor model to a NODDIlike multicompartment model, gathering the advantages of both techniques in a unified model called NODDISH. NODDISH is able to obtain accurate microstructural information, such as the intracellular, extracellular, and CSF volume fractions, in voxels presenting any kind of fiber distribution, from crossing to fanning. Unlike the classical NODDI, NODDISH permits to reconstruct at the same time also the fODF in each voxel, making it suitable for tractography and connectomic studies. Thanks to the SMT, NODDISH elegant mathematical framework enables an extremely fast computation of the volume fractions and fast fODF reconstruction. NODDISH results were compared to those provided by NODDI and CSD, respectively, for the volume fractions and the fODF, demonstrating competitive results in both respects. We also tested NODDISH performance when only a limited number of samples are provided for the fitting. In addition, robustness to the reduction in the number of samples was determined in order to establish the clinical exploitability of the proposed model.
The manuscript is organized as follows. Section 2 revisits the FORECAST model and introduces the NODDISH model, as well as the simulated and invivo data used for the experiments. Section 3 shows the comparison of NODDISH with NODDI in the simulated and invivo experiments, as well as the robustness of the NODDISH with respect to the number of samples. Section 4 discusses our findings, and Section 5 summarizes the conclusions of our work.
2 Methods
Assuming that in the considered voxel the fiber population is homogeneous, the diffusion signal can be modeled as
(1) 
where is the diffusion signal of a fiber aligned along v and is the fODF. Since is probability density function it must satisfy the following conditions
(2)  
(3) 
In diffusion MRI, the fODF is generally modeled using Spherical Harmonics (see Appendix A) [1], [21], [25]. Replacing the SH equation into equation (1) we obtain
(4) 
The choice of the single fiber signal will define the resulting signal basis. In this section we will consider two cases. In the first case is modeled as an axially symmetric Gaussian function (FORECAST), while in the second it results from the linear combination of the contribution arising from a three compartments model (NODDISH).
2.1 FORECAST model
FORECAST considers to be an axially symmetric tensor, with parallel diffusivity and perpendicular diffusivity
(5) 
where is the bvalue at the given diffusion time , and is the diffusion tensor aligned along v
(6) 
Replacing equation (5) in equation (4), we can calculate as
(7) 
Using the FunkHecke Theorem [27] it is possible to solve the integral in equation (7) (see Appendix B for the complete derivation) obtaining the final FORECAST model
(8) 
where is the solution of , with the Legendre polynomial of order .
The estimation of the fODF coefficients requires the estimation of and in each voxel. In order to estimate these diffusivities, [21] proposed to exploit the average of the diffusion signal for a given bvalue, which is given by the element of the FORECAST basis
(9) 
Equation (9) exploits the fact that the coefficient is equal to , as detailed in Appendix B.0.1, in order to fulfill the requirements of equation (2). If the signal is acquired in more than one shell (e.g. using more than one bvalue), it is possible to use equation (9) to estimate and accurately.
2.2 The NODDISH model
The NODDISH model considers the single fiber diffusion signal as the weighted sum of three tissue compartments: intracellular , extracellular , and cerebrospinal fluid
(10) 
where , , and are the relative volume fractions in the tissue, with the constraint . The compartment is modeled as a tensor with , as in [2], the extracellular compartment as a tensor with a certain parallel diffusivity and [11], and the CSF compartment as an isotropic tensor with mm/s
(11) 
Replacing equation (10) into equation (4) (see Appendix C) leads to the final NODDISH model formulation
(12) 
in which the cerebrospinal fluid compartment volume fraction , influences only the first coefficient of the basis (the one referring to the isotropic SH), while intra and extra cellular compartments affect all the coefficients.
Using the SMT [21], [25], [24], the first element of the basis at and represents the mean value of the diffusion signal at a given bvalue
(13) 
The estimation of the volume fractions is performed as follows: we initially create a dictionary of 383 combinations of possible volume fractions. The dictionary was created by sampling in a nonuniform manner the volume fractions space. In particular, we considered the number of combinations of and to be proportional to . Using this method the dictionary presents a higher variability of intracellular and extracellular volume fractions where the CSF volume fraction is low, and a lower variability where the CSF volume fraction is high. For each combinations of these volume fractions we calculated the estimated using equation (13), for each bvalue. In order to select the best combination, we evaluated the mean of the diffusion signal at each bvalue as with N the number of samples on each shell and thus selected the set of volume fractions which better approximates the average diffusion signal. The evaluation of equation (13) is extremely fast, and even if we test all the 383 combinations of volume fractions, the total estimation time is of the order of fraction of second per voxel (see Section 3.4 for computational time details).
2.3 Basis coefficients estimation
The estimation of the basis coefficients for both FORECAST and NODDISH is performed via convex optimization using CVXOPT^{1}^{1}1http://cvxopt.org allowing to enforce the positivity constraint as in equation (3) on a spherical grid composed of 181 points. In detail, our objective is to find
(14) 
subject to
Yc  (15)  
(16) 
where is the measured diffusion signal, c is the coefficients vector, M the NODDISH or FORECAST basis matrix, and Y is the SH basis matrix. We implemented both algorithms under the Diffusion Imaging in Python (DIPY) [28], [29] software library^{2}^{2}2http://dipy.org and the code is available on request. The CSD response function was fixed for all the brain and modeled using the average eigenvalues of the diffusion tensor in voxels presenting fractional anisotropy higher than 0.8. We adopted the same strategy to estimate NODDISH parameter. In simulated data we set mm/s and mm/s. SH basis of order 8 was chosen for both techniques.
NODDI results were calculated using the standard NODDI toolbox^{3}^{3}3http://mig.cs.ucl.ac.uk/index.php?n=Tutorial.NODDImatlab. In NODDI model [2] intracellular volume fraction is weighted by a factor in order to be compared to NODDISH [30]. In our simulations this normalization led to more accurate estimation for the NODDI model.
2.4 Simulated dataset
To generate a model of white matter, we simulated the diffusion signal as the weighted sum of an intracellular and an extracellular contribution. The intracellular part was modeled as a cylinder with diameter using the Multiple Correlation Function framework [31], [32], [30], while the extracellular part was modeled as a symmetric Gaussian. Both the cylinder and the Gaussian share the same mm/s, while the extracellular depends on the extracellular volume fraction as detailed in Section 2.2.
To model different cases of dispersion and crossing, we initially extract random directions from the Kent distribution. Kent distribution [33], [30] is able to model symmetric probability density functions on the sphere. In particular Kent distribution is modeled as
(17) 
where is the concentration parameter which rules the degree of orientation dispersion and is the parameter which models the anisotropy of such a dispersion. The parameter is bounded between 0 and , with meaning perfectly isotropic dispersion and completely anisotropic dispersion, respectively. The variable represents the central direction of the pdf, while the two orthogonal unit vectors and define the orientation of the anisotropic dispersion. The function is defined as
(18) 
with the Gamma function and the modified Bessel function. Figure 1 shows some examples of the Kent distribution^{4}^{4}4Kent distribution has only one lobe. We added the antipodal symmetric lobe for visualization purposes. , given the parameters and .
,  ,  , 
,  ,  , 
To model dispersion and crossing we extract random directions from the Kent distribution. For each of these directions, we simulate the white matter signal as the sum of a cylinder and a Gaussian aligned in the direction . The total signal is than calculated as
(19) 
where the directions , the intracellular volume fraction and the extracellular volume fraction are given. In our simulations we consider the following parameters: , , three different rotations of the Kent distribution on its axis, and 11 different orientations of the Kent distribution for a total of 297 combinations of parameters only for the Kent distribution. We tested also different combinations of intracellular and extracellular volume fractions, with and . We added Rician noise with a Signal to Noise Ratio (SNR) equal to 20 (10 different instances per combination of parameters). The total number of voxels generated with such a scheme is 26730.
With this model, we can simulate any kind of fiber fanning, but not fiber crossings. In order to overcome this limitation, we considered two Kent distributions, both with and but centered on two different directions and (see Figure 2). The total signal equation in this case becomes
(20) 
Crossing angles of 90, 60, and 45 degrees were considered, and, as in the previous case, the orientation of the crossing was changed in 11 different directions. The same volume fractions combination of the “single fiber” dataset was used for the crossing. Rician noise at was added to the signal in ten different instances. The total number of voxels, in this case, was 2970.
All the simulations were generated using the Human Connectome Project (HCP) [34] sampling scheme, which will be presented in the next subsection.
2.5 Invivo dataset
In order to test the NODDISH model invivo we considered one subject of the HCP diffusion MRI dataset. HCP acquisition scheme is composed of 18 volumes sampled at bvalue 0 s/mm, plus 90 volumes at bvalue 1000, 2000, and 3000 s/mm, for a total of 288 volumes. Each of the 270 volumes which composed the three shell are acquired using an unique gradient direction. HCP pulse separation time ms and pulse width ms corresponding to an effective diffusion time ms. Each volume consists of voxels with a resolution of 1.25mm.
HCP scheme was created using an incrementally optimal sampling [35]. This means that considering only the first samples, still provides an almost optimal angular coverage. This property was exploited for reducing the number of qsamples as detailed in Section 3.3. In particular, two subsets of and samples per shell were considered and the corresponding signal was compared to that obtained using the fully sampled data ( samples). In the same section, we also considered the case where only the first two shells with bvalue 1000 and 2000 s/mm are provided to the model. The results of the NODDISH fitting were evaluated by calculating the mean square error (MSE) between the original HCP signal and the NODDISH signal for all the samples, also those which were not used for the fitting.
3 Results
3.1 Results of simulated data
ISOTROPIC FANNING  ANISTROPIC FANNING 

Figure 3 (left column) shows the estimation of the intracellular volume fraction for NODDISH and NODDI in the case of isotropic dispersion (), compared to the ground truth (GT) values (black circles). NODDISH is able to provide the most accurate estimation of in case of highly dispersed fibers (), while NODDI in case of . For NODDISH, the estimation error increases while decreasing the dispersion ( and ). Both models tend to slightly overestimate in all conditions. The average errors for NODDISH are 3.4, 3.1, and 1.6 for , , and respectively. The same errors for NODDI are 3.4, 2.2, and 3.7. The standard deviation (plus sign in the graphs) of NODDI is slightly lower than NODDISH, with an average standard deviation of 2.1 and 2.7 for NODDI and SHADOW, respectively. In case of anisotropic dispersion (Figure 3, right column) NODDISH estimation is closer to the GT value with respect to the case of isotropic dispersion. NODDI model provides better estimation of for , but progressively underestimates the intracellular volume fraction as the anisotropy of the dispersion decreases. NODDISH, on the contrary, presents an inverse pattern with its best estimation for case, while the largest error is reached for . The average errors for the , , and are respectively 2.8, 2.2, and 1.5 for NODDISH and 2.6, 3.3, and 3.8 for NODDI. The average standard deviations in this case are 1.8 and 2.5 for NODDI and NODDISH, respectively.
CROSSING FIBERS  

Figure 4 shows the estimation of for the crossing fibers dataset. In this case, the estimation errors of both, NODDISH and NODDI, increase with respect to the fiber fanning dataset. Although estimating the intracellular volume fraction for crossing fibers represents a more difficult task with respect to the single fiber case, NODDISH estimation results slightly more accurate with respect to NODDI. In this dataset, average errors across crossing angles are 7.58 and 6.24 for NODDI and NODDISH, respectively. The amplitude of the crossing angle does not seem to have an impact on the estimation for both models. The average standard deviations in the crossing fiber case are 2.8 and 3.9 for NODDI and NODDISH, respectively.
CSD  FORECAST 
NODDISH  AVERAGE 
In order to asses the precision of the fODF, the angular error (AE) between the fODF peaks and the ground truth mean directions of the simulated fibers was determined. Figure 5 shows the AE estimated for CSD, FORECAST, and NODDISH fODFs. In particular, the AE was estimated as the average of the angular distance of the two ground truth directions with the closest fODF peak. CSD was able to retrieve the principal directions of diffusion almost perfectly in the 90 and 60 degrees cases () while its AE increased in the most difficult case, the 45 degrees crossing, with an average AE of 7.9 and the highest standard deviation. For both, FORECAST and NODDISH, the AE was less than 5 in the 90 and 60 degrees crossing cases. In the 45 degrees crossing the error was lower for low values and increased progressively with the neurite density. CSD performed better than FORECAST and NODDISH in voxels featuring high restriction. However, it is important to stress that CSD results depend on the choice of the response function. In our experiments, we purposely selected the response function to be very close to the intracellular signal in order to obtain highresolution fODFs. Both NODDISH and FORECAST are able to adapt its response function in every voxel, which is a clear advantage in the case of real tissues.
3.2 Results on HCP data
NODDISH three compartments model is designed to capture the heterogeneity of the brain tissues. However, there are no guarantees that NODDISH estimation of neurites density and volume fraction corresponds to the underlying tissue composition, as it is the case for all models. Therefore, results on invivo data should be taken with care since we lack a histological ground truth to validate these data. Due to the lack of ground truth, the performance of NODDISH in invivo data was assessed by comparison with NODDI that was used as a benchmark.
Figure 6 shows the intracellular and CSF volume fractions estimated using NODDI and NODDISH, respectively. The maps obtained with the two techniques are very similar. NODDISH intracellular volume fraction appears to be slightly more contrasted with respect to NODDI, in particular in single fiber areas such as the corpus callosum and the corticospinal tract. NODDISH CSF map presents more voxels with higher partial volume of CSF near the cortex, but in general, the agreement between the two maps is extremely high. For both models, CSF fraction maps present low, but nonzero, values in the white matter. Although the presence of free water between the axons is not likely, these values could be explained by an intrinsic error of the three compartment model used by both techniques.
INTRACELLULAR  CSF  
NODDI 


NODDISH 

Figure 7 shows the MSE calculated between the normalized diffusion signal and the estimated diffusion signal of NODDI and NODDISH. It is worth to mention that in this picture we use a very narrow range of values in order to emphasize the contrast and that both techniques were able to estimate the diffusion signal accurately. NODDI MSE appears to be higher in white matter, and in particular in single fiber areas, with respect to gray matter or CSF where it is close to zero. NODDISH MSE presents a more uniform pattern in the white matter, gray matter, and CSF. The highest MSE values for NODDISH were found in the corpus callosum and in the basal nuclei.
NODDI MSE  NODDISH MSE 

Figure 8 (left) highlights the regions of the brain in which NODDISH features a lower MSE than NODDI. Almost all the white matter appears to be white, which seems to indicate that NODDISH is able to reconstruct the diffusion signal more accurately. These results could be explained by the fact that the Watson distribution is not able to model most of the brain fiber configurations, while NODDISH can model any kind of fiber arrangement including fanning and crossing. In most of the gray matter and CSF NODDI fits the diffusion signal slightly better than NODDISH, although the absolute value of the difference of the MSE is lower than that obtained in white matter (Figure 8, right).
NODDI MSE NODDISH MSE  NODDI MSE  NODDISH MSE 

ROI  CSD 
FORECAST  NODDISH 
Figure 9 shows the fODFs obtained by CSD, FORECAST, and NODDISH on a region of interest (ROI) of the HCP brain (top left) which includes crossings, fannings, a part of the corpus callosum, and the right ventricle. There are not substantial differences between the three fODFs. The main advantage of using NODDISH to calculate the fODF is that it is possible to directly filter out the fODFs which do not belong to white matter. In Figure 9 (bottom right) NODDISH fODFs, for instance, only the voxels with were considered for display. In order to do the same with CSD fODFs, it would be necessary to integrate the CSD results with those obtained by a multicompartmental model, such as NODDI.
3.3 Robustness with respect to the number of samples
HCP acquisition scheme consists of 270 gradients plus 18 volumes acquired at bvalue = 0 s/mm. Such acquisition is far from being considered adapt to clinical practice, which normally consists of 60 diffusion weighted images. In order to test the ability of NODDISH to fit more practical acquisition schemes, we subsampled the HCP data. Figure 10 shows the mean signal value of the HCP data for each bvalue, while keeping all the gradients (90 directions), 60 directions, and only 30 directions, respectively. As expected, the diffusion signal decays for increasing bvalues. In regions of fast diffusion, such as the CSF, the diffusion signal is almost zero even at bvalue = 1000 s/mm. At bvalue = 2000 s/mm it is possible to see some contrast between white matter (high diffusion restriction) and gray matter regions, which is more pronounced at bvalue = 3000 s/mm.
90 directions  60 directions  30 directions  






At a first glance, it is impossible to distinguish between the mean value images obtained using a different number of samples. This because the mean difference in a white matter region between the mean signal values at bvalue = 1000 s/mm derived from the 90 and the 30 samples cases is only of . The mean signal over each shell is necessary for NODDISH to retrieve the volume fractions, and the robustness of the mean signal to subsampling reflects the robustness of the volume fractions estimation.
Figure 11 shows the variation of the intracellular volume fraction calculated using NODDISH when changing the number of directions per shell and the number of shells (bvalue = [1000, 2000, 3000] s/mm or bvalue = [1000,2000] s/mm ). As expected, the intracellular volume fraction does not change significantly with the number of samples when all the shells are kept (=3000 s/mm). On the contrary, excluding all the samples at bvalue = 3000 s/mm affects the estimation of , in particular if we consider only 30 samples per shell. The map in this case results more noisy. Nevertheless, the overall quality of the image does not degrade too much with respect the fully sampled image. This result is extremely promising, considering that we move from a 270 directions acquisition with s/mm to a 60 directions acquisition with s/mm (more than 76 samples less).
90 directions  60 directions  30 directions  
s/mm 

s/mm 

In order to quantify the error that is introduced while decreasing the number of samples, we calculated the MSE between the reconstructed signal obtained after fitting NODDISH to the subsampled dataset and the full diffusion signal dataset. It is worth to mention that calculating the MSE in such a way clearly advantages the NODDISH model which fits the 90 directions, s/mm data because the MSE is evaluated on the same points used for the fitting. On the contrary, when only 30 directions with s/mm are used for the fitting, the NODDISH model requires to blindly estimate 210 points in order to reconstruct the complete diffusion signal. Figure 12 shows the MSE for the reconstructed signal for the different sampling schemes. As for the maps, it is challenging to perceive any difference.
90 directions  60 directions  30 directions  
s/mm 

s/mm 

In order to better highlight the changes, we calculated the normalized histogram of the MSE image in voxels presenting high level of intracellular volume fraction (). Figure 13 shows the histograms of the MSEs obtained from the different sampling schemes. As it possible to see, 90 directions, s/mm presents the highest number of voxels with low MSE, followed by the 60 directions, s/mm, and the 90 directions, s/mm. The worst schemes are the ones using only 30 directions per shell. The high MSE of the 30 directions scheme is probably due to the low angular resolution which makes it less suited to capture crossing fibers signal.
3.4 Computational time
NODDISH volume fractions estimation and fODF are substantially similar to NODDI and CSD results, respectively. However, to obtain the same information provided by NODDISH it is necessary to run both NODDI and CSD on the same data. Table 1 shows the computational time for CSD, FORECAST, NODDI, and NODDISH obtained on the 2970 voxels of the simulated crossing dataset. The computational time was computed using a single core on a laptop equipped with an Intel®Core™i76500 CPU @2.50GH. DIPY CSD implementation is extremely fast, followed by NODDISH which is six times slower. NODDI is the slowest model (1.5 seconds per voxel). FORECAST is the slowest model for what concerns fODF estimation. Compared to the default NODDI toolbox, the calculation of the volume fractions for NODDISH is 237 times faster.
It is worth mentioning that using AMICO [17] it is possible to fit the NODDI model extremely fast, less than 10 minutes for a full brain, according to the original paper [17]. However, AMICO is not a model in itself but a convex optimization framework and does not solve the theoretical limitations of NODDI in modeling only isotropic fanning fibers. However, if the goal is only to estimate the volume fractions and not the fODF, only 1.30s are necessary for the NODDISH to fit 2970 voxels. This makes NODDISH even faster than NODDI implemented using AMICO. A complete comparison of the performance of the NODDISH with respect to AMICO is out the scope of this paper and will be the aim of future research.
2970 voxels (s)  seconds per voxel  

CSD  3.13  0.001 
FORECAST  61.05  0.021 
NODDI  4419.45  1.488 
NODDISH  18.63  0.007 
NODDISH (no fODF)  1.30  0.0004 
4 Discussion
NODDISH embeds both microstructural parameters estimation and fODF reconstruction in the same mathematical framework. We demonstrated that NODDISH was able to perform on par with NODDI for what concerns the volume fraction estimation, and provided comparable performance with respect to CSD for the fODF estimation. Considering CSD and NODDI individually, it is impossible to obtain both the volume fractions and the fODF at the same time as in the case of NODDISH. The possibility to obtain all these information using a single model can greatly simplify the processing pipeline for diffusion MRI data. Another important insight which can be derived from NODDISH is that using the SMT, the estimation of the volume fractions depends only on the mean value of the signal. Using this feature, the NODDISH parameters fitting is not only fast, but also requires few sampling points (60 directions spread on two shells) in order to obtain good results. NODDISH efficiency and speed make it suitable to be used in clinical application. We also showed that 30 samples are enough to provide a good approximation of the average signal for each shell. Acquisition schemes which aim at detecting the volume fractions should, therefore, increase the number of bvalues, using only a moderate number of samples per shell.
As mentioned in Section 3.2, there is no guarantee that the volume fractions estimated using compartmental models are linked to real changes in the tissue composition in each voxel. The fitting procedure always provides a result, regardless of the correctness of the model itself. Low MSE, as the one presented in Figure 7, can be considered only as a general indication of the goodness of the model, and not the proof that the underlying tissue architecture corresponds to the actual model parameters. In this work, our goal was to compare the NODDISH with the classical NODDI model. NODDI exploits several assumptions for reducing the number of parameters to be fit that have been shown to be wrong for the human white matter. In particular the fact that the parallel diffusivity is set to mm/s has been shown to be wrong in [36], [37], [26], and [25]. Another oversimplification used in NODDI is the direct dependence of the perpendicular diffusivity to the intracellular volume fraction which is generally not true [37]. In our implementation, NODDISH shares the same formulation of NODDI model and, therefore, some of its limitations. More accurate results would be obtained by making voxelwise adaptive as in [24], [26], [36], and [37]. In this case, the only modification needed is the creation of new entries in the volume fractions dictionary with different . For what concerns the microstructural parameters, [26] Bayesian estimation of the volume fraction and water diffusivity seems to provide more reliable results with respect to NODDI. As it is the case of NODDISH, [26] exploit the mean value of the signal, as well as other rotation invariant features, in order to estimate such parameters. This approach could be easily integrated into NODDISH potentially improving the resulting fODF and MSE.
In this paper, we considered only the HCP acquisition scheme and subsampling from it. Some model assumptions may not be valid in the case of different acquisitions. For example, we used a “stick” to model the intracellular compartment. Stick signal presents no decay in the perpendicular plane. While this model can be realistic in the case of small axon diameters at low bvalues, it could not be the best choice for larger axons (e.g. 4 ) at higher bvalues (in the order of at least s/mm). Another limitation is the fact that the HCP scheme uses a single diffusion time. In the case of multiple diffusion times, the Gaussian assumption for the extracellular compartment may not be true [38], [39], and a different model should be used instead. Thanks to the flexibility of the NODDISH these diffusion time dependent model could be easily integrated into its framework.
5 Conclusions
NODDISH is a flexible model that allows obtaining both tissues volume fraction parameters and the fODF at the same time. Using NODDISH, the signal representation is analytically derived from the convolution of a NODDI inspired single fiber response with the fiber ODF modeled using SH. In this work, we explicitly derived a single fiber response function containing three compartments (isotropic, intracellular and extracellular compartments, respectively). The result of the convolution forms a basis in the diffusion signal space. The basis first element can be used to estimate the volume fraction of the threecompartments model. The current NODDISH model can be viewed as an extension of generic distribution of fibers of previously proposed compartmental models, such as NODDI, with an efficient parameters estimation procedure that is efficient and robust to fanning and crossing fibers. The processing pipeline of big datasets, such as the HCP, can be drastically improved using NODDISH. Results show that NODDISH can also be used for dataset presenting a limited number of samples, which makes it feasible for processing clinical data.
6 Acknowledgements
Data were provided by the Human Connectome Project, WUMinn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.
Appendix A Real Symmetric Spherical Harmonics
The real symmetric SH basis is defined as
(21) 
where is the general SH basis, written as
(22) 
with the polar representation of u, and the associated Legendre Polynomial. A possible way to represent is to consider its SH series expansion
(23) 
Appendix B Analytical derivation of the FORECAST model
In the FORECAST model the diffusion signal can be calculated as
(24) 
with . The tensor can calculated as
(25) 
assuming that the two perpendicular eigenvalues are equal. We can now calculate as
(26) 
Using the FunkHecke Theorem [42], [27] stating that, given a unit vector u, an integral on the sphere can be converted to a one dimensional integral
(27) 
assuming continuous on and any SH of order , and with the Legendre polynomial of degree .
Replacing it into FORECAST equation the problem became:
(28) 
The only step left is to solve . This integral does not have a general solution, but we can solve it for each single Legendre polynomial of a certain order. The first eight even Legendre polynomials are

=

=

=

=

=
Considering , we can now calculate each integral necessary to integrate each Legendre polynomials
(29) 
When is indeterminate but the following limits hold true
(30) 
We can define as the value of the Legendre polynomials when , which for the first 8 Legendre polynomials corresponds to

=

=

=

=

=
The FORECAST signal equation can then be written as
(31) 
where N [0,8].
b.0.1 Special case b=0
In the case, the signal equation can be rewritten as
(32) 
This because for and zero otherwise, and . In order to have a normalized signal with , must be equal to . This is also the same condition which we need to verify equation (2):
(33) 
this because
(34) 
where is equal to one when and and zero elsewhere. This equation is valid only for .
Appendix C Analytical derivation of the NODDISH
References
 [1] Tournier, J.D., Calamante, F., Connelly, A.: Robust determination of the fibre orientation distribution in diffusion mri: nonnegativity constrained superresolved spherical deconvolution. NeuroImage 35(4) (2007) 1459–1472
 [2] Zhang, H., T., S., WheelerKingshott, C., Alexander, D.: NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61(4) (2012) 1000–1016
 [3] Basser, P.J., Mattiello, J., LeBihan, D.: Estimation of the effective selfdiffusion tensor from the nmr spin echo. Journal of Magnetic Resonance, Series B 103(3) (1994) 247–254
 [4] Pierpaoli, C., Jezzard, P., Basser, P.J., Barnett, A., Di Chiro, G.: Diffusion tensor mr imaging of the human brain. Radiology 201(3) (1996) 637–648
 [5] Descoteaux, M. In: High Angular Resolution Diffusion Imaging (HARDI). John Wiley & Sons, Inc. (2015)
 [6] Tournier, J.D., Calamante, F., Gadian, D.G., Connelly, A.: Direct estimation of the fiber orientation density function from diffusionweighted mri data using spherical deconvolution. NeuroImage 23(3) (2004) 1176–1185
 [7] Acqua, F.D., Rizzo, G., Scifo, P., Clarke, R.A., Scotti, G., Fazio, F.: A modelbased deconvolution approach to solve fiber crossing in diffusionweighted mr imaging. Biomedical Engineering, IEEE Transactions on 54(3) (2007) 462–472
 [8] Dell’Acqua, F., Scifo, P., Rizzo, G., Catani, M., Simmons, A., Scotti, G., Fazio, F.: A modified damped richardson–lucy algorithm to reduce isotropic background effects in spherical deconvolution. NeuroImage 49(2) (2010) 1446–1458
 [9] Descoteaux, M., Deriche, R., Knosche, T.R., Anwander, A.: Deterministic and probabilistic tractography based on complex fibre orientation distributions. IEEE Transactions on Medical Imaging 28(2) (2009) 269–286
 [10] Assaf, Y., Basser, P.J.: Composite hindered and restricted model of diffusion (charmed) {MR} imaging of the human brain. NeuroImage 27(1) (2005) 48 – 58
 [11] Alexander, D.C., Hubbard, P.L., Hall, M.G., Moore, E.A., Ptito, M., Parker, G.J., Dyrby, T.B.: Orientationally invariant indices of axon diameter and density from diffusion MRI. NeuroImage 52(4) (2010) 1374 – 1389
 [12] Zhang, H., Hubbard, P.L., Parker, G.J., Alexander, D.C.: Axon diameter mapping in the presence of orientation dispersion with diffusion {MRI}. NeuroImage 56(3) (2011) 1301 – 1315
 [13] Aboitiz, F., Scheibel, A.B., Fisher, R.S., Zaidel, E.: Fiber composition of the human corpus callosum. Brain research 598(1) (1992) 143–153
 [14] Liewald, D., Miller, R., Logothetis, N., Wagner, H.J., SchÃ¼z, A.: Distribution of axon diameters in cortical white matter: an electronmicroscopic study on three human brains and a macaque. Biological Cybernetics 108(5) (2014) 541–557
 [15] Nilsson, M., Alerstam, E., Wirestam, R., Sta, F., Brockstedt, S., Lätt, J., et al.: Evaluating the accuracy and precision of a twocompartment kärger model using monte carlo simulations. Journal of Magnetic Resonance 206(1) (2010) 59–67
 [16] Nilsson, M., van Westen, D., Ståhlberg, F., Sundgren, P.C., Lätt, J.: The role of tissue microstructure and water exchange in biophysical modelling of diffusion in white matter. Magnetic Resonance Materials in Physics, Biology and Medicine 26(4) (2013) 345–370
 [17] Daducci, A., CanalesRodríguez, E.J., Zhang, H., Dyrby, T.B., Alexander, D.C., Thiran, J.P.: Accelerated microstructure imaging via convex optimization (amico) from diffusion mri data. NeuroImage 105 (2015) 32–44
 [18] Tariq, M., Schneider, T., Alexander, D.C., WheelerKingshott, C.A.G., Zhang, H.: Bingham–noddi: Mapping anisotropic orientation dispersion of neurites using diffusion mri. NeuroImage 133 (2016) 207–223
 [19] Kaden, E., Knösche, T.R., Anwander, A.: Parametric spherical deconvolution: inferring anatomical connectivity using diffusion mr imaging. NeuroImage 37(2) (2007) 474–488
 [20] Sotiropoulos, S.N., Behrens, T.E., Jbabdi, S.: Ball and rackets: inferring fiber fanning from diffusionweighted mri. NeuroImage 60(2) (2012) 1412–1425
 [21] Anderson, A.W.: Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magnetic Resonance in Medicine 54(5) (2005) 1194–1206
 [22] Jespersen, S.N., Kroenke, C.D., Ãstergaard, L., Ackerman, J.J., Yablonskiy, D.A.: Modeling dendrite density from magnetic resonance diffusion measurements. NeuroImage 34(4) (2007) 1473 – 1486
 [23] Jespersen, S.N., Bjarkam, C.R., Nyengaard, J.R., Chakravarty, M.M., Hansen, B., Vosegaard, T., Ãstergaard, L., Yablonskiy, D., Nielsen, N.C., VestergaardPoulsen, P.: Neurite density from magnetic resonance diffusion measurements at ultrahigh field: Comparison with light microscopy and electron microscopy. NeuroImage 49(1) (2010) 205 – 216
 [24] Kaden, E., Kelm, N.D., Carson, R.P., Does, M.D., Alexander, D.C.: Multicompartment microscopic diffusion imaging. NeuroImage 139 (2016) 346 – 359
 [25] Kaden, E., Kruggel, F., Alexander, D.C.: Quantitative mapping of the peraxon diffusion coefficients in brain white matter. Magnetic Resonance in Medicine 75(4) (2016) 1752–1763
 [26] Reisert, M., Kellner, E., Dhital, B., Hennig, J., Kiselev, V.G.: Disentangling micro from mesostructure by diffusion mri: A bayesian approach. NeuroImage (2016) –
 [27] Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R.: Regularized, fast, and robust analytical qball imaging. Magnetic Resonance in Medicine 58(3) (2007) 497–510
 [28] Garyfallidis, E., Brett, M., Amirbekian, B., Nguyen, C., Yeh, F.C., Halchenko, Y., NimmoSmith, I.: Dipy–a novel software library for diffusion mr and tractography. In: 17th annual meeting of the organization for human brain mapping. (2011) 1–5
 [29] Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., Van Der Walt, S., Descoteaux, M., NimmoSmith, I., Contributors, D.: Dipy, a library for the analysis of diffusion mri data. Frontiers in neuroinformatics 8 (2014)
 [30] Zucchelli, M., Brusini, L., Méndez, C.A., Daducci, A., Granziera, C., Menegaz, G.: What lies beneath? diffusion eapbased study of brain tissue microstructure. Medical image analysis 32 (2016) 145–156
 [31] Grebenkov, D.S.: Analytical solution for restricted diffusion in circular and spherical layers under inhomogeneous magnetic fields. The Journal of Chemical Physics 128(13) (2008)
 [32] Özarslan, E., Shemesh, N., Basser, P.J.: A general framework to quantify the effect of restricted diffusion on the nmr signal with applications to double pulsed field gradient nmr experiments. The Journal of Chemical Physics 130(10) (2009) –
 [33] Kent, J.T.: The fisherbingham distribution on the sphere. Journal of the Royal Statistical Society. Series B (Methodological) 44(1) (1982) 71–80
 [34] Sotiropoulos, S.N., Jbabdi, S., Xu, J., Andersson, J.L., Moeller, S., Auerbach, E.J., Glasser, M.F., Hernandez, M., Sapiro, G., Jenkinson, M., et al.: Advances in diffusion MRI acquisition and processing in the human connectome project. NeuroImage 80 (2013) 125–143
 [35] Caruyer, E., Lenglet, C., Sapiro, G., Deriche, R.: Design of multishell sampling schemes with uniform coverage in diffusion mri. Magnetic resonance in medicine 69(6) (2013) 1534–1540
 [36] Jelescu, I.O., Veraart, J., Fieremans, E., Novikov, D.S.: Degeneracy in model parameter estimation for multicompartmental diffusion in neuronal tissue. NMR in Biomedicine 29(1) (2016) 33–47 NBM150204.R2.
 [37] Novikov, D.S., Veraart, J., Jelescu, I.O., Fieremans, E.: Mapping orientational and microstructural metrics of neuronal integrity with in vivo diffusion mri. arXiv preprint arXiv:1609.09144 (2016)
 [38] Fieremans, E., Burcaw, L.M., Lee, H.H., Lemberskiy, G., Veraart, J., Novikov, D.S.: In vivo observation and biophysical interpretation of timedependent diffusion in human white matter. NeuroImage 129 (2016) 414–427
 [39] De Santis, S., Jones, D.K., Roebroeck, A.: Including diffusion time dependence in the extraaxonal space improves in vivo estimates of axonal diameter and density in human white matter. NeuroImage 130 (2016) 91–103
 [40] Yolcu, C., Özarslan, E.: Diffusionweighted magnetic resonance signal for general gradient waveforms: Multiple correlation function framework, path integrals, and parallels between them. In: Visualization and Processing of Higher Order Descriptors for MultiValued Data. Springer (2015) 3–19
 [41] Zucchelli, M., Azfali, M., Yolcu, C., Westin, C.F., Menegaz, G., Özarslan, E.: The confinement tensor model improves characterization of diffusionweighted magnetic resonance data with varied timing parameters. (2016)
 [42] Andrews, G.E., Askey, R., Roy, R.: Special functions, volume 71 of encyclopedia of mathematics and its applications (1999)