Neurite orientation dispersion and density imaging of the healthy cervical spinal cord in vivo

Here we present the application of neurite orientation dispersion and density imaging (NODDI) to the healthy spinal cord in vivo. NODDI provides maps such as the intra-neurite tissue volume fraction (vin), the orientation dispersion index (ODI) and the isotropic volume fraction (viso), and here we investigate their potential for spinal cord imaging. We scanned five healthy volunteers, four of whom twice, on a 3T MRI system with a ZOOM-EPI sequence. In accordance to the published NODDI protocol, multiple b-shells were acquired at cervical level and both NODDI and diffusion tensor imaging (DTI) metrics were obtained and analysed to: i) characterise differences in grey and white matter (GM/WM); ii) assess the scan-rescan reproducibility of NODDI; iii) investigate the relationship between NODDI and DTI; and iv) compare the quality of fit of NODDI and DTI. Our results demonstrated that: i) anatomical features can be identified in NODDI maps, such as clear contrast between GM and WM in ODI; ii) the variabilities of vin and ODI are comparable to that of DTI and are driven by biological differences between subjects for ODI, have similar contribution from measurement errors and biological variation for vin, whereas viso shows higher variability, driven by measurement errors; iii) NODDI identifies potential sources contributing to DTI indices, as in the brain; and iv) NODDI outperforms DTI in terms of quality of fit. In conclusion, this work shows that NODDI is a useful model for in vivo diffusion MRI of the spinal cord, providing metrics closely related to tissue microstructure, in line with findings in the brain.


Introduction
Neurite orientation dispersion and density imaging (NODDI) (Zhang et al., 2012) is a model-based diffusion-weighted (DW) MRI technique that allows the quantification of specific microstructural features directly related to neuronal morphology. Unlike model-free techniques, such as diffusion spectrum imaging (DSI) (Wedeen et al., 2005) or diffusion propagator imaging (Descoteaux et al., 2011), which do not make any particular assumption about the local tissue microstructure, NODDI relies on the formulation of a geometric model that aims to capture the salient features of neuronal microarchitecture. This formulation tries to overcome the main limitation of phenomenological models, such as diffusion tensor imaging (DTI) (Basser et al., 1994), which are sensitive to changes in the local microstructure but only provide unspecific, surrogate information. NODDI parametrises the signal as a function of biophysically meaningful indices. In particular, the NODDI model assumes that water protons in neuronal tissue can be considered as belonging to three different pools: i) free water, modelling isotropic diffusion such as in areas contaminated by CSF; ii) restricted water within dispersed sticks, modelling dendrites and axons; and iii) anisotropically hindered water, modelling diffusion within glial cells, neuronal cell bodies and the extracellular environment. In the NODDI framework, DW data are acquired and parametric maps describing the properties of the compartments within which water pools diffuse are obtained fitting the model to the data. Such parametric maps represent indices such as neurite density and neurite orientation dispersion. The former estimates the fraction of axons and dendrites within tissue. The latter quantifies how parallel neurites are to each other. Low orientation dispersion is indicative of coherent organisation, since the orientation of single neurite elements does not deviate much from the mean overall orientation. On the other hand, high orientation dispersion occurs when neurites are dispersed in space and their orientations vary considerably from each other. In practice, NODDI can be performed in a clinically feasible time and it was recently applied with encouraging results to brain tumour (Wen et al., 2014), multiple sclerosis (MS) (Magnollay et al., 2014;Schneider et al., 2014), in vivo g-ratio estimation (Campbell et al., 2014), focal cortical dysplasia (Winston et al., 2014), NeuroImage 111 (2015) [590][591][592][593][594][595][596][597][598][599][600][601] neurofibromatosis (Billiet et al., 2014), neonatal encephalopathy (Lally et al., 2014) and healthy newborn brain (Kunz et al., 2014).
NODDI, originally developed for the brain, may potentially be a valuable and useful imaging technique also for spinal cord applications. Several diseases are known to alter the normal structure of the healthy spinal cord, such as MS (DeLuca et al., 2004;Lukas et al., 2013;Mottershead et al., 2003), spinal cord injury (Cohen-Adad et al., 2011;Tator and Fehlings, 1991), amyotrophic lateral sclerosis (Cohen-Adad et al., 2013;Sasaki et al., 1992) and others. Nonetheless, to date little investigation of neuronal morphology alteration in the pathological human spinal cord has been carried out in vivo. Although some applications of more advanced DW MRI techniques than DTI have been reported (Cohen-Adad et al., 2011;Duval et al., 2014;Farrell et al., 2008;Rangwala et al., 2013;Schneider et al., 2011), the majority of studies still rely on DTI, due to several technical challenges Wheeler-Kingshott et al., 2014). The spinal cord is a relatively small structure compared to the brain. Nonetheless, high resolution is required for precise localisation of grey and white matter (GM/WM) (Mohammadi et al., 2013). Physiological and instrumental artifacts may add undesired distortions and biases to the data Verma and Cohen-Adad, 2014). Moreover, acquisitions are often made longer due to the employment of cardiac gating in order to reduce the effects of cerebrospinal fluid (CSF) pulsation and physiological noise (Wheeler-Kingshott et al., 2002a).
The practicality of NODDI and its encouraging results in brain studies make its potential application to the spinal cord appealing. For instance, neurite density estimates in WM may be employed to characterise axonal loss or to provide new insights about the pathological mechanisms underlying spinal cord atrophy in MS, which is strongly associated to clinical disability (Lukas et al., 2013). The quantification of neurite orientation dispersion may instead be of interest at the level of nerve roots or to assess the integrity of neuronal processes. Neurite orientation dispersion has been found to be a non-negligible feature at the typical voxel scale even in coherent areas such as the human corpus callosum (Budde and Annese, 2013), positively affecting the performance of DW MRI models when accounted for (Ferizi et al., 2013). Therefore, we speculate that it may be an important feature and a potential biomarker also in another organised region such as the spinal cord.
In this work, we present the first application of NODDI to the spinal cord in vivo. As a first exploratory step, the published NODDI diffusion encoding protocol developed for brain imaging was employed to acquire data at cervical level from five healthy volunteers, four of whom were scanned twice. Results from NODDI analysis were compared to those obtained with standard DTI, routinely employed in DW MRI studies in the spinal cord . We investigated: i) NODDI metrics, characterised by region-of-interest (ROI) analysis, which focussed on the differences between GM and WM; ii) the reproducibility of NODDI in GM and WM, assessed in the subjects who underwent the second scan; iii) the relationships linking NODDI and DTI indices, supported by the results from computer simulations; and iv) the goodness of fit of both NODDI and DTI models on the acquired data.

