NODDI-SH: a computational efficient NODDI extension for fODF estimation in diffusion MRI

NODDI-SH: a computational efficient NODDI extension for fODF estimation in diffusion MRI

Mauro Zucchelli Department of Computer Science, University of Verona, Italy    Maxime Descoteaux Sherbrooke Connectivity Imaging Lab, Computer Science, Université de Sherbrooke, Canada    Gloria Menegaz Department of Computer Science, University of Verona, Italy

Diffusion Magnetic Resonance Imaging (DMRI) is the only non-invasive 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 NODDI-SH, is derived by convolving the single fiber response with the fODF in each voxel. NODDI-SH 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 in-vivo 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 NODDI-SH 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 NODDI-SH 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 ill-conditioned problem which can be affected by noise or error in the estimation of the response function. Constrained Spherical Deconvolution (CSD) [1] added non-negative 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 multi-compartment 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 in-between 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 non-linear 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 two-compartmental 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 b-value 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 NODDI-like multi-compartment model, gathering the advantages of both techniques in a unified model called NODDI-SH. NODDI-SH 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, NODDI-SH 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, NODDI-SH elegant mathematical framework enables an extremely fast computation of the volume fractions and fast fODF reconstruction. NODDI-SH 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 NODDI-SH 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 NODDI-SH model, as well as the simulated and in-vivo data used for the experiments. Section 3 shows the comparison of NODDI-SH with NODDI in the simulated and in-vivo experiments, as well as the robustness of the NODDI-SH 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


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


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


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 (NODDI-SH).

2.1 FORECAST model

FORECAST considers to be an axially symmetric tensor, with parallel diffusivity and perpendicular diffusivity


where is the b-value at the given diffusion time , and is the diffusion tensor aligned along v


Replacing equation (5) in equation (4), we can calculate as


Using the Funk-Hecke Theorem [27] it is possible to solve the integral in equation (7) (see Appendix B for the complete derivation) obtaining the final FORECAST model


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 b-value, which is given by the element of the FORECAST basis


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 b-value), it is possible to use equation (9) to estimate and accurately.

2.2 The NODDI-SH model

The NODDI-SH model considers the single fiber diffusion signal as the weighted sum of three tissue compartments: intracellular , extracellular , and cerebrospinal fluid


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


Replacing equation (10) into equation (4) (see Appendix C) leads to the final NODDI-SH model formulation


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 b-value


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 non-uniform 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 b-value. In order to select the best combination, we evaluated the mean of the diffusion signal at each b-value 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 NODDI-SH is performed via convex optimization using CVXOPT111 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


subject to

Yc (15)

where is the measured diffusion signal, c is the coefficients vector, M the NODDI-SH 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 library222 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 NODDI-SH 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 toolbox333 In NODDI model [2] intracellular volume fraction is weighted by a factor in order to be compared to NODDI-SH [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


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


with the Gamma function and the modified Bessel function. Figure 1 shows some examples of the Kent distribution444Kent distribution has only one lobe. We added the antipodal symmetric lobe for visualization purposes. , given the parameters and .

, , ,
, , ,
Figure 1: Instances of Kent distribution 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


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

Figure 2: Instances of crossings generated with two Kent distributions with and .

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 In-vivo dataset

In order to test the NODDI-SH model in-vivo we considered one subject of the HCP diffusion MRI dataset. HCP acquisition scheme is composed of 18 volumes sampled at b-value 0 s/mm, plus 90 volumes at b-value 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 q-samples 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 b-value 1000 and 2000 s/mm are provided to the model. The results of the NODDI-SH fitting were evaluated by calculating the mean square error (MSE) between the original HCP signal and the NODDI-SH signal for all the samples, also those which were not used for the fitting.

3 Results

3.1 Results of simulated data

Figure 3: Estimation of the intracellular volume fraction of the simulated data for NODDI-SH (top row) and NODDI (bottom row) in the case of isotropic fanning () and anisotropic fanning.

Figure 3 (left column) shows the estimation of the intracellular volume fraction for NODDI-SH and NODDI in the case of isotropic dispersion (), compared to the ground truth (GT) values (black circles). NODDI-SH is able to provide the most accurate estimation of in case of highly dispersed fibers (), while NODDI in case of . For NODDI-SH, the estimation error increases while decreasing the dispersion ( and ). Both models tend to slightly overestimate in all conditions. The average errors for NODDI-SH 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 NODDI-SH, 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) NODDI-SH 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. NODDI-SH, 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 NODDI-SH 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 NODDI-SH, respectively.

