NODDI-DTI: extracting neurite orientation and dispersion parameters from a diffusion tensor

NODDI-DTI is a simplification of the NODDI model that, when its underlying assumptions are met, allows extraction of biophysical parameters from standard single-shell DTI data, permitting biophysical analysis of DTI datasets. In contrast to the NODDI signal model, the NODDI-DTI signal model contains no CSF compartment, restricting its application to voxels with very low CSF partial-volume contamination. This simplification allowed derivation of simple analytical relations between parameters representing axon density, v and dispersion, τ, and DTI invariants (MD and FA) through moment expansion of the NODDI-DTI signal model. These relations formally allow extraction of biophysical parameters from DTI data. It should be emphasised that the NODDI-DTI model inherits the strong assumptions of the NODDI model, and so is nominally restricted to analysis of healthy brains. NODDI-DTI parameter estimates were computed by applying the proposed analytical relations to DTI parameters extracted from the first shell of data, and compared to parameters extracted by fitting the NODDI-DTI model to (i) both shells (recommended) and (ii) the first shell (as for DTI) of data in the white matter of three different in vivo diffusion datasets. NODDI-DTI parameters estimated from DTI and NODDI-DTI parameters estimated by fitting the model to the first shell of data gave similar errors compared to two-shell NODDI-DTI estimates. The NODDI-DTI method gave unphysical parameter estimates in a small percentage of voxels, reflecting voxelwise DTI estimation error or NODDI-DTI model invalidity. In the course of evaluating the NODDI-DTI model it was found that diffusional kurtosis strongly biased DTI-based MD values, and so a novel heuristic correction requiring only DTI data was derived and used to account for this bias. Our results demonstrate that NODDI-DTI is a promising model and technique to interpret restricted datasets acquired for DTI analysis with greater biophysical specificity, though its limitations must be borne in mind.


