Estimating axial diffusivity in the NODDI model

To estimate microstructure-related parameters from diffusion MRI data, biophysical models make strong, simplifying assumptions about the underlying tissue. The extent to which many of these assumptions are valid remains an open research question. This study was inspired by the disparity between the estimated intra-axonal axial diffusivity from literature and that typically assumed by the Neurite Orientation Dispersion and Density Imaging (NODDI) model (d∥ = 1.7 μm2/ms). We first demonstrate how changing the assumed axial diffusivity results in considerably different NODDI parameter estimates. Second, we illustrate the ability to estimate axial diffusivity as a free parameter of the model using high b-value data and an adapted NODDI framework. Using both simulated and in vivo data we investigate the impact of fitting to either real-valued or magnitude data, with Gaussian and Rician noise characteristics respectively, and what happens if we get the noise assumptions wrong in this high b-value and thus low SNR regime. Our results from real-valued human data estimate intra-axonal axial diffusivities of ~ 2 – 2.5 μm2/ms, in line with current literature. Crucially, our results demonstrate the importance of accounting for both a rectified noise floor and/or a signal offset to avoid biased parameter estimates when dealing with low SNR data.


Introduction
In diffusion MRI, biophysical models aim to relate macroscopic diffusion signals to microscopic, biologically meaningful tissue parameters such as fibre orientation, dispersion or diameter. The primary challenge for biophysical models is being able to describe the complex tissue microstructure in only a handful of parameters that can be estimated reliably from the diffusion signal. Consequently, the model must make strong, simplifying assumptions about the underlying architecture and/or diffusion properties of the tissue. Some parameters are often constrained or set to a single value which is either taken from the literature or decided by some other means (e.g. by fitting the model multiple times with different values for some assumed parameter, and taking the result which maximises the quality of the fit across voxels or specimens (Grussu et al., 2017 ;Guerreroid et al., 2019). Though these assumptions make the estimation of otherwise degenerate parameters tractable, if the modelling constraints are inaccurate, the estimation of the remaining model parameters will be biased (Jelescu et al., 2016;Howard et al., 2019).
Neurite Orientation Dispersion and Density Imaging (NODDI) (Zhang et al., 2012) is a commonly used biophysical model in diffusion MRI which, due to its clinically feasible scan times requirements, has been reported in numerous studies of patient populations (Adluru et al., 2014;Timmers et al., 2016;Schneider et al., 2017;Hagiwara et al., 2019;Taoka et al., 2020). NODDI is a variant of the white matter 'standard model' (Novikov et al., 2018b) in which specific assumptions about the tissue diffusivity allow for voxelwise estimation of fibre dispersion and neurite 'density' or signal fraction. One assumption is that the intra-axonal axial diffusivity, i.e. the diffusion of water molecules inside the axon as they travel along the primary axis or orientation, is a fixed, global value known a priori and typically set to 1.7 μm 2 /ms. An axial diffusivity of 1.7 μm 2 /ms has been previously justified in adult white matter by minimising the model residuals when fitting to in vivo data from multiple subjects (Guerreroid et al., 2019). However, optimising over the residuals may be sub-optimal for a degenerate model like NODDI, where multiple parameter sets can produce the same diffusion signal. This degeneracy likely explains the broad range of reported axial diffusivities (~ 1 -2.5 μm 2 /ms) with similar residual error. Though these axial diffusivities may produce a similar residual error, the estimated NODDI parameters from each fit will be different. Furthermore, there now exists a body of work in which the intra-axonal axial diffusivity is generally estimated to be higher, typically in the range of ~ 2 -2.5 μm 2 /ms (Dhital et al., 2019;Jelescu and Budde, 2017;Jelescu et al., 2020;Kaden et al., 2016;Kunz et al., 2018;Lampinen et al., 2020;McKinnon et al., 2018;Nilsson et al., 2021;Ramanna et al., 2020;Zheng et al., 2017). Notably, these studies use different methods to achieve compartment specific selectivity (including high b-value data (McKinnon et al., 2018;Veraart et al., 2019), gadolinium injection (Kunz et al., 2018) and planar filtering (Dhital et al., 2019)), as well as different microstructure models and/or parameter constraints to estimate axial diffusivity. Inspired by the difference between the reported and assumed intra-axonal axial diffusivity, this study first aims to explore how changing the predefined axial diffusivity in the NODDI model affected the remaining estimated parameters. A second aim of this study was to investigate whether, by utilising high b-value data, we could simultaneously estimate axial diffusivity within the NODDI framework. Here it is useful to estimate axial diffusivity within the NODDI framework (rather than e.g, from the voxels with highest fractional anisotropy (Basser et al., 1994)) as it facilitates estimation of both axial diffusivity and fibre dispersion on a voxelwise basis. Were the fibres instead assumed to be coherently aligned, the axial diffusivity estimates would be biased (Dhital et al., 2019;Howard et al., 2019). Crucial to our second aim was the use of ultra-high b-value data, where it can be assumed that the higher-diffusivity extra-axonal water is eliminated such that only the intra-axonal compartment contributes signal (Jensen et al., 2016;Kleban et al., 2020;Lundell et al., 2021;McKinnon et al., 2018;McKinnon and Jensen, 2019;McKinnon et al., 2017;Novikov et al., 2018a;Veraart et al., 2019), thus overcoming known degeneracies between diffusion characteristics of the intra-and extra-axonal space (Howard et al., 2019). However, high b-value data also posed several challenges. In particular, we found that in this low SNR regime, the model had to account for both the rectified noise floor and/or a signal offset to avoid parameter degeneracy and bias. We demonstrate how, with appropriate modifications, the NODDI framework can be applied to both magnitude and real-valued data to estimate axial diffusivities in line with previous literature.

Methods
This paper is organised as follows. First we demonstrate how changing the assumed intraaxonal axial diffusivity affects the output of the 'standard NODDI model' (described below), highlighting the importance of the discrepancy between the assumed axial diffusivity of NODDI and many estimates of axial diffusivity found in literature. We then investigate how a 'modified NODDI model' can be used to estimate both the axial diffusivity and ODI concurrently from high b-value data. Specifically, this modified NODDI model (i) considers the intra-axonal compartment only and (ii) explicitly accounts for Rician noise floor and a signal offset. Finally, the modified model was applied to in vivo, human data, where we explore how noise can bias model estimates in this high b-value and thus low SNR regime.

NODDI output sensitivity to the assumed axial diffusivity
To evaluate how the NODDI output changed with respect to the assumed axial diffusivity, the standard NODDI model (Zhang et al., 2012) was applied to the diffusion data with d ∥ = 1.7, 2.3 or 3 μm 2 /ms. Briefly, the NODDI model consists of a Watson-like fibre orientation distribution which is convolved with three compartments that are typically associated with the CSF, extra-axonal and intra-axonal space. The first compartment has isotropic, free diffusion, the second compartment describes tensor-like diffusion, and the third compartment describes stick-like diffusion. The model fitting involves 5 parameters being estimated (the intra-axonal signal fraction, the isotropic signal fraction, the fibre orientation and orientation dispersion index, [f in , f iso , θ, ϕ, ODI]) whilst assuming some fixed, global axial diffusivity which may differ from d ∥ = 1.7 μm 2 /ms, that the radial diffusivity of the tensor-like compartment is given by the tortuosity model d ⊥ = d ∥ (1 − f in ) and that the diffusivity of the isotropic compartment is that of free diffusion d iso = 3 μm 2 /ms.
Here we utilised preprocessed T1-weighted and diffusion-weighted data for the first 10 subjects of the WU-Minn Human Connectom Project (HCP); for details of the acquisition protocol and preprocessing pipeline, please see Sotiropoulos et al., 2013;Van Essen et al., 2013). Briefly, the diffusion-weighted data included 90 gradient directions each at b-values of b = 1, 2 and 3 ms/μm 2 and 18 interspersed volumes with negligible diffusion weighting. The distortion corrected ("pre-processed ") b ~ 0 ms/μm 2 data were linearly registered to each subject's T1-weighted structural scan (FLIRT (Jenkinson et al., 2002;Jenkinson and Smith, 2001)), and the T1 non-linearly registered to the MNI standard space (FNIRT (Andersson et al., 2007;Woolrich et al., 2009)). The NODDI fitting (Zhang et al., 2012) was performed in subject space using the cuDIMOT framework (Hernandez-Fernandez et al., 2019) for GPU acceleration with Rician noise modelling.