Figure 4: Estimation of the intracellular volume fraction of the simulated data for NODDI-SH (left) and NODDI (right) in the case of crossing fibers.

Figure 4 shows the estimation of for the crossing fibers dataset. In this case, the estimation errors of both, NODDI-SH 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, NODDI-SH 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 NODDI-SH, 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 NODDI-SH, respectively.

       CSD        FORECAST
       NODDI-SH        AVERAGE
Figure 5: Angular error (AE) corresponding to the fODF calculated with CSD, FORECAST, and NODDI-SH. The average represents the mean of the AE at 45, 60, and 90 degrees for each models.

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 NODDI-SH 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 NODDI-SH, 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 NODDI-SH 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 high-resolution fODFs. Both NODDI-SH 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

NODDI-SH three compartments model is designed to capture the heterogeneity of the brain tissues. However, there are no guarantees that NODDI-SH estimation of neurites density and volume fraction corresponds to the underlying tissue composition, as it is the case for all models. Therefore, results on in-vivo 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 NODDI-SH in in-vivo 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 NODDI-SH, respectively. The maps obtained with the two techniques are very similar. NODDI-SH 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. NODDI-SH 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 non-zero, 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.




Figure 6: Intracellular volume fraction (left) and CSF volume fraction (right) estimation of a coronal slice of one HCP brain calculated using NODDI (top row) and NODDI-SH (bottom row).

Figure 7 shows the MSE calculated between the normalized diffusion signal and the estimated diffusion signal of NODDI and NODDI-SH. 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. NODDI-SH MSE presents a more uniform pattern in the white matter, gray matter, and CSF. The highest MSE values for NODDI-SH were found in the corpus callosum and in the basal nuclei.

Figure 7: NODDI MSE (left) and NODDI-SH MSE (right) of a coronal slice of one HCP brain.

Figure 8 (left) highlights the regions of the brain in which NODDI-SH features a lower MSE than NODDI. Almost all the white matter appears to be white, which seems to indicate that NODDI-SH 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 NODDI-SH 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 NODDI-SH, although the absolute value of the difference of the MSE is lower than that obtained in white matter (Figure 8, right).

Figure 8: Regions of the brain in which NODDI-SH obtain a lower MSE with respect to NODDI (left), and absolute value of the difference of the MSE (right).
Figure 9: CSD, FORECAST, and NODDI-SH fODFs, calculated on region of interest taken from a coronal slice of the HCP data.

Figure 9 shows the fODFs obtained by CSD, FORECAST, and NODDI-SH 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 NODDI-SH 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) NODDI-SH 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 multi-compartmental 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 b-value = 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 NODDI-SH to fit more practical acquisition schemes, we subsampled the HCP data. Figure 10 shows the mean signal value of the HCP data for each b-value, while keeping all the gradients (90 directions), 60 directions, and only 30 directions, respectively. As expected, the diffusion signal decays for increasing b-values. In regions of fast diffusion, such as the CSF, the diffusion signal is almost zero even at b-value = 1000 s/mm. At b-value = 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 b-value = 3000 s/mm.

90 directions 60 directions 30 directions