Introduction
The white matter (WM) of the human brain consists of dense bundles of neuronal axons connecting its functional areas. Neural circuits thus formed allow these areas to work together as a coherent entity. Changes in WM impact these neural circuits, and are thus the subject of studies investigating pathology [1,2,3], and cognition and learning [4,5].
Diffusion tensor imaging (DTI) [6] is, at present, the most commonly used method to observe WM changes in-vivo [4,5,7]. This is because DTI is simply implemented and time efficient while allowing robust extraction of complementary parameters (e.g. 'fractional anisotropy' (FA) and 'mean diffusivity' (MD) [8]) sensitive to microstructural WM changes [9], even in clinical contexts (see e.g. References [1,3]). Despite its microstructural sensitivity, the model underlying DTI (gaussian anisotropic diffusion [6]) is unspecific to biological changes. Numerous studies show MD and FA change in white matter (e.g. due to learning a new skill [4] or the pathology of Alzheimer's disease [2]), but cannot, in the absence of further information, distinguish e.g. changes in axon density from changes in axon arrangement.
In order to extract parameters of direct neurobiological relevance from diffusion MRI, we need biophysical models [10,11,9]. The majority of biophysical models (including the model introduced below) are 'multicompartment models'. Such models assume voxelwise diffusion contrast arises from linear combination of diffusion signals from distinguishable water compartments. Numerous multicompartment models have been proposed (see e.g. References [12,13,14,15,16,17,18,19,20,21]), but the complexity and lack of robustness of most of these models hinder their routine use in neuroscientific and clinical studies.
The NODDI (neurite orientation dispersion and density imaging) model [17] is a multicompartment model allowing robust and time-efficient extraction of maps of parameters representing neurite (in WM: axon) density and dispersion, and represents a trade-off between complexity, robustness, and acquisition-time duration. Robustness is achieved by fixing the values of several model parameters from earlier models [22,18], reducing the number of fitted parameters. As a result, the amount of data required to invert the model is reduced, giving acquisition-time durations approaching those available in clinical settings [17]. NODDI is thus gaining popularity in diffusion application studies [23,24,15,25,26,20,27], though the potential of fixed parameters to lead to bias in the fitted parameters has been a source of criticism [15,28,21,29].
Herein, we investigate a simplification of the NODDI model in which the CSF compartment (see below) is removed, and prove that there exist explicit relations between the remaining parameters of this model and MD and FA. We call this restricted version of the NODDI model 'NODDI-DTI' because it allows extraction of neurite orientation and dispersion parameters using MD and FA from DTI. The existence of these relations explains previously observed correlations between NODDI parameters and MD and FA [17,30,25,31,32]. Below, we investigate how accurately NODDI-DTI parameters extracted from DTI parameters in WM correspond to parameters fitted using the full NODDI-DTI model, and characterise limitations of the model and method.

Theory
The NODDI signal model supposes three compartments: intraneurite water, extraneurite water, and free water [17]. The biophysical parameters fitted in the model are neurite density (volume fraction of the intraneurite compartment), ν; a measure of neurite dispersion, κ; a vector giving the main neurite orientation; and a volume fraction accounting for partial-volume effects with free water (nominally CSF) [33,34,17]. An important fixed parameter is the intrinsic diffusivity of the intraneurite compartment, d = 1.7 × 10 −3 mm 2 s −1 [17]. The primary neurite orientation [17] is formally equivalent to the principal eigenvector of the diffusion tensor (DT) (see Appendix A), as previously observed empirically [35].
The model which we investigate here, which we call 'NODDI-DTI', is a reduced form of the NODDI model with no CSF volume fraction; this model has been previously observed to give reasonable estimates of ν and κ from single-shell data [36]. For ease of computation in the following we use τ instead of κ as our measure of dispersion, where [15,17,18] and erfi is the imaginary error function. τ ranges from 1/3 (isotropically distributed neurites) to 1 (perfectly aligned neurites)-increasing τ corresponds to increasing neurite alignment-and is the average of cos 2 (ψ) over the neurite distribution, where ψ is the angle between a given neurite and the main neurite orientation [15]. By expanding the NODDI-DTI signal model in moments, one can derive a corresponding DT [18]. As shown in Appendix A, appropriate combination of the eigenvalues of this DT allows expression of ν and τ in terms of MD and FA of this DT: These relations are exemplified in Figure 1, and we note that Equation (2) has been previously derived independently [37]. Eqs. (2) and (3) demonstrate a one-to-one mapping from (MD,FA) to (ν, τ ), demonstrating that, formally, NODDI-DTI parameters can be extracted from DTI data. We can predict domains within which MD and FA should lie if the NODDI-DTI model provides a valid representation: substituting ν ∈ [0, 1] and τ ∈ [1/3, 1] into Eqs. (2) and (3) gives the domains Exceptions will arise in practice either from DT measurement error or the invalidity of the NODDI-DTI model as a representation in a given voxel. Estimation of a quantitative DT (i.e. the first moment of the diffusion signal) is nontrivial, as immediately apparent from examination of Figure 1: the main bulk of the MD values lies below the prediction of Eq. (2). A large part of this bias is due to neglecting higher order moments in estimating MD [38]. In order to reduce this bias, we define the heuristically corrected MD, where λ i is the ith eigenvalue of the measured DT and δ ij is the Kronecker delta. Eq. (5), derived in Appendix B, pragmatically assumes that: only the first higher moment, diffusional kurtosis [39], contributes; the square of the apparent diffusion coefficient is uncorrelated with the apparent diffusional kurtosis; the mean diffusional kurtosis can be taken to be unity (approximately true over much healthy human brain WM [39,40,41]); and the effect of diffusional kurtosis on each individual eigenvalue is negligible. Figure 1 shows much closer agreement between the prediction of Eq. (2) and the experimental data when using MD h ; further justification of this correction is given below. Substituting Eq. (5) into Eq. (2) gives the relation used in the following to compute ν from DT invariants: The effect of failing to correct for diffusional kurtosis is much less pronounced for FA [38], and preliminary experiments (data not shown) showed that applying diffusional kurtosis correction to only MD in Eq. (3) resulted in a modest increase in the number of unphysical τ parameter estimates. This latter observation can be explained using Eq. (4): whenever heuristic diffusional kurtosis correction leads to overestimation of MD, the upper bound for allowed FA values is artificially decreased, potentially leading to unphysical τ estimates. We therefore apply no correction to Eq. (3).

Data collection and preprocessing
All data were collected by scanning healthy volunteers in a MAGNETOM Tim Trio 3 T MRI system (Siemens AG, Healthcare Sector, Erlangen, Germany) as part of a study approved by a local ethics committee. In each case informed written consent was obtained prior to scanning.
In all cases (subjects 1-3) subject motion, eddy currents, and susceptibility distortions were corrected for using the ACID toolbox [45]; for details see References [46,40,47,48]. Additionally for the multiband data (subjects 1 and 2), after the above corrections were applied, the corrected data from the two phase encoding directions were summed for use in subsequent analysis. This final step was unnecessary for subject 3.

Parameter estimation and comparison
Parameters were estimated only in WM voxels determined to be largely unaffected by CSF or grey matter partial volume effects. This determination was made by thresholding at 50 % probability a WM probability map obtained by segmenting the first b = 0 image of each respective dataset in SPM 12 [49].
The ACID toolbox was used to compute FA, MD, and the eigenvalues of the DT from the low b-value shell of each dataset, and home-written SPM scripts were then used to generate ν and τ using Eqs. (6) and (3), respectively, from these DT parameters; we refer to these as the 'NODDI-DTI method' results. For subjects 1 and 2, the ACID toolbox was also used to simultaneously estimate the diffusion and kurtosis tensors [40], giving silver standard mean diffusivity estimates (MD DKI ) less biased by the effects of diffusional kurtosis [38] allowing evaluation of the validity of Eq. 5. MD DKI was not computed for subject 3 as the second b-value was deemed too large to allow accurate estimation of the kurtosis tensor [39].
NODDI-DTI silver standard results were obtained by fitting both shells of data using the NODDI toolbox [17,50] with CSF volume fraction fixed at zero, followed by conversion of κ into τ using a home-written SPM script implementing Eq. (1). We refer to these results as the 'two-shell NODDI-toolbox fitted' results, and they represent a silverstandard because two shells of data are sufficient to make inversion of multi-compartment signal models well-posed [51].
In order to investigate the magnitude of the differences between NODDI-DTI methodand the NODDI-toolbox fitted-results, fits were also made using the NODDI toolbox of the subset of the diffusion data used for the DTI fitting. The designation 'one-shell NODDI-toolbox fitted' distinguishes these results from the NODDI-DTI method-and two-shell NODDI-toolbox fitted-results. The NODDI toolbox has been shown previously to give reasonable results fitting single-shell data when CSF compartment fraction is fixed at zero [36].

Results
Before examining the results of NODDI-DTI parameter extraction, we examine the heuristic correction of MD for diffusional kurtosis (Eq. (5)) used to calculate ν. Figure 1 and the BA plots in Figure 2 show that the heuristically corrected MD, MD h , is less biased compared to the uncorrected MD, justifying use of the correction. Numerical values (mean ± one standard deviation) of the differences in Figure 2 are for subject 1: 0.117 ± 0.040 (MD), 0.025 ± 0.046 (MD h ); and for subject 2: 0.119 ± 0.041 (MD), 0.029 ± 0.041 (MD h ).
The similarity of parameters estimated using the NODDI-DTI method to the silver standard results can be seen in Figure 3, which shows parameter maps computed with each method, along with maps showing the differences between the parameter estimates. Differences between the parameter estimates are further presented in several complementary ways: BA plots in Figure 4 show general behaviour, plots of the means and standard deviations of the differences in Figure 5 compare this general behaviour across subjects, and the series of slices in Figures 6 and 7 show the behaviour of NODDI-DTI method parameter estimates throughout the brain. Subjects 1 and 2, measured using the same protocol, show the best agreement, but all three subjects behave similarly, demonstrating robustness of the NODDI-DTI method.
The magnitudes of the means and standard deviations of the differences between oneand two-shell NODDI-toolbox fitting are shown in Figure 5. One-shell NODDI-toolbox fitting gave stable fits in this case because the CSF compartment fraction was fixed at zero [36]. The NODDI-DTI method showed smaller mean differences and one-shell NODDI-toolbox fitting showed smaller standard deviations of the differences. Overall, however, both methods of parameter estimation were comparable, further demonstrating the validity of the NODDI-DTI method.

Discussion
This work has demonstrated that, with caveats to be discussed below, parameters of potential neurobiological relevance can be extracted from DTI parameters using the NODDI-DTI relations, Equations (3) and (6). The improved interpretability thus gained would not only be applicable to future DTI studies, but also to existing DTI studies.
As an example of using NODDI-DTI to reinterpret existing DTI studies, we apply the method superficially to the study of Scholz, et al. [4], which demonstrated a statistically significant FA increase in WM "underlying the intraparietal sulcus" after participants learned to juggle [4]. Assuming no concomitant change in MD (as suggested by the authors not reporting any significant change), this FA increase could be interpreted, using Eq. (3), as an increase in τ , i.e. an increase in alignment of neuronal axons in this area with training. This result is much more specific than a change in FA, though we note that quantitative analysis would require reanalysis of the original data. A follow up study recording multi-shell data, combined with proper mechanistic analysis of the WM plasticity mechanisms, would allow investigation of this effect in more detail.
Unfortunately, the NODDI-DTI method did not always give physically plausible parameter estimates for the datasets studied herein. We posit four overlapping explanations for these unphysical estimates, related to NODDI-DTI assumptions.
Assumption 1: DTI parameters can be accurately estimated from the diffusion signal. Errors in DTI parameter estimation will lead to errors in parameters estimated using Eqs. (3) and (6), potentially giving unphysical parameter estimates. The b-values used likely resulted in overestimates of FA in regions of high anisotropy due to poor estimates of the low DT eigenvalues [8,53], explaining why many of the unphysical τ estimates were found in the highly anisotropic corpus callosum (Figure 7). This hypothesis is bolstered by the results of subject 3: a lower b-value increased the proportion of unphysical estimates, especially in highly anisotropic regions.
Assumption 2: CSF can be ignored in voxels with a high probability of being WM. NODDI-DTI could give unphysical parameter estimates whenever a voxel contains a significant amount of CSF: the high diffusivity of CSF [17] can take MD h outside the limits of NODDI-DTI (Eqs. (4)). Figures 6 and 7 show that many of the voxels where NODDI-DTI gave unphysical parameter estimates are close to the edge of the WM mask, in line with partial volume effects being important. The larger voxel size used for scanning subject 3 means that partial volume effects are more prominent, partially explaining the poorer ν performance in this case. Because CSF volume fraction was fixed at zero in our NODDI-toolbox fits, residual partial volume effects may also have affected these parameter estimates.
Assumption 3: MD can be heuristically corrected for diffusional kurtosis bias. Figure 1 shows that diffusional kurtosis affects our estimates of MD; such effects could take MD estimates out of the range of applicability for NODDI-DTI. While our heuristic correction (Eq. (5)) substantially mitigates this issue, it does not eliminate it. This is evident in the ν Bland-Altman plots (Figure 4), where the mean and standard deviation of the differences visibly vary with the mean ν estimate, implying (via Eq. (6)) residual correlation between the errors and the corrected MD.
Residual bias in the region around the corpus callosum of the ν difference maps in Figure 6 can be explained by incomplete diffusional kurtosis correction. Here the mean diffusional kurtosis is greater than unity [39], and the low DT eigenvalues are poorly estimated (see Assumption 1), meaning that the heuristic correction is not entirely valid.
Assumption 4: the NODDI-DTI model is a valid representation of the diffusion signal. Assumptions underlying the NODDI-DTI signal model can be considered overly restrictive: the model cannot formally represent WM voxels containing perpendicularly crossing fibre bundles [54], and it is a subject of controversy as to whether intra-and extraneurite intrinsic diffusivities can be taken to be equal [28]. With respect to the latter point, NODDI-DTI will be unrepresentative whenever a voxel contains a significant amount of iron (e.g. through partial voluming with iron rich grey matter nuclei): susceptibility effects lower diffusivities artificially [55], and so lower d, assumed to be a fixed value in NODDI-DTI. In voxels where the NODDI-DTI signal model is unrepresentative, the NODDI-DTI method may give unphysical parameter estimates. Such failures will not be immediately apparent in the NODDI-toolbox fitted parameters because constraints in the fitting procedure mean parameters outside the physical range can never be returned, regardless of the model's biological implausibility in a given voxel.
The greater number of unphysical NODDI-DTI method parameter estimates for ν as compared to τ can be explained by ν estimation being more sensitive to partial volume and diffusional kurtosis effects. This is borne out by the locations of the failures ( Figure 6): mainly either close to the edge of the WM mask (implying partial volume effects), or in regions of high anisotropy (implying residual diffusional kurtosis effects).
Pathology could further undermine the assumptions underlying NODDI-DTI: pathological processes can lead to free water located far from CSF compartments [56], can affect mean kurtosis values [57], and could affect the 'true' value of d [15].
NODDI-DTI could be improved and made more appropriate for clinical studies through investigation of the following points. Unphysical parameter estimates could be pragmatically eliminated by constraining DT fitting using Eqs. (4) (appropriately corrected using Eq. (5)). Estimates of CSF volume fraction could be incorporated into NODDI-DTI without requiring extra data acquisition using the free water elimination method [56,34,58]. Known values of mean diffusional kurtosis in the brain [39,40,41] could be used to construct mean diffusional kurtosis Bayesian priors [59,51], or diffusional kurtosis corrected MD and FA could be measured directly using time-efficient methods [60,61,62].
We finish by providing practical recommendations for a minimal NODDI-DTI acquisi-tion scheme. Results at b = 1000 mm 2 s −1 were reasonable, and we would recommend this as a lower b-value bound; b = 800 mm 2 s −1 unfortunately gave many unphysical parameter estimates in the corpus callosum. An upper bound on b-value comes from ensuring diffusional kurtosis does not constitute the majority of the diffusion contrast. Eq. (B.3) shows that (assuming the apparent diffusional kurtosis is approximately unity [39,40,41]) choosing b 6/d ≈ 3500 mm 2 s −1 means that the DT dominates diffusion contrast, giving an upper b-value bound. High resolution acquisitions which maintain good signal to noise ratio [53] are recommended to reduce partial volume effects. Accurate DT estimation requires measurement of at least 30 distinct diffusion directions [6]; we recommend at least this number for application of NODDI-DTI, although the lowest number of orientations tested was 60.

Conclusions
We have estimated biophysical parameters representing neurite density and dispersion with reasonable accuracy from diffusion tensor parameters extracted from single-shell diffusion data. Heuristic kurtosis correction of MD was necessary to remove diffusional kurtosis bias; use of corrections such as that derived here could improve other analyses of single-shell diffusion data requiring quantitative MD estimates.
NODDI-DTI has the potential to open up two new opportunities: (a) more specific neurobiological interpretation of observed microstructural changes in DTI data (including interpretation of existing datasets), and (b) simple and time efficient estimation of biophysical parameters from smaller diffusion datasets, despite limitations due to the underlying model and difficulties estimating accurate diffusion tensors.

Appendix A.1. Diffusion tensor of the NODDI-DTI signal model
We begin derivation of the NODDI-DTI relations by deriving the DT arising from the NODDI-DTI signal model. This derivation is similar to that of Reference [18], where the DT of a precursor to the NODDI-DTI model [22] was derived.
The normalised signal arising from the NODDI signal model can be written [17] S = ν p(κ, µ, n) exp{−bd q t n n t q} d n where the first term represents the intraneurite water compartment with diffusivity d parallel to the neurite and zero perpendicular to it; the second term represents the extraneurite water compartment; arrows denote normalised vectors; · t denotes transposition; q is the diffusion gradient vector; ν represents neurite density; and D ec = d p(κ, µ, n) n n t + (1 − ν)(1 − n n t ) d n, the DT of the extraneurite compartment. The form of the extraneurite DT arises from assuming that: the diffusivity of the extraneurite space in the absence of neurites is equal to the intraneurite diffusivity along the direction of the neurite [17], the neurites reduce the diffusivity in a long-time-limit tortuous manner [17], and extracellular water is in fast exchange among all neurite orientations [21]. The probability density is a Watson distribution giving the distribution of neurites about main orientation µ with dispersion parameter κ [17]. Isotropically distributed neurites correspond to κ = 0, neurites perfectly aligned along µ correspond to κ → ∞. = (1 − ν) 2 d + (2 − ν)νd p(κ, µ, n) n n t d n . (A.8) The integral appearing on the right-hand side of Eq. (A.8) is given by [18] p(κ, µ, n) n n t d n = τ µ µ t + where τ is defined in Eq. (1). Inserting Eq. (A.9) into Eq. (A.8) gives from which, by inspection, the largest eigenvalue (corresponding to an eigenvector colinear with the main neurite orientation) is 11) and the other two eigenvalues are degenerate (with respective eigenvectors arbitrarily defined in the plane perpendicular to µ): Because λ 1 ≥ λ 2 , when the primary eigenvector of D NODDI−DTI is well-defined (i.e. when D NODDI−DTI is not isotropic), this eigenvector is formally equivalent to the main neurite orientation, as observed empirically [35]. (A.14) Solving this quadratic equation for ν one obtains where the sign ambiguity is resolved by recalling ν ≤ 1 to give Eq. (2).
which can be rearranged into the quadratic equation for which the solutions are: Eq. (A. 19), reveals that all we must do to express τ in terms of MD and FA is (i) express τ in terms of MD and λ 2 and then (ii) substitute Eq. (A.19) into the resulting expression. Part (i) is achieved by substituting Eq. (2) into Eq. (A.12), then simplifying to give which, after rearranging for τ , reveals We can now perform part (ii): inserting Eq. (A.19) into Eq. (A.21) and simplifying gives the result The sign ambiguity is resolved by recalling that τ ≥ 1/3 [15] and that both MD and FA are nonnegative, resulting in Eq. (3). Eq. (3) is ill-defined at MD = d (the denominator of the second term goes to zero); we classify values at this point as 'unphysical' unless FA is also zero. This latter situation corresponds to the complete absence of fibres (as confirmed by inserting MD = d into Eq. (2)), and so τ is taken to equal its isotropic value, 1/3.

