Joint Multi-Fiber NODDI Parameter Estimation and Tractography Using the Unscented Information Filter

Tracing white matter fiber bundles is an integral part of analyzing brain connectivity. An accurate estimate of the underlying tissue parameters is also paramount in several neuroscience applications. In this work, we propose to use a joint fiber model estimation and tractography algorithm that uses the NODDI (neurite orientation dispersion diffusion imaging) model to estimate fiber orientation dispersion consistently and smoothly along the fiber tracts along with estimating the intracellular and extracellular volume fractions from the diffusion signal. While the NODDI model has been used in earlier works to estimate the microstructural parameters at each voxel independently, for the first time, we propose to integrate it into a tractography framework. We extend this framework to estimate the NODDI parameters for two crossing fibers, which is imperative to trace fiber bundles through crossings as well as to estimate the microstructural parameters for each fiber bundle separately. We propose to use the unscented information filter (UIF) to accurately estimate the model parameters and perform tractography. The proposed approach has significant computational performance improvements as well as numerical robustness over the unscented Kalman filter (UKF). Our method not only estimates the confidence in the estimated parameters via the covariance matrix, but also provides the Fisher-information matrix of the state variables (model parameters), which can be quite useful to measure model complexity. Results from in-vivo human brain data sets demonstrate the ability of our algorithm to trace through crossing fiber regions, while estimating orientation dispersion and other biophysical model parameters in a consistent manner along the tracts.


INTRODUCTION
Diffusion MRI (dMRI) is a non-invasive technique to study the microstructure of brain tissue. However, we need a mathematical model to interpret the diffusion weighted signal to study the microstructure of white matter fibers. Broadly, such models fall into two categories: parametric and nonparametric. The simplest parametric model is the diffusion tensor model, which describes a Gaussian estimate of the strength and diffusion orientation at each voxel (Basser et al., 1994;Alexander et al., 2002). While robust, this model can be inadequate in cases of mixed fiber presence or more complex orientations (Frank, 2002). To handle more complex diffusion patterns, parametric models have been introduced including mixtures of tensors and directional functions (Alexander et al., 2001;Anderson, 2005;Kreher et al., 2005;Parker and Alexander, 2005;Peled et al., 2006;Kaden et al., 2007;Zhang et al., 2007;Rathi et al., 2008). Several techniques attempt to reconstruct pathways based on these models. In this case, tractography is done by following the principal diffusion direction(s).
With multifiber models, care must be taken in fitting the model parameters to the recorded signal consistently, such that, the correlation in diffusion along the tract is accounted for while estimating the parameters (Malcolm et al., 2009a). As demonstrated in this paper, estimation based on information from previous estimates aids in this process.
Nonparametric techniques, unlike estimating a fixed number of fibers in parametric models, estimate an oriented distribution function (ODF) describing an arbitrary configuration of fibers. For this estimation, Tuch (2004) introduced Q-ball imaging to numerically compute the ODF via the Funk-Radon transform, and later, spherical harmonics were used to simplify the computation with an analytic form (Anderson, 2005;Hess et al., 2006;Ozarslan et al., 2006;Descoteaux et al., 2007). Another method to estimate the fiber ODF is to assume a model for the signal response of a single-fiber and use spherical deconvolution (Jansons and Alexander, 2003;Tournier et al., 2004;Jian and Vemuri, 2007;Kaden et al., 2007) to obtain a much sharper orientation profile. A good review of both parametric and nonparametric models and diffusion MRI in general can be found in Assemlal et al. (2011).
Recently, a model of diffusion was proposed by Zhang et al. (2012) called Neurite Orientation Dispersion Diffusion Imaging or NODDI, which accounted for the dispersion in orientation of the axonal fibers. In this work, the authors proposed an algorithm to estimate the parameters of the NODDI model assuming the existence of a single fiber pathway at each voxel. However, it is well known that most of the white matter tissue in the brain has crossing fibers (Tuch et al., 2002), which must be taken into account for proper analysis of the estimated microstructural parameters. Further, the correlation in diffusion (and consequently the estimated NODDI parameters) is not taken into account in the work of Zhang et al. (2012). Thus, estimating the orientation dispersion in the presence of multiple fiber crossing would be quite useful in analyzing the geometrical structure of white matter in healthy as well as diseased subjects.

