Fingerprinting Orientation Distribution Functions in Diffusion MRI detects smaller crossing angles

TL;DR

This paper introduces ODF-Fingerprinting, a novel method that improves detection of small crossing angles in diffusion MRI, enhancing the accuracy of white matter tractography and brain connectivity analysis.

Contribution

The paper presents a fingerprinting-based approach that uses the entire ODF shape for better fiber direction detection, especially for small crossing angles, outperforming existing methods.

Findings

Improved detection of small crossing angles in noisy ODFs.

Enhanced tractography accuracy with better fiber bundle delineation.

Maintained fiber direction precision despite noise.

Abstract

Diffusion tractography is routinely used to study white matter architecture and brain connectivity in vivo. A key step for successful tractography of neuronal tracts is the correct identification of tract directions in each voxel. Here we propose a fingerprinting-based methodology to identify these fiber directions in Orientation Distribution Functions, dubbed ODF-Fingerprinting (ODF-FP). In ODF-FP, fiber configurations are selected based on the similarity between measured ODFs and elements in a pre-computed library. In noisy ODFs, the library matching algorithm penalizes the more complex fiber configurations. ODF simulations and analysis of bootstrapped partial and whole-brain in vivo datasets show that the ODF-FP approach improves the detection of fiber pairs with small crossing angles while maintaining fiber direction precision, which leads to better tractography results. Rather than focusing on the ODF maxima, the ODF-FP approach uses the whole ODF shape to infer fiber directions to improve the detection of fiber bundles with small crossing angle. The resulting fiber directions aid tractography algorithms in accurately displaying neuronal tracts and calculating brain connectivity.

Equations

$$S(\vec{q})=f_W{e}^{-\overrightarrow{q^2}{D}_W(\Delta -\frac{\delta }{3})}$$

$$+{\Sigma }_{j=1}^N{f}_j{e}^{-{\vec{q}}^T{\bar{D}}_j\vec{q}(\Delta -\frac{\delta }{3})}$$

$$l_m=ar{g}_{l}max(L_{ODF}\cdot OD{F}_m^T).$$

$$l_m=ar{g}_{l}max\left(log(L_{ODF}\cdot OD{F}_m^T)-\frac{n_{par}}{4n}{\sigma }_n\right)$$

Full text

In ODF-FP, fiber configurations are selected based on the similarity between measured ODFs and elements in a pre-computed library. In noisy ODFs, the library matching algorithm penalizes the more complex fiber configurations.

ODF simulations and analysis of bootstrapped partial and whole-brain in vivo datasets show that the ODF-FP approach improves the detection of fiber pairs with small crossing angles while maintaining fiber direction precision, which leads to better tractography results.

Rather than focusing on the ODF maxima, the ODF-FP approach uses the whole ODF shape to infer fiber directions to improve the detection of fiber bundles with small crossing angle. The resulting fiber directions aid tractography algorithms in accurately displaying neuronal tracts and calculating brain connectivity.

*Keywords: Diffusion MRI, Orientation Distribution Function, Fingerprinting, Fiber Tractography, Fiber Identification, Crossing angle, Radial Diffusion Spectrum Imaging, multi-shell Q-ball imaging

††journal: NeuroImageAbbreviations: ADC, Apparent Diffusion Coefficient; AIC, Akaike Information Criterion; CHARMED, Composite Hindered and Restricted Diffusion Model; CI, 95% Confidence Intervals; CSF, CerebroSpinal Fluid; dODF, diffusion Orientation Distribution Function; DSI, Diffusion Spectrum Imaging; DTI, Diffusion Tensor Imaging; DWI, Diffusion Weighted Imaging; EPI, Echo Planar Imaging; FA, Fractional Anisotropy; fODF, fiber Orientation Distribution Function; GQI, Generalized Q-Space Sampling; HARDI, High Angular Resolution Diffusion Imaging; HCP, Human Connectome Protocol; ODF, Orientation Distribution Function; ODF-FP, Orientation Distribution Function Fingerprinting; QA, Quantitative Anisotropy; qODF, q-ball Orientation Distribution Function; RDSI, Radial Diffusion Spectrum Imaging; ROI, Region of Interest; SNR, Signal to Noise Ratio;

1 Introduction