Appendix B. Heuristic correction of MD for diffusional kurtosis (Eq. (5))
When diffusional kurtosis and higher order moments are zero, the normalised diffusion signal S is related to the apparent diffusivity, D app , by [64] log(S) = −bD app , (B.1) giving where · denotes averaging over all diffusion directions. The complicated microstructure of white matter requires higher order moments to represent the diffusion signal [63,39,38]. To the order of the diffusional kurtosis the normalised diffusion signal is where K app is the apparent diffusional kurtosis [63]. The effective mean diffusivity MD eff derived from this signal (as per Eq. (B.2)) would be which differs from the true MD by a term which we call 'diffusional kurtosis bias'. While good estimates of unbiased MD can be obtained from multi-b-value data [38], such extra data is not available for most DTI acquisitions, and so we derive and use an heuristic correction to mitigate diffusional kurtosis bias.
Defining the covariance of D app 2 and K app : where MK = K app , the mean kurtosis, we can write Eq. (B.4) in the form We pragmatically assert that cov (D app 2 , K app ) = 0, i.e. assume that the apparent diffusivity and apparent diffusional kurtosis are uncorrelated. This assertion results in which is further simplified by assuming MK = 1 (true in much healthy WM [39,40,41]): To compute the average in Eq. (B.9), we express D app in components of the DT, D, and orientation vector, q, i.e. [39] As we integrate over all q on the sphere, we can freely choose the basis of q. We thus choose the diagonal basis of D, simplifying Eq. (B.10) to: where λ i is the ith eigenvalue of D and is independent of orientation. Averages over the products of the components q i , evaluate to q 2 i q 2 k = (1 + 2δ ik )/15, where δ ik is the Kronecker delta, thus Inserting Eq. (B.12) into Eq. (B.9) and rearranging gives the heuristically corrected MD: which becomes independent of diffusional kurtosis upon assuming the measured (diffusional kurtosis biased) eigenvalues can be substituted for the 'true' eigenvalues, giving Eq. (5). : Absolute value maps of the difference between τ computed using the NODDI-DTI method and fitting two-shells with the NODDI toolbox. Data from all three subjects are shown, and slice numbers are given for each row (slice) and column (subject). The extent of the colour scale at the top right shows the windowing for all slices. Blue denotes voxels where the NODDI-DTI method gave an unphysical parameter estimate. The large number of unphysical τ values in the corpus callosum of subject 3 is discussed in the text.