1 Introduction
Among the many uses of diffusion MRI, using tractography methods to depict the underlying white matter fiber tracts of tissues may be the most important. Procedures of identifying tracts are referred as to fiber tracking. For clinical practice, fiber tracking provides potential benefits for presurgical planning (Chung et al., 2011). In neuroscience, understanding the anatomical connection of the brain is an important component of the connectome (Sporns et al., 2005).
Many technologies have been developed for fiber tracking in recent years. Among these methods, DiST (short for Diffusion Direction Smoothing and Tracking) (Wong et al., 2016) is one of the most prominent (Kang and Li, 2016; Schwartzman et al., 2016; Lazar et al., 2016). DiST
is composed of three major steps: In Step 1, voxelwise diffusion directions are estimated using diffusionweighted signals; In Step 2, the estimated diffusion directions obtained in Step 1 are smoothed over space; Finally in Step 3, the smoothed diffusion directions are taken as the inputs of a fiber tracking algorithm that determines if some voxels construct a fiber.
Although DiST has many appealing features, it has some limitations (Kang and Li, 2016; Lazar et al., 2016; Schwartzman et al., 2016). A major limitation may be that the separate steps of DiST
make it difficult to properly account for statistical uncertainties. In light of this, we propose a Bayesian hierarchical approach which allows valid statistical inference for fiber tracking. First, we assume that the logarithm signals follow a normal distribution, simplifying the model by avoiding the challenging Rician distribution
(Wong et al., 2016). Also, we induce spatial smoothness using a random field for spatially dependent positive definite matrices instead of the optimizationbased smoothing procedure of DiST. This avoids the optimization issue raised by Schwartzman (Schwartzman et al., 2016).In the rest of the paper, we first introduce our method, referred as to SpDiST (short for Spatial Diffusion Direction Smoothing and Tracking). To demonstrate our proposal, we use both real data analysis and a simulation study. The real data analysis demonstrates that our proposal provides a valid and efficient means to quantify the uncertainties of fiber tracking. Moreover, the simulation study shows that our proposal produces accurate estimation. Finally, we conclude with a discussion.
2 Method: SpDiST
2.1 Spatial Tensor Model
In this section, we introduce the spatial tensor model based on directed acyclic graph autoregression for positive definite matrices. The diffusion MRI has
measurements at voxel , denoted as . The measurements are used to estimate the diffusion tensor for voxel . is a positive definite matrix interpreted as covariance matrix of a local Brownian motion, indicating the local tensor direction. The goal is to use the measurements to obtain tensor direction information from .The noiseless signal intensity can be expressed in terms of (Mori, 2007) as
In this expression, , , and are nondiffusion weighted intensity, scale parameter, and
unitnorm gradient vector, respectively. A detailed explanation of these three quantities can be found in Soares et al.
(Soares et al., 2013). Given ,can be understood as the probability intensity of the Gaussian motion when measuring at direction
. For statistical modeling, , , and can simply be understood as fixed and known values.The observations are noisy realizations of . The Rician distribution (Wong et al., 2016) is reasonable for modeling but it causes computational issues. Here, we assume that the noise is a multiplier to the and it follows a lognormal distribution. The model is
(1) 
where
is the noise following a meanzero normal distribution with variance
.To induce spatial smoothness, an image is treated as a directed graph whose nodes are voxels and whose directed edges are from node to nodes in . Following Datta et al. (Datta et al., 2017), we use the directed acyclic graph (directed and no loops) to construct to , leading to a valid joint density function of . In particular, we assume that the conditional mean of is the average of its neighboring tensors, denoted as , where is a set containing neighboring voxel indices of voxel , and is the set size.
In a directed acyclic graph, we have at least one voxel whose is an empty set. For is an empty set, we assume that follows a Wishart distribution with mean matrix
. Otherwise, conditional on , we assume that follows a Wishart distribution with mean matrix and degrees of freedom . The model is(2)  
In Equation (2), to preserve the designed mean realizations, we parameterize the Wishart distribution for to have
. The probability density function is
where is the matrix dimension and is the multivariate gamma function.
Here, we give an approach to construct a directed acyclic graph. For an image, we construct an undirected graph whose voxels are nodes, and the neighboring nodes are connected. We order the voxels by their coordinates, i.e., for a 2D image on a xy axis, we first order the voxels according to their coordinates of the yaxis, then next we order the voxels according to their coordinates of the xaxis. For each edge of the undirected graph, we modify the undirected edge to a directed edge which is from the node with a smaller rank to a node with a larger rank. The modified graph is a directed acyclic graph whose edges connect neighboring voxels. In Figure 1, we give an example describing how a directed acyclic graph for a image is constructed.
2.1.1 MCMC Algorithm
We use MCMC for model fitting. We give and . The primary challenge in the MCMC algorithm is to sample the posterior of . Since the prior of is not conjugate, we sample it using singlesite MetropolisHastings sampling with Wishart distribution as the proposal distribution. The algorithm is described below:
 Candidate Generation:

Generate a candidate sample using ;
 Acceptance Rate:

Calculate the acceptance rate , where
(3) where . and are the density functions of Wishart distribution and normal distribution, respectively.
 Decision:

Generate . If , accept .
The acceptance rate can be tuned by the degrees of freedom , where smaller leads to smaller acceptance rate. We tune to make the acceptance rate around .
We use MetropolisHastings algorithm with lognormal random walk as proposal distribution to update the degrees of freedom and use Gibbs sampling to update based on its posterior .
2.2 Probabilistic Fiber Tracking Algorithm
We collect the MCMC samples of , denoted as
. For each sample, we compute the principal eigenvector of
, denoted as . For each posterior draw, we use as inputs of a fiber tracking algorithm. In this paper, we continue to use the Fiber Assignment by Continuous Tracking (FACT) (Mori et al., 1999), following Wong et al. (Wong et al., 2016). The algorithm can be stated as
Initialization: Starting from seed voxels;

Recursive: Starting with voxel , we search neighboring voxels and compute the two angles: is the angle between the two tensor directions ( and ) and is the angle between the current tensor () and betweenvoxel direction (). See Figure 2 for details. We move to the voxels with and . If there are multiple voxels statisfying this condition, we move to all the voxels and treat each voxel as a current voxel for next iteration.;

