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

TL;DR

This paper introduces a probabilistic fiber tracking method for diffusion MRI using a directed acyclic graph auto-regressive model of positive definite matrices, enhancing uncertainty quantification in tissue connection mapping.

Contribution

It proposes a novel directed acyclic graph auto-regressive model for positive definite matrices and applies it to probabilistic fiber tracking in diffusion MRI, addressing uncertainty quantification.

Findings

Demonstrates effectiveness through real data analysis.

Shows improved uncertainty quantification in fiber tracking.

Validates approach with numerical studies.

Abstract

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.

Tables

Table 1. Table 1 : Summary of simulation results based on Metric 1 and Metric 2 . The mean estimates by averaging over 50 50 50 replications and the associated standard errors (in parentheses) are summarized.

Metric	Noise ($\tau$)	Least Square	SpDiST	DiST
$d_1$	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)
$d_2$	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)

Equations

$${\bar{S}}_{mv}=S_{0v}\exp (-b{\mathbf{g}}_m^T{\mathbf{A}}_v{\mathbf{g}}_m).$$

$$\log S_{mv}=\log S_{0v}-b{\mathbf{g}}_m^T{\mathbf{A}}_v{\mathbf{g}}_m+{ϵ}_{mv},{ϵ}_{mv}\stackrel{iid}{\sim }\mathcal{N}(0,{\sigma }^2),$$

$${\mathbf{A}}_v|{\mathbf{A}}_u,u\in N(v)$$

$${\mathbf{A}}_v$$

$$f(\mathbf{X})=\frac{|\mathbf{X}{|}^{(k-p-1)/2}{e}^{-\mathrm{tr}([\mathbf{V}/k{]}^{-1}\mathbf{X})/2}}{2^{\frac{kp}{2}}|\mathbf{V}/k{|}^{k/2}{\Gamma }_p(\frac{k}{2})},$$

$$\mathcal{L}({\mathbf{A}}_v^{\ast }|.)$$

$$\mathcal{W}({\mathbf{A}}_v^{\ast }|{\bar{\mathbf{A}}}_v^{\ast },q)\prod _{u:v\in N(u)}\mathcal{W}({\mathbf{A}}_u|{\bar{\mathbf{A}}}_u^{\ast },q),$$

$$\log S_{mv}^{(r)}=\log {\bar{S}}_{mv}+E_{mv},E_{mv}\sim \mathcal{N}(0,{\tau }^2),$$

$$\arg \underset{{\mathbf{A}}_v}{\min }\sum _m||\log S_{mv}-\log S_{0v}-b{\mathbf{g}}_m^T{\mathbf{A}}_v{\mathbf{g}}_m|{|}^2.$$

$$d_1({\mathbf{m}}_v,{\hat{\mathbf{m}}}_v)=\text{arccos}(|{\mathbf{m}}_v^T{\hat{\mathbf{m}}}_v|),$$

$$d_2({\hat{\mathbf{m}}}_v,{\hat{\mathbf{m}}}_u)=|\text{arccos}(|{\hat{\mathbf{m}}}_v^T{\hat{\mathbf{m}}}_u|)-\text{arccos}(|{\mathbf{m}}_v^T{\mathbf{m}}_u|)|,$$

Full text

Keywords: Directed Acyclic Graph Auto-Regressive Model, Diffusion MRI, Fiber Tracking, Positive Definite Matrices

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 $m\in\{1,2,...,M\}$ measurements at voxel $v\in\{1,2,...,n\}$, denoted as $S_{mv}\in\mathbb{R}^{+}$. The measurements $S_{mv}$ are used to estimate the diffusion tensor $\bm{A}_{v}$ for voxel $v$. $\bm{A}_{v}$ is a $3\times 3$ positive definite matrix interpreted as covariance matrix of a local Brownian motion, indicating the local tensor direction. The goal is to use the measurements $S_{mv}$ to obtain tensor direction information from $\bm{A}_{v}$.

The noiseless signal intensity $\bar{S}_{mv}$ can be expressed in terms of $\bm{A}_{v}$ (Mori, 2007) as

[TABLE]

In this expression, $S_{0v}$, $b$, and $\bm{g}_{m}$ are non-diffusion weighted intensity, scale parameter, and $3\times 1$ unit-norm gradient vector, respectively. A detailed explanation of these three quantities can be found in Soares et al. (Soares et al., 2013). Given $\bm{A}_{v}$, $\bar{S}_{mv}$ can be understood as the probability intensity of the Gaussian motion when measuring at direction $\bm{g}_{v}$. For statistical modeling, $S_{0v}$, $b$, and $\bm{g}_{m}$ can simply be understood as fixed and known values.