Modified NODDI for high b-value data
The NODDI model was modified for high b-value data where we assume the entirety of the diffusion signal can be attributed to the intra-axonal compartment (Jensen et al., 2016;Kleban et al., 2020;Lundell et al., 2021;McKinnon et al., 2018;McKinnon and Jensen, 2019;McKinnon et al., 2017;Novikov et al., 2018a;Veraart et al., 2019). Here, the diffusion signal along gradient direction g is given by the convolution of the fibre orientation distribution, which was assumed to be a Watson distribution (Mardia and Jupp, 2000;Zhang et al., 2012), and a fibre response function for stick-like fibres with Gaussian axial diffusion and no radial diffusion: The integrand is over x ∈ S 2 ; f in being the signal fraction of the intra-axonal compartment; μ(θ, ϕ), the direction of the fibre; d ∥ , the intra-axonal axial diffusivity; x, a unit vector on the sphere S 2 ; b, the b-value and S 0 , the non-diffusion weighted signal. C W is the normalising constant, where, C W = ∫ S 2 exp κ μ ⊤ x 2 dx = 4π 1 F 1 1/2; 3/2; κμμ ⊤ , and 1 F 1 (α; β; X) is the confluent hypergeometric function of the first kind with a matrix argument X. The dispersion parameter κ is typically rewritten in terms of the orientation dispersion index, ODI = 2/π arctan(1/κ) (Zhang et al., 2012) which ranges from 0 to 1, representing perfectly aligned and isotropic fibre distributions respectively.
The integral in Eq. 1 can also be written in terms of the confluent hypergeometric function of the first kind with matrix argument X = κμμ Tbd ∥ gg T . Consequently, Eq. 1 becomes (Sotiropoulos et al., 2012): , g ≈ f in S 0 1 F 1 1/2; 3/2; κμμ ⊤ − bd gg ⊤ 1 F 1 1/2; 3/2; κμμ ⊤ . (3) By taking the ratio of S b, g /S b (rather than S b,g /S 0) we can derive an analytic solution where the model is independent of both f in and S 0 , and described by only four parameters S(ODI, d ∥ , θ, ϕ). S b represents the powder averaged signal, i.e. the average signal across all diffusion gradients for a given b-value, where at high b-value, for stick-like diffusion convolved with an arbitrary fibre orientation distribution Kaden et al. (2016), Dividing S b,g (Eq. 3) by S b (Eq. 4), the diffusion signal has an analytic form given by: S b, g ≈ S b 2 bd π erf bd 1 F 1 1/2; 3/2; κμμ ⊤ − bd gg ⊤ 1 F 1 1/2; 3/2; κμμ ⊤ .
Here the model depends on only two free parameters, d ∥ and ODI (or κ). During fitting, S b was calculated for each shell in turn.
In our initial investigations we aimed to simultaneously estimate the three parameters d ∥ , κ and F in = f in S 0 on a voxelwise basis according to Eq. 3. However, upon closer examination, the parameters were found to be degenerate, leading to an overestimation of d ∥ . Here we used simulated data (shown later) to investigate the model parameter degeneracy and bias in relation to the noise characteristics of the data. Crucially, we found that the parameter degeneracies and bias could be overcome if we fitted to high SNR data, and explicitly model both a signal offset and the Rician noise floor. The presence of a signal offset in the in vivo data was supported by the presence of non-zero signal in the ventricles at high b-value. Consequently, we accounted for a signal offset i.e. some background signal that was independent of diffusion weighting, where Y being the data, S the diffusion signal and c the offset, which is sometimes referred to as a dot compartment. Then in magnitude data, we also accounted for the rectified Rician noise floor using Koay's inversion technique Koay and Basser (2006): ϵ being the Rician scaling parameter, which is equivalent to the standard deviation of Gaussian noise for complex data. Typically, ϵ is estimated a priori from noise estimation methods . However, as shown below, a priori estimation might not be ideal as even a small misestimation of ϵ could lead to large parameter biases. An alternative way of circumventing Rician noise floor effects is to consider real-valued data with Gaussian noise characteristics (Fan et al., 2020). Using complex data, the signal phase is first removed from each voxel after which the real component of the signal is extracted (Fan et al., 2020). In this study, the model was fitted to both real valued data (where c was estimated) and magnitude data (where both c and ϵ were estimated) from the same subjects.
As even a small misestimation of c or ϵ can lead to a large parameter bias, both were estimated as parameters of the model that were fitted voxelwise during model optimisation.