Result: Sequences of voxels constructing fibers.
Since we apply the algorithm for each posterior draw, the algorithm returns possible fibers. We summarize distinct patterns from the outputs and calculate the associated probability for pattern defined as , where is the frequency of the pattern . This procedure is known as probabilistic fiber tracking and quantifies the uncertainties of fiber tracking result.
3 Real Data Application
In this section, we use a real data example (Dryden et al., 2009, Section 6) to demonstrate our proposal. In particular, we focus on uncertainty quantification. The real data has voxels and measurements. A detailed description can be found in Dryden et al. (Dryden et al., 2009, Section 6). We sample MCMC samples after discarding samples as burnin and thin the MCMC chain by retaining every iterations of the chain.
Since it is more efficient to visualize tensor directions in a 2D environment and the image is 2D, we focus on the first two dimensions of and compute the corresponding principal eigenvector . To quantify the uncertainties of tensor direction estimation in each voxel, we overlay the MCMC samples on a map (Figure 3). In Figure 3, the voxels with heterogeneous directions have large uncertainties. Otherwise, there are small uncertainties.
Figure 3 only provides voxelwise uncertainties but our fullyBayesian approach can propagate spatial uncertainties through to uncertainty in fiber tracking. In this way, the MCMCbased SpDiST also provides a probabilistic approach to quantifying the uncertainties of fiber tracking. To have a concise and representative illustration, we focus on the region in the orange box of Figure 3. We apply the FACT algorithm as described in Section 2.2. In light of the conventions in setting the threshold (Chung et al., 2011), we consider ranging from to . There are two distinct patterns (Pattern A and Pattern B) for (Figure 4
) as dominating the posterior probability of the tract. These two tracts differ only by how far the tract continues vertically in column 18. The probabilities of each pattern vary by different thresholds
.Kang and Li (Kang and Li, 2016, Section 3) show that the FACT algorithm hinges on the tuning parameter . It requires a sensitivity analysis to explore the impact of . Here, we give a sensitivity analysis. We apply the FACT algorithm with and . Since there are only two distinct patterns, we report the probabilities of Pattern B with different thresholds (Figure 5). We find that the result is sensitive to the choice of unless it is ranging from to .
4 Numerical Study
4.1 Data Description
In this section, we use synthetic diffusionweighted signals in Wong et al. (Wong et al., 2016, S6) and further modify them for our numerical study. In total, we have voxels where the three digits represent the dimension of axis, axis, and axis, respectively. The underlying tensors and fibers from the synthetic signals are displayed in Figure 6. A comprehensive description of example data generation can be found in Wong et al. (Wong et al., 2016, S6) (i.e., generating model, parameters, true tensor directions, etc.). Here, we give a brief description. The fibers are essentially arcs with the center point at right/left bottom points. For voxels composing fibers, its principal eigenvector is tangent to the arc. The noiseless signal in the example data is given as , a reparameterized model of Model 1 (Wong et al., 2016).
To mimic lowquality images with signal noise, we further add noise on the log scale simulated from a meanzero normal distribution with standard deviation
. That is, the simulated data for each replication () is(4) 
where is the simulated signals for each replication (), is the logarithm signal from the example data, and is simulated noise.
4.2 Simulation Details
We construct as described before. We use the posterior mean estimate of SpDiST to compare with the estimates of alternatives. We compute posterior mean of based on MCMC samples after samples as burnin. In comparison, we compare our method to DiST. In addition, we also compare our method to a nonspatial method: the least squares method (Niethammer et al., 2006). The least squares method (Niethammer et al., 2006) is to estimate via
For DiST, the estimates are the principal eigenvectors. To compare to DiST, for SpDiST, we compute the principal eigenvectors of the posterior means of diffusion tensor. For comparison, we also compute the principal eigenvector of the diffusion tensor estimate of the least squares method.
4.3 Results
To quantify the performance of the three methods, we introduce two metrics. For voxels with fiber directions, we use Metric 1
a metric measuring acute angle between true tensor direction and estimated . A small indicates that the fiber direction is estimated accurately. We also introduce Metric 2 measuring the difference between true betweenneighbor angle and estimated betweenneighbor angle:
where are neighbors. A small leads to an accurate decision if two voxels belong to the same fiber.
We summarize the results in Table 1, including the mean estimates by averaging over
replications and the associated standard errors (in parentheses). From the result, we find that
SpDiST and DiST have an overall better performance in comparison to the nonspatial method. From Table 1, the SpDiST is more robust to noise, which may motivate a study on the robustness of tensor direction estimates based on different parameterization. However, the noise may have little effect on Metric 2, leading to the same fiber tracking results. Although the SpDiST and DiST have similar performance, however, the MCMCbased SpDiST provides a means to quantify the uncertainties of fiber tracking, unlike DiST.Metric  Noise () 

SpDiST  DiST  