Materials and methods
In this section we provide a description of the NODDI model and implementation, of the in vivo data acquisition and of the steps followed for the analysis of the fitted metrics.

The NODDI model and implementation
NODDI consists of a multi-compartment representation of the DW MRI signal in each voxel (the NODDI signal model), and a related acquisition protocol, optimised according to the experiment design framework of Alexander (2008).
The NODDI model represents the signal in each voxel for a pulsedgradient spin-echo (PGSE) experiment (Stejskal and Tanner, 1965) as the sum of the contribution from three non-exchanging tissue compartments (isotropic, intra and extra-neurite), and is written as: In Eq. (1), A 0 stands for the non-DW signal; v iso and A iso stand for the relaxation-weighted voxel volume fraction of the isotropic compartment, e.g. the isotropic volume fraction, and its associated signal decay; v in and A in stand for the intra-neurite volume fraction (namely the volume fraction of non-isotropic tissue occupied by neurites, also relaxation-weighted) and its relative signal decay; and A en stands for the characteristic signal decay of the extra-neurite tissue compartment, including cell bodies (Zhang et al., 2012). Recently, v in has also been referred to as neurite density index (NDI) (Billiet et al., 2014;Lally et al., 2014;Magnollay et al., 2014;Schneider et al., 2014), but here we shall consistently employ the abbreviation v in , to be consistent with Eq. (1). Also, it should be noted that in the original work that introduced NODDI (Zhang et al., 2012), the intra-neurite and extra-neurite compartments were respectively named as "intracellular" and "extracellular", and that the volume fraction of the former was indicated as v ic , rather than v in . Here, we employ the new nomenclature since it summarises with greater precision the physical sources of the diffusion-weighted signal that the compartments account for.

Subjects
For this study, we recruited five healthy volunteers (3 females, median age 34 years, range 25-47), here referred to as subjects S1 to S5. All volunteers provided informed written consent and the experimental sessions were approved by the local research ethics committee.

MRI acquisition
Subjects were scanned on a 3 T Philips Achieva scanner, equipped with maximum gradient strength of 65 mT m − 1 , employing a 16channel neurovascular receive-only RF coil. Subjects S1 to S4 were scanned twice, with the second acquisition being performed within 8 months of the first one and consisting of the same acquisition protocol of the first scan. The acquisition protocol for all scans relied on a reduced field-of-view (FOV), cardiac gated PGSE ZOOM-EPI sequence (Symms et al., 2000;Wheeler-Kingshott et al., 2002b) with outer volume suppression (Wilm et al., 2007). A peripheral pulse oximeter was employed for monitoring the cardiac cycle and triggering cardiac gating. We implemented the diffusion encoding scheme according to the published NODDI protocol (Zhang et al., 2012), acquiring sequentially two high angular resolution diffusion imaging (HARDI) (Tuch et al., 2002) shells. The first shell consisted of thirty measurements at b = 711 s mm −2 , whereas the second shell consisted of sixty measurements at b = 2855 s mm −2 . Six non-DW measurements were also acquired, interleaved with the DW ones (three b = 0 images for each shell).
Scans were performed axial-oblique, i.e. perpendicular to the cord longitudinal direction, which was carefully aligned with the sliceselection direction (z) on a sagittal localiser. Twelve slices were acquired with the following parameters: TR = 4 RR repeats, TE = 65.50 ms (minimum TE achievable for a diffusion-weighting strength of b = 2855 s mm − 2 ), reduced FOV of 64 × 48 mm 2 , SENSE factor of 1.5 in the anterior-posterior direction, resolution of 1 × 1 × 5 mm 3 , and triggering delay for cardiac gating of 150 ms. Different b-values were achieved by varying the gradient strength, while fixing Δ = 32.20ms and δ = 20.50ms for both shells. This corresponded to an effective diffusion time of τ d = Δ − δ/3 = 25.37ms. The total acquisition time, depending on subjects' heart rate, was approximately TA = 35 min, including survey scans and coil calibrations.

Pre-processing of the in vivo data
Pre-processing of the in vivo data was performed before running any further analysis, in order to reduce the effects of motion-related artifacts. In this work, we corrected for motion with a slice-wise linear registration, implemented using FSL FLIRT (Greve and Fischl, 2009;Jenkinson et al., 2002;Jenkinson and Smith, 2001) and described in detail in Appendix A. This approach has shown the highest benefits in DW MRI of the spinal cord (Mohammadi et al., 2013). Briefly, for each slice, the first acquired image, which was non-DW (b = 0), was chosen as the registration target. All following non-DW images were then registered to the target and the estimated registration transformations stored. Lastly, DW images, which were acquired in-between the b = 0 ones, were also warped to the target employing the transformations previously estimated, so that each DW image was warped employing the transformation of the closest preceding non-DW one.