The observations $S_{mv}$ are noisy realizations of $\bar{S}_{mv}$. The Rician distribution (Wong et al., 2016) is reasonable for modeling $S_{mv}$ but it causes computational issues. Here, we assume that the noise is a multiplier to the $\bar{S}_{mv}$ and it follows a lognormal distribution. The model is

[TABLE]

where $\epsilon_{mv}$ is the noise following a mean-zero normal distribution with variance $\sigma^{2}$.

To induce spatial smoothness, an image is treated as a directed graph whose nodes are voxels and whose directed edges are from node $v$ to nodes in $N(v)$. Following Datta et al. (Datta et al., 2017), we use the directed acyclic graph (directed and no loops) to construct to $N(v)$, leading to a valid joint density function of $[\bm{A}_{1},\bm{A}_{2},...,\bm{A}_{n}]$. In particular, we assume that the conditional mean of $\bm{A}_{v}$ is the average of its neighboring tensors, denoted as $\mathbb{E}[\bm{A}_{v}|\bm{A}_{u},u\in N(v)]=\frac{1}{|N(v)|}\sum_{u\in N(v)}\bm{A}_{u}$, where $N(v)$ is a set containing neighboring voxel indices of voxel $v$, and $|N(v)|$ is the set size.

In a directed acyclic graph, we have at least one voxel $v$ whose $N(v)$ is an empty set. For $N(v)$ is an empty set, we assume that $\bm{A}_{v}$ follows a Wishart distribution with mean matrix $\bm{I}$ and degrees of freedom $k$. Otherwise, conditional on $\bm{A}_{u},u\in N(v)$, we assume that $\bm{A}_{v}$ follows a Wishart distribution with mean matrix $\bar{\bm{A}}_{v}=\frac{1}{|N(v)|}\sum_{u\in N(v)}\bm{A}_{u}$ and degrees of freedom $k$. The model is

[TABLE]

In Equation (2), to preserve the designed mean realizations, we parameterize the Wishart distribution for $\bm{X}\sim\mathcal{W}(\bm{V},k)$ to have $\mathbb{E}\bm{X}=\bm{V}$. The probability density function is

[TABLE]

where $p$ is the matrix dimension and $\Gamma_{p}({\frac{k}{2}})$ 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 $5\times 5$ image is constructed.

2.1.1 MCMC Algorithm

We use MCMC for model fitting. We give $k\sim\mathcal{U}(3,50)$ and $\sigma^{-2}\sim\mathcal{GA}(0.01,0.01)$. The primary challenge in the MCMC algorithm is to sample the posterior of $\bm{A}_{v}$. Since the prior of $\bm{A}_{v}$ is not conjugate, we sample it using single-site Metropolis-Hastings sampling with Wishart distribution $\mathcal{W}(\bm{A}_{v}^{\prime}|\bm{A}_{v},q)$ as the proposal distribution. The algorithm is described below:

Candidate Generation:

Generate a candidate sample $\bm{A}_{v}^{\prime}$ using $\bm{A}_{v}^{\prime}\sim\mathcal{W}(\bm{A}_{v}^{\prime}|\bm{A}_{v},q)$;

Acceptance Rate:

Calculate the acceptance rate $r(\bm{A}_{v}^{\prime},\bm{A}_{v})=\frac{\mathcal{L}(\bm{A}_{v}^{\prime}|.)\mathcal{W}(\bm{A}_{v}|\bm{A}_{v}^{\prime},q)}{\mathcal{L}(\bm{A}_{v}|.)\mathcal{W}(\bm{A}_{v}^{\prime}|\bm{A}_{v},q)}$, where

[TABLE]

where $\bar{\bm{A}}_{u}^{*}=\frac{\sum_{u\in N(u)/v}\bm{A}_{u}+\bm{A}_{v}^{*}}{|N(v)|}$. $\mathcal{W}(.|\bm{A},\nu)$ and $\mathcal{N}(.|\mu,\sigma^{2})$ are the density functions of Wishart distribution and normal distribution, respectively.

Decision:

Generate $u\sim\mathcal{U}(0,1)$. If $u<r(\bm{A}_{v}^{\prime},\bm{A}_{v})$, accept $\bm{A}_{v}^{\prime}$.

The acceptance rate can be tuned by the degrees of freedom $q$, where smaller $q$ leads to smaller acceptance rate. We tune $q$ to make the acceptance rate around $0.4$.