Diffusion weighted MRI (DWI, (Le Bihan et al., 1986)) non-invasively captures the complex microstructure of the brain. The angular dependence of DWIs sensitivity to water molecule motion powers techniques such as Diffusion Tensor Imaging (DTI, (Basser and Pierpaoli, 1996; Mori et al., 2002)), Q-ball imaging (Tuch et al., 2002; Tuch, 2004) and Diffusion Spectrum Imaging (DSI, (Wedeen et al., 2005, 2008; Baete et al., 2016; Baete and Boada, 2018)) to identify fiber bundle directions in each voxel. This fiber bundle information, collected throughout the whole brain, is the input for tractography algorithms which produce representations of long-range axonal structure (Pierpaoli et al., 1996; Basser et al., 2000; Bammer et al., 2005; Fernandez-Miranda et al., 2012; Shin et al., 2012). While the anatomical accuracy of these representations is a subject of much investigation (Knosche et al., 2015; Schilling et al., 2016), the ability to non-destructively obtain brain structural connectivity information has led to the adoption of tractography for use in neuroscience and clinical applications (Fernandez-Miranda et al., 2012; Shin et al., 2012; Jbabdi et al., 2015; Galantucci et al., 2016; Mitra et al., 2016).

High Angular Resolution Diffusion Imaging (HARDI) methods, such as multi-shell Q-ball and DSI, capture the complex intra-voxel crossings (Wedeen et al., 2012; Fernandez-Miranda et al., 2012) in Orientation Distribution Functions (ODFs). Direct calculation of these ODFs (diffusion ODF or dODF) requires a sufficiently dense sampling of diffusion weightings and directions as in Cartesian (Wedeen et al., 2005) and Radial DSI (Baete et al., 2016), typically distributed on several shells. From single shell acquisitions, q-ball ODFs (qODF) can be estimated using a spherical Radon transform (Tuch, 2004). Both these ODFs can be transformed to fiber ODFs (fODF) by spherical deconvolution with an estimated Fiber Response Function (Tournier et al., 2008; Jeurissen et al., 2014; Dhollander et al., 2016) - the prototypical expected response of a single fiber.

A key step in employing ODFs in tractography algorithms, whether they are dODFs, qODFs or fODFs, is the correct identification of fiber directions in each voxel. To this end, many algorithms have been proposed ranging from simple finite difference methods (Descoteaux et al., 2007; Frey et al., 2008) to more complex numerical optimization solutions combining gradient ascent (Berman et al., 2008) or Newton-Rhapson (Tournier et al., 2004) algorithms with appropriate thresholds on fiber proximity. Other methods transform ODFs to a constrained polynomial basis to aid the numerical identification of stationary points on the ODFs surface (Aganj et al., 2010). Yet another approach uses Bayesian estimation to fit models to the diffusion data in each voxel (probabilistic estimation, (Behrens et al., 2007)).

All these fiber direction identification approaches are however flawed due to the intrinsic ODF peak width ((Barnett, 2009; Jensen and Helpern, 2016), Fig. 1a). This ODF peak width (Barnett, 2009) and the limited angular resolution of the acquisition make it difficult to accurately estimate the directions of fibers crossing at shallow angles ((Kuo et al., 2008; Jeurissen et al., 2013), Fig. 1b). Most proposed methods indeed fail to detect crossing angles less than 40∘(Kuo et al., 2008; Jeurissen et al., 2013; Tournier et al., 2008; Descoteaux et al., 2007; Daducci et al., 2014; Wilkins et al., 2015). Even deconvolving the ODFs with a Fiber Response Function fails to reliably detect crossing angles smaller than 30∘(Tournier et al., 2008; Jeurissen et al., 2014). Failing to identify all bundle directions prohibits tractography algorithms from correctly following fibers in areas with crossing bundles.

By focusing on the maxima of the ODF to identify fiber directions, the methods listed above ignore the information captured in the shape of the ODF. For example, two fiber bundles crossing at an angle smaller than the intrinsic ODF peak width of the reconstruction give rise to one single ODF maximum. The flattened shape of the peak, however, reveals the presence of two fibers (Fig. 1b).

Here we propose a new approach to fiber bundle identification inspired by key concepts first introduced in MR Fingerprinting (Ma et al., 2013; Cloos et al., 2016). Instead of a dictionary with spin evolutions at different $T_{1}$ and $T_{2}$ relaxation times, we generate a library of ODF-fingerprints and identify the fiber directions of ODFs by assessing the similarity between the measured data and the elements in our library (Fig 1e). This ODF-Fingerprinting (ODF-FP) approach utilizes the whole ODF shape to infer fiber directions rather than more narrowly focusing on the ODF maxima.