Segmentation
For each subject and scan, whole cord segmentation was carried out on the mean b = 0 volume evaluated after motion correction. A semiautomatic active surface method (Horsfield et al., 2010) implemented in Jim (http://www.xinapse.com/home.php) was employed to obtain the cord outline ("cord finder" tool) and to create a binary mask (fitting cord mask) of the whole cord. Salient parameters adopted for the outlining were: nominal cord diameter of 8 mm, number of shape coefficients of 25, and order of longitudinal variation of 7. These parameters control the radial expansion of the surface describing the spinal cord and constrain its spatial smoothness, as described in detail elsewhere (Horsfield et al., 2010). The fitting cord mask was then eroded slice-by-slice to limit CSF contamination and cropped to the 6 central slices (whole-cord mask).
Afterwards, manual grey matter (GM) outlining was carried out on the average DW volume obtained according to the procedure described in Kearney et al. (2014). The GM and WM masks were obtained directly from the manually drawn outlines and were limited to voxels within the whole-cord mask.

Model fitting
The NODDI and DTI models were fitted within the fitting cord mask of all scans. The NODDI MATLAB (The MathWorks, Inc., Natick, Massachusetts, USA) Toolbox (http://nitrc.org/projects/noddi_toolbox) was employed for NODDI, whereas in-house MATLAB code was employed to fit DTI. In order to fit the NODDI model, the two diffusivities representing the diffusion coefficient of the isotropic compartment (d iso ) and the intrinsic diffusivity of the intra-neurite compartments (d || ) were fixed as in Zhang et al. (2012) to d iso = 3.00μm 2 ms −1 and d || = 1.70μm 2 ms −1 , which are the values commonly employed in literature for the free diffusivity of water particles in CSF and neural tissue in vivo at body temperature. Both NODDI and DTI were fitted to the whole double-shell data set. Furthermore, DTI was also fitted to a reduced set of measurements, corresponding to the shell acquired at b = 711 s mm −2 . This was performed in order to evaluate the diffusion tensor in a regime where the contribution of non-Gaussian diffusion is small (Clark and Le Bihan, 2000;Farrell et al., 2008;Frøhlich et al., 2006). In practice, the DTI metrics obtained from the b = 711 s mm −2 shell were used to study the reproducibility and the relationship between DTI and NODDI, whereas DTI fitting performed to the whole double-shell set of measurements aimed to assess the overall goodness of fit of the model.
The fitting was carried out maximising the likelihood for a Rician noise model (Gudbjartsson and Patz, 1995) and the same objective function routines were employed to fit DTI and NODDI. The signal level at b = 0 was estimated as the mean of the non-DW measurements, as by default in the NODDI toolbox. Practically, the optimal parameters were estimated with a gradient descent algorithm, initialised by a grid search, which minimised the opposite of the log-likelihood. The parameter σ controlling the spread of the Rician distribution was always fixed before running the fitting and it was estimated voxel-by-voxel as the standard deviation of the non-DW measurements. In the freely available NODDI toolbox, the estimate of σ obtained as such is by default further scaled by a factor of 100, practically increasing the apparent SNR and pulling the noise model towards a Gaussian regime. However, in this work, we did not scale σ, since simulations proved that such a scaling can bias the fitted metrics if true Rician noise is added to the measurements at SNR levels of 10 or lower.
The following voxel-wise maps were obtained. For NODDI: the isotropic volume fraction (v iso ), the intra-neurite volume fraction (v in ) and the orientation dispersion index (ODI). Furthermore, we also calculated the effective volume fraction of the intra-neurite (restricted) com- For DTI, we evaluated fractional anisotropy (FA), axial, radial and mean diffusivity (AD, RD and MD respectively).

Analysis
Our analysis aimed to achieve four main goals: 1. characterisation of the fitted metrics; 2. investigation of the reproducibility of NODDI metrics; 3. investigation of the relationship between NODDI and DTI metrics; 4. characterisation of the quality of fit of NODDI and DTI.

Characterisation of the fitted metrics
We visually inspected the fitted metrics in all subjects and scans and characterised regional variation calculating the medians of each metric within each ROI (GM, WM and whole-cord). We also quantified the contrast (C) and the contrast-to-noise ratio (CNR) between GM and WM for all NODDI and DTI metrics as done in other studies (Gringel et al., 2009). We employed the relations Above, symbols μ WM , μ GM , σ WM 2 and σ GM 2 respectively indicate the mean value of a metric within the WM and the GM masks, and the variance of the same metric within the WM and the GM masks. C and CNR were calculated only for the first scan of the subjects. For DTI, medians, C and CNR were calculated from the metrics derived from the b = 711 s mm −2 shell.

Reproducibility
The reproducibility of all NODDI metrics was investigated in GM, in WM and at whole-cord level separately, and the investigation aimed to characterise the variability of each metric within these three ROIs across subjects and scans. In particular, two main aspects were studied: i) the total variability of each metric; and ii) the within-subject and between-subject contributions (Bartlett and Frost, 2008), in order to estimate the fraction of variability respectively due to measurement errors and biological differences in the recruited cohort. The statistical analysis was performed in MATLAB and was carried out consistently for all ROIs and metrics in subjects S1 to S4.
In practice, for each ROI we calculated a percentage coefficient of variation (CV) and the intraclass correlation coefficient (ICC), as explained in Appendix B. CV relates the total variability of a metric to the average value of the same metric observed across subjects and scans, and is expressed in percentage points. On the other hand, ICC ranges from 0 to 1 and estimates the fraction of the total variability due to biological differences between subjects. It follows that the fraction of the total variability due to measurements errors (i.e. within-subject variability) can be estimated as 1 − ICC. For comparison, the two statistics CV and ICC were also calculated for DTI indices obtained from the b = 711 s mm −2 shell.