Figure 10: Mean value of the diffusion signal for each b-values. The mean diffusion signal is calculated using the full dataset (90 directions) or keeping only 60 and 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 b-value = 1000 s/mm derived from the 90 and the 30 samples cases is only of . The mean signal over each shell is necessary for NODDI-SH 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 NODDI-SH when changing the number of directions per shell and the number of shells (b-value = [1000, 2000, 3000] s/mm or b-value = [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 b-value = 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



Figure 11: Intracellular volume fraction calculated using NODDI-SH considering only a subsample of the HCP gradients for the fitting.

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 NODDI-SH 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 NODDI-SH 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 NODDI-SH 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



Figure 12: MSEs calculated using NODDI-SH considering only a subsample of the HCP gradients for the fitting.

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.

Figure 13: Normalized histogram of the MSE for the different sampling schemes calculated only in voxel with .

3.4 Computational time

NODDI-SH volume fractions estimation and fODF are substantially similar to NODDI and CSD results, respectively. However, to obtain the same information provided by NODDI-SH it is necessary to run both NODDI and CSD on the same data. Table 1 shows the computational time for CSD, FORECAST, NODDI, and NODDI-SH 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™i7-6500 CPU @2.50GH. DIPY CSD implementation is extremely fast, followed by NODDI-SH 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 NODDI-SH 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 NODDI-SH to fit 2970 voxels. This makes NODDI-SH even faster than NODDI implemented using AMICO. A complete comparison of the performance of the NODDI-SH 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
NODDI-SH 18.63 0.007
NODDI-SH (no fODF) 1.30 0.0004
Table 1: Computational times of CSD, FORECAST, NODDI, and NODDI-SH.

4 Discussion

NODDI-SH embeds both microstructural parameters estimation and fODF reconstruction in the same mathematical framework. We demonstrated that NODDI-SH 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 NODDI-SH. 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 NODDI-SH is that using the SMT, the estimation of the volume fractions depends only on the mean value of the signal. Using this feature, the NODDI-SH parameters fitting is not only fast, but also requires few sampling points (60 directions spread on two shells) in order to obtain good results. NODDI-SH 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 b-values, 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 NODDI-SH 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 over-simplification 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, NODDI-SH shares the same formulation of NODDI model and, therefore, some of its limitations. More accurate results would be obtained by making voxel-wise 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 NODDI-SH, [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 NODDI-SH potentially improving the resulting fODF and MSE.

In this paper, we considered only the HCP acquisition scheme and sub-sampling 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 b-values, it could not be the best choice for larger axons (e.g. 4 ) at higher b-values (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 NODDI-SH these diffusion time dependent model could be easily integrated into its framework.

5 Conclusions

NODDI-SH is a flexible model that allows obtaining both tissues volume fraction parameters and the fODF at the same time. Using NODDI-SH, 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, intra-cellular and extra-cellular 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 three-compartments model. The current NODDI-SH 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 NODDI-SH. Results show that NODDI-SH can also be used for dataset presenting a limited number of samples, which makes it feasible for processing clinical data.

Future works will include the extension of the NODDI-SH to other models [40], [41] and the estimation of the fiber dispersion from the fODF.

6 Acknowledgements

Data were provided by the Human Connectome Project, WU-Minn 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


where is the general SH basis, written as


with the polar representation of u, and the associated Legendre Polynomial. A possible way to represent is to consider its SH series expansion


Appendix B Analytical derivation of the FORECAST model

In the FORECAST model the diffusion signal can be calculated as


with . The tensor can calculated as


assuming that the two perpendicular eigenvalues are equal. We can now calculate as


Using the Funk-Hecke Theorem [42], [27] stating that, given a unit vector u, an integral on the sphere can be converted to a one dimensional integral


assuming continuous on and any SH of order , and with the Legendre polynomial of degree .

Replacing it into FORECAST equation the problem became:


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


When is indeterminate but the following limits hold true


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


where N [0,8].

b.0.1 Special case b=0

In the case, the signal equation can be rewritten as


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):


this because


where is equal to one when and and zero elsewhere. This equation is valid only for .

Appendix C Analytical derivation of the NODDI-SH

The three-compartments NODDI-SH model can be derived equation (10) into equation (4):


solving each integral as a single FORECAST basis (see Appendix B). The final formulation of the NODDI-SH model can be rearranged as



  • [1] Tournier, J.D., Calamante, F., Connelly, A.: Robust determination of the fibre orientation distribution in diffusion mri: non-negativity constrained super-resolved spherical deconvolution. NeuroImage 35(4) (2007) 1459–1472
  • [2] Zhang, H., T., S., Wheeler-Kingshott, 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 self-diffusion 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 diffusion-weighted 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 model-based deconvolution approach to solve fiber crossing in diffusion-weighted 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 electron-microscopic 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 two-compartment 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., Canales-Rodrí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., Wheeler-Kingshott, 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 diffusion-weighted 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., Vestergaard-Poulsen, 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.: Multi-compartment microscopic diffusion imaging. NeuroImage 139 (2016) 346 – 359
  • [25] Kaden, E., Kruggel, F., Alexander, D.C.: Quantitative mapping of the per-axon 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 q-ball imaging. Magnetic Resonance in Medicine 58(3) (2007) 497–510
  • [28] Garyfallidis, E., Brett, M., Amirbekian, B., Nguyen, C., Yeh, F.C., Halchenko, Y., Nimmo-Smith, 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., Nimmo-Smith, 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 eap-based 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 fisher-bingham 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 multi-compartmental diffusion in neuronal tissue. NMR in Biomedicine 29(1) (2016) 33–47 NBM-15-0204.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 time-dependent 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 extra-axonal space improves in vivo estimates of axonal diameter and density in human white matter. NeuroImage 130 (2016) 91–103
  • [40] Yolcu, C., Özarslan, E.: Diffusion-weighted 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 Multi-Valued 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 diffusion-weighted 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)
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description