Probabilistic Diffusion MRI Fiber Tracking Using a Directed Acyclic Graph Auto-Regressive Model of Positive Definite Matrices

06/15/2019 ∙ by Zhou Lan, et al. ∙ 0

Diffusion MRI is a neuroimaging technique measuring the anatomical structure of tissues. Using diffusion MRI to construct the connections of tissues, known as fiber tracking, is one of the most important uses of diffusion MRI. Many techniques are available recently but few properly quantify statistical uncertainties. In this paper, we propose a directed acyclic graph auto-regressive model of positive definite matrices and apply a probabilistic fiber tracking algorithm. We use both real data analysis and numerical studies to demonstrate our proposal.



There are no comments yet.


page 8

page 9

page 10

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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, voxel-wise diffusion directions are estimated using diffusion-weighted 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 optimization-based 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 auto-regression 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 non-diffusion weighted intensity, scale parameter, and

unit-norm 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



is the noise following a mean-zero 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

and degrees of freedom

. Otherwise, conditional on , we assume that follows a Wishart distribution with mean matrix and degrees of freedom . The model is


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 x-y axis, we first order the voxels according to their coordinates of the y-axis, then next we order the voxels according to their coordinates of the x-axis. 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.

Figure 1: The construction of a directed acyclic graph based on an undirected graph. The left panel is the undirected graph of a image. The right pannel is the corresponding directed acyclic graph after modifying the edges.

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 single-site Metropolis-Hastings 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


where . and are the density functions of Wishart distribution and normal distribution, respectively.


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 Metropolis-Hastings algorithm with log-normal 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 between-voxel 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.

Figure 2: Two angles between two voxels. This figure is modified from Chung et al. (Chung et al., 2011, Figure 3).

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 burn-in 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: For each voxel, the MCMC samples of are overlaid on the its location on a map. For each voxel, we plot . The voxels with heterogeneous directions have large uncertainties. Otherwise, there are small uncertainties.

Figure 3 only provides voxel-wise uncertainties but our fully-Bayesian approach can propagate spatial uncertainties through to uncertainty in fiber tracking. In this way, the MCMC-based 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


Figure 4: The consecutive orange arrows construct a fiber. Two patterns are identified, where the left pattern and the right pattern are denoted as Pattern A and Pattern B, respectively.

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 .

Figure 5: The probability of Pattern B varies with different threshold . When is ranging from to the probability is insensitive to .

4 Numerical Study

4.1 Data Description

In this section, we use synthetic diffusion-weighted 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).

Figure 6: The tensor directions (left panel) and underlying fibers (right panel) of the example data .

To mimic low-quality images with signal noise, we further add noise on the log scale simulated from a mean-zero normal distribution with standard deviation

. That is, the simulated data for each replication () is


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 burn-in. In comparison, we compare our method to DiST. In addition, we also compare our method to a non-spatial 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 between-neighbor angle and estimated between-neighbor 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 non-spatial 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 MCMC-based SpDiST provides a means to quantify the uncertainties of fiber tracking, unlike DiST.

Metric Noise ()
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)
Table 1: Summary of simulation results based on Metric 1 and Metric 2. The mean estimates by averaging over replications and the associated standard errors (in parentheses) are summarized.

5 Discussion

In the numerical study, we find that DiST and SpDiST have similar performances. However, the MCMC-based 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.


  • 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 auto-regressive (DAGAR) models. arXiv preprint arXiv:1704.07848 .
  • Dryden et al. (2009) Dryden, I. L., Koloydenko, A., and Zhou, D. (2009). Non-Euclidean 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). Three-dimensional 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.