Relationship between NODDI and DTI indices
We investigated the experimental relationship between NODDI metrics v in and ODI, and DTI indices, with scatter plots of v in and ODI, colour-coded according to the values of FA, AD, RD and MD in turn. The analysis reveals how key microstructural features such as neurite density and orientation dispersion contribute to the patterns of DTI indices in the spinal cord. We employed DTI metrics obtained from the b = 711 s mm −2 shell and we focussed on voxels with negligible isotropic volume fraction (v iso b 0.05) to match as much as possible the assumption v iso = 0 made to perform our simulations (in line with the underlying NODDI model hypothesis according to which the contribution of the isotropic compartment in areas not contaminated by CSF should be low).
In order to support the relationships observed in vivo, we also ran computer simulations and evaluated the theoretical patterns of FA, AD, RD and MD as functions of v in and ODI that would be observed in a tissue perfectly matching the assumptions of the NODDI model at b = 711 s mm −2 , as described in detail in Appendix C. In practice, we sampled metrics v in and ODI in a uniform grid of 64 × 64 values in the interval [0.05; 0.95] × [0.005; 0.5]. We visually matched the results from simulations and the scatter plots of the in vivo data as follows. Values of v in and ODI from the five subjects were quantised and mapped to the same grid employed for simulations. Then, values of DTI metrics from different voxels whose pairs (v in , ODI) were mapped to the same location in the discrete grid were averaged, before colour-coding the scatter plots.

Quality of fit of NODDI and DTI
We investigated the quality of fit of both NODDI and DTI models on the full set of DW measurements of the five subjects (rescans omitted). The Bayesian Information Criterion (BIC) (Schwarz, 1978) was employed to compare the quality of fit of the two models, as in Ferizi et al. (2014), Ferizi et al. (2013), and Panagiotaki et al. (2012). We calculated in all voxels within the fitting mask the percentage relative difference between BIC values of NODDI and DTI, with respect to those of DTI, as with BIC NODDI and BIC DTI indicating respectively BIC values for NODDI and DTI.

In vivo data acquisition
Illustrative raw images and ROIs relative to several slices of the first scan of one of the subjects (S2) are reported in Fig. 1. In the first column, the mean b = 0 image and the outline of the masks are shown. In the second column, DW images obtained at b = 711 s mm −2 and for a gradient direction almost perpendicular to the cord longitudinal axis are displayed. In the third column, we report DW images obtained at b = 711 s mm −2 but for a diffusion encoding gradient nearly along the cord axis. Lastly, in the fourth and fifth columns, we show DW images as in the second ad third columns but for b = 2855 s mm −2 . The first column shows that the fitting cord mask (in yellow) is likely to contain voxels with CSF partial volume. Secondly, the visual inspection of the DW images reveals that diffusion weighting with gradients nearly along the cord axis causes higher attenuation than diffusion weighting with gradients almost orthogonal to the same direction. Lastly, it can be noticed that at b = 2855 s mm −2 , gradients nearly along the cord axis cause the MRI signal to decay almost completely to the noise floor.

Analysis
3.2.1. Characterisation of the fitted metrics Fig. 2 shows illustrative voxel-wise NODDI metrics and DTI FA for a slice from both scans of subjects S2 and S3. The specific cases reported in the figure, consistent with the metrics trend observed in all subjects, show a good anatomical correspondence between scan and rescan.
The maps reveal a number of facts. Firstly, v iso is maximum on the cord boundaries, where CSF partial volume is likely to occur. Moreover, although this metric is often close to 0, v iso is close to 0.1 in a significant proportion of voxels. Secondly, v in is close to 1 on the boundary of the fitting mask, where voxels are likely to be characterised by partial volume with CSF. In those voxels, v iso is close to 1 and the effective volume fraction of the restricted compartment (v r ) is therefore very low. Also, on visual inspection, v r shows slightly reduced contrast between GM and WM areas compared to v in . Lastly, ODI shows a good contrast between GM and WM, comparable to that seen in FA maps. Fig. 3 shows the medians of the fitted metrics within each ROI (GM, WM and whole-cord) for the five volunteers (rescans omitted), and Table 1 summarises the median and the range of these five values. Firstly, the figure confirms that v iso is non-negligible on average in WM (median across the five subjects of 0.12, whereas it is 0.004 in GM). Secondly, v in shows more contrast between GM and WM compared to v r : the medians across the five subjects are 0.49 in GM and 0.57 in WM for v in and 0.45 in GM 0.49 in WM for v r . Next, the ODI contrast between GM and WM seen in Fig. 2 corresponds to values of this metric approximately three times higher in GM than in WM (median of ODI across the five volunteers of 0.086 in GM and 0.027 in WM). Lastly, Fig. 3 shows that as far as DTI metrics are concerned, FA is higher in WM compared to GM (median across the five subjects of 0.57 and 0.80 in GM and WM respectively), AD and MD are higher in WM than in GM (medians across the five subjects of 1.60μm 2 ms −1 and 0.92μm 2 ms −1 in GM and of 2.16μm 2 ms −1 and 0.97μm 2 ms −1 in WM for AD and MD respectively) and RD is higher in GM compared to WM (median across the five subjects of 0.54μm 2 ms −1 in GM and 0.36μm 2 ms −1 in WM). Fig. 4 shows the contrast C and the contrast-to-noise CNR between GM and WM for the five subjects, reported as box plots. NODDI metrics v iso and especially ODI provide higher values of C than DTI indices, whereas v in and v r provide similar values of such index of contrast to those of DTI MD. The CNR of NODDI metrics is comparable to that of DTI. In particular, ODI shows the highest CNR among all NODDI indices, and such value is slightly lower than the CNR of DTI FA and AD. The latter metric shows the highest CNR among all indices.

