A framework for optimal whole-sample histological quantification of neurite orientation dispersion in the human spinal cord

BACKGROUND
The complexity of fibre distributions in tissues is an important microstructural feature, now measurable in vivo by magnetic resonance imaging (MRI) through orientation dispersion (OD) indices. OD metrics have gained popularity for the characterisation of neurite morphology, but they still lack systematic validation. This paper demonstrates a framework for whole-sample histological quantification of OD in spinal cord specimens, potentially useful for validating MRI-derived OD estimates.


NEW METHOD
Our methodological framework is based on (i) sagittal sectioning; (ii) Palmgren's silver staining; (iii) structure tensor (ST) analysis; (iv) directional statistics. Novel elements are the data-driven optimisation of the spatial scale of ST analysis, and a new multivariate, weighted directional statistical approach for anisotropy-informed quantification of OD.


RESULTS
Palmgren's silver staining of sagittal spinal cord sections provides robust visualisation of neuronal elements, enabling OD quantification. The choice of spatial scale of ST analysis influences OD values, and weighted directional statistics provide OD maps with high contrast-to-noise. Segmentation of neurites prior to OD quantification is recommended.


COMPARISON WITH EXISTING METHODS
Our framework can potentially provide OD even in demyelinating diseases, where myelin-based histology is not suitable. As compared to conventional univariate approaches, our multivariate weighted directional statistics improve the contrast-to-noise of OD maps and more accurately describe the distribution of ST metrics.


CONCLUSIONS
Our framework enables practical whole-specimen characterisation of OD in the spinal cord. We recommend tuning the scale of ST analysis for optimal OD quantification, as well as neurite segmentation and weighted directional statistics, of which examples are provided herein.


Introduction
Microstructural magnetic resonance imaging (MRI) methods are a novel, emerging family of quantitative techniques that enable the non-invasive and clinically feasible mapping of tissue microstructure. Examples of such microstructural methods are diffusion MRI (dMRI) (Assaf et al., 2013;Le Bihan, 2003), quantitative susceptibility mapping (qSM) (Haacke et al., 2015) and quantitative magnetisation transfer (qMT) (Henkelman et al., 1993). dMRI exploits the ceaseless Brownian motion of water molecules in tissues as a unique probe of cell density and morphology. qSM methods estimate the local magnetic susceptibility from the MRI signal phase, which is related to iron concentration in tissues as well as to the local fibre orientation (Liu, 2010;Wharton and Bowtell, 2012). qMT provides indices of macromolecular content, specific to myelin in the central nervous system and also influenced by the primary fibre orientations in white matter (Pampel et al., 2015). Although different microstructural methods rely on different biophysical sources of contrast, their common goal is the estimation of features at a scale smaller than that of the image resolution via mathematical modelling of the MRI signal. This has the potential of revealing exquisite information about microstructural damage in early stages of severe conditions (Tur et al., 2016). Interestingly, the strengths and the complementary insights that each microstructural method provide can also be combined in more general multi-modal approaches, as for example in myelin g-ratio mapping (Stikov et al., 2015).
A number of microstructural imaging methods estimate properties related to the local fibre orientation (Descoteaux et al., 2011;Jeurissen et al., 2011;Liu, 2010;Pampel et al., 2015;Tournier et al., 2007;Wharton and Bowtell, 2012). This enables the mapping of features such as the dispersion of axon and dendrite orientations at the MRI volume element (voxel) scale (Kaden et al., 2016;Sotiropoulos et al., 2012;Tariq et al., 2016;Zhang et al., 2012). Neurite orientation dispersion (OD) is a relevant characteristic even in white matter areas where axons are well aligned among each other, and has been shown to be affected in neurological conditions, as in the case of multiple sclerosis (MS) (Schneider et al., 2014). Fibre OD in coherent white matter is believed to be a consequence of the underlying fine foliate structure (Budde and Annese, 2013), while regional variation in cortical grey matter reveals differences in terms of laminar organisation (Budde and Annese, 2013). The characterisation of neuronal fibre OD in the healthy and diseased brain may cast new light on variations of brain complexity due to pathology, maturation (Chang et al., 2015) and normal aging (Billiet et al., 2015). For instance, it is known that variations of dendritic complexity are per se potentially capable of causing cognitive deterioration, beyond neuronal cell loss due to neurodegeneration (Dickstein et al., 2007). Furthermore, OD mapping is also appealing in other areas of the nervous system, such as the spinal cord (Grussu et al., 2015), or even in other tissues whose fibrous structure can be investigated with MRI, as for example skeletal muscles (Damon et al., 2002), tendons (Fechete et al., 2005), heart (Hsu et al., 1998;Lombaert et al., 2012) and prostate (Bourne et al., 2012;Panagiotaki et al., 2015).
Histological validation is essential in microstructural imaging research (Jespersen et al., 2010;Jespersen et al., 2012;Wang et al., 2014), in order to confirm the specificity of the metrics and to rule out the influence of confounding factors not accounted for by the models. Recently, structure tensor (ST) analysis (Bigun and Granlund, 1987) has become popular as a robust algorithm for the estimation of fibre orientations from optical histological images, thus providing foundation for the validation of OD indices from MRI. Several studies have demonstrated that neurite orientations obtained from dMRI faithfully replicate ST-derived measures. For instance, 2D ST analysis of the rat brain (Budde and Frank, 2012) and of the human cortex (Seehaus et al., 2015) has shown good correspondence between ST estimates of neurite orientations and dMRI indices of dominant diffusion direction. Recently, 3D ST analysis was employed to relate confocal microscopy data of the rhesus macaque (Khan et al., 2015) and of the squirrel monkey (Schilling et al., 2016) brains to dMRI, with encouraging results.
In this work, we demonstrate a methodological framework for whole-sample histological quantification of OD in the human spinal cord. The procedure provides histology-based measures of OD for the validation of dMRI methods. It consists of: i) histological processing and Palmgren's silver staining; ii) optical imaging; iii) ST analysis; iv) directional statistical analysis for the quantification of neurite OD from ST orientations. We demonstrate the framework on two specimens of post mortem non-pathological human spinal cord, obtaining maps of OD. Moreover, we study in detail the implementation of the ST calculation and the directional statistical analysis. Specifically, we quantify how the spatial scale of ST analysis affects the estimation of OD, effectively demonstrating a data-driven approach for ST analysis optimisation. Also, we assess the potential utility of supporting ST-based OD quantification with neurite segmentation, and compare the performance of unexplored weighted directional statistical models and conventional univariate directional approaches for the evaluation of OD indices.

Materials and methods
This section describes the implemented methodological framework and the tissue specimens employed for its demonstration. All experimental sessions were approved by a local research ethics committee and followed appropriate consent as per Human Tissue Authority guidelines.

Tissue specimens
For the demonstration of our methodological framework, we studied two specimens of formalin-fixed, non-pathological human spinal cord, named samples 1 and 2, derived from the Oxford Brain Bank. Sample 1 measured 3.1 mm and was obtained from a 66 years old male subject at upper thoracic level. Sample 2 measured 2.1 mm and was obtained from a 67 years old female subject at upper lumbar level.

Description of the methodological framework
Our methodological framework consists of a number of steps, which are illustrated in Fig. 1 and summarised below.
1. The first step involves histological procedures to obtain sections of histological material. The sections are obtained with unconventional sagittal cuts, which enable the visualisation of the directions along which fibres run. The histological material is subsequently stained with the Palmgren's silver staining method (Palmgren, 1960), demonstrating neuronal and axonal elements. 2. The second step is optical imaging. 3. The third step is the estimation of axon and dendrite orientations from optical images with ST analysis. This estimation relies on the integration of neighbourhood image information. The scale at which information is integrated needs to be optimised. 4. The fourth step is the calculation of the dispersion of ST orientations using the tools of directional statistics.
For the practical demonstration of the framework as presented in this work, we proceeded as follows.

Histological procedures
Histology in spinal cord tissue is routinely performed axially. Here, we performed cuts sagittally in order to demonstrate the directions along which neuronal processes run. Specifically, 10 m-thick formalin-fixed paraffin embedded sections were used. Sections were sliced with a microtome from a surface exposed by a preliminary mid-sagittal cut. Subsequently, the derived histological material was impregnated with the Palmgren's silver staining method (Palmgren, 1960), optimised for optical imaging (DeLuca et al., 2004;DeLuca et al., 2013).

Optical imaging
Two stained sections underwent optical imaging (one for each specimen). High resolution optical images were acquired in digital format with an Aperio slide scanner (ScanScope AT Turbo, 400 × magnification). Digital images were analysed at a pixel resolution of 1 m × 1 m.

ST analysis
We performed ST analysis on the grey-scale converted optical images with the aim of estimating axon and dendrite orientations, similarly to other studies (Budde and Annese, 2013).
The ST of an image at each pixel can be evaluated directly from the image intensity gradient. It contains information about the local dominant orientation and about the strength of the underlying local structure. Specifically, in this work, we derived the conventional ST pixel-wise maps of local dominant orientation Â (x, y) and anisotropy index ˛(x, y) at each image pixel (x, y), as explained in detail in Supplementary material S1. The index ˛ at pixel (x, y) ranges in [0; 1] and quantifies the certainty of measuring the orientation Â at the same pixel location, given the local image structure (Bigun and Granlund, 1987).
For the practical evaluation of the ST, we employed in-house MATLAB code (The MathWorks, Inc., Natick, Massachusetts, USA) that has been freely released in GitHub (github.com/fragrussu/ StructureTensorToolbox). The code relies on the convolution of the image intensity with Gaussian filters, as this is a standard approach (Budde and Annese, 2013). The spread of the filters was controlled by a parameter here referred to as .
Since the optimal value of is unknown, we sampled a range of potential values for this parameter, performing ST analysis for 19 unique values of . These were the 15 evenly space values in the range [0.1; 6.0] m, as well as ∈ 7, 9, 12, 15 m. Image areas characterised by the presence of large vessels and high concentration of neuronal cell bodies were manually masked out for further processing. Moreover, we produced the conventional hue-saturation-value (HSV) representation of the ST (Budde and Annese, 2013), encoding Â in the hue channel, ˛ in the saturation channel and the image intensity in the value channel.

Directional statistics
We studied pixel-wise ST metrics with the aim of estimating neurite OD. This was achieved quantifying the spread of Â (x, y) within image patches of size equal to that of a plausible post mortem dMRI voxel (200 m × 200 m), using directional statistics.
The distribution of Â (x, y) in each patch can be modelled according to univariate probability density functions defined over the unit circle, as for example the Von Mises distribution (Budde and Frank, 2012). While this approach is common and capable of providing useful OD metrics, it has the main limitation of treating equally all orientations, including those coming from areas with weak local structure. The ST always provides a dominant orientation even in absence of underlying structure, as its eigenvalues would most likely never equal each other due to noise. As a consequence of this, the indiscriminate inclusion of meaningless random orientations from areas with weak local structure (i.e. low values of the anisotropy index ˛) may affect the estimation of ST-derived OD.
To address this limitation, approaches that model jointly the distribution of the local dominant orientation as well as its corresponding anisotropy may be adopted. This alternative methods would benefit from the additional layer of information provided by the anisotropy index ˛, allowing the evaluation of indices of spread of angular random variables that account for variable reliability of the observation of the angular variables themselves (Brunsdon and Charlton, 2006).
Such approaches would introduce weighted directional distributions. Weighting of directional probability densities has been used, for instance, to describe astrophysical data, as in the case of the study of cosmic rays (Linsley, 1975;Mardia and Edwards, 1982). The application of weighted directional statistics for OD quantification has not been explored yet, and here we hypothesise that they could provide a more robust description of the underlying STderived metrics, as compared to conventional univariate angular distributions.
In this work, we propose to employ weighted directional statistics to quantify the OD of ST-derived orientation Â (x, y) within image patches. Building upon (Bigun and Granlund, 1987), we employ the anisotropy index ˛(x, y) as a measure of reliability of Â (x, y), and handle both as random variables. We move from conventional univariate directional statistical models, and employ the measure of reliability ˛ as a weighting factor.
In this work, we map OD according to three of the most common directional statistical models: the Watson, the Von Mises and the Gaussian distributions (Mardia and Jupp, 2009). Moreover, we employ their weighted counterparts, which rely on the anisotropy index ˛ as a measure of certainty of measuring Â (the weighted-Watson, weighted-Von Mises and weighted-Gaussian distributions), as our anisotropy weighting approach is general and can be potentially applied to any probability density of Â. With this rich set of distributions, we evaluate on the one hand the performance of the most common directional models for OD mapping, while on the other hand we test generally the benefits of our innovative, anisotropy-weighted approach as compared to conventional univariate models.
The analytical expressions of all six distributions, as well as the expressions of their maximum likelihood parameter estimates and the values of the maximised (log-)likelihood have been reported in Supplementary material S2.

Experiments
The following experiments were performed to evaluate: i) the effect of the choice of the spatial scale of ST analysis; ii) the impact of supporting OD quantification with neurite segmentation; iii) the performance of weighted directional statistical models, as opposed to conventional univariate approaches.
In practice, we adopted the circular variance (CV) as a metric of OD for the Watson and Von Mises models, as well as for their weighted counterparts (expressions of CV reported in Supplementary material S2). CV is a common index quantifying the spread of angular variables (Mardia and Jupp, 2009). It ranges from 0 to 1 and increasing values imply increasing OD. For the Gaussian and weighted-Gaussian models instead, the metric of OD was the spread of the distribution (SD, corresponding to parameter in Supplementary material S2). SD has in practice been used in other studies (Budde and Frank, 2012), and increasing SD implies increasing OD. For all models, we studied the values of CV/SD in grey and white matter and the grey/white matter contrast-to-noise ratio (CNR) (Grussu et al., 2015) of CV/SD, as well as the values of the log-likelihood functions at maximum, as the scale parameter varies.
For each value of and for each sample, we followed the following steps.

1) We evaluated patch-wise OD (CV for Watson/weighted-Watson, Von Mises/weighted-Von Mises distributions; SD for
Gaussian/weighted-Gaussian distributions) and the value of the log-likelihood at maximum (log L max ) for all six statistical models. Patches of 200 m × 200 m were employed. 2) We calculated the CNR of CV/SD between grey and white matter as (1) In Eq.
(1),m GM/WM and s 2 GM/WM are respectively the mean and the variance of CV/SD within a region of interest (ROI) in grey/white matter. The ROIs were drawn manually on the histological images downsampled at the resolution of the OD (CV or SD) map. The ROIs had the same size for grey and white matter.
3) We evaluated patch-wise distributions of Â and ˛ via calculation of their normalised histograms (128 bins). 4) We ranked the six models (Watson, weighted-Watson, Von Mises, weighted-Von Mises, Gaussian, weighted-Gaussian) in each patch depending on how well they describe the underlying distribution of ST metrics ˛ and Â. For this purpose, we sorted patch-by-patch the models depending on the value of the loglikelihood function at maximum (log L max ), such that the higher log L max , the higher the position in the ranking.
In order to evaluate the effect of including non-tissue pixels for directional statistics, steps 1-4 were also repeated considering only pixels within neuronal elements, while excluding the extra-neurite space. This was done after segmenting neuronal elements on the Palmgren's silver images (k-means algorithm with k-means++ initialisation (Arthur and Vassilvitskii, 2007), 4 clusters).
All image processing was carried out using in-house MATLAB functions and scripts, freely released under GitHub (github.com/ fragrussu/StructureTensorToolbox). The code was run on a 16-core Linux machine (64 GB memory, 3 GHz processor speed).

Data acquisition: histological procedures and optical imaging
Fig. 2 shows digital images obtained from two sagittal sections, which were impregnated with the Palmgren's silver method. The figure also shows magnified details of size equal to the image patches used for quantitative analysis, from both grey and white matter. In both spinal cord samples, sagittal sectioning followed by Palmgren's silver staining enables the effective visualisation of neuronal and axonal processes (i.e. axons in white matter and dendrites in grey matter), which are stained in dark brown and black onto a lighter background.  ment between ST orientations (mapped onto the hue channel of the images) and the dominant neurite orientation. Nevertheless, as becomes larger and larger, bigger parts of the patch show spurious orientations (for example: blue-violet bits, as compared to red, apparent for = 5.549 m in white matter). The plotted patches show that OD occurs in white matter due to phenomena such as axon undulation, and also highlight the complexity of dendritic arborisations in grey matter.

ST analysis
For very low , the HSV encoding is well saturated, while higher provides less saturated images. Also, very low values such as = 0.100 m provide grainy orientation maps, which instead become smoother for higher values of . Nevertheless, as high as about 5 m and more leads to the presence of blobs characterised by orientations that appear to differ from the local neurite orientation, and whose size increases as increases. Similar trends were also observed for sample 1.

Experiments
3.3.1. OD patch-wise map Fig. 4 shows the patch-wise OD maps (CV for Watson, Von Mises, weighted-Watson, weighted-Von Mises models; SD for the Gaussian and weighted-Gaussian models) for sample 2 as varies. Trends shown in the figure agree well with those measured in sample 1. The superior panel illustrates maps obtained when all pixels within a patch are included for the calculation of CV/SD. The maps showed in the inferior panel, instead, were calculated including only pixels segmented as part of a neurite.
For all statistical models and both post processing approaches, the values of CV and SD vary strongly as varies. For instance, in white matter, all models provide values of OD that decrease rapidly as increases from 0.100 m up to about 2.207 m, and then increase again as increases further. Both CV and SD show high grey/white matter contrast for the smallest value of (0.100 m). The contrast decreases slightly as increases further, and is substantially degraded as approaches the highest values sampled in this study. The OD maps obtained after neurite segmentation resemble those obtained when all pixels in a patch are considered for analysis, although values of CV and SD are qualitatively higher in the latter case, as compared to the former. Lastly, we observe that the weighted distributions provide OD maps that appear smoother than those from the conventional, univariate counterparts. Fig. 5 shows the median value of the OD metrics CV and SD within grey and white matter for the two samples. CV/SD are shown in two cases: when all pixels in a patch (circles) and when only pixels belonging to neurites (crosses) are considered for OD mapping. For all spatial scales, all statistical models and all tissue types (grey/white matter), values of CV/SD obtained including all pixels in a patch are higher than values obtained considering only pixels within neurites. This effect is more pronounced for low to intermediate spatial scales, corresponding to values of up to roughly 5 m. Fig. 6 shows the distribution of ST orientation Â, as well as the corresponding maximum-likelihood Watson/weighted-Watson, Von Mises/weighted-Von Mises, Gaussian/weighted-Gaussian fits, for a white matter patch from sample 2. The distributions are reported for different values of .

Distributions of ST orientation
Overall, ST analysis identifies the dominant, underlying neurite orientation in the image plane. Small, intermediate and high values of Â appear to capture similar dominant orientations. When information from all pixels is used and is as low as 0.100 m, the measured distribution of ST orientations shows several spurious peaks, that cannot be captured by the fits. As grows, the distributions of ST Â get smoother, and all models capture the dominant direction. Finally, for growing even more and surpassing a value of about 5 m, the distribution of Â becomes again less smooth and characterised by a number of spurious peaks. Lastly, we note that excluding non-neurite pixels provides in general similar trends as considering the whole image patch for analysis, although fewer spurious peaks are seen at low . Fig. 7 plots representative distributions of ST anisotropy index ˛ (normalised histogram, in grey), as well as the maximum-likelihood weighted-Watson, weighted-Von Mises and weighted-Gaussian fits (in blue, violet and cyan respectively, as marginal distributions of ˛), for five different choices of . Notice that no information about the conventional Watson, Von Mises and Gaussian models is displayed, since they do not model ˛.

Distributions of ST anisotropy
In both post-processing approaches (inclusion or not of nonneurite pixels), the distributions show that the probability of observing a certain value of ˛ increases as ˛ increases for low and For higher values of (on the order of the tens of m), the distribution of ˛ becomes more complex and multimodal, although a trend for increasing probabilities densities as ˛ increases is still seen. The weighted distributions capture the monotonic increase of the normalised histogram for intermediate values of (for instance, 0.943 m and 2.207 m), but cannot model multimodal peaks as those observed at = 12 m. Lastly, we report that for very low (0.100 m), the distribution of r esembles a Dirac-delta concentrated in ˛ = 1, whose shape cannot be fully captured by the weighted-Watson, weighted-Von Mises and weighted-Gaussian fits. Fig. 8 shows the grey/white matter CNR of the OD metrics CV/SD as a function of . The CNR was evaluated when all pixels in a patch are included for the calculation of CV and also when neurite pixels only are considered. The CNR is relative to CV for Watson/weighted-Watson and Von Mises/weighted-Von Mises models, while it is relative to SD for the Gaussian and weighted-Gaussian models.

Grey/white matter CNR of OD
In all cases, CNR is maximum when is minimum (0.100 m). As increases, CNR decreases, reaching a plateaux for ranging approximately in [1; 3] m, and then monotonically decreases, degrading to values smaller than 1 for of the order of 10 m. The inclusion of neurite-only pixels for the calculation of CV/SD does not enhance systematically the CNR, as compared to the case when all pixels in a patch are included. Furthermore, the plots demonstrate that the CNR of CV/SD is always higher when the metric is evaluated according to the anisotropy-weighted distributions, as compared to the conventional univariate distributions (Watson vs weighted-Watson; Von Mises vs weighted-Von Mises; Gaussian vs weighted-Gaussian). Fig. 9 plots the patch-wise ranking of the directional models employed to map OD. The figure shows results for sample 2, which are consistent with those from sample 1. Also, the ranking is displayed for the two post-processing cases (inclusion or not of background pixels from the extra-neurite space), which provide similar trends.

Ranking of directional models
It can be seen that for low to intermediate spatial scales the weighted distributions always rank higher than their univariate counterparts in the vast majority of patches. For high spatial scales instead (as =12 m), the weighted and conventional distributions show similar quality of fit, as they share a similar number of patches where one is better than the other.
Lastly, the figure also demonstrates that the Gaussian and Von Mises distributions provide better quality of fit as compared to the Watson distribution.

Discussion
This paper demonstrates a methodological framework for whole-sample quantification of neurite OD in post mortem specimens of human spinal cord. The framework consists of histological procedures (dehydration, paraffin-embedding, unconventional sagittal sectioning on a microtome), Palmgren's silver staining, optical imaging, ST analysis of the optical images in digital format and directional statistics of ST metrics. Two specimens of post mortem spinal cord tissue were employed for the demonstration of the framework. The acquired digital images were also analysed in detail to quantify: i) the effect of the choice of the spatial scale of ST analysis; ii) the impact of supporting OD quantification with neurite segmentation; iii) the performance of weighted directional statistical models to conventional univariate approaches for the evaluation of OD.

Feasibility of the framework
We demonstrated that sagittal sectioning and Palmgren's silver staining are excellent methods for whole-sample ST-based quantification of neurite OD in the spinal cord. Palmgren's silver staining was previously proven to demonstrate neuronal and axonal elements throughout the length of the cord more consistently than neurofilament immunostains, which are susceptible to suboptimal antigenicity due to formalin fixation and influenced by neuronal/axonal phosphorylation status (DeLuca et al., 2004). Moreover, unlike myelin-based methods such as myelin-specific Fig. 6. Effect of the choice of on the distribution of ST orientation Â within a representative white matter patch from sample 2. In panel A, top-left, the patch is pictured. In panel B, to the right, the distributions of ST orientations obtained at various are illustrated ( increases from bottom to top; notice that only one of the two antipodal peaks has been shown for the Watson/weighted-Watson distributions), as well as the fits from the six distributions (marginal distribution in Â is shown for the anisotropy-weighted models). Panel B has two columns: plots in the left column are obtained using orientations from all pixels within a patch; plots in the right column are obtained including orientations from neurite pixels only. In panel C, bottom-left, the polar histograms of ST orientation Â is also reported. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 7. Effect of the choice of on the distribution of ST anisotropy index ˛ within the same patch pictured in Fig. 6. The normalised histogram of ˛ is displayed using the black colour, whilst the maximum-likelihood fit of the weighted-Watson, weighted-Von Mises and weighted-Gaussian distributions (marginal distributions of ˛) are respectively shown in blue, violet and cyan. Notice that no fits for the conventional univariate models are plotted, since they do not account for ˛. In the superior row, distributions are obtained using values of ˛ from all pixels in a patch. In the inferior row, only values of ˛ from neurite pixels are used. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) immunohistochemistry or polarised light imaging (Axer et al., 2011), it can be applied reliably in demyelinating diseases (DeLuca et al., 2004;DeLuca et al., 2013).
From high-resolution digital images of the silver stained sections, the local neurite orientation can be estimated with ST analysis, a popular gradient-based method. ST provides pixel-wise maps of local orientation and anisotropy, from which neurite OD indices can be readily calculated using the tools of directional statistics, as also reported in other studies (Budde and Annese, 2013;Budde and Frank, 2012). This provides potential for the systematic histological validation of dMRI OD metrics, as for example NODDI ODI, a potential biomarker in the brain (Zhang et al., 2012) and in the spinal cord (Grussu et al., 2015).

Choice of the spatial scale of ST analysis
ST analysis was performed and characterised for various choices of the parameter .
Visual inspection of the conventional HSV encoding of the ST in both grey and white matter demonstrates that the technique recovers the underlying dominant neurite orientation at all scales. However, different values of provide tensors and tensor-derived statistics that also differ notably from each other.
Very low provides grainy HSV encodings of the ST, characterised by high anisotropy, as suggested by the high saturation of the HSV images. This is due to the presence of spurious peaks in ST orientation distributions evaluated at the patch level, which are not described well by smooth distributions such as Watson/weighted-Watson, Von Mises/weighted-Von Mises, Gaussian/weighted-Gaussian. Nevertheless, low allows the detection of differences in terms of neurite morphology between grey and white matter, since CNR of OD metrics as high as 6 in sample 2 was measured for = 0.100 m. We speculate that very low values of may define a spatial scale that is even smaller than the characteristic length of the morphological features whose orientation needs to be estimated. In the limit case of → 0 the Gaussian kernels employed in ST analysis collapse to a Dirac delta ı (x, y) (Supplementary material S1). In such a limit condition, the local orientation is simply that of the direction orthogonal to the image gradient, whose calculation is notoriously prone to noise amplification (Jähne, 2005), possibly explaining the high number of spurious orientations and the high values of the OD metrics observed at very low .
As grows up to roughly 3 m, smoother HSV encoding of the ST is produced in both grey and white matter. This fact is associated with the observation of smooth distributions of orientation and anisotropy, which are captured by the Watson/weighted-Watson, Von Mises/weighted-Von Mises, Gaussian/weighted-Gaussian models. The grey/white matter CNR of CV/SD decreases slightly, but is well above 1 in both specimens.
As increases even more, the patch-wise distributions of ST metrics become again less smooth, due to the re-emergence of spurious peaks in the orientation distribution, and the directional models studied in this paper provide poorer descriptions of ST metrics. The grey/white matter CNR of CV/SD reduces to values well below 1 for both statistical models, as in the case of as high as about 5 m in sample 1 and as high as 10 m in sample 2. This suggests that when the scale at which image information is integrated is too coarse, differences in terms of neurite orientation between the two tissue types can be blurred and lost.
Finally, we note that the values of CV/SD depend themselves on the choice of . This is not surprising, since pre-convolution for the calculation of the ST (Budde and Annese, 2013), which has the desirable property of providing robustness to noise, has been demonstrated analytically to alter the anisotropy of the ST . Nevertheless, the approach proposed in this work reli-ably detects differences in term OD between tissue types, and can be employed to study regional or group-level variations of OD, as for the validation of dMRI. In future, more sophisticated noise reduction methods could be explored, and the impact of the choice of the denoising algorithm on the accuracy of CV/SD could be evaluated.

Role of neurite segmentation
Our results demonstrate that background noise, i.e. voxels not including neurites, could affect some of the metrics. Trends in terms of CV, SD, CV/SD CNR, log-likelihood and distribution of orientation/anisotropy as a function of obtained after neurite segmentation are indeed similar and consistent with those obtained when the segmentation was not performed. Nevertheless, the values of OD obtained without neurite segmentation are always higher than values obtained when OD is calculated from neuritesegmented pixels only. Further exploratory tests on numerical phantoms confirm that the higher CV/SD is due to increased variability that follows the inclusion of random orientations from background pixels of non-stained areas (Supplementary material S3). We remark that in ST analysis it is common to consider images as a continuous medium , without pre/postprocessing steps aiming to identify cellular or neuronal structures. Although this approach seems valid for capturing differences in neurite orientation complexity in healthy tissue, dependence of the actual OD values on the amount of non-stained background should be handled with care. In particular, increased betweenaxonal space due to axonal loss, known to characterise diseases such as MS (DeLuca et al., 2004), may cause artefactual increases of OD values in pathological tissue as compared to non-pathological areas. Hence, careful design of further pre/post-processing steps for ST analysis in pathology appear desirable, especially in studies focussing on the effect of pathology on OD.

Comparison between directional models
In this work we have compared six different models describing the patch-wise distribution of ST metrics. These models belong to two groups: the first was a set of conventional probability density functions often employed in literature to describe fibre orientation distributions, i.e. the Watson distribution (Mardia and Jupp, 2009;Zhang et al., 2012), the Von Mises distribution (Budde and Frank, 2012) and the Gaussian distribution (Budde and Annese, 2013). The second group was instead obtained adding anisotropy-informed weighting to the previous three models, providing the weighted-Watson, weighted-Von Mises and weighted-Gaussian models. This latter set of models describes ST anisotropy ˛ and ST orientation Â jointly, in a bivariate fashion, with the former metric effectively acting as a measure of certainty of Â. In the conventional univariate framework, the spread of the sample of orientations can be seen as an overall measures of uncertainty at the image patch level. On the other hand, with the innovative anisotropy-weighted approach, a further level of uncertainty is modelled. Due to its dependency on the anisotropy index ˛, the weighted distributions can also capture uncertainty at the level of the individual pixel, down-weighting unreliable observations of the orientation Â coming from areas with no structure.
We compared all models studying the grey/white matter CNR of their patch-wise OD maps (CV for Watson, Von Mises and weighted counterparts; SD for Gaussian and weighted Gaussian). We also ranked the models depending on the values of their corresponding log-likelihood functions at maximum, indicative of how well they describe patch-wise samples of ˛ and Â. Our results demonstrate that the anisotropy-weighted approach performs better than the conventional univariate approach for small to intermediate spatial scales. In those cases, the weighted version of each model (i.e. weighted-Watson, weighted-Von Mises, weighted-Gaussian) is ranked always at higher position than its non-weighted counterpart (Watson, Von Mises and Gaussian) in the vast majority of voxels. This also corresponds to an increase of the grey/white matter CNR of the OD map, more evident for the Watson and Von Mises models, thus enabling the detection of potentially subtler variations of OD of one of the two tissue types due to pathology. Our experiments also demonstrate that at high spatial scale (i.e. as big as 12 m and beyond), the weighted models perform similarly to their non-weighted, conventional counterparts. However, it has already pointed out that such scales are not likely to be used in practice, as they are too coarse for the type of data studied in this paper. In practice, at those scale, the differences between grey and white matter are essentially blurred and lost, as also suggested by the very low CNR shown by all models.
Finally, our model ranking shows that the Von Mises and Gaussian distributions rank at higher position than the Watson distribution. This is due to the fact that the two distributions exhibit a single peak, while the Watson distribution has antipodal symmetry and is hence characterised by two peaks separated by (180 • ), i.e. any two orientations Â and Â ± are modelled with the same probability density value. In this work, we forced angles to be defined within the range [0; ] (i.e. in each patch there is a single peak), as two fibre orientations separated by are considered the same. Nevertheless, it should be noted that the apparent poorer performance of the Watson distribution comes with a much higher flexibility as compared to the Von Mises and Gaussian distributions. The Watson model can handle any angle without the need of implementing any check or rephasing prior to directional statistics, which is instead essential for a correct employment of the Von Mises and especially of the Gaussian models.

Potential extensions of the anisotropy-weighted approach
In this paper, we have analysed 2D optical images. As a result of this, one angular variable sufficed to parametrise neurite orientations in the image plane. Nevertheless, the directional statistical models studied in this work could be generalised to more dimensions. For example, 3D ST analysis could provide pixel-wise azimuthal and polar angles that could be modelled with a Watson/Bingham distribution defined over the unit sphere (Mardia and Jupp, 2009;Tariq et al., 2016;Zhang et al., 2012), as well as anisotropy indices (Khan et al., 2015) that could be employed as weighting terms in a trivariate weighted-Watson/Bingham formalism.

Potential applications of the framework
Our framework offers a number of potential applications. Firstly, it could be employed alone for the investigation of the alteration of neurite orientation distribution complexity in neurological conditions. Histopathology studies of the spinal cord in diseases such as MS or spinal cord injury have historically focussed mainly on structural/axonal integrity and on conventional immunohistochemical labelling (DeLuca et al., 2004;Noristani et al., 2015). Our methodological framework offers new opportunities, such as the possibility of investigating how pathology alters the complexity of dendritic trees. Alterations in dendrite OD could potentially provide new indices that closely reflect impaired transmission of sensory or motor information, as OD has been related directly to functional connectivity in the nervous system (Nazeri et al., 2015).
Secondly, the framework could also be a valid tool for validating MRI-based OD mapping in the healthy and diseased spinal cord (Grussu et al., 2015), as well as to inform MRI signal models. For example, application to ex vivo specimens may allow the establishment of direct links between the diffusion MRI signal and the underlying fibre orientation distribution and dispersion, which could be employed for histology-informed MRI protocol optimisation. Also, the framework would enable the practical validation of recent finding suggesting that OD is affected in diseases such as MS (Schneider et al., 2014). For instance, the in-plane dispersion as measured by techniques such as NODDI (Zhang et al., 2012) could be compared directly to the dispersion indices obtained from the Watson, Von Mises or Gaussian models, as long as the location from which the MRI signal originates is identified accurately in the histological material.

Limitations
The framework shown in this article is inherently two dimensional, and relies on the assumption that sagittal sectioning is adequate for the complete characterisation of neurite orientations, thus neglecting contribution of through-plane fibres. Although this assumption may be reasonable in white matter areas with negligible amount of collateral axons, it may not hold true in grey matter or at the interface between white and grey matter, due to the presence of intricate dendritic trees. Nevertheless, the variation of neurite orientations in our 2D images was sufficient to characterise the different OD of grey and white matter, as proven by grey/white matter CNR of the metrics CV/SD well above 1. Moreover, the acquisition of 2D images for ST analysis, as opposed as 3D data (as for example confocal microscopy (Khan et al., 2015;Schilling et al., 2016)), has the advantage of keeping the information to be processed to amounts that enable whole-sample characterisation with commonly available computational resources. Our framework enables the practical characterisation of neurite OD in areas of several millimetres in size, as for example focal or widespread tissue damage. 3D approaches could provide complementary and potentially more specific information, but are likely to be employed to study much smaller portions of tissue.
Another limitation of the study is the specificity of the chosen Palmgren's silver method. While Palmgren's silver staining consistently demonstrates axons in black/dark brown, it also stains blood vessels and neuropil in lighter shades of brown. In this work, manual preprocessing was necessary to remove the biggest vessels and neuronal somas before quantitative analysis. Nevertheless, Palmgren's silver method is technically robust and overcomes the limitations imposed by formalin-fixed post mortem spinal cord material on neurofilament immunohistochemistry, namely compromised antigenicity and phosphorylation variability. Future improvement to the presented framework could be achieved introducing further pre/post-processing steps that automatic detect and remove such unspecific stained areas.
Finally, k-means clustering of the pixel RGB values was used to segment neurites with the aim of assessing how the inclusion of image areas that do not contain axons or dendrites affects the estimation of OD. Although such a segmentation could potentially be refined using morphological operators or shape descriptors, careful visual inspection suggests that it sufficed for the explorative purposes of this study.

Conclusions
In this paper, we demonstrate a methodological framework for whole-sample histological quantification of OD in the human spinal cord. We show that Palmgren's silver staining of spinal cord material sectioned sagittally provides unique visualisation of neuronal and axonal elements, enabling quantitative OD mapping via ST analysis.
Our experiments also suggest that segmentation of neurites prior to the quantification of OD from ST-derived orientations may be a useful step, especially for the characterisation of OD in pathology. Also, we believe that the spatial scale of ST analysis needs to be tuned for robust OD quantification. For this purpose, we have demonstrated that the between-tissue CNR of the OD indices is a useful metric that can suggest optimal choices of the spatial scale of analysis. In view of our data, we recommend spatial scales that do not exceed roughly 3 m, in terms of the parameter employed in this work.
Lastly, we also conclude that the quantification of neurite OD benefits from the information conveyed by ST anisotropy, which should be accounted for together with ST-derived local orientation. For this purpose, weighted directional statistical distributions such as the weighted-Watson, weighted-Von Mises or weighted-Gaussian models studied in this paper could be used.

Funding sources
This work was funded by the UCL Grand Challenges scheme (an initiative of the UCL School of Life and Medical Sciences, UCLH/UCL Biomedical Research Centre and the Specialist Biomedical Research Centres at Moorfields/UCL and Great Ormond Street/UCL) and by the H2020-EU.3.1 CDS-QUAMRI grant (ref: 634541), that support FG. TS is an employee of Philips HealthTech (Guildford, Surrey, UK). DCA receives funding from the Engineering and Physical Sciences Research Council (EPSRC, grants G007748, I027084 and M020533). CGWK and DCA receive funding from the H2020-EU.3.1 (CDS-QUAMRI grant, ref: 634541). CGWK also receives funding from the EPSRC (EP/I027084/1). GCD is supported by the NIHR BRC, Oxford. The NMR Research Unit is supported by a programme grant from the UK MS Society (grant 892/08) and receives generous support from the Department of Health's National Institute for Health Research Biomedical Research Centres (BRC, Capital Project Number R&D03/10/RAG0449). Funding sources had no involvement in the design of the study, in the collection and analysis of the data and in the interpretation of the results.

Conflicts of interest
Authors have no conflicts of interest to disclose with regard to the matter of the work.