In this work, we show that the ODF-FP method not only identifies smaller crossing angles more accurately, it also improves the performance at larger crossing angles. To this end, we discuss the different facets of the ODF-FP algorithm for fiber direction identification and look at the methods’ performance in a simulated phantom containing crossing fiber bundles, at angular precision in individual simulated crossing fiber ODFs, at the performance in a multi-resolution Human Connectome Protocol (HCP) dataset, at reproducibility and noise sensitivity in in vivo bootstrapped datasets and at the impact on tractography results in in vivo whole brain RDSI acquisitions.

2 Methods

2.1 ODF-Fingerprinting

In ODF fingerprinting, measured ODFs are matched to a library of ODF-elements (Fig. 1e). This library is generated by simulating diffusion weighted signals for a wide range of possible fiber combinations. Here we modeled each voxel as a composition of a water component $f_{W}$ and $N$ fibers, each with a volume fraction $f_{j}$ and a cylindrical diffusion tensor $\bar{D_{j}}$ ($\lambda_{2}\,=\,\lambda_{3}$). The simulated diffusion weighted signals are calculated as (Alexander et al., 2001; Tuch et al., 2002; Wilkins et al., 2015)

[TABLE]

with $\Delta$ and $\delta$ the diffusion time and the gradient duration, $D_{W}$ the water component diffusion coefficient and $\vec{q}$ the q-space wave vector. For each fiber combination, the simulated diffusion weighted signals are then reconstructed to form ODFs. These ODFs, expressed as vectors of numerical values are the entries of the ODF-library.

Matching of a measured $ODF_{m}$ is done by searching the ODF-dictionary for the ODF-entry with the best agreement. That is, for each $ODF_{m}$, we can find the index $l_{m}$ in the library matrix $L_{ODF}$ of the ODF-entry with the largest dot-product (Ma et al., 2013; Cloos et al., 2016):

[TABLE]

Although this straightforward library matching algorithm works, we found that it tends to favor multiple fiber configurations in noisy ODFs (Baete et al., 2018). Therefore, a penalty term was added to Eq. 2, which weighs complex fiber configurations in the library by a factor proportional to the noise estimate $\sigma_{n}$ of the input diffusion data, thus adapting the penalty to the dataset SNR. This approach is similar to model selection approaches such as the Akaike Information Criterion (AIC). The library matching algorithm then becomes

[TABLE]

with $n_{par}$ a measure of library element complexity, here $n_{par}\,=\,1+5N$ ($N$ the number of fibers). $\sigma_{n}$ is estimated using a linear minimum mean square error estimator for the variance of the noise of the diffusion weighted signals (Aja-Fernandez et al., 2008).

The size of the library is reduced by rotating the maximum values of both the library elements and $ODF_{m}$ to the Z-axis before matching. Note that the ODF fingerprinting approach differs from earlier methods where dictionaries of ODF elements were learned from the acquired dataset using iterative compressed sensing methods (Merlet et al., 2013; Bilgic et al., 2012; Gupta and Awate, 2017; Sun et al., 2013; Ye et al., 2012; Awate and Dibella, 2013) and individual ODFs are reconstructed as a combination of a sparse subset of these learned elements. In contrast, in ODF-FP the library elements are generated independently of the acquired dataset based on a simple diffusion model (Eq. 1). The use of the diffusion model also distinguishes ODF-FP from methods where the library elements are generated from Monte-Carlo simulations of diffusing water in specific microstructural configurations (Rensonnet et al., 2019).

The ODF-libraries used in this work simulate fiber combinations with up to 3 fibers. The main fiber is aligned along the Z-axis and subsequent fibers are sampled on a 642-point tessellation of the unit sphere (angular resolution of 2.3∘). Other parameters are: fiber FA ranging from 0.3 to 1.0 (in vivo) and 0.2 to 0.6 (simulations) in steps of 0.1 (2 fibers; 0.16 for 3 fiber libraries); fiber bundle volume ranging from 0 to (100%-water component) in steps of 10% (2 fibers; 0.16 for 3 fiber libraries); 10% water component with an ADC of 0.9 mm2/s (in vivo) and 1.0 mm2/s (simulations); minimum angle between fibers is 20∘in the 3 fiber libraries. For the in-house datasets, diffusion signals are sampled on a Radial Diffusion Spectrum Imaging (RDSI, (Baete et al., 2016)) grid with 236 q-space samples on four shells (250, 1000, 2250, 4000 s$/$mm2 (simulations) and 200, 1500, 2750, 4000 s$/$mm2 (in vivo)) and ODFs are generated with the RDSI-reconstruction. In the HCP dataset, diffusion signals are sampled on 256 q-space points on three shells (1000, 3000, 5000 s$/$mm2) and ODFs are generated with Generalized Q-Space Imaging (GQI,(Yeh et al., 2010; Yeh and Tseng, 2011)). The libraries contain 15,366 (2 fibers, simulations), 46,091 (2 fibers, in vivo) and 4,753,266, (3 fibers in vivo) elements and are created in 45 s, 51 s and 6 h 7 min 13 s respectively on a standard high-end laptop (Dell Precision 5510, Quad core Intel Xeon E3-1505).

The ODF-Fingerprinting method is compared to peak identification using local maximum search (DSIStudio111http://dsi-studio.labsolver.org (Yeh et al., 2010), compiled from source on Nov 21st, 2018; a version of the algorithm in Matlab (Mathworks) is also used), Newton search along a set of specified directions (MRtrix3222https://github.com/MRtrix3, v3.0_RC3, compiled from source on Jan 27th, 2019,sh2peaks, default parameters), multi-shell multi-tissue constrained spherical deconvolution (CSD, MRtrix3, dwi2fod msmt_csd(Jeurissen et al., 2014) and sh2peaks, default parameters, unsupervised estimation of response functions using dwi2response dhollander (Dhollander et al., 2016))and probabilistic estimation (FSL333https://fsl.fmrib.ox.ac.uk/fsl (Jenkinson et al., 2012), bedpostx (Behrens et al., 2007), v5.0.9). For bedpostx the following parameters were used: deconvolution with sticks and a range of diffusivities, constrained non-linear fitting, 250 burn-in points, 500 jumps, automatic relevance determination.

The matching process itself takes 2min14 (2 fiber library, 1h43min2s for the 3 fiber library) on a standard high-end laptop (Dell Precision 5510, Quad core Intel Xeon E3-1505) to reconstruct a full-brain 2.5 mm isotropic acquisition with a 46,091-element library (4,753,266 elements in the 3 fiber library). This compares to 12 s, 26 s, 3 min 24 s and 5 h 11 min 25 s for peak finding on the same dataset with local maximum search, Newton search, CSD and probabilistic estimation approaches, respectively.

Source code (Matlab) for the ODF-Fingerprinting approach is available for download at https://bitbucket.org/sbaete/odffingerprinting.

2.2 Computer simulations

Performance of fiber direction identification is evaluated using a composite hindered and restricted diffusion model (CHARMED, (Assaf and Basser, 2005)) as implemented in the Phantomas-software444http://www.emmanuelcaruyer.com/phantomas.php (Caruyer et al., 2014). The simulated volumetric fiber configuration was previously used in the HARDI reconstruction challenge at the ISBI 2013 conference.

In addition, ODFs of crossing fiber bundles are simulated using an in-house simulation. Diffusion weighted samples are generated as above (Eq. 1), with random fiber directions, on a RDSI q-space grid on four shells (250, 1000, 2250, 4000 s$/mm^{2}$. Intra-voxel fiber orientation dispersion of 20∘(Jelescu and Budde, 2017) is added in some simulations and Rician noise is added where necessary (SNR of $b_{0}$ indicated). After ODF-reconstruction, fiber directions are identified with ODF-Fingerprinting and the other methods described above. Results are binned based on the simulated crossing angle to calculate number of identified fibers, detected crossing angular error ($|\alpha_{found}-\alpha_{simulated}|$), angular precision and dispersion of the fiber directions (Kuo et al., 2008) and the detected crossing angle. The angular dispersion is a measure of angular accuracy as it assesses the uncertainty of mapping the fiber orientations (Kuo et al., 2008). If more than two fiber directions are found, the directions closest to the simulated directions are used for further calculations.

2.3 In Vivo acquisitions

In vivo DSI acquisitions of healthy volunteers are acquired on a 3T clinical scanner (Prisma, Siemens, Erlangen; 64ch head coil; 80 mT/m; Twice Refocused Spin Echo EPI sequence; RDSI q-space sampling (Baete et al., 2016), 4 shells, 200, 1500, 2750, 4000 mm2/s; TR = 2000 (bootstrap) / 8500 (whole brain) ms, TE = 86 ms, 10 (bootstrap) / 50 (whole brain) slices, field of view 200$\times$200 mm2, 2.5$\times$2.5$\times$2.5 mm3 resolution, partial Fourier 5/8). For each volunteer either whole brain or limited coverage datasets (five replications) for bootstrapping are acquired. For image correction an extra $b_{0}$-dataset is acquired with PA phase encoding rather than AP phase encoding. A $T_{1}$-weighted gradient-echo sequence (MPRAGE) serves as a reference for image registration (TR = 2300 ms, TE = 2.87 ms, 192 slices, 1$\times$1$\times$1 mm3 resolution, TI = 900/1000 ms, 5 min 03 s) The protocol was approved by an Institutional Review Board.

A high resolution preprocessed in vivo DWI acquisition was provided by the Human Connectome Project (HCP) consortium led by Washington University, University of Minnesota, and Oxford University. We used a single subject from the MGH datasets (3T Siemens Skyra System; 64ch head coil; 1000, 3000, 5000 mm2/s, 256 q-space volumes, TR = 8800 , TE = 57 ms, 96 slices, field of view 210$\times$210 mm2, 1.5$\times$1.5$\times$1.5 mm3 resolution, partial Fourier 5/8, GRAPPA 3; healthy volunteer).

2.4 DWI processing

Post-processing of datasets is performed offline. In-house images are denoised (Veraart et al., 2016) and corrected for susceptibility, eddy currents and subject motion using topup and eddy (FSL (Jenkinson et al., 2012)). The preprocessed(Glasser et al., 2013) HCP images were corrected for gradient non-linearity, motion (FreeSurfer) and eddy currents (FSL(Jenkinson et al., 2012) eddy). RDSI reconstructions, incorporating variable sample density correction, or GQI(Yeh et al., 2010; Yeh and Tseng, 2011) reconstructions (HCP dataset) are performed using custom-made software (Matlab, Mathworks) and displayed with DSIStudio (Yeh et al., 2010). Bootstrapped DSI datasets (500) are generated from five original datasets with limited brain coverage using repetition bootknife sampling (Cohen-Adad et al., 2011). From the bootstrapped datasets reproducibility and noise sensitivity metrics are calculated: the number of fibers identified and 95% confidence intervals (CI), coherence $\kappa$ (Jones, 2004) and Quantitative Anisotropy (QA) values of the first and second fibers.

Multi-resolution HCP datasets were created as follows. The 1.5 mm isotropic dataset was down sampled (MRtrix3, mrresize) to a 3 mm isotropic resolution such that each voxel in the low resolution (LR, 3 mm isotropic) dataset corresponds to 8 high resolution (HR, 1.5 mm isotropic) voxels. Hence, for each 3 mm isotropic voxel we compared the identified fiber directions relative to the fibers found in the 8 corresponding HR voxels. From this comparison we calculated the number of correctly (true positive) and wrongly (false positive) identified fibers and the number of missed fibers (false negative).

Fiber tracts are generated with a deterministic tracking algorithm ((Yeh et al., 2013), implemented in DSIStudio (Yeh et al., 2010), parameters as suggested in (Fernandez-Miranda et al., 2012)) in both simulated and in vivo datasets. Tracts originate from random seeding points uniformly distributed in user-defined seeding regions, propagate along the most prominent fiber direction with a step size of 2 mm, are smoothed with 20% of the previous direction and halt when the turning angle $>$ 60∘or QA $<$ 0.25. Resulting fiber tracts shorter than 30 mm and longer than 200 mm are discarded until a predetermined number of fiber tracts ($1\times 10^{4}$) is created or a maximum number of seed points is reached ($1\times 10^{6}$). While the choice of these tractography parameters influences the final results, relative differences in the tracts resulting from different peak identification methods remained constant as long as the parameters were consistent for all methods.

3 Results

3.1 Simulation Results

Fig. 2 compares the performance of the fiber identification of ODF-FP to the local maximum search algorithm in a diffusion dataset simulated with the Phantomas-software (Caruyer et al., 2014). ODF-FP better identifies small crossing angles (Fig. 2b,c,e,f, yellow arrows) even though the diffusion dataset is simulated with a diffusion model different from the model we used to generate the library (composite hindered and restricted diffusion model (Assaf and Basser, 2005) vs. sum of diffusion tensors (Alexander et al., 2001; Tuch et al., 2002; Wilkins et al., 2015)). ODF-FP also correctly identifies more fibers (Fig. S1a,e), produces less false negatives (missed fibers, Fig. S1c,g) and reduces the angular error of identified fiber directions (Fig. S1d,h). Detection of smaller crossing angles leads to improved fiber tractography in individual simulated fiber bundles (Fig. 3).

The results of individual crossing fiber ODF simulations are summarized in Fig. 4, 5 and S2. These simulations, with their inherent reference, allow detailed evaluation of the fiber identification performance as a function of the simulated crossing angle at different SNR levels (no added noise in the left column of Fig. 4, 5 and S2, SNR = 50 in the middle column and SNR = 25 in the right column). For SNR = 25, the b-value was limited to 3000 mm/s2 as the SNR of the b = 4000 mm/s2-shell was too low for processing. The different search strategies correctly identify most larger crossing angle pairs (i.e. two fibers are identified) when no noise is added but fail clearly when the crossing angles become smaller or noise is added (Fig. 4a,b,c). In particular, the ODF maximum search methods (local maximum and Newton search) detect too many fiber bundles in noisy data while CSD does not find all fibers; the performance of the probabilistic method seems to depend on the noise level. ODF-FP in contrast successfully finds just two fibers for a range of crossing angles (50 - 90∘) for noisy data.

Looking at the crossing angular error (Fig. 4d,e,f), the ODF maximum search methods lead to disproportionately large errors when estimating the angle of shallow crossings in the presence of noise. The probabilistic, CSD and ODF-FP methods, on the other hand, show errors that scale with the simulated crossing angle when the crossing angle is small (as only one fiber is identified Fig. 4b,c). The crossing angular error then drops for larger crossing angles as more fiber pairs are identified. Compared to CSD, ODF-FP is able to detect smaller crossing angle pairs (40∘  vs. 50∘, or 50∘  vs. 70∘  when the data is noisy; Fig. 4a-c). In more detail, the ODF maximum search methods and the probabilistic method tend to underestimate large crossing angles (Fig. S2a,b,c), while CSD and ODF-FP mostly correctly estimate these large crossing angles ($>\,$50∘, Fig. S2a,b,c). At smaller crossing angles ODF-FP, CSD and the probabilistic method do not find all of the crossing fiber pairs. These methods however also do not generate incorrect estimates in noisy datasets in contrast to the ODF maximum search methods (Fig. S2b,c).

Fig. 5 and S2d-i plot the angular precision and dispersion of the identified fiber directions. While Fig. 5 plots the values including all identified fiber combinations, the results in Fig. S2d-i are split by the number of fiber directions found for each ODF ($\Box$ and $\bigcirc$ when two fibers are found vs. $\Diamond$ when one fiber is found). Angular precision and dispersion of the ODF maximum, CSD and ODF-FP search methods is similar when no noise is added (Fig. 5a,d) with ODF-FP showing a higher precision and less dispersion. In general, the angular precision is high at large and small crossing angles when either a single fiber or the average of two crossing fibers is found. The angular precision and dispersion of the probabilistic method is lower (higher value for the precision) since this method does not identify all of the fibers. When noise is added, ODF-FP and CSD retain most of their angular precision and dispersion, while the ODF maxima search methods loose precision due to erroneous detection of fibers (Fig. 5b,c,e,f); the probabilistic method seems to perform better at lower SNR. Angular precision and dispersion split by number of fiber directions (Fig. S2d-i) largely follow the same trend as above.

Intra-voxel fiber orientation dispersion(Jelescu and Budde, 2017) can complicate the identification of fiber directions. Fig. 4a,d and 5a,d explore the implications of a 20∘(Jelescu and Budde, 2017) fiber dispersion (smaller dot symbols). These figures show a general decrease in precision and increase in crossing angular error of all methods. Overall, the impact of fiber dispersion is minimal.

3.2 In Vivo Results

A closer look at a few selected in vivo ODFs (Fig. 6) illustrates the differences in fiber directions identified by the different algorithms. ODF maximum search methods reliably detect strong peaks in the ODFs, but tend to miss more subtle smaller peaks (Fig. 6, ODF 3). Probabilistic estimation, on the other hand, seems to read too much into minor ODF peaks (Fig. 6, ODF 1,4). In contrast, ODF-FP and CSD identify the desired peaks whilst refraining from detecting non-existing fiber bundles (Fig. 6, ODF 3).

Evaluation of in vivo fiber direction identification can be performed by comparing fiber directions found in down sampled DWI with those present in high resolution DWI. The latter thus form an internal reference which is otherwise absent in in vivo data. Fig. 7 illustrates the setup of such an experiment; the directions found in a LR voxel (yellow arrows) are compared to those found in the corresponding HR voxels (yellow square in the panel on the right). In the indicated voxel, ODF-FP, CSD and the probabilistic method successfully identify the second fiber bundle. However, the probabilistic method and CSD also identify a number of false positive fibers. Maps of the number of correctly (true positive) and wrongly (false positive) identified fibers (Fig. 8, second and third row) confirm these findings over the whole volume.

Analysis of 300 bootstrapped RDSI datasets (Fig. 9) shows that ODF-FP identifies crossing fibers where expected (Fig. 9, #fib-column on the left) while CSD finds them to a lesser extent, ODF maximum search methods miss fibers and probabilistic estimation identifies an artificially high number of fibers. In addition, CSD and ODF-FP do not identify as many fibers in areas where no fibers are expected such as the Cerebrospinal fluid (CSF). This leads to different patterns in the reproducibility and noise sensitivity statistics (95% confidence intervals (CI) and coherence $\kappa$, Fig. 9). QA-maps further illustrate these altered patterns (Fig. 9, $QA$-columns).

When fibers are identified with ODF-FP, the 95% CI and $\kappa$ of the fibers are similar to the reproducibility values found for ODF maximum search methods and CSD, even for the second identified fiber ($CI_{2}$, $\kappa_{2}$-columns of Fig. 9). The 95% CI and $\kappa$ maps further show the noise sensitivity of the probabilistic method with higher 95% CI and lower $\kappa$-values.

Fiber tractography on a whole brain in vivo dataset shows that the fiber directions identified by ODF-FP (Fig. 10) allow the algorithm to improve results (Fig. 11). Fiber bundles generated with ODF-FP input probe the expected anatomical extent of the tracts in contrast to the fiber bundles based on ODF maximum search methods and CSD, in particular for the Corticospinal tract. The tendency of probabilistic methods to derive erroneous fiber directions (e.g. Fig. 6, ODF 1,4) produces more spurious tracts (Corticospinal tract, Arcuate Fasciculus and Optic Radiation in particular).

4 Discussion

As tractography based representations of the brain gain importance in clinical and neuroscientific applications, so grows the desire to resolve evermore detailed and complex neuronal pathways. One of the last remaining methodological challenges on the path towards such high-fidelity tractography representations is the identification of multiple intra-voxel fiber crossings. In particular, when dealing with small crossing angles. Although great progress has already been made (Baete et al., 2016; Descoteaux et al., 2007; Frey et al., 2008; Berman et al., 2008; Tournier et al., 2004; Aganj et al., 2010), these methods still fail to reliably detect crossing angles less than 40∘(Kuo et al., 2008; Jeurissen et al., 2013; Tournier et al., 2008; Descoteaux et al., 2009; Daducci et al., 2014; Wilkins et al., 2015). In this work, we have sought to overcome these limitations by changing the paradigm for fiber direction identification from a search for maxima on the surface of the ODF to an assessment of similarity relative to a dictionary of pregenerated ODFs.

The application of key concepts from fingerprinting to the ODF based fiber direction identification task (Fig. 1e) improves the detection of fiber pairs crossing at small angles as shown in simulation results (Fig. 2,4). This is achieved while maintaining angular precision of fiber directions over the whole range of crossing angles (Fig. 5a,b,c S2d,e,f). In vivo bootstrap analysis shows that ODF-FP detects crossing fiber pairs where anatomically expected (Fig. 9) while not over-identifying fiber bundles in areas where no fiber bundles are expected such as in the CSF (Fig. 9, # fibers). Furthermore, the fiber directions identified with ODF-FP are reproducible over the bootstrapped datasets (Fig. 9, CI and $\kappa$) and reproduce the internal reference of an in vivo multi-resolution HCP dataset (Fig. 7-8). Consequently, the improved fiber detection results in increased adherence of fiber tractography to the underlying simulated microstructure (Fig. 3). Although no gold standard is available for in vivo tractography, the data suggest that the fingerprinting based method to improve tractography results (Fig. 10, 11).

The performance of the ODF-FP methods compares favorably to ODF maximum search methods (local maximum and Newton search), CSD and probabilistic methods. These methods generally underestimate the crossing angle (Fig. S2a) or, at smaller angles, do not detect the crossing fiber pair but rather a single fiber (Fig. 4a). While the former situation is suboptimal, the latter biases the tractography (Fig. 3). Both problems originate in the focus on the maxima on the ODF surface in combination with the intrinsic ODF peak width (Barnett, 2009; Jensen and Helpern, 2016). The ODF maximum search methods perform better when attempting to identify fiber pairs with crossing angles between 40∘  and 55∘  in noisy data (Fig. 4e,f), though still overestimate the crossing angle (Fig. S2b,c). This improved detection does come at the cost of erroneous detection of non-existent fiber bundles (Fig. 4b,c) which might confuse tractography algorithms. The ODF-FP algorithm more precisely determines (smaller) crossing angles. Nevertheless, no approach to fiber identification from ODFs is perfect.

The ODF-library in this work is generated with the simple diffusion model described in Eq. 1. Results are consistent when this library is applied to simulations using the same generative diffusion model (Fig. 4, 5, S2), to ODFs simulated with a different diffusion model (CHARMED, Phantomas, Fig. 2, 3) and to in vivo multi-shell HCP (Fig. 7,8) and RDSI datasets (Fig. 6, 9, 10 and 11). Results are also consistent when multi-shell sampling and GQI ODF-reconstructions are used (Fig. 7-8)) rather than RDSI sampling and reconstruction (Fig. 2-6,9-11). This illustrates that the ODF-Fingerprinting approach can be used regardless of q-space sampling and ODF-reconstruction method.

The diffusion phenomenon in vivo is more complex than described by Eq. 1 (Novikov et al., 2016). Careful consideration of non-Gaussian diffusion in intra- and extra-axonal space, the axon diameter and dispersion Ghosh et al. (2016) of the fiber bundles leads to models such as CHARMED (Assaf and Basser, 2005), ActiveAx (Alexander et al., 2010), Neurite Orientation Dispersion and Density Imaging (Zhang et al., 2012), and White Matter Tract Integrity (Fieremans et al., 2011). Each of these models, with an appropriate range of tissue parameters, can be used to generate more complex dictionaries for the ODF-fingerprinting method, including fiber configurations such as branching and kissing fibers. By adhering closer to in vivo microstructure, these complex dictionaries may eventually be better suited for ODF-FP and additional model parameters associated with the diffusion models may be simultaneously estimated. By removing the need for direct parameter fitting, we then adhere to the philosophy behind quantitative MR fingerprinting (Ma et al., 2013; Cloos et al., 2016).

The computational bottleneck in the ODF fingerprinting method is the matching algorithm (Eq. 3). Fortunately, the central dot-product in Eq. 3, $L_{ODF}\cdot ODF_{m}^{T}$, can be efficiently implemented as a matrix multiplication. Another consideration is the dictionary size. Simulating two crossing fibers with random orientations on a 642 point grid (321 options due to symmetry) gives 321x320 = 102,720 possibilities. This number has to be multiplied with the desired number of fiber parameter combinations. In the ODF-library used here for the in vivo reconstructions for instance, 144 parameter combinations were considered leading to a total of $~{}14.8\times 10^{6}$ possibilities. A major reduction in dictionary size is possible by rotating the maximum value of the ODF-traces to e.g. the Z-axis before matching, reducing the library size by a factor of 321 to $~{}46.1\times 10^{3}$ in our example for 2 fibers. Further reductions may be possible by rotating the ODF to align the second maximum when present. The ODF-FP approach as currently implemented is slower than the ODF maximum search methods (2 min 14 for 2 fiber and 1 h 43 min 2 s for 3 fiber dictionaries compared to 12 s and 26 s for local maximum search and Newton search in a full-brain acquisition); though fast enough for use in general processing pipelines. Acceleration is possible by porting the implementation from Matlab to more general-purpose programming languages. Furthermore, since the ODF-fingerprinting method operates on ODFs, it is easily incorporated in existing postprocessing pipelines.

Tractography algorithms guided by the more accurate fiber detection of ODF-FP adhere better to the underlying tissue microstructure (Fig. 11), thus enhancing the utility of tractography representations in clinical and neuroscientific applications. In neurosurgery, precise fiber bundle delineation informs surgical decisions (Fernandez-Miranda et al., 2012; Shin et al., 2012), while structural brain connectivity analysis is leveraged in neuroscientific research (Jbabdi et al., 2015; Galantucci et al., 2016; Mitra et al., 2016). Improved tractography will thus aid these endeavors.

A limitation of the ODF-FP method is the assumption of a diffusion model. This diffusion model, independent of which model was selected, may be insufficiently generalizable to encompass diseased tissue. Also, in the ODF-library, depending on the model used and number of grid points, the entries may be non-unique. This is expected to be more of a problem when working with more complex diffusion models with multiple parameters. A more complete validation of the ODF-FP and other fiber direction identification methods can be performed using hollow fiber phantoms (Guise et al., 2016; Hubbard et al., 2015; Pathak et al., 2017).

In conclusion, we demonstrated a novel method for fiber direction identification from ODFs based on key concepts first introduced in MR Fingerprinting. In ODF-Fingerprinting, fiber configurations are selected based on the similarity of the shape of measured ODFs with pre-computed library elements. This approach improves detection of fiber pairs with small crossing angles while maintaining fiber direction precision. The resulting, more precise, fiber directions aid fiber tracking algorithms in more accurately calculating brain connectivity for clinical and neuroscientific applications.

5 Acknowledgements

This project is supported in part by the National Institutes of Health (NIH, R01CA111996, R01NS082436 and R01MH00380). Some of the data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

6 References