The final model & optimisation-Combining
Eq. 5 with either Eq. 6 (realvalued data) or 7 (magnitude data) and assuming a known fibre orientation μ = V 1 (the primary eigenvector of the diffusion tensor (Basser et al., 1994)), the final model was dependent on only three or four parameters: the orientation dispersion index ODI, the axial diffusivity d ∥ , the signal offset c, and, in magnitude data only, the noise floor parameter ϵ. μm 2 /ms where the diffusivity in free water at 37°C is ~ 3 -3.1 μm 2 /ms. Due to the observed possible degeneracy between the two noise parameters c and ϵ, when both were estimated in magnitude data, they were constrained to be 50 -150% of c and ϵ as estimated from high b-value data in the ventricles (a Rician distribution was fitted to b = 17.8 ms/μm 2 data, where c = the non-centrality parameter and ϵ = the scale parameter). When only one noise parameter was estimated (in real-valued data or simulations with c = 0), the noise parameter was constrained such that c ∈ [0, 0.5 · S 0 ] or ϵ ∈ [0, 0.5 · S 0 ]. The powder averaged signal S b was calculated for each shell in turn, where the data Y were for first corrected for any signal offset c or noise floor ϵ according to Eq. 6 or 7 (i.e. Y was converted to S prior to signal averaging). The model was optimised using the Metropolis Hastings (MH) algorithm which afforded estimation of each parameter's posterior distribution. The initial parameters for MH were found by grid search.
Our investigations using simulated data demonstrated how higher SNR leads to more precise estimates of ODI and d ∥ . Consequently, for in vivo data, the model was fitted to the concatenated signal across many (N) voxels to boost SNR. Here, the signal from each voxel was first rotated such that the primary eigenvector of the diffusion tensor was aligned and the secondary eigenvector was randomly orientated, since the Watson distribution describes symmetric dispersion about the primary fibre orientation. As NODDI assumes a single-fibre population per voxel, we fitted to high FA voxels from the corpus callosum which, to enforce some data consistency, were selected to have similar S 0 S 0 ± 10% . Instead of concatenating the signal, we could have averaged the signal across voxels (as we did using simulated data below). However, this would require 1) the diffusion signal from each voxel to first be rotated so that the primary fibre orientations align, and 2) the rotated signal to be resampled along consistent gradient orientations g. Consequently, signal concatenation was deemed preferable to signal averaging across voxels, as our results would not be biased by interpolation effects which can introduce smoothing and effect noise properties in non-trivial ways.

Simulated data
Data were simulated with a known ground truth to investigate the precision and accuracy of the parameter estimates. The S 0 , b-values and gradient directions were chosen to mimic the in vivo data below. The ground truth parameter values were d ∥ = 2.2, ODI = 0.03, f in = 0.6 and c = 10 unless otherwise stated. The SNR of a single voxel was defined as SNR vox = S 0 /σ = 16.5, which is similar to the lower bound of the SNR in the in vivo data used later in this study, where σ represents the standard deviation of Gaussian noise in the complex data. As in in vivo data, the model was fitted to data from N simulated voxels. The SNR across many voxels was then approximated as SNR SNR vox × N.
Note, for in vivo data we fitted the model to the concatenated signal across voxels, whilst for simulated data we fitted to the averaged signal. In simulated data, we fixed the fibre orientation and could thus compute the average signal across voxels without interpolation. Furthermore, by averaging the signal, we fitted to fewer gradient directions and minimised computation time.

In vivo human data
The modified NODDI model was applied to pre-existing data from 6 healthy participants where both magnitude and real-valued diffusion images were reconstructed from the complex MR data.
Diffusion MRI data were previously acquired, reconstructed and pre-processed according to Tian et al. (2022) and Fan et al. (2020). Briefly, complex diffusion data were acquired on the 3T MGH Connectom scanner (Siemens Healthcare, Erlangen, Germany) using gradients up to 300mT/m (Jones et al., 2018;McNab et al., 2013) and a pulsed gradient, spin echo EPI sequence (Setsompop et al., 2012): TR/TE = 4000/77 ms, 2mm isotropic resolution, 17 b-values from 0 to 17. 8 ms/μm 2 , two different diffusion times δ/Δ = 8/19 ms and δ/Δ = 8/49 ms. The complex images were first corrected for background phase contamination after which the real component of the diffusion images was extracted and the imaginary part discarded (Fan et al., 2020). The magnitude images were also extracted. The data were pre-processed using tools from FSL and bespoke code. The real-valued and magnitude data were separately corrected for gradient nonlinearities as well as susceptibility and eddy current distortions 2003;Smith et al., 2004). Here we fitted the modified NODDI model to data from 6 subjects with Δ = 49 ms, b = 6.75, 9.85 and 13.5 ms/μm 2 and 64 gradient directions per shell. The b = 0 ms/μm 2 (30 repeats) and b = 0.95 ms/μm 2 images with 32 gradient directions were used for fitting the diffusion tensor model only (Basser et al., 1994). The modified NODDI model was then fitted to selected voxels from the corpus callosum where the assumptions of the NODDI model (e.g. white matter voxels with a single fibre population) are most valid.
T1-weighted structural images were also acquired and here used for white matter segmentation (Zhang et al., 2001) and to register a corpus callosum mask from MNI standard space to subject space (Andersson et al., 2007;Woolrich et al., 2009). Fig. 2 shows example preprocessed diffusion images. In Fig. 2 right we see how the images retain good signal at very high b-values. Figs. 2 left shows the distribution of high b-value (b = 17.8 ms/μm 2 ) signal from voxels in the ventricles for a single subject. The noise from the magnitude images follows a Rician distribution, whilst that from real-valued data is Gaussian distributed as expected. The Rician non-centrality parameter and Gaussian mean are both non-zero and indicative of a positive signal offset. The data SNR, here taken to be S 0 /σ, where σ was taken to be the Guassian standard deviation or Rician scaling parameter, was found to be SNR ∈ [17, 30] with a median SNR of 21.