Our Contributions
In this work, we propose to use the unscented information filter (UIF) based framework to perform joint NODDI parameter estimation and tractography. Earlier works have used the unscented Kalman filter (UKF) with a parametric multi-tensor model or a non-parametric spherical harmonic representation to do simultaneous model estimation and tractography (Malcolm et al., 2010;Baumgartner et al., 2012). The UIF has several advantages over the UKF as has been noted in Lee (2008), namely, it propagates the Fisher-information matrix as opposed to the state-covariance matrix leading to a significant reduction in computational load, while providing a more robust estimation of the estimated parameters. Further, the Fisher-information is calculated for each of the model parameters along the tract, providing the variance in the estimation of the model parameters.
Existing methods either use a region-of-interest (ROI) or an atlas within which the NODDI parameters are estimated and analyzed or follow the principal orientation to perform tractography. Thus, the fit is performed independently at each voxel disregarding the correlation of diffusion along the fiber path. In the UIF framework, we perform model estimation and tractography simultaneously to trace fiber tracts. The methodology assumes a "time varying" Gaussian distribution of the model parameters as the algorithm moves from one location to the next. This allows for consistent estimation of the NODDI parameters, while allowing for smooth estimation of the fiber tracts. In the current work, we demonstrate the performance of our method in the case of single fiber and 2fiber NODDI model. To the best of our knowledge, this is a first extension of the NODDI model to the multi-fiber case within a tractography framework. Thus, the ability to trace multiple crossing fibers while robustly estimating the NODDI model parameters is one of the major contributions of this work. We demonstrate our technique on in-vivo human data set from the human connectome project (HCP) data set. We expect, that the proposed method will be quite useful to study the tissue microstructure in several disorders.

APPROACH
The main idea of our approach is to trace the local fiber orientations employing the three compartment NODDI model (intra-cellular, extra-cellular, and isotropic free-water) and using the estimation at previous positions to guide estimation at the current position. In a recursive fashion, the information filter estimates the model at the current position, moves a step in the most consistent direction, and then begins estimation again. In this case, the model parameters, which form the state of the filter, are assumed to have a Gaussian distribution, whose mean and covariance change at each step. Thus, the space of possible solutions (among the infinitely many possibilities) is greatly reduced leading to better accuracy in resolving individual orientations and yielding inherently smooth tracts despite the presence of noise and uncertainty. This is in complete contrast to existing methods, where the model parameters at each voxel are estimated from an arbitrary initialization or a large number of random initializations to finally obtain a solution that best fits the data. This could potentially lead to vastly different solution even in neighboring voxels since the space of possible solutions is large. In the proposed method, since each iteration begins with a nearoptimal solution provided by the previous estimation, and since the step size is kept to be very small, the UIF quickly converges to the optimal solution and many local minima are naturally avoided.

Modeling Local Fiber Orientations
The three compartment NODDI model consists of an intra-cellular, an extra-cellular and an isotropic compartment (Zhang et al., 2012;Daducci et al., 2015). The normalized signal E can be written as: (1) where E ic and V ic are the normalized signal and volume fractions of the intra-cellular compartment; E ec is the normalized signal corresponding to the extra-cellular compartment; and E iso and V iso are the normalized signal and volume fractions of the isotropic compartment. The space bounded by the membrane of neurites is called the intra-cellular compartment. In this paper, we adopt a model where neurites are modeled as a set of cylinders of zero radius capturing the unrestricted nature of diffusion along the neurites and restricted diffusion perpendicular to them, i.e., is the orientation dispersion function around vector n; and d is intrinsic diffusivity along the cylinder. As in Zhang et al. (2012), the watson distribution is used to model the orientation dispersion function f : where m is the mean orientation about which the dispersion is symmetric and κ is the concentration parameter, which determines the extent of orientation dispersion along the mean orientation; and M is the confluent hypergeometric function. The orientation dispersion index is given by OD = 2 π arctan 1 κ . The space around the neurites is called the extra-cellular compartment. In this compartment, water diffusion is assumed to be hindered and is modeled as a Gaussian anisotropic tensor, with where D(n) is a cylindrically symmetric diffusion tensor with principal orientation n, and d and d ⊥ are the coefficients of diffusion parallel and perpendicular to n. Both these parameters can expressed in terms of the intra-cellular fraction V ic and the concentration parameter κ (see Zhang et al., 2012 for more details). The isotropic compartment is modeled with an isotropic Gaussian diffusion with diffusivity set to d iso .