0.1  0.09(0.007)  0.08 (0.006)  0.08(0.010)  
0.5  0.20(0.04)  0.12(0.010)  0.19(0.020)  
0.1  0.06 (0.005)  0.06(0.014)  0.06(0.012)  
0.5  0.24(0.05)  0.08(0.010)  0.09(0.015) 
5 Discussion
In the numerical study, we find that DiST and SpDiST have similar performances. However, the MCMCbased SpDiST provides a probabilistic means to quantify the result of fiber tracking. This provides some potentially important information for neuroscientists to understand brain anatomical connection. Furthermore, we also give a sensitivity analysis to the tuning parameter , addressing the issue raised by Kang and Li (Kang and Li, 2016).
Although the current methodologies might be sufficient for preliminary fiber tracking, there are still several issues. One problem is that the current methods focus on developing an imaging processing tool but not on scientifically and statistically explaining the outcomes (Lazar et al., 2016). However, proposing a statistical approach which characterizes factors affecting the outcomes might be critical in further studies, providing more insightful information in neuroscience. However, this is challenging because to incorporate covariates in the model and to properly combine the model to a fiber tracking algorithm are not straightforward. Another issue is crossing fibers. That is the single tensor model (Mori et al., 1999) fails to account for voxels where there are multiple fibers. Although it is assumed that increasing the resolution of the image may handle this issue, Schilling et al. (Schilling et al., 2017) give an unexpected result that increasing the resolution is not a solution. This needs to be rigorously studied with close interdisciplinary collaboration.
References
 Chung et al. (2011) Chung, H.W., Chou, M.C., and Chen, C.Y. (2011). Principles and limitations of computational algorithms in clinical diffusion tensor mr tractography. American Journal of Neuroradiology 32, 3–13.
 Datta et al. (2017) Datta, A., Banerjee, S., and Hodges, J. S. (2017). Spatial disease mapping using directed acyclic graph autoregressive (DAGAR) models. arXiv preprint arXiv:1704.07848 .
 Dryden et al. (2009) Dryden, I. L., Koloydenko, A., and Zhou, D. (2009). NonEuclidean statistics for covariance matrices, with applications to diffusion tensor imaging. The Annals of Applied Statistics pages 1102–1123.
 Kang and Li (2016) Kang, J. and Li, L. (2016). Discussion of “fiber direction estimation in diffusion MRI”. The Annals of Applied Statistics 10, 1162.
 Lazar et al. (2016) Lazar, N. A. et al. (2016). Discussion of “fiber direction estimation in diffusion MRI”. The Annals of Applied Statistics 10, 1160–1161.
 Mori (2007) Mori, S. (2007). Introduction to diffusion tensor imaging. Elsevier.
 Mori et al. (1999) Mori, S., Crain, B. J., Chacko, V. P., and Van Zijl, P. C. (1999). Threedimensional tracking of axonal projections in the brain by magnetic resonance imaging. Annals of Neurology: Official Journal of the American Neurological Association and the Child Neurology Society 45, 265–269.
 Niethammer et al. (2006) Niethammer, M., Estepar, R. S. J., Bouix, S., Shenton, M., and Westin, C.F. (2006). On diffusion tensor estimation. In 2006 International Conference of the IEEE Engineering in Medicine and Biology Society, pages 2622–2625. IEEE.
 Schilling et al. (2017) Schilling, K., Gao, Y., Janve, V., Stepniewska, I., Landman, B. A., and Anderson, A. W. (2017). Can increased spatial resolution solve the crossing fiber problem for diffusion MRI? NMR in Biomedicine 30, e3787.
 Schwartzman et al. (2016) Schwartzman, A. et al. (2016). Discussion of “fiber direction estimation in diffusion MRI”. The Annals of Applied Statistics 10, 1157–1159.
 Soares et al. (2013) Soares, J. M., Marques, P., Alves, V., and Sousa, N. (2013). A hitchhiker’s guide to diffusion tensor imaging. Frontiers in neuroscience 7,.
 Sporns et al. (2005) Sporns, O., Tononi, G., and Kötter, R. (2005). The human connectome: A structural description of the human brain. PLoS Computational Biology 1, e42.
 Wong et al. (2016) Wong, R. K., Lee, T. C., Paul, D., Peng, J., ADNI, et al. (2016). Fiber direction estimation, smoothing and tracking in diffusion MRI. The Annals of Applied Statistics 10, 1137–1156.
Comments
There are no comments yet.