Changes in the NODDI output for different axial diffusivity
Figs. 3 and 4 show how the parameters of NODDI (Zhang et al., 2012) change when the assumed axial diffusivity is set to d ∥ = 1.7 μm 2 /ms -as is typically assumed -or 2.3 μm 2 /ms, or 3 μm 2 /ms, the latter being the approximate diffusivity of free water. For completeness we show results for the full brain, whilst acknowledging that NODDI parameter estimates from the grey matter are highly challenging to interpret as here the NODDI assumptions of non-exchanging compartments and negligible soma contribution are violated (Jelescu et al., 2022;Olesen et al., 2022;Palombo et al., 2020). when the assumed axial diffusivity is increased from 1.7 to 3 μm 2 /ms. Interestingly, when d ∥ = 3 μm 2 /ms, both the extra-axonal signal fraction and radial diffusivity (d ⊥ ) show two distinct distributions associated with the white and grey matter. As d ∥ is increased, the ODI is seen to increase in both the grey and white matter, and the signal fraction associated with isotropic diffusion, f iso is reduced close to zero across most of the brain and particularly in the white matter where we would not typically expect to find isotropic, free diffusion.
In Fig. 4 we see the NODDI parameter maps with axial diffusivity d ∥ = 3 μm 2 /ms as a percentage of equivalent maps for d ∥ = 1.7 μm 2 /ms. Since d ∥ = 3 μm 2 /ms represents the upper bound of water diffusion in vivo, this comparison represents the maximum difference we may obtain when increasing the assumed axial diffusivity from d ∥ = 1.7 μm 2 /ms. When d ∥ = 3 μm 2 /ms, the isotropic signal fraction and extra-axonal signal fraction decrease on average to 30 -50% and 50 -60% respectively of their value when d ∥ = 1.7 μm 2 /ms.
Concurrently, the signal fraction associated with the intra-axonal compartment and the ODI increase on average to ~ 150% and ~ 120%. Here we do not see a global, step change, but rather one which varies across the tissue. In particular, the ODI is substantially increased in the corticospinal tract as well as by 200 -250% in areas of the optic radiation and corpus callosum. We see a large ODI change across much of the corpus callosum, though not at the midline along the left-right axis, a known region of increased fibre dispersion.

Modified NODDI model for high b-value data: simulated data
A second aim of this study was to investigate whether, by applying NODDI to high b-value data, it was possible to also estimate axial diffusivity. Crucially, at high b-value, signal contributions from extra-axonal water are assumed negligible and thus the observed signal can be attributed to the intra-axonal compartment, here described by diffusion along dispersed sticks. To examine the accuracy and precision with which this high b-value 'modified' model could simultaneously estimate axial diffusivity and ODI, we began by examining simulated data, with known ground truth parameters.
3.2.1. Parameter distributions and model fits-To investigate degeneracy, Fig. 5 shows the estimated parameter distributions (left) and model fit (right) for high SNR data (N = 100, SNR ~ 165). The model was fitted to data simulated with both a Gaussian and Rician noise model, to mimic real-valued and magnitude diffusion data. In both datasets, the model appeared to fit the data well, where the residuals between the data and the predicted signal were small.
Upon inspection (Fig. 2), the in vivo data used in this study appeared to contain a signal offset (i.e the complex noise was not zero-mean), that was subsequently included as a parameter of the model (c). Fig. 5 shows how, in real-valued simulated data, the signal offset can be estimated with ease: it is not correlated to the other model parameters and does not appear to affect their estimation. In magnitude data, the signal offset is highly correlated with the noise floor parameter ϵ, causing difficulties in parameter fitting where multiple minima may exist. In comparison, for data with Gaussian noise (meaning that ϵ = 0) the parameter distributions are approximately Gaussian and close to the ground truth values, demonstrating that axial diffusivity and ODI can be estimated reliably and without degeneracy from real-valued data. Supplementary Fig. 1 shows similar plots for magnitude data without a signal offset (c = 0) where parameter estimation is again improved. Fig. 5 examines parameter estimation in relatively high SNR data. As most in vivo data have an SNR < 165, Fig.   6 shows how the precision and accuracy of the model parameters vary as a function of SNR. To produce data of varying SNR, here the model was fitted to the signal averaged over 1-1000 voxels, producing an SNR range of ~ 16.5 -520. The fitting was repeated 20 times per SNR. At low SNR, the parameter estimation is highly sensitive to noise and the parameters are degenerate, as indicated by a large spread (standard deviation) in the parameter distributions (see also green box). The precision of the model parameters increases with SNR and the parameter degeneracies are largely overcome when SNR ≳ 100.