Reproducibility
The reproducibility analysis aims to quantify the total variability associated to the fitted metrics, and estimates how biological variation in the cohort and measurement errors contribute to this total variability. For this purpose, the indices CV (quantifying the total variability of a metric with respect to its mean across subjects and scans) and ICC (estimating the fraction of the total variability due to biological variation) were calculated and their values are reported in Tables 2 and 3.
The CV for NODDI parameters is the highest for v iso (140% in GM, 41% in WM and 69% in the whole-cord ROIs) and is below 10% for v in and v r in all ROIs. Furthermore, CV for ODI is equal to 44% in GM, 7% in WM and whole-cord ROIs. The ICC is bigger than 0.5 for v iso only in GM, and is just above 0.5 for v in in all ROIs (ICC of 0.54, 0.62, 0.54 in GM, WM and whole-cord ROIs respectively) and it is no less than 0.66 for v r and ODI Example of voxel-wise metrics for a slice from the two scans performed on subjects S2 and S3. Top to bottom: first row refers to subject S2, first scan, slice z = 7; second row refers to subject S2, second scan, slice z = 7; third row refers to subject S3, first scan, slice z = 9; fourth row refers to subject S3, second scan, slice z = 9. Left to right: first column shows the mean b = 0 image with overlaying fitting mask (yellow), WM mask (light blue) and GM mask (red); second column shows NODDI isotropic volume fraction (v iso ); third column shows NODDI intra-neurite volume fraction (v in ); fourth column shows NODDI effective volume fraction of the restricted compartment (v r ); fifth column shows NODDI orientation dispersion index The CV for DTI metrics is smaller than 10% for FA in the WM and whole-cord ROIs and for AD and MD in all ROIs. CV is equal to 13% in GM for FA and it is equal to 16% in GM, 19% in WM and 17% in the whole-cord ROI for RD. Lastly, the ICC statistic is well above 0.50 for all DTI metrics in all ROIs. For instance, it was maximum for FA in GM (0.95) and minimum for FA in WM and whole-cord ROIs (0.75). Fig. 5 shows the relationships between DTI metrics, evaluated at low b-value, and NODDI indices v in and ODI, for negligible v iso . We report the theoretical relationships for SNR → ∞ and for SNR = 10 respectively to the left and in the central column, whereas we plot the relationships observed in vivo to the right.

Relationship between NODDI and DTI indices
It is apparent from the synthetic data with SNR → ∞ that several different combinations of v in and ODI, simulating different cytoarchitectures of the neural tissue, produce the same value of a DTI metric such as FA. For instance, a value of FA close to 0.80, typical of WM, can potentially be explained by a broad range of combinations of v in and ODI (see Fig. 4.A), such as for v in = 0.65 and ODI = 0.06 (relatively high neurite density but modest neurite coherence) or for v in = 0.55 and ODI = 0.02 (lower neurite density but neurites well parallel to each other). Furthermore, an increase of FA can be independently caused by an increase in v in or a decrease in ODI. Also, AD decreases with increasing ODI, and has little dependence on v in for combinations of ODI and v in that are likely to be measured in the spinal cord in vivo (e.g. ODI smaller than 0.25 and v in greater than 0.25). As far as RD is concerned, increasing values of v in imply a reduction of RD, given the lower amount of diffusion occurring in the extra-neurite compartment, whereas an increase of ODI contributes to an increased RD, if v in is not too low, e.g. higher than 0.2. Lastly, MD shows little dependence on ODI in general, and decreases for increasing v in , as a result of increased contribution from the restricted compartment. For completeness, it is worth noticing that the simulation suggests that DTI can underestimate the intrinsic diffusivities of the substrate if restricted diffusion dominates; this is the scenario for very high v in values, which however was not observed in vivo.
Results from simulations at a SNR level of 10 are reported in the central plots of Fig. 4. They represent the relationship between DTI indices and metrics v in and ODI that could be observed if the imaged tissue matched exactly the NODDI model assumptions and the parameters employed to run the simulations, for a SNR level plausible in the spinal cord in vivo. Hence, they can be considered as the ground truth, which  , v r in C), and ODI in D). Plots in the bottom row refer to DTI metrics: FA in E); AD in F); RD in G); and MD in H). For each ROI, different markers stand for different subjects and represent the median of a metric within that particular ROI. Colours encode subjects: black stands for S1, red stands for S2, green stands for S3, blue stands for S4 and light blue stands for S5.

Table 1
Medians and ranges (within round brackets) across the five subjects of NODDI and DTI metrics for all ROIs. The medians reported here are the medians of the five data points (black, red, green, blue and light blue) shown in Fig. 3 the experimental scatter plots should be compared to. The presence of Rician noise does not prevent the possibility of identifying the same general patterns seen at SNR → ∞, although contours are not as precisely delineated as before.
In vivo scatter plots, reported to the right hand side of Fig. 4, replicate the trends observed for simulations at SNR = 10 well. In the observed range, FA increases for increasing v in and decreasing ODI. AD depends more on ODI rather than v in , especially if the former metric is very low, and appears slightly higher than values provided by simulations. RD increases for decreasing v in and for increasing ODI. Lastly, although the patterns exhibited by MD are less clear than in simulations at SNR = 10, MD decreases as v in increases, especially if ODI is not too low. Also, MD appears slightly higher than the theoretical values obtained from synthetic data. Overall, in vivo scatter plots confirm results from simulations since the same value of all DTI metrics can be observed in voxels characterised by combinations of v in and ODI considerably different to each other.

Quality of fit of NODDI and DTI
In this section, we report the results from the BIC analysis. BIC is a useful statistics that quantifies how well a model fits to some data while controlling for model complexity (number of model parameters), and here we employed it to compare NODDI and DTI. For this purpose, we calculated the relative percentage difference between BIC values of the NODDI and DTI models, here referred to as δBIC. Fig. 6 shows voxel-wise δBIC maps in slices 4 to 9 of the first scan of all subjects, within the fitting mask. In Fig. 6, voxels with negative δBIC favour NODDI, positive DTI, since we always obtained positive BIC values for both models and since lower BIC implies better fitting. It is clear from the illustration that the number of voxels with negative δBIC (supporting NODDI) is overwhelmingly higher than those with positive δBIC (supporting DTI): a few voxels with δBIC N 0 are orange and yellow in slices 4, 5, 7, 8 and 9 of subject S2 and in slices 4, 5 and 6 of S5. The figure also demonstrates that voxels with extremely negative δBIC values, as low as roughly −70% (dark blue and violet voxels), are located in the cord boundary, as shown in slice 4 and 8 of subject S1, slice 4 of subject S2 or in slice 7 of subject S5. Results were consistent with other slices 1 to 3 and 10 to 12.

Discussion
In this work we have studied for the first time the application of NODDI to the spinal cord in vivo. NODDI provides metrics that map directly neurite architecture and that could potentially be important biomarkers for spinal cord conditions, and here we have shown the feasibility of applying the published technique to the healthy cervical cord in vivo and in a clinical setting.
Our work reveals that several known anatomical features can be identified in NODDI metrics, but also suggests that further investigation is required to better relate them to the real microstructural characteristics of the spinal cord, and in relation to other possible diffusion DW MRI models employed in brain applications, such as "ball and stick" (Behrens et al., 2003), CHARMED (Assaf et al., 2004), AxCaliber (Assaf et al., 2008) and others. The trend for our DTI indices is consistent with previous findings (Vedantam et al., 2013;Wheeler-Kingshott et al., 2002a;Xu et al., 2013) and in agreement with tract-specific measures in the lateral and dorsal columns of the cervical cord WM (Smith et al., 2010).
NODDI ODI shows the highest contrast between GM and WM among all fitted metrics, surpassing DTI indices. ODI in GM is roughly three times bigger than in WM, reflecting known existing differences Fig. 4. To the left, in a), box plot of the contrast between GM and WM for NODDI and DTI metrics. To the right, in b), box plot of corresponding CNR values between GM and WM. The plots are obtained studying data from the first scan of all subjects, omitting rescans.  between the two tissue types in terms of neurite architecture and confirming qualitative findings in the brain (Zhang et al., 2012). On the other hand, volume fractions v in and v r show less contrast between GM and WM compared to ODI, similar to that of DTI MD. Though no clear quantitative comparison between GM and WM was presented in Zhang et al. (2012), results reported here seem to suggest that in the spinal cord v in and v r are more homogeneous between GM and WM compared to the brain. This could simply be a partial volume effect due to the small number of voxels in each tissue type (GM/WM), which increases the chance of contamination of the indices. On the other hand, since v in is designed to map the density of both axons (prominent in WM) and dendrites (mainly found in GM), it may be possible that the reduced contrast between the two tissue types is indicative of the presence of either axons in GM or dendrites in WM, which means a reduced heterogeneity between the two tissue types in the spinal cord compared to the brain in terms of neurite density, possibly also contributing to the relatively high FA in GM. For instance, it is known from anatomy that GM strands, separated by interwoven nerve fibres, invade the WM lateral funiculus, especially at cervical level (formatio reticularis). Furthermore, axons originating from white funiculi and spinal nerve roots are known to be present in spinal GM: as an example, nerve fibre bundles are found in laminae I, II and IV of the posterior grey horns. The calculation of CNR demonstrates that the CNR levels of NODDI are similar to those of DTI. Also, although ODI has the highest contrast among all metrics, its CNR is surpassed by FA and AD. This may be an effect of the higher variability of NODDI indices compared to those of DTI, which also affects the reproducibility scores. The quantification of the The in vivo scatter plots were instead obtained omitting rescans and studying voxels within the whole-cord mask of subjects S1 to S5 where v iso b 0.05. DTI metrics used for colour-coding were obtained from the b = 711 s mm −2 shell.
contrast and of the CNR therefore points towards the potential of NODDI for spinal cord applications, since the technique clearly differentiates the microstructural differences of GM and WM in terms of neurite orientation dispersion. However, future optimisation of the technique for the spinal cord should help to reduce the variability of the metrics, and hence to enhance the CNR measured in this work.
Our analysis also shows that v iso is maximum in the cord boundaries, where CSF contamination is likely. Potential CSF contamination is also suggested in those areas by v in and v r , which are respectively very high and very low. Moreover, when v iso is very high, namely v iso → 1, then v in is also usually high and differs considerably from v r . However, in such cases v in is poorly defined as the tissue signal is very low. We also observe that an amount of v iso close to 0.1 is seen in WM. Residual CSF partial volume effects in correspondence of the anterior median fissure and the central canal may contribute in voxels within the neural tissue, but it is likely that v iso explains the fraction of signal decay that can not be associated directly with restricted diffusion, i.e. relatively isotropic diffusion but slower than diffusion in CSF. We speculate that such relatively isotropic contributions to the DW signal in WM may partly arise from the water pool within the biggest axons. A tiny but nonnegligible percentage of axons is known to be characterised by large diameters (Feirabend et al., 2002;Häggqvist, 1936;Wesselink et al., 1998) with the tails of the fibre diameter distribution at cervical level extending up to 15 μm and beyond (Makino et al., 1996). These axons, even if not numerous, can make a substantial contribution to the DW signal, given their large size and therefore relatively large volume fraction . Since in this study we employed a diffusion time of 25.37 ms, which would imply a radial mean square diffusion displacement of L d ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2D τ d p ¼ 10:5μm, assuming a value for the free diffusivity of water particles in neural tissue at body temperature of D = 2.16μm 2 ms −1 (median AD in WM from Table 1), it can be hypothesised that spins diffusing within largest axons did not fully explore the boundaries of the compartment making diffusion appear relatively free. Moreover, the presence of collateral WM fibres (Lundell et al., 2011) or heterogeneity in terms of glial cell composition between GM and WM (Liuzzi and Miller, 1987) is not captured explicitly in the NODDI model, and may have had an impact on the observed patterns of the metrics. A preliminary analysis included as supplementary material suggests that collateral fibres in WM may be associated to a slight increase of isotropic volume fraction and neurite orientation dispersion, although a bigger sample size is needed for a definite conclusion. Radially oriented astrocytes, characteristic of the WM glial cell population in the spinal cord of mammals (Liuzzi and Miller, 1987), may contribute in a similar way, being their fine processes also collateral to the direction of WM fibre bundles.
The reproducibility analysis quantifies the general trend in terms of within and between-subject variability. Results suggest that metrics v in , v r and ODI are characterised by a total variability across subjects and scans of the same order of DTI metrics, whereas the total variability of v iso is considerably greater. Nonetheless, the total variability of ODI is slightly higher than that of v in and v r . Moreover, the calculation of the ICC shows that the total variability is driven by different sources for v iso , v in , v r and ODI. An ICC greater than 0.5 suggests that the variability due to biological differences is bigger than that of measurement errors, whereas the opposite holds if ICC is smaller than 0.5. Here, some of the metrics are characterised by ICC values much greater than 0.5, whilst in other cases ICC was well below 0.5. Biological differences are predominantly captured by v r and ODI, whereas measurement errors dominate for v iso . Similar contribution between these two sources of variability is instead seen for v in . The worst reproducibility scores are observed for v iso , which is not surprising given that it is in general very low, and negligible in GM ( Table 1 shows that very little contribution of the isotropic compartment is seen in GM, causing v iso to be poorly defined). For instance, a CV of 140% in GM, caused by a tiny median v iso across subjects and scan, smaller than the total variability, suggests that big cohorts of subjects may be required to detect variations of v iso caused by pathology. v iso is intrinsically highly variable in the healthy GM and values ranging from 0 to about 0.1 appear all as biologically plausible. On the other hand, for all DTI indices, the total variability is always associated almost entirely to biological variation. We speculate that the higher amount of variability and dependence on measurement errors of NODDI metrics compared to DTI could be due to the fact that NODDI non-linear signal model is more susceptible to noise and to signal distortions arising from the unavoidable magnetic field inhomogeneities that hamper spinal cord MRI Summers et al., 2014;Verma and Cohen-Adad, 2014). Such distortions were likely to be different between scans and rescans, because of different shimming and cord alignment with the main field. Nonetheless, overall the analysis shows that the reproducibility of the three metrics v in , v r and ODI is sufficient to warrant future efforts to optimise the acquisition for spinal cord tissue characteristics as well as the employment of NODDI in its current form. Our simulations allow the visualisation of the relationships between NODDI metrics v in and ODI and DTI indices, in the case of low diffusion weighting strength, i.e. in a regime where non-Gaussian contributions are small (Clark and Le Bihan, 2000;Farrell et al., 2008;Frøhlich et al., 2006). The colour-coded scatter plots presented here agree with the trends obtained at SNR = 10 and even at SNR → ∞, and show that different configurations of parameters v in and ODI can potentially lead to similar values of each DTI index, hence replicating findings in the brain (Zhang et al., 2012). Therefore, as long as metrics v in and ODI are believed to provide physical quantities mapping directly neurite morphology, NODDI appears to disentangle specific factors contributing to the exhibited patterns of DTI indices in the spinal cord in vivo.
The last piece of analysis aims to compare the quality of fit of NODDI and DTI. For this purpose, we compare voxel-wise BIC maps calculated for both models. The results show clearly that NODDI does fit the DW signal better than DTI in the vast majority of voxels, since the index δBIC, i.e. the percentage relative difference between BIC values of NODDI and DTI with respect to those of DTI, is positive in only a few scattered voxels. Moreover, the results reveal that voxels where NODDI outperforms DTI in terms of goodness of fit (e.g. voxels with δBIC values as negative as −70% or more) are located in the cord border, where CSF contamination is likely, thus suggesting that multicompartment models can provide a better description of the signal in areas with elevated partial volume effects.
Despite the achievements of this work, our approach presents a number of limitations. The most significant is the fact that the diffusion encoding protocol employed for NODDI as published in (Zhang et al., 2012) is not optimised for the spinal cord. For instance, single fibre orientation optimisation approaches similar to that followed in Schneider et al. (2010) or multi-band strategies as in Wen et al. (2014) may be adopted in order to reduce dramatically the acquisition time. This would facilitate the inclusion of NODDI in studies where the scan time is limited, which at present may not be straightforward given a total acquisition time of roughly 35 min. Also, different priors for the intrinsic diffusivities in the signal model may be used. As an example, the prior d iso = 3.00μm 2 ms −1 for the diffusion coefficient of the isotropic compartment may not be optimal in areas where CSF flow related effects or vascular components are not negligible, such as in those portions of tissue occupied by the branches of the anterior spinal artery and vein, which pass into the anterior median fissure. In fact, in those areas an apparent diffusivity greater than that of free water at 37°C may be needed to explain isotropic signal decay, as observed for instance in Baron and Beaulieu (2014) for the brain ventricles in a PGSE experiment with diffusion time of 40 ms. Preliminary analysis not included in this manuscript found that MD of voxels in areas within the CSF and likely to be affected by partial volume with large vessel, is on average roughly equal to 4.30μm 2 ms −1 . Exploratory tests demonstrated that employing priors for d iso bigger than 3.00μm 2 ms −1 and up to 4.30μm 2 ms −1 , may help to reduce the value of fitted v iso . Nevertheless, since this is the first attempt to apply the full brain NODDI protocol to the spinal cord in vivo, we adhered to the same diffusivity priors adopted in Zhang et al. (2012) for consistency. Another limitation of the work is that the reproducibility of DTI, as reported here, may be underestimated, since its signal model was fitted to a smaller number of measurements than NODDI and to a b-shell that was not acquired employing the shortest echo time achievable on the scanner at that particular b-value. Here, the same echo time was set for both b-shells to achieve the same T2weighting for the three tissue compartments. This choice though caused the SNR of the lowest shell to be sub-optimal.
A further limitation of this study is related to the motion correction strategy. In the current implementation, motion corrupting the DW images was corrected employing transformations evaluated among non-DW images, due to the challenge of registering reliably images obtained at very high b-values with those obtained at b = 0, as reported in other spinal cord studies performed at high diffusion weighting (Farrell et al., 2008) and as evident from the observation of Fig. 1. Although in the literature registration transformations have been estimated directly for the DW volumes (Mohammadi et al., 2013;Xu et al., 2013) or at least for the mean DW volume (Cohen-Adad et al., 2011), it should be noted that in such studies the diffusion weighting strength was much lower than here: the maximum b-value was equal to 500 s mm −2 in Mohammadi et al. (2013), to 800 s mm −2 in Xu et al. (2013) and to 1000 s mm −2 in Cohen- Adad et al. (2011). Our procedure implicitly assumes that negligible motion has occurred during the acquisition of the DW volumes between two subsequent b = 0 ones, and visual inspection of the uncorrected and corrected data proved that this hypothesis described our data reasonably well. Nevertheless, improvements to the NODDI analysis pipeline could be achieved with a more precise modelling of the patterns of motion occurring during the acquisition of the DW images.
Lastly, the a priori choice for the intrinsic diffusivity of the neural tissue, i.e. d || = 1.70μm 2 ms −1 , may not be optimal in the presence of pathological phenomena such as inflammation and degeneration. In these cases, different values for d || may be necessary, and data-driven approaches could be followed to estimate suitable d || priors, such as analysis of DTI-derived AD in WM. More complex fitting procedures relying on Markov chain Monte Carlo (MCMC) methods may also be employed, in order to fit d || and to obtain distributions of plausible d || values.

Conclusion
In summary, here we have studied for the first time multi-shell DW MRI data of the healthy cervical cord in vivo with NODDI, and related the results to those provided by routine DTI. The main achievement of this work is to demonstrate that NODDI can be applied successfully in the cord on a standard clinical system, and that the technique replicates key findings in the brain, as reported in Zhang et al. (2012). In particular, NODDI fits the acquired data better than DTI and neurite orientation dispersion differentiates clearly GM and WM appearing as a key factor determining the level of diffusion anisotropy. This result is in line with recent findings showing that fibre orientation dispersion is an important feature at the voxel scale even in highly organised areas such as the human corpus callosum (Budde and Annese, 2013;Ferizi et al., 2013). Moreover, the reproducibility of NODDI, quantified in both GM and WM, is seen to lie in a range that would allow its application in studies involving larger cohorts of subjects, since the CV is well below 10% for metrics v in and v r in all ROIs and for ODI in WM. The higher CV in GM for ODI, driven by biological differences between subjects (ICC of 0.86), implies though that NODDI may require larger sample sizes for studies where GM neurite orientation dispersion is key to the findings. In conclusion, NODDI is a feasible alternative to DTI for in vivo spinal cord imaging, and has major advantages in measuring indices specific to neurite morphology able to disentangle the main factors contributing to DTI-derived anisotropy. Future work is warranted to optimise the technique to account for spinal cord anatomy. For instance, single fibre direction optimisation of the acquisition may lead to a shorter diffusion-weighting protocol. Furthermore, the histological correlates of NODDI indices should be investigated to confirm the specificity of the metrics and the model assumptions in the presence of pathology, by means of analysis of DW MRI scans of ex vivo spinal cord samples.
(EPSRC) (grants EP/I027084/01, EP/G007748, EP/L022680/1). Also, the group receives generous support from the Department of Health's National Institute for Health Research Biomedical Research Centres (BRC), which funds the infrastructure (Capital Project Number R&D03/ 10/RAG0449). The authors are thankful to Dr Daniel R. Altmann for statistical advice and to the volunteers who took part in this study.