Estimating the NODDI Parameters
As described in Zhang et al. (2012), we fix the following parameters within the NODDI model: d = 1.7 × 10 −3 mm 2 s −1 and d iso = 3 × 10 −3 mm 2 s −1 . Consequently, for a single fiber NODDI model, the free parameters are: V ic , κ, m and V iso . For the two-fiber NODDI model, we used the following formulation: In this case, the free parameters to be estimated are: V ic1 , κ 1 , m 1 , V ic2 , κ 2 , m 2 , and V iso . Given the measured signal at a particular voxel, we want to estimate the above model parameters that best explain the signal. We propose to achieve this using the unscented information filter (UIF), which is a recursive non-linear least squares estimator, giving the maximum likelihood estimate of the model parameters. Further, the Fisher information computed by the UIF provides a lower bound on the precision with which the model parameters can be estimated given the data. This statistic can be summarized into a single number as the estimated uncertainity at each point, by computing the matrix norm of the estimated covariance matrix. This uncertainity measure can be quite useful in removing unlikely fibers (with high uncertainity) from the tractography obtained.
Note that, such information is typically not available using other model estimation methods. As in streamline tractography, we treat the fiber as the trajectory of a particle which we trace out. At each step, we examine the measured signal at that position, use that measurement to update our model parameters within the filter, and propagate forward in the most consistent direction. To use a state-space formulation for estimating the model parameters, we need the following application-specific definition of four filter components: For our state (x), we directly use the parameters of the NODDI model (in the case of 2-fibers): For the state transition function s, we assume identity dynamics since the local fiber configuration does not undergo drastic change as it moves from one location to the next, when the step size is kept very small. The predicted signal is computed using: y = h[x] given by Equations (1) or (2) depending on the model used (1 fiber or 2-fiber NODDI) and our actual measurement is the actual signal interpolated directly on the diffusion weighted images at the current position.
It is important to note that we chose the unscented Information filter for its low computational complexity compared to the UKF. In the UKF formulation, one not only estimates the state, but also the state covariance matrix P, whose dimension is k × k, where k is the number of diffusion weighted gradients. At each step, the covariance matrix P is updated, which involves inversion of a large matrix that depends on k. This makes the filter computationally very expensive, especially in cases like the human connectome data, which has about 270 measurements. In contrast, in the UIF filter, only the Fisher-Information matrix is computed and propagated, which involves inversion of a matrix whose maximal dimension is substantially small, i.e., the length of the state vector. Thus, in the case of the two-fiber NODDI model, one only needs to invert the FIGURE 1 | Visualization of intracelluar volume fraction, fiber dispersion, isotropic volume fraction, data fitting error and uncertainity in parameter estimation using 1-fiber noddi model in the arcuate fasciculus. The background slice is the single tensor FA map. information matrix of size 11 × 11. This significantly increases the computational speed of the algorithm.
We use the unscented information filter with constraints, where the final solution is projected onto the physiological range of the parameters using a quadratic programming problem. For example, V iso , V ic are constrained to lie between [0, 1], while κ is always set to be positive. For a more thorough treatment of the UKF and UIF filters with constraints (see Lee, 2008;Malcolm et al., 2009bMalcolm et al., , 2010. In Appendix, we describe the mathematical equations for the prediction and measurement update steps of the UIF filter.

The tractography Algorithm
The UIF filter provides the maximum a-posteriori estimate of the model parameters, given the signal at each location. We embed this into a tractography framework, as described in detail in Malcolm et al. (2010). Briefly, we begin by initializing the UIF filter at each seed point by using the coarse-grid search method FIGURE 2 | Corticospinal tract showing orientation dispersion from two different views. In the sagittal view, the background is FA, while in the coronal view the 1-fiber NODDI dispersion map (Daducci et al., 2015) is shown in the background. as described in Zhang et al. (2012). Subsequently, we run the UIF filter at each seed point, which provides an estimate of the model parameters, including the principal orientation(s) of the fiber bundle. A small step (using a fixed pre-determined step size) is taken along the direction of the principal fiber orientation that is most consistent as compared to the previous estimate. At sub-voxel locations, the signal is interpolated using an isotropic Gaussian kernel with its width (variance) given by the smallest voxel length in any direction. This interpolated signal is then used as the acquired measurement (y), within the UIF filter to estimate the model parameters as well as the Fisher information and the covariance matrix. In a loop, the fiber is then traced until a termination criteria is reached. In the present case, we used generalized fraction anisotropy (GFA) threshold of 0.08 and κ of 0.06 as the stopping criteria.

EXPERIMENTS
We tested the proposed algorithm by tracing several fiber bundles from the in-vivo human data set obtained from HCP (Van Essen et al., 2013). All b-values of {1000, 2000, 3000}s/mm 2 were used to perform whole brain tractography using single fiber and twofiber NODDI models.
On a 2.4 GHz processor with 16 cores, the UKF filter (with 2fiber NODDI model) required about 72 hours of computational time for a single subject whole-brain tractography with 1 seed per voxel on the HCP data set with 270 gradient directions. With the same set of parameters, the UIF-based whole brain tractography required about 30 hours, an improvement by more than a factor of 2. Whole brain tractography was performed using with the following parameter settings: expected rate of change of orientation q m = 0.002, rate of change of κ (dispersion), q κ = 0.015 and expected rate of change of intracellular and isotropic volume fraction was set to q ic = q iso = 0.0007. The parameter r s is akin to a regularization parameter that should be set based on the expected noise level in the data. In the case of HCP data, we set it to 0.02, based on the noise level in the data. Thus, for higher noise in the data, r s should be increased, which implicitly implies that the algorithm will trust the model more than the data. On the other hand, for high SNR data, r s should be reduced to allow the tractography method to trust the data more than the expected model. For all the experiments, we set the step length parameter to 0.5 mm.
The following figures show the traced fiber bundles extracted using the white matter query language (WMQL) (Wassermann FIGURE 4 | Corticospinal tract traced using the 2-fiber NODDI model. In the background is a slice of orientation dispersion obtained from 1-fiber NODDI model of Daducci et al. (2015). Several NODDI specific measures are shown along the tracts along with the data fitting error, which is below 3% in most cases. et al., 2013). In particular, Figure 1 shows the arcuate fasciculus traced using the 1-fiber NODDI model. Estimates of several diffusion measures of interest, such as, intracellular volume fraction, orientation dispersion, isotropic volume fraction, normalized mean squared error (NMSE) in fitting the data and uncertainity in the estimated parameters are shown along the tract with the standard single diffusion tensor based FA map shown in the background.
In Figure 2, the corticospinal tract is shown as traced using the 1-fiber NODDI model. In the coronal view, the background is a slice of fiber orientation obtained using the method in Daducci et al. (2015). Note the similarity in the obtained measurements using both the methods, i.e., high dispersion in the CSF and gray matter areas and low in the deep white matter regions. Figure 3 shows results for the 2-fiber NODDI model obtained using the UIF filter. The cortico-spinal tract (shown in color) FIGURE 5 | Estimated uncertainty in the model parameters in the CST (left) and IOFF (right) fiber bundles. Fibers with high uncertainty (likely false positives), as marked by white arrows, can be easily removed by thresholding the estimated uncertainty.
FIGURE 6 | Comparision of results between 1-fiber NODDI model tractography(red) and 2-fiber NODDI model tractography(green) for the IOFF fibers. and the superior-longitudinal fasiculus II (SLF-II) intersect in the centrum-semiovale region. As can be seen, the proposed algorithm is able to trace fibers through crossing regions, which is not possible using the 1-fiber NODDI model. Also, the lateral fibers of the corticospinal tract (CST) that go to the hand and face area are missing in the 1-fiber NODDI model, but can be nicely traced using the 2-fiber NODDI model.
Results shown in Figure 4 demonstrate how the proposed method can be used to trace the CST to the hand and face motor areas. The estimated model parameters, such as intracellular volume fraction (for both the fibers of the 2-fiber NODDI model), orientation dispersion and the isotropic volume fraction are also shown. Note once again, that we show the orientation dispersion and intracellular volume fraction for the second fiber as well, although the tracts were obtained by following the primary fibers. The results demonstrate the smooth and robust estimation of FIGURE 7 | Arcuate Fasciculus traced using 1-fiber NODDI model. Background is the b = 0 image, where very bright regions indicate CSF areas. As seen, higher error in data fitting occurs only in the CSF areas, which is extremely noisy with free isotropic diffusion.
all the model parameters, including the ones for the second fiber.
The UIF filter can also estimate the uncertainty in the estimation of the model parameters. This information can be used to removed unlikely fibers or false positives from the tractography, which offers a powerful way to automatically detect such fibers and remove them. One such example is shown in Figure 5 as pointed by the white arrows. Note that, to the best of our knowledge, only the UIF (and UKF) based tractography methods allow an inherent way to detect unlikely fibers. Figure 6 shows a comparison of the traced inferior occipitofrontal fibers (IOFF) as traced by the 1-fiber (red) and 2-fiber (green) NODDI model. The 2-fiber NODDI manages to trace and connect a different part of gray matter region that is missed by the 1-fiber NODDI model. In fact the 1-fiber tracts only trace the medial portion of the lateral occipital cortex, whereas the 2-fiber NODDI tracts cover most of the lateral-occipital cortex as labeled by Freesurfer. Thus, the 2-fiber NODDI potentially provides a better estimate of fiber connectivity.

DISCUSSION AND LIMITATIONS
In this work, we applied a new computationally efficient and numerically robust unscented information filtering framework for joint estimation of NODDI parameters and tractography. The proposed UIF filter provides a 2-fold improvement in processing speed, while computing the uncertainty in the estimated parameters, which is generally lacking in most existing methods. The method allows for tracing crossing fibers while accounting for the correlation in diffusion along the tracts. To the best of our knowledge, this is the first time, multi-fiber NODDI model has been used to estimate fiber dispersion index and intracellular-extracellular volume fractions for each of the crossing fibers separately, within a tractography framework. This can be quite useful in understanding the geometric properties of the white matter as it is traced along the tracts. The proposed algorithm is an open-source software and can be downloaded from https://github.com/pnlbwh/ukftractography, with the option of using the NODDI model as one of the possible choices.
We should however note that, the fiber dispersion dispersion (OD) index computed using the NODDI model is different than the one obtained using the methods in Savadjiev et al. (2010Savadjiev et al. ( , 2012. In particular, the latter are computed from the fiber tracts or tensor fields at a macroscopic level, whereas the fiber dispersion obtained from NODDI is inherently microscopic. While there could be regions where they agree (some regions of white matter), yet in the gray matter the OD measures from NODDI is very different than the fiber dispersion computed using the method in Savadjiev et al. (2012).
Nevertheless, the proposed method has a few limitations. First, we assume equal volume fraction for each of the crossing fibers in our 2-fiber NODDI model, which may not be accurate. Second, the model fit to the data in the CSF areas is poor due to high noise in the data itself, as seen in Figure 7. However, the error in most white and gray matter areas is quite low, i.e., less than 2%. Another limitation of the current implementation of the proposed method is its inability to trace more than 2 fiber crossings. While an extension to trace 3 fiber model is straightforward, it can be done in areas which are a-priori known to have 3 fiber crossings. However, we believe that using a 3fiber model for tracing all fibers would result in over-fitting of the data. Yet, we should note that the a majority of the white matter voxels have two crossings, and a very small region has more than 2 fiber crossings. Thus, the proposed method can be applied in most parts of the brain to trace fibers and estimate the dispersion index of each fiber separately, which is a significant improvement over the existing single-fiber based model of Zhang et al. (2012).
Another limitation, as is the case with most tractography algorithms, is the directional nature of the estimated tracts. For example, the tracts seeded in region A and reaching B, may not be obtained if seeding was done in region B. However, a typical way this particular problem is addressed is by seeding the whole brain and extracting tracts of interest from the whole brain tractography as has been done here and in most works using deterministic tractography.
Nevertheless, we believe that the present method allows to estimate parameters of the NODDI model along fiber tracts and allows to trace fibers through crossing regions. This could be useful in neuroscience studies to detect changes in white matter structure due to disease.

AUTHOR CONTRIBUTIONS
YR: Conception, Design, Manuscript writing and data analysis. PR: Implementation, paper writing, data analysis.