Parameter estimation as a function of SNR-
Although an SNR of ≳ 100 is often not realised in vivo (where the in vivo data in this study have an SNR of ~ 20), Fig. 6 motivates parameter estimation through signal averaging or, as is done for the in vivo data in this study, the concatenation of signal across voxels.
In magnitude data, the axial diffusivities appear biased and overestimated with respect to the ground truth. This may be related to the calculation of the spherical mean whilst correcting for the rectified noise floor S b = Re Y b 2 − ϵ 2 − c, where the positivity constraint can lead to an overestimation of S b . In comparison, parameter accuracy is increased in real-valued data, where the mean axial diffusivity is generally more aligned with the ground truth value.
Together, Figs. 5 and 6 motivate the use of real-valued data for more reliable and accurate estimation of axial diffusivity and orientation dispersion via the modified NODDI model.

Assuming incorrect noise characteristics biases parameter estimates
-In low SNR regimes, assuming the wrong noise distribution (Rician/Gaussian), or neglecting the presence of a signal offset, could have a considerable effect on the estimated model parameters. Consequently, Fig. 7 uses simulated data to examine the effect of getting either of the noise parameters c and ϵ wrong. We examine both what happens when we neglect the Rician noise floor and signal offset completely (c, ϵ = 0), but also what happens when we get them only slightly incorrect. The latter could illustrate a situation where the noise parameters are set to global, predefined values, rather than estimating them on a voxelwise basis to account for local parameter fluctuations. In Figure 7a,b, data were simulated with a signal offset c = 10 and the model was fitted assuming some fixed offset.
In Figure 7c,d, magnitude data were simulated with c = 0 and fitted for a fixed noise floor.
Note, when the assumed noise floor = 0, this is equivalent to assuming Gaussian noise (Fig.  7d). In Fig. 7 a,c, even a small error in the noise parameter leads to biased estimates of axial diffusivity (overestimation) and ODI. This is concurrent with increased parameter degeneracy, as indicated by the high standard deviation in the parameter distributions (bottom). In Fig. 7c, the ground truth noise floor indicates the standard deviation of the complex Gaussian noise. This value is equivalent to what would be typically measured using de-noising methods and input to Koay inversion method to account for Rician bias. In contrast, the noise floor value which produces an axial diffusivity and ODI closest to the ground truth values is higher. This difference is likely related to Koay's inversion method being only well suited to data with SNR > 2, where the diffusion signal along many gradients at high b-value will likely not meet this requirement. Again, these results point to potential pitfalls when using Koay's inversion method to correct for Rician bias using a priori estimates of the noise. In Fig. 7b, where c is fixed but ϵ estimated in magnitude data, a small error in the signal offset c does not have a large effect on the other parameters. This is likely due to the correlation between c and ϵ (Fig. 5), where ϵ can be adjusted to account for the error in c.
Under the wrong noise model (Fig. 7d), the ODI and axial diffusivity estimates are both highly biased (overestimated) and correlated (degenerate), even though the SNR of the data is high (SNR ~ 165). When considering the opposite situation i.e. assuming a Rician noise model when fitting to real-valued data, the Rician noise floor ϵ tended to zero as expected, and parameter estimation was otherwise largely unaffected.

Modified NODDI model for high b-value data: in vivo data
The modified NODDI model was then fitted to high b-value in vivo data from 6 healthy subjects. Here we fitted to both magnitude and real-valued images reconstructed from the same complex data. To increase the SNR, data was concatenated across N voxels from the corpus callosum, chosen to have similar S 0 and high FA. The results are shown in Fig. 8, where the colours indicate data from different subjects.
As in simulated data (Fig. 6), the spread of the ODI/diffusivity standard deviation decreases as a function of N where in single voxel data the parameters appear degenerate, but can be estimated with increased precision for larger N. Note, whereas in simulated data, each voxel was simulated with identical ODI and axial diffusivity, in in vivo data, each voxel is likely to exhibit slightly different diffusivity and fibre orientation distribution. Thus, concatenating signal across many voxels has both the advantage of increasing SNR, and some disadvantage, where the estimated axial diffusivity and ODI are likely some average across the voxel population. The latter was minimised by selecting voxels from within the corpus callosum with similar S 0 and FA. Nonetheless, we generally obtain fairly similar axial diffusivity values when fitting to either N=50 or 100 voxels, giving some confidence to these results.
Both magnitude and real-valued data estimate somewhat similar axial diffusivities and ODI, though the spread of axial diffusivities across subjects is higher for the magnitude data. The mean intra-axonal axial diffusivity across subjects for N = 100 was d ∥ = 2.26 μm 2 /ms and d ∥ = 2.42 μm 2 /ms for real-valued and magnitude data respectively. Crucially, these values are much higher than the d ∥ = 1.7 μm 2 /ms typically assumed when fitting NODDI, suggesting that many current studies may be producing biased NODDI parameter estimates.