We use Metropolis-Hastings algorithm with log-normal random walk as proposal distribution to update the degrees of freedom $q$ and use Gibbs sampling to update $\sigma^{2}$ based on its posterior $[\sigma^{-2}|.]\sim\mathcal{GA}(Mn/2+0.01,\sum_{m,v}(\log S_{mv}-\log S_{0v}+b\bm{g}_{m}^{T}\bm{A}_{v}\bm{g}_{m})^{2}/2+0.01)$.

2.2 Probabilistic Fiber Tracking Algorithm

We collect the $T$ MCMC samples of $\bm{A}_{v}$, denoted as $\{\bm{A}_{v}^{(t)}:t=1,2,...,T\}$. For each sample, we compute the principal eigenvector of $\bm{A}_{v}^{(t)}$, denoted as $\bm{m}_{v}^{(t)}$. For each posterior draw, we use $\bm{m}_{v}^{(t)}$ 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 $u$, we search neighboring voxels and compute the two angles: $\delta_{uv}=\text{arccos}\left(\frac{\bm{m}_{v}^{T}\bm{m}_{u}}{|\bm{m}_{v}||\bm{m}_{u}|}\right)$ is the angle between the two tensor directions ($\bm{m}_{u}$ and $\bm{m}_{u}$) and $\theta_{uv}=\text{arccos}\left(\frac{\bm{m}_{u}^{T}\bm{l}_{u,v}}{|\bm{m}_{u}||\bm{l}_{u,v}|}\right)$ is the angle between the current tensor ($\bm{m}_{u}$) and between-voxel direction ($\bm{l}_{u,v}$). See Figure 2 for details. We move to the voxels with $\theta<C$ and $\delta<C$. 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 $T$ possible fibers. We summarize $K$ distinct patterns from the outputs and calculate the associated probability for pattern $k\in\{1,2,...,K\}$ defined as $\frac{T_{k}}{T}$, where $T_{k}$ is the frequency of the pattern $k$. 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 $50\times 20$ voxels and $M=15$ measurements. A detailed description can be found in Dryden et al. (Dryden et al., 2009, Section 6). We sample $2000$ MCMC samples after discarding $3000$ samples as burn-in and thin the MCMC chain by retaining every $100$ 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 $\bm{A}_{v}$ and compute the corresponding principal eigenvector $\bm{m}_{v}$. To quantify the uncertainties of tensor direction estimation in each voxel, we overlay the MCMC samples on a $50\times 20$ map (Figure 3). In Figure 3, 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 $C$ (Chung et al., 2011), we consider $C$ ranging from $18^{o}$ to $28^{o}$. There are two distinct patterns (Pattern A and Pattern B) for $C\in[18^{o},28^{o}]$ (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 $C$.

Kang and Li (Kang and Li, 2016, Section 3) show that the FACT algorithm hinges on the tuning parameter $C$. It requires a sensitivity analysis to explore the impact of $C$. Here, we give a sensitivity analysis. We apply the FACT algorithm with $C=18+0.01\times s$ and $s=\{0,1,2,...,1000\}$. Since there are only two distinct patterns, we report the probabilities of Pattern B with different thresholds $C$ (Figure 5). We find that the result is sensitive to the choice of $C$ unless it is ranging from $24^{o}$ to $28^{o}$.

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 $8\times 7\times 2$ voxels where the three digits represent the dimension of $x$-axis, $y$-axis, and $z$-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 $\bm{m}_{v}$ is tangent to the arc. The noiseless signal in the example data is given as $\bar{S}_{mv}=S_{0v}\exp[-b(\bm{g}_{m}^{T}\bm{m}_{v})^{2}]$, a reparameterized model of Model 1 (Wong et al., 2016).

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 $\tau=0.1,0.5$. That is, the simulated data for each replication ($r=50$) is

[TABLE]

where $\log S_{mv}^{(r)}$ is the simulated signals for each replication ($r$), $\log S_{mv}$ is the logarithm signal from the example data, and $E_{mv}$ is simulated noise.

4.2 Simulation Details

We construct $N(v)$ as described before. We use the posterior mean estimate of SpDiST to compare with the estimates of alternatives. We compute posterior mean of $\bm{A}_{v}$ based on $2000$ MCMC samples after $3000$ 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 $\bm{A}_{v}$ via

[TABLE]

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

[TABLE]

a metric measuring acute angle between true tensor direction $\bm{m}_{v}$ and estimated $\hat{\bm{m}}_{v}$. A small $d_{1}(\bm{m}_{v},\hat{\bm{m}}_{v})$ 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:

[TABLE]

where $u,v$ are neighbors. A small $d_{2}(\bm{m}_{v},\hat{\bm{m}}_{u})$ 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 $50$ 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.

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 $C$, 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.