Discussion
To describe the complex tissue microstructure in only a handful of modelling parameters, biophysical diffusion models make strong, simplifying assumptions about the underlying tissue architecture. The extent to which many of these assumptions are valid in both healthy and diseased tissue remains an open research question. This study was inspired by the disparity between the estimated intra-axonal axial diffusivity from literature (~ 2 -2.5 μm 2 /ms (Dhital et al., 2019;Jelescu and Budde, 2017;Jelescu et al., 2020;Kaden et al., 2016;Kunz et al., 2018;Lampinen et al., 2020;McKinnon et al., 2018;Nilsson et al., 2021;Zheng et al., 2017)) and that typically assumed by the NODDI model and often utilised in simulations (d ∥ = 1.7 μm 2 /ms). We first demonstrate the extent to which the NODDI output is dependent on the assumed axial diffusivity. Second, we illustrate how the NODDI framework can be adapted for high b-value data and overcome known parameter degeneracies (Jelescu et al., 2016) to estimate axial diffusivity as a free parameter of the model. Crucially, by utilising the NODDI framework and high b-value data, we can forgo modelling of the extra-axonal space and attribute the diffusion attenuation to only two parameters, the intra-axonal axial diffusivity and fibre dispersion, which are simultaneously estimated on a voxelwise basis. Were we to estimate axial diffusivity without accounting for dispersion, our axial diffusivity estimates would be biased (Howard et al., 2019).
Our results cannot ascertain the "correct " axial diffusivity for the NODDI model, nor do the changes in the NODDI outputs reported here necessarily negate group differences in NODDI parameters. Instead they challenge the interpretation of the NODDI outputs as accurate, biophysical parameters, rather than biased indices which are dependent on many modelling assumptions, including the input axial diffusivity. In Fig. 4 -when d ∥ = 3 μm 2 /ms (the most extreme case) is compared to d ∥ = 1.7 μm 2 /ms -we frequently see the parameter estimates rise or fall by ~ 50%, where in some cases the parameter may be as little as 20% or as much as 250% of its former value. Notably, both the ODI and the signal fractions associated with each of the three compartments change. With higher axial diffusivities, a substantially larger fraction of the white matter signal is associated with the intra-axonal compartment and the associated isotropic compartment is reduced, which may be more in line with our microstructural expectations. This is coupled with an overall ~ 20% increase in the orientation dispersion index, with regions of the corpus callosum and optic tract sometimes doubling. Here, more signal attenuation perpendicular to the fibre is explained by the interplay of intra-axonal diffusion and fibre orientation dispersion, rather than the radial diffusivity associated with the extra-axonal compartment via the tortuosity model.
The change in estimated ODI invites future ODI validation against microscopy gold standards. However, validation is complicated by having to select an appropriate axial diffusivity, a choice that is especially challenging in postmortem tissue where the diffusivities are considerably reduced when compared to their in vivo values. In previous work, Schilling et al. (2017) found the NODDI ODI to correlate well with dispersion values derived from histology (r=0.66), though the ODI was consistently higher than that from histology. In comparison, Grussu et al. (2017) compared the NODDI outputs to histological data from the human spinal cord to find an approximate one-to-one mapping (r=0.84 in controls, r=0.64 in MS cases). The discrepancy between the results from the two studies may indeed be related to the choice of assumed axial diffusivity. In postmortem spinal cord, Grussu et al. assume an axial diffusivity of 1.5 μm 2 /ms which was found to optimise the fitting across all samples. Schilling et al. imaged postmortem squirrel monkey brain, though the assumed axial diffusivity wasn't reported.
Though Figs. 3 and 4 assume a single, global diffusivity across the brain and subjects, there is likely great value in being able to account for both between-subject and across-brain variations in intra-axonal axial diffusivity. Recent, promising work by Ramanna et al. (2020) and Nilsson et al. (2021) utilise triple diffusion encoding and data with multiple b-tensor encodings respectively to estimate intra-axonal axial diffusivity on a voxelwise basis. Ramanna et al. estimate axial diffusivities across the white matter of three subjects to be 2.24 ± 0.18, whilst Nilsson et al. demonstrate particularly high axial diffusivities of d ∥ ~ 2.7 μm 2 /ms in the corticospinal tract. Studies of intra-axonal axial diffusivity are typically limited to analysing data from relatively few, healthy subjects, making characterisation of normative values or between brain variations challenging. This challenge is likely further exacerbated when considering either development, ageing or pathology which likely either directly or indirectly alter the observed axial diffusivity. Indeed, although the values of axial diffusivity presented here should not be over-interpreted, our results do potentially indicate some across subject variation in intra-axonal axial diffusivity across a series of six healthy participants. Together, these arguments challenge the assumption of a fixed, predefined axial diffusivity and instead advocate for the co-estimation of both architectural features and diffusion characteristics on a voxel-wise or subject-wise basis.
A second aim of the study was to demonstrate how by modifying the NODDI model for high b-value data, the framework could be used to simultaneously estimate axial diffusivity and orientation dispersion. Though the NODDI model here explicitly considers only macroscopic fibre orientation, it is likely that microscopic fibre dispersion, including but not limited to axon undulations, also affects intra-axonal axial diffusivity. Indeed Andersson et al. (2020) and Lee et al. (2020) combine Mote Carlo simulations of diffusion with realistic axon morphologies from 3D X-ray nano-holotomography and electron microscopy data respectively, to demonstrate how complex axon morphology, including axon microdispersion, caliber variations, and the presence of mitochondria, all likely influence the observed diffusivities of the tissue. Andersson et al. explicitly show a reduction in the apparent axial diffusivity as a function of the complexity in morphological complexity, with micro-dispersion a key contributing factor.
One of the primary limitations of the NODDI framework is the description of the fibre orientation distribution as that of a single fibre population with symmetric dispersion. Consequently, we would not recommend applying the model to voxels that contain multiple fibre populations, grey matter or otherwise complex microstructure that violate these (or any other) model assumptions. Furthermore, in the current study, we assume the fibre orientation can be described by the primary eigenvector of the diffusion tensor. This allows us to concatenate the signal across multiple voxels to boost SNR. However, were the model applied to very high SNR data, therefore facilitating voxelwise model fitting, the model could be relaxed to also estimate the primary fibre orientation in single fibre voxels. Future work should consider replacing the Watson with a Bingham distribution to account for dispersion asymmetry, modelling multiple fibre populations per voxel and applying the model to high SNR data to achieve voxelwise estimates of axial diffusivity. An alternative to estimating orientation dispersion as a model parameter is to either assume isotropic dispersion (e.g. when fitting to grey matter voxels (Kroenke et al., 2004)) or integrating over all possible orientations and fitting to the powder averaged signal which is independent of the fibre orientation distribution (Kaden et al., 2016;Lundell et al., 2021;Nilsson et al., 2021;Ramanna et al., 2020). The powder averaged signal of a diffusion tensor was recently combined with data from double diffusion encoding magnetic resonance spectroscopy to confirm the assumption of negligible extra-axonal signal at high b-value (b ~ 7.2 ms/μm 2 ) and estimate d ∥,water = 1.95 μm 2 /ms in the white matter for Δ = 45ms, with lower axial diffusivities reported in the grey matter Lundell et al. (2021).
In addition, the NODDI framework assumes each axon can be described as having identical, stick-like, time-independent diffusion characteristics, where any observed signal kurtosis is attributed to the interplay of Gaussian axial diffusivity and fibre dispersion. Any deviation from this assumption, e.g. diffusion time dependence or sensitivity to axon diameter at high b-value, the presence of non-identical diffusion properties leading to a range of diffusion tensors, or other forms of intra-compartmental kurtosis (e.g. structural disorder), will likely bias our results (Henriques et al., 2021;Jespersen et al., 2019;Veraart et al., 2020).
Here we performed basic simulations of the diffusion signal from a single cylinder and the acquisition scheme used in this study (data not shown), where we only began to see substantial deviation from the S ∝ b −0.5 behaviour expected for stick-like diffusion (Veraart et al., 2020) for axons with radii r ≳ 5μm. In this study, the modified NODDI model was applied to voxels from the corpus callosum where microscopy data suggests that the majority of axon radii are ≲ 1μm, with very few axons with r ≳ 3μm (Veraart et al., 2020).
Consequently, we expect there to be minimal influence of axon diameter on our results, though caution should be taken if applying the model to data with shorter diffusion times or higher b-values, or areas of the brain with larger diameter axons, such as the spinal cord or corticospinal tract.
Nonetheless, it is reassuring how, even within this restrictive framework, we estimated axial diffusivities from real-valued data in the range d ∥ ~ 2 -2.5μm 2 /ms, with an across subject mean of 2.26 μm 2 /ms, in line with current literature (Dhital et al., 2019;Jelescu and Budde, 2017;Jelescu et al., 2020;Kaden et al., 2016;Kunz et al., 2018;Lampinen et al., 2020;McKinnon et al., 2018;Nilsson et al., 2021;Ramanna et al., 2020;Zheng et al., 2017). In real-valued data, the axial diffusivity across subjects varied by ~ 0.5μm 2 /ms. Though this variation is larger than that reported in other studies (Dhital et al., 2019;Ramanna et al., 2020), it is similar to the range of estimated axial diffusivities observed in simulated data due to sensitivity to different noise realisations (Fig. 6). Consequently, we would be cautious of interpreting this as true, subject specific variation, rather than expected parameter uncertainty given the model and the data. The variation in estimated axial diffusivities in magnitude data was larger, where simulations demonstrated the fit to be less reliable due to estimation of both the noise offset and noise floor.
A final limitation of the general NODDI model (Zhang et al., 2012) is that it assumes the axial diffusivity of the intra-and extra-axonal compartment to be equal, and that the extra-axonal radial diffusion is linked to the axial diffusivity via the tortuosity model. Various studies now provide evidence for different intra-and extra-axonal axial diffusivities (Dhital et al., 2018;Jespersen et al., 2018;Kunz et al., 2018;Skinner et al., 2017;Szczepankiewicz et al., 2015;Veraart et al., 2018). Notably, when the assumptions of the NODDI model were relaxed to also estimate intra-axonal d ∥,in , extra-axonal d ∥,ex and extra-axonal d ⊥,ex (NODDIDA), Jelescu et al. (2016) found two possible sets of solutions: one with d ∥,in > d ∥,ex , the other with d ∥,in < d ∥,ex . Consequently, it is challenging to ascertain a "correct " axial diffusivity to represent both compartments when fitting NODDI to low b-value data.
A third aim of this study was to examine the effect of getting the noise characteristics of the data wrong in low SNR data. This is of particular importance as many studies aim to acquire data at high b-value, high spatial resolution, or short scan times -all of which are typically low SNR regimes. When working with our high b-value data, both the presence of a rectified noise floor or signal offset had to be carefully considered to avoid parameter degeneracy and bias. Fig. 7 demonstrates the effect of assuming Gaussian noise in magnitude data, as may be often naively done. When the rectified Rician noise floor is not accounted for, the parameter estimates are both biased and correlated (degenerate). This challenges the idea that it is valid to naively apply biophysical models that assume Gaussian noise to low SNR magnitude images. Careful consideration should be taken when transferring models to data with different noise characteristics to what they have been designed and tested on.
Finally, Fig. 7 further advocates for estimating noise parameters within the model, rather than as fixed, predefined properties that are e.g. first measured from background voxels and assumed constant across the brain. There are various methods to remove the Rician noise bias that include but are not limited to a) utilising real-valued rather than magnitude images (Fan et al., 2020;Tian et al., 2022), b) approximating the Rician signal as Gaussian via e.g. the Koay inversion technique (Koay and Basser, 2006) and, as was not explicitly done in this study, c) applying a Rician noise model during optimisation. Note, the latter requires knowledge of the fibre orientation distribution and so cannot be applied in models fitted to the spherical mean, or powder averaged signal, which explicitly aim to circumvent modelling of the fibre orientation distribution (Fan et al., 2020). Further, as the Koay inversion technique is an approximation that breaks down for low SNR data (which likely causes the discrepancy between the ground truth ϵ and that estimated from simulated data in Figures 6,7), this study suggests that it may be preferable to use Rician noise modelling rather than Koay's inversion where possible in low SNR data.
The data used in this study is openly available via http://www.humanconnectomeproject.org/ and (Tian et al., 2022). At https://git.fmrib.ox.ac.uk/amyh/noddi-axialdiffusivity we provide a cuDIMOT implementation of NODDI (Hernandez-Fernandez et al., 2019) where the assumed diffusivities d ∥ and d iso are user-defined at runtime, and MATLAB scripts to implement the modified NODDI model. Interesting avenues for future work include: applying the model to data with both high b-value and high SNR (e.g. from animal models) to explore spatial variations in axial diffusivity; the relationship between axial diffusivity and features of the tissue architecture such as axon diameter; whether intraaxonal axial diffusivity is affected by tissue degradation or disease; the relationship between axial diffusivity and brain development; as well as the estimation of axial diffusivity in postmortem tissue, both in situ and after perfusion or immersion fixation. Furthermore, the model could be extended to estimate multiple fibre populations per voxel and applied brain-wide to provide maps of axial diffusivity across the brain. Finally, co-registered MRI and microscopy data should be used to validate the orientation dispersion estimates reported here. This will likely require the acquisition of a bespoke postmortem dataset which combines data from multiple shells at ultra-high b-values (accounting for the reduced diffusivity of fixed postmortem tissue) with corresponding microscopy imaging of the white matter fibres.

Conclusion
This study focuses on the assumption of a fixed axial diffusivity in the NODDI model in diffusion MRI. We first demonstrate a considerable dependency of the NODDI parameters on the assumed axial diffusivity, which challenges the interpretation of the NODDI parameters as biophysical metrics, rather than biased indices which are dependent on the modelling assumptions. Second, we demonstrate how the axial diffusivity could be estimated within the NODDI framework by utilising high b-value data. In this challenging, low SNR regime, we demonstrate the importance of locally estimating a rectified noise floor and signal offset i.e. background signal independent of diffusion weighting, both of which could otherwise result in parameter degeneracy or bias. Our results show an axial diffusivity of d ∥ ~ 2 -2.5 μm 2 /ms in real-valued in vivo human data, which is both in line with current literature and substantially above that typically assumed by NODDI (d ∥ = 1.7 μm 2 /ms). This motivates the use of more advanced diffusion acquisitions and/or modelling that can resolve parameter degeneracies whilst making minimal assumptions about the tissue microstructure. The modified NODDI model for co-estimation axial diffusivity and orientation dispersion in high b-value data. The fibre response function (FRF) and fibre orientation distribution (FOD) can be first convolved and then multiplied by the non-attenuated diffusion signal associated with the intra-axonal compartment (F = f in · S 0 ) to produce the diffusion signal S. A signal offset c, and (in magnitude data only) rectified noise floor ϵ, were then added to the signal to produce the real or magnitude data, Y as required. Here we use the powder averaged signal S b to avoid estimating F as a parameter of the model. See the main text for full definition of parameters. Left: Example magnitude and real-valued diffusion images from a single subject. Here we show the mean signal across all gradient directions. We observe a considerable rectified noise floor in the magnitude images, which appears greatly reduced in the real valued data. Right: The distribution of signal from voxels in the ventricles at b = 17.8 ms/μm 2 , where we assume the signal to be purely noise. In both cases, the signal is not zero-mean: the magnitude data follows a Rician distribution with the non-centrality parameter v = 13.7, and σ = 8.4; the real-valued data has Gaussian noise with a positive offset, μ = 10.4, σ = 9.2.
Consequently a signal offset and, in the case of magnitude data, rectified noise floor were estimated as parameters of the model. The dependence of NODDI parameters on the assumed axial diffusivity. The NODDI model (Zhang et al., 2012) was fitted to 10 subjects from the HCP dataset with various assumed axial diffusivities d ∥ = 1.7, 2.3, 3 μm 2 /ms. The parameter maps were co-registered after which the average map across subjects was calculated. The distribution of these parameters is shown for all voxels in the brain (top) and for the white matter (middle, dashed) and grey matter (bottom, dotted) separately. f ansio describes the signal fraction of the anisotropic compartment, where f ansio = 1 − f iso . d ⊥ is in units of μm 2 /ms. Estimated NODDI parameter maps with axial diffusivity d ∥ = 3 μm 2 /ms as a percentage of equivalent maps for d ∥ = 1.7 μm 2 /ms, as is typically assumed. The blue shows where parameters decrease, and the red where they increase, as we change d ∥ from 1.7 to 3 μm 2 .
The signal fraction associated with isotropic diffusivity and the extra-axonal compartment is substantially reduced, as is the radial diffusivity. The intra-axonal signal fraction is largely increased as is the ODI, though to a broadly lesser extent. a) The extra-axonal signal fraction is reduced substantially in the cerebellum. b,c) The ODI increased to < 250% in areas of the optic radiation and corpus callosum. Left: Parameter distributions output from MCMC. Each data point in the density plots represents a combination of parameters that fits the signal equally well. Grey dashed lines represent ground truth values. Right: Model fit (line) to simulated data (dots) and the associated error (prediction-data). g is the gradient direction, μ the fibre orientation and b the b-value in ms/μm 2 . The diffusivities are given in units of μm 2 /ms. The signal is averaged over 100 voxels, each with SNR = 16.5. Note, as the noise floor parameter is an approximation, we do not plot its 'ground truth' value. The precision and accuracy of parameter estimates as a function of SNR. The modified NODDI model was fitted to simulated data with a known ground truth (grey dashed line). In high SNR real-valued data, the estimated parameters are estimated more precisely (smaller standard deviation) and tend towards the ground truth values. In magnitude data, the axial diffusivity tends to be overestimated, even at high SNR. The green boxes show the output from MCMC at the SNR of a single voxel and 100 voxels. We see how low SNR leads to parameter degeneracy and biased parameter estimates, both of which are overcome with increased SNR. The diffusivities are given in units of μm 2 /ms. Assuming inaccurate noise parameters biases the estimated axial diffusivity and ODI. Top: data were simulated with c = 10 and fitted assuming some fixed offset. Bottom: Simulated data with offset c = 0 was fitted assuming a fixed noise floor ϵ. The plots show the mean (top) and standard deviation (std, bottom) of the parameter distribution output from MCMC. The model was fitted to high SNR data with N = 100 and the diffusivities are given in units of μm 2 /ms. Fitting the modified NODDI model to both magnitude and real-valued in vivo data. The model was fitted to the concatenated signal across N=1, 25, 50 or 100 voxels. This corresponds to an approximate SNR given on the x-axis (lower bound, assuming SNR vox = 16.5). Plotted are the mean and standard deviation (std) of the parameter distributions output from MCMC. Each colour represents data from a single subject. The diffusivities are given in units of μm 2 /ms.