Cortico-Cortical,Cortico-Striatal,and Cortico-Thalamic White Matter Fiber Tracts Generated in the Macaque Brain via Dynamic Programming

Abstract

Probabilistic methods have the potential to generate multiple and complex white matter fiber tracts in diffusion tensor imaging (DTI). Here, a method based on dynamic programming (DP) is introduced to reconstruct fibers pathways whose complex anatomical structures cannot be resolved beyond the resolution of standard DTI data. DP is based on optimizing a sequentially additive cost function derived from a Gaussian diffusion model whose covariance is defined by the diffusion tensor. DP is used to determine the optimal path between initial and terminal nodes by efficiently searching over all paths, connecting the nodes, and choosing the path in which the total probability is maximized. An ex vivo high-resolution scan of a macaque hemi-brain is used to demonstrate the advantages and limitations of DP. DP can generate fiber bundles between distant cortical areas (superior longitudinal fasciculi, arcuate fasciculus, uncinate fasciculus, and fronto-occipital fasciculus), neighboring cortical areas (dorsal and ventral banks of the principal sulcus), as well as cortical projections to the hippocampal formation (cingulum bundle), neostriatum (motor cortical projections to the putamen), thalamus (subcortical bundle), and hippocampal formation projections to the mammillary bodies via the fornix. Validation is established either by comparison with in vivo intracellular transport of horseradish peroxidase in another macaque monkey or by comparison with atlases. DP is able to generate known pathways, including crossing and kissing tracts. Thus, DP has the potential to enhance neuroimaging studies of cortical connectivity.

Introduction

The advent of diffusion tensor imaging (DTI) technology is enabling clinicians and neuroscientists to examine normal and impaired connectivity in the mammalian brain (Basser et al., 2000; Conturo et al., 1999; Mori et al., 1999; Pace and Pierpaoli, 1963; Parker et al., 2002a; Poupon et al., 2000). DTI encapsulates the directionality of water diffusion along fibers constituting tracts or bundles between these structures (Basser et al., 1994a; Le Bihan et al., 1989; Pierpaoli et al., 1996; Xue et al., 1999). The diffusion tensor, D, is a 3×3 tensor, that is, a symmetric matrix whose eigenvectors and eigenvalues define an ellipsoid that can describe the diffusive characteristics of water at each voxel in the image despite being at a lower resolution than that of the individual axons in the tracts (Basser et al., 1994a; Beaulieu, 2002). The principal eigenvector associated with the largest eigenvalue represents the direction of fastest diffusion of water in the volume. Application of these essential properties of water diffusion have led to the development of a wide range of methods for elucidating connectivity of gray matter structures and segmentation of white matter (WM) in the past decade (Basser and Jones, 2002; Johansen-Berg and Rushworth, 2009; Le Bihan et al., 2001; Mori and van Zijl, 2002; Ramnani et al., 2004; Tournier et al., 2011).

Fiber tracking methods can be generally classified as either deterministic (Basser et al., 2000; Conturo et al., 1999; Hageman et al., 2009; Mori et al., 1999; Poupon et al., 2000) or probabilistic (Acar and Yörük, 2009; Behrens et al., 2003; Friman and Westin, 2005; Jbabdi et al., 2007; Jonasson et al., 2005; Sherbondy et al., 2008). However, both approaches have limitations such as inability to discriminate crossing fibers, low signal-to-noise ratio (SNR), and poor tracking in regions of low anisotropy (Jbabdi and Johansen-Berg, 2011; Lazar and Alexander, 2003; Le Bihan et al., 2006; Moldrich et al., 2010; Mukherjee et al., 2008; Tournier et al., 2002). In particular, deterministic methods rely on tracking the principal eigenvector and thus do not continue tracking into regions in which anisotropy is typically low, and the principal eigenvector is a less reliable indication of fiber orientation. More recently, a new class of probabilistic methods has emerged that can overcome some of these shortcomings. Posing the problem as an optimization problem enables computation of a “shortest path” between chosen initial and terminal points that globally minimizes a sequentially additive energy constraint defined by the tensor in the spirit of the classical Djikstra's algorithm (Everts et al., 2009; Fout et al., 2005; Iturria-Medina et al., 2007; Lal, 2004; Lifshits et al., 2009; Merhof et al., 2006a, 2006b; Poynton et al., 2005; Zalesky, 2008; Zalesky and Fornito, 2009). These assign a probability distribution to the local orientation of fibers at each voxel, and use path finding methods to compute the optimal path between two regions.

This article focuses on the use of dynamic programming (DP) as the optimization procedure by implementing the algorithm originally developed by Lal (2004) and Poynton et al. (2005) in the widely used MRIStudio software (Jiang et al., 2006). In particular, the probabilistic law for diffusion in the Lal algorithm is modified to ensure that optimal paths through isotropic tensors are not affected by the magnitude of the tensor and then, a modification of cost per unit length as a connectivity metric, originally proposed by Poynton et al. (2005), is introduced to facilitate pruning of fibers that are not anatomically valid. An ex vivo macaque hemi-brain is used for validation, and the advantages of this approach as an alternative to using phantom data are discussed.

Methods

Probabilistic labeling of paths

By defining a node as the voxel center and a path as an initial node x ₀ and a set of direction vectors \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ v_0 , v_1 , \cdots , v_{N - 1} \} $$ \end{document} connecting nodes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ x_1 , \ldots , x_N \} $$ \end{document} , a probability distribution can be induced on fiber tract orientation at each node, enabling the computation of the probability of a path as the product of the probabilities of the transitions between subsequent nodes along the path. A voxel j is in the neighborhood of voxel i if j is immediately adjacent to i in the 26-connected sense. There is a direct transition from node x_i to node x_j , representing the centers of voxels i and j, respectively, if and only if j is in the neighborhood of i. The state space of all nodes should be defined by appropriate thresholds on scalar quantities such as fractional anisotropy, ignoring voxels below a specified value. DP can then be used to compute the maximal probability path linking two nodes in the state space and the basic algorithm extended to compute K nonintersecting paths between a set of initial and terminal nodes.

Since evidence suggests that the directionality of the diffusion depends on the orientation of axonal fibers, it is reasonable to assume that the orientation of fibers follows the same Gaussian distribution as that of the diffusion of water molecules (Alexander et al., 2000). With this assumption, the problem of tracking fibers is reduced to the problem of computing the path between two nodes in a graph that minimizes a cost function determined by the eigenvalues of the covariance representation of the quadratic form in a sequentially additive quadratic cost. More sophisticated probability models characterizing diffusion have recently been developed (Friman and Westin, 2005; Sherbondy et al., 2008; Tuch, 2004) and can be incorporated into this method as long as the probability distribution for the diffusion at a particular voxel remains locally defined.

It is essential that the probability of Gaussian diffusion over unit time be the same for isotropic diffusion; so, the diffusion tensor at voxel i, D_i , is normalized by its trace Tr(D_i ), that is, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\bar{D}_i = D_i / { \rm Tr} ( D_i )$$ \end{document} (Frandsen et al., 2007). Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( \lambda_1^i , \lambda_2^i , \lambda_3^i )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( e_1^i , e_2^i , e_3^i )$$ \end{document} be, respectively, eigenvalues and unit eigenvector of D_i such that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$D_i = V_i^t L_iV_i$$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$L_i = { \rm diag} ( \lambda_1^i , \lambda_2^i , \lambda_3^i )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$V_i = [ e_1^i , e_2^i , e_3^i ]$$ \end{document} . If \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\bar{L}_i = { \rm diag} ( \bar{ \lambda}_1^i , \bar{ \lambda}_2^i, \bar{ \lambda}_3^i )$$ \end{document} where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\bar{ \lambda}_k^i = \lambda_k^i / ( \lambda_1^i + \lambda_2^i + \lambda_3^i ) = \lambda_k^i / { \rm Tr} ( D_i )$$ \end{document} , it can be shown that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \rm det} ( \bar{D}_i ) = { \rm det} ( \bar{L}_i )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\bar{D}_{ i}^{ - 1} = V_{ i} \bar{L}_{ i}^{ - 1}V_{ i}^t$$ \end{document} .

We need to define the probability associated with a transition between connected nodes x_i and x _i+1 to have a Gaussian distribution with covariance matrix \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\bar{D}_i$$ \end{document} . Then, the probability of an N-length path π_N(x ₀, x_N ) is given by the product of the quadratic form of the individual directional vectors d_j between nodes along the path as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}P [ \pi_N ( x_0 , x_N ) ] = \prod_{i = 0}^{N - 1} \frac {1} {( 2 \pi)^{3/2} \mid \bar {D} _i \mid ^ {1/2}} \exp \left( - \frac {d_j ^t \bar {D}_i^{- 1} d_j} {2} \right) \tag{1} \end{align*} \end{document}

that implicitly assumes diffusion between any adjacent voxel irrespective of length in a fixed time of τ=0.5 (Alexander et al., 2000).

By considering logarithms, the maximum probability N-length path π* (x ₀, x_N ) is given by the minimization of the negative log likelihood \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{\pi}^* ( x_0 , x_N ) = & {{\rm argmin} \atop \pi \in P ( x_0 , x_N )} \mathop \sum_{i = 0}^{N - 1} \left\{\mid \mid d_j \mid \mid^2 \left( \frac {\langle \bar{d}_j , e_1^i \rangle^2}{\bar{\lambda}_1^i} + \frac {\langle \bar{d}_j , e_2^i \rangle^2}{\bar{\lambda}_2^i} + \frac {\langle \bar{d}_j , e_3^i \rangle^2}{\bar{\lambda}_3^i} \right) + \ln (\bar {\lambda}_1^i \bar{\lambda}_2^i \bar{\lambda}_3^i ) + 3 \ln (2 \pi) \right\} \tag{2} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\bar{{d}}_j = {d}_j / \mid \mid {d}_j \mid \mid$$ \end{document} is the unit directional vector, and P(x ₀, x_N ) is the set of all paths between nodes x ₀ and x_N ; see Lal (2004) for proof. Note that the constant term is retained to ensure that the transition cost is positive. The path \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\pi^* = \{ x_0 , x_1 , \ldots , x_N \} $$ \end{document} that minimizes this function is the optimal path.

Dynamic programming

The problem of determining the optimal path reduces to one of computing the minimum energy path linking initial and terminal nodes in the state space. Let S be the finite collection of nodes of size ||S||=M and define c^k (i, j) as the cost of the transition from \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i \in S$$ \end{document} to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j \in S$$ \end{document} at time k. If the cost is additive over the length of the path, and the optimal path between two points is assumed to pass through no more than N nodes, then there may be as many as M^N paths linking two nodes, assuming the most complex case in which all nodes are connected to each other, that is, have a valency of M. Carrying out a brute force search would require considering the cost associated with each of these paths, and such an approach would very quickly become computationally impractical as N and M become large. DP overcomes this problem, as it reduces the complexity of the search to order of NM ² by taking advantage of the fact that the cost function being minimized is sequentially additive. Note that while most practical implementations simplify the graph such that the connectedness or valency (V) of each node is less than M, for example, 6 or 26, the brute force search's performance will be in the order of MV^N , which is still computationally impractical as N gets large.

There are many ways of implementing DP. The present implementation is described in Algorithm 1 next. Let H_x be the set of nodes {n} such that a direct transition between x and n exists, the cost c^k (x_i , x_j ) for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x_j \in H_{x{_i}}$$ \end{document} following the probabilistic model described by Equation (2), and c^k (x_i , x_j )=∞ for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x_j \,\notin H_{x{_i}}$$ \end{document} . The optimal N-length cost J₀ (s) from s to t is given by the final step of the algorithm evaluated at i=s. Note that the algorithm permits the degenerate move from a node i to itself with cost c^k (i, i)=0 for all i, thus allowing paths with length less than N to be included.

Algorithm 1. Initialize: J_k (i)=∞ i≠t, J_k (t)=0 for all k; S_N ={t};

For every k=N−1 to 0,

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$S_k = \{ i \mid i \in H_j , j \in S_{k + 1} \} $$ \end{document}

;

set c^k (i, j) from Equation (2),

and

set

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$J_k ( i ) = \min\limits_{j \in \{ S_{k + 1} \cap H_i \} } \{ c^k( i , j ) + J_{k + 1} ( j ) \} , i \in S_k$$ \end{document}

end

Here, the state space S_k is defined as the subset of nodes that can be reached from the initial node in k steps with a finite cost. This further optimizes the search pattern by only considering the restricted state spaces S_k ⊂ S at each iteration. Thus, many of the M² paths at each iteration can be ignored without compromising the condition of optimality.

It follows that if the cost of the optimal k length paths from the initial node to all nodes in the state space is known, then the optimal k+1 length paths to all nodes in the state space can be found by examining at most M ² paths. So, since k is iterated from 0 to N−1, DP computes the optimal path by considering NM ² or fewer paths; instead of the M^N paths, a brute force search would be required.

Computing a single optimal path between two regions has little practical use, as any two brain structures are linked by a bundle comprising multiple fibers. Thus, there is a need to find a most probable set of distinct fibers connecting the two regions. This is achieved by altering the state space as follows. Given an initial set of nodes S and a terminal set of nodes T, introduce two “dummy nodes” s and t into the state space, such that there is a transition from s to every node in S in the first time step with cost 0, and there is a transition from every node in T to node t with cost 0, at all time steps. Hence, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$c^0 ( s , n ) = 0 , \forall n \in S$$ \end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$c^k ( n , t ) = 0 , \forall n \in T , \forall k$$ \end{document} .

The path formed by removing the first and last arcs from the optimal path between nodes s and t is the optimal path between the sets of nodes S and T; see Lal (2004) for proof. So, the introduction of the “dummy nodes” into the state space allows for the computation of the optimal path between an initial set of nodes and a terminal set of nodes. So, to extend the search process to find K distinct paths between two sets of nodes, the DP algorithm is performed via a modification in the search space described in Algorithm 2 next.

One iteration of the DP algorithm gives the optimal path \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\pi^* = \{ s = x_{ - 1} , x_0 , \cdots , x_{N - 1} , x_N = t \} $$ \end{document} . Let H_x represent the neighborhood of node x. For all nodes x_i , 0≤i≤N−1 along the optimal path π*, set H_{x_i} =Ø. This step cuts off all the connections from the nodes along the path to all its neighboring nodes, effectively removing the nodes from the search space. The next iteration of the DP algorithm directly follows, which will find the most probable path in the modified search space, a path that is completely distinct from the previous paths found.

Algorithm 2. Determine the optimal path from every node in S to x_N , and let the optimal path from x ₋₁ to x_N be p ₁. Let F be a list of paths, initially empty.

Initialize p=p ₁, k′=1 and F[1]=p.

While an alternative to p exists and k′≤K

For every

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x \in \{ x_0 , \cdots , x_{N - 1} \} $$ \end{document}

update H_{x_i} =Ø

end

set p equal to the optimal path from x ₋₁ to x_N , k′=k′+1 and F[k′]=p.

end

The number of fibers K is a free parameter that can be empirically used. This can be done by generating a large number of fibers and selecting the first K generated fibers based on when the fibers begin to deviate from the bundle's true geometry and location. Unfortunately, this information may not be known a priori. If one chooses K sufficiently small, it is much less likely that a fiber generated will deviate to a great extent from the most optimal fiber in the bundle.

Connectivity metric

Following Lal (2004), Poynton et al. (2005) proposed a measure of the fiber bundle containing K fibers via \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$ \frac { 1 } { K } \sum \nolimits_ { k \prime = 1 } ^K J ( F [ k^ \prime ] ) / N_f$$ \end{document} , where J(F[k′]) is the cost of a fiber, F[k′], computed from Equation (2) and N_f is the number of nodes along F[k′]. However, a more meaningful measure can be obtained from J(F[k′])/L(F[k′]), where L(F[k′]) is the length of the fiber, that can be construed as a connectivity metric (Jackowski et al., 2004; Marigonda and Orlandi, 2011; Merhof et al., 2006b; Parker et al., 2002b; Pichon et al., 2005).

Data acquisition and processing via MRIStudio

Following Oishi et al. (2011), one rhesus monkey (Macaca mulatta) brain was perfusion fixed using 4% paraformaldehyde and scanned using a 4.7T Bruker scanner. A three-dimensional multiple spin-echo diffusion tensor sequence was used to acquire a set of diffusion-weighted images in seven linearly independent directions (Zhang et al., 2003). Parameters were as follows: echo time (TE)=32.5 ms, TR=0.7 sec, field of views (FOV)=80×58×60 mm (zero-filled to a data matrix=186×256×192 with nominal resolution 0.312×0.312×0.312 mm³), and b=1000 sec/mm². Co-registered T2-weighted images with the same FOV were also acquired with the fast spin-echo sequence via TE=15 ms and TR=1 sec (Huang et al., 2006). A left-sided hemi-brain volume was then created with the data resampled to 0.2×0.2×0.2 mm³ and zero-filled matrix size 240×480×320. The dataset was transferred to a workstation and processed using MRIStudio (Jiang et al., 2006). After inspection, the six independent elements of the 3×3 diffusion tensor were calculated for each pixel by using the multivariate linear-fitting method. After diagonalization, three eigenvalues and three eigenvectors were obtained (Basser et al., 1994b) and used to calculate the fractional anisotropy. DP was implemented in MRIStudio to take advantage of the editing features such as Boolean operations developed for fibers generated by the Fiber Assignment by Continuous Tracking (FACT) algorithm. Three-dimensional screenshots of fibers were taken in CAWorks (www.cis.jhu.edu/software/caworks).

Anatomical regions of interest

Anatomical regions of interest (ROIs) were manually segmented in the T2 volume using the Schmahmann and Pandya (SP) atlas (Schmahmann and Pandya, 2006). The T2 modality was chosen, as it had a similar contrast as the atlas. It should be noted that electrophysiologic stimulation was not used to confirm anatomical correspondence and that with the exception of the putamen, it was necessary to include WM subjacent to the gray matter structure.

(1) Walker's area 46 encompasses the WM lying subjacent to the principal sulcal cortex (area 46) along its entire rostro-caudal extent.

(2) Dorsal Cingulum Bundle (CB) consists of the CB at the level of the cingulate motor area just rostral to the intersection of the central sulcus with the medial face of the hemisphere. The CB at this level is surrounded by the cortex of area 23, 30, and 29 medially and by cortical WM laterally.

(3) Ventral CB is situated in the WM protrusion that underlies the most caudal portion of the presubiculum and the transitional cortices (areas 29 and 30) forming the caudomedial lobule of Goldman–Rakic et al. (1984).

(4) Mediodorsal Thalamus is located in the ventral portion of the subcortical bundle (SB) just caudal to the decussation of the anterior commissure. This portion of SB forms a portion of the boundary between the rostral globus pallidus and the body of the caudate nucleus.

(5) Principal Sulcus Dorsal Bank is located in the cortex lining the dorsal bank of the principal sulcus (area 46), encompassing the middle third of the rostral-caudal extent of area 46 and stretching from the dorsal border with area 9 along the bank of the sulcus approximately two-third the distance to the depth of the sulcus.

(6) Principal Sulcus Ventral Bank is situated in the middle third of the rostrocaudal extent of the principal sulcal cortex (area 46) and stretches from the ventral lip of the sulcus to approximately two-third the distance to the depth of the sulcus.

(7) Frontal Superior Longitudinal Fasciculus II-Version A (SLFII-FrontalA) begins midway along the rostro-caudal extent of the principal sulcus and includes the cortex of the dorsal bank of the principal sulcus (area 46) and cortex on the lateral two-third of the superior frontal gyrus (area 9d). More caudally, the ROI extends into the cortex surrounding the superior spur of the arcuate sulcus (8Ad) and onto the adjacent cortex on the dorsal convexity (6d). Far caudally, this ROI is limited to the depths and medial bank of the superior arcuate sulcus (6d).

(8) Parietal Superior Longitudinal Fasciculus II-Version A (SLFII-ParietalA) is located in the cortex, forming the caudal half of the lateral bank of the intraparietal sulcus (area POa) and adjacent ventral convexity (areas PG and Opt).

(9) Frontal Superior Longitudinal Fasciculus II-Version B (SLFII-FrontalB) includes the nearly the entire dorsal bank of the principal sulcus (area 46), excluding only the rostral quarter. The ROI encroaches onto the adjoining dorsal convexity (area 9d), but only slightly. More caudally situated cortices (8Ad and 6d) were not included in order to generate kissing or crossing fibers with the arcuate fasciculus (AF; see below).

(10) Parietal Superior Longitudinal Fasciculus II-Version B (SLFII-ParietalB) includes nearly the entire lateral bank of the intraparietal sulcus (area POa), with only the rostral quarter not included, and the most caudal region of the cortex on the inferior parietal lobule (areas PG and Opt). Near the caudal tip of the lateral fissure, there is marginal encroachment of the ROI onto the cortex of the convexity (PFG). This ROI was intentionally made larger than that of SLFII-ParietalA to try to generate kissing or crossing fibers.

(11) Frontal Superior Longitudinal Fasciculus III (SLFIII-Frontal) is situated in the cortex forming the ventral bank of the inferior spur of the arcuate sulcus (areas 44, 6v) and in the caudal extension of premotor area 6v on the precentral gyrus.

(12) Parietal Superior Longitudinal Fasciculus III (SLFIII-Parietal) is located in the rostral parietal cortex just ventral to the rostral tip of the intraparietal sulcus (area PF) and encompassing the cortex between the intraparietal sulcus and the lateral fissure (areas PFG and PFop).

(13) Frontal Arcuate Fasciculus (AF-Frontal) consists of the cortex forming the ventral (area 8A) and dorsal banks (area 6d) of the superior spur of the arcuate sulcus. Note that the area 46/9 components of this ROI were intentionally excluded to generate kissing or crossing fibers.

(14) Temporal Arcuate Fasciculus (AF-Temporal) includes the cortex lining the dorsal bank of the superior temporal sulcus (area TPO) and caudal portions of the superior temporal gyrus (area Tpt) adjacent to TPO.

(15) Frontal Occipital Fasciculus (FOF-Frontal) extends rostrally from the approximate midpoint of the principal sulcus, where the ROI occupies the dorsal bank of the principal sulcus (46d). It extends onto the adjacent dorsal convexity (9d), encroaching on the adjacent medial cortex (9m). More caudally, the ROI includes the cortex surrounding the superior spur of the arcuate sulcus and bordering convexities (8Ad and 6d), as well as 46d. The most caudal portion of the ROI is present at the level where the principal sulcus ends and includes only the dorsal bank of the superior spur of the arcuate sulcus (6d).

(16) Frontal Occipital Fasciculus (FOF-PO) is located in the medial occipital lobe in the cortex lining the medial bank of the intraparietal sulcus (PO) and extending through the adjacent WM to include the cortex on the medial wall of the hemisphere (PGm). The most caudal origin of the FO (PO) was chosen as the second ROI in order to generate the full longitudinal extent of the FOF (as illustrated by Case 17 of SP).

(17) Fornix-Mammillary Bodies (MMB) includes the MMB of the thalamus.

(18) Fornix-Fimbria is located in the caudal most portion of the fimbria near the splenium of the corpus callosum. The fimbria at this level is situated just medial to the caudal extent of the lateral ventricle and to the caudal most, elongated portion of the caudate nucleus.

(19) Frontal Uncinate Fasciculus (UF-Frontal) is limited to a finger of WM lying subjacent to area 14 of the rostral orbitofrontal cortex in one coronal section (95) corresponding roughly to Section 20 of SP.

(20) Temporal Uncinate Fasciculus (UF-Temporal) also consists of a slender finger of WM located in the rostral temporal lobe between the amygdala and the lateral temporal cortices TE1, TE2 and extending ventrally to lie just subjacent to areas TL and 35. This ROI is restricted to one coronal section (164) that corresponds to Section 57 of SP.

(21) SP Motor Case 24 is situated in cortex lining the frontal operculum, or dorsal lip of the lateral fissure, in the precentral gyrus (areas 1 and 2).

(22) SP Motor Case 25 is located in the face representation cortex in the ventral precentral gyrus (area 4).

(23) SP Motor Case 26 is located in the hand motor area of the precentral gyrus (area 4).

(24) SP Motor Case 27 is situated in the motor trunk representation of the precentral gyrus (area 4).

(25) SP Motor Case 28 includes the dorsal portion of the precentral gyrus (area 4), corresponding to the foot motor area.

(26) SP Motor Case 29 is located on the medial surface of the superior frontal gyrus, just dorsal to the cingulate sulcus, in the rostral supplementary motor area associated with the face representation (area 6m).

(27) Putamen includes the entire rostro-caudal extent of the putamen. It extends from the rostral-most portion located ventrolateral to the anterior limb of the internal capsule (IC) to its most caudal portion lateral to the globus pallidus and dorsal to the tail of the caudate nucleus.

Delineation of fiber tracts with intracellular transport of horseradish peroxidase

A comparison with postmortem tracing data was performed via visual examination of histologically processed sections from a different Rhesus macaque monkey. This monkey had horseradish peroxidase (HRP) pellets implanted in the dorsal and ventral banks of the principal sulcus in order to identify fiber tracts connecting Walker's area 46 with other brain regions. Details of the surgical procedure, perfusion, and histological processing have been previously described (Selemon and Goldman-Rakic, 1985) (case 14). Eight HRP pellets were implanted in this monkey, resulting in a large injection site that encompassed the caudal two-third of the Walker's area 46 (Brodmann areas 9, 10) and encroached caudally on adjacent Brodmann areas 8 and medially on area 6 (see Figs. 8 and 9 in Selemon and Goldman–Rakic, 1985). This monkey also received an injection of tritiated amino acids in the posterior parietal cortex.

Results

Except where stated, the following results were obtained with a threshold of 0.4 for FA, which is consistent with Oishi et al. (2011).

Single transition cost

The cost of a transition from one node to seven possible adjacent nodes for different types of diffusion tensors is both tensor and length dependent (Table 1). This cost assessment is consistent with the Gaussian diffusion model first formulated by Alexander et al. (2000).

Table 1.

Transition Costs from One Node to Seven Adjacent Nodes via the Black Directional Vector (Top Row) Based on the Tensor Shape and Orientation at the Start Node (First Column)


1.9353	10.6853	10.4853	11.9353	20.6853	11.9353	21.9353
2.6735	6.2449	11.2449	7.6735	16.2449	12.6735	17.6735
3.1629	4.8296	11.4962	6.4962	14.8296	13.1629	16.4962
3.8363	3.8363	11.6140	6.0585	13.8663	13.8363	16.0585
6.3103	6.3103	10.6853	3.1853	16.3103	16.3103	13.1853
8.1021	8.1021	7.1018	10.3523	8.4723	8.4723	4.5790

Caveats due to high anisotropy bias, resolution, and ROI proximity

The number of fibers (K) should be as large as possible to resolve problems that arise due to fibers with high anisotropy bias and to closely aligned fiber bundles. These are illustrated by generation of SLFII using SLFII-FrontalA and SLFII-ParietalA ROIs with K=20 (Fig. 1A) in a comparison to generation of this same fiber bundle with K=200 (Fig. 1B). When K=20, three distinct paths are generated: one through the CB, one along the expected SLFII, and one through the IC and the SB. The deviant bundles through the CB and IC probably arise due to the strong projections from both frontal and parietal cortices to the CB and IC (Selemon and Goldman-Rakic, 1988). At this voxel resolution, DP is confounded by the close proximity of these frontal and parietal bundles in the CB and the IC. It should be noted, however, that the lateral bundle, that is, the one which courses through the SLFII, has lower cost (indicated by the lighter color) than the other two. At larger values of K such as 200 (Fig. 1B), the lower cost or equivalently higher probability SLFII bundle is partially masked by the presence of a fourth bundle, not seen with the fewer fibers, which appears to be a part of the AF and exhibits higher anisotropy bias than the one through the CB. A strategy which utilizes the connectivity metric-based color code to identify fibers for editing and then employs the AND/NOT editing feature in MRIStudio to prune these anatomically inconsistent fibers results in a representation of the SLFII that is consistent with the classic Schmahmann and Pandya (2006) anatomical tracings (Fig. 1C). Thus, several factors can generate incorrect fibers: the bias to high anisotropy, the tendency to pass from one to the other of two closely aligned fiber bundles, and the tendency to prefer shortest paths. To compensate for these factors, it is important to generate as many fibers as possible (K=200) and to use the color-coded connectivity metric to prune incorrect fibers.

FIG. 1.

Parietal Superior Longitudinal Fasciculus II-Version A (SLFII-ParietalA) to Frontal Superior Longitudinal Fasciculus II-Version A (SLFII-FrontalA): (A) K=20 fibers; (B) K=200 fibers; (C) K=67, after pruning. Central Sulcus (CS), arcuate fasciculus (AF), cingulum bundle (CB) and subcortical bundle (SB) [internal capsule (IC)] are marked as reference points. Coronal views for 1C shown in last row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

The need for “way stations” in long association bundles

Some fiber bundles that exhibit a long traverse coupled with an extreme change of direction may require a “way station.” This is illustrated by the CB, a fiber bundle that connects many cortical areas with the hippocampal formation. Only connectivity between the dorsolateral prefrontal cortex (Walker's area 46) and the hippocampal formation was examined. In order to override the tendency for DP to prefer the shortest path between these two ROIs and thereby avoid the most caudal portion of the CB that curves around the splenium of the corpus callous, a third ROI, the dorsal CB was used as a way station (Fig. 2).

FIG. 2.

Walker's area 46 to CB Dorsal with K=49 and thence to CB Ventral with K=68. CS and the superior temporal gyrus (STG) are marked as reference points. Coronal views for Walker's area 46 to CB Dorsal shown in middle row and for CB Dorsal to CB Ventral in bottom row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

Comparison with in vivo tract tracing

Coronal sections through a macaque brain in which HRP had been placed in the dorsolateral prefrontal cortex (Walker's area 46) and anterogradely transported into the fiber bundles of the CB and SB were compared with DP-generated fibers (Fig. 3). HRP-labeled fibers were found in the fiber bundles of the SB in the anterior limb of the IC (Fig. 3A, B) and in the CB (Fig. 3C, D).

FIG. 3.

Automated fiber tracking of the SB and CB in comparison to the same tracts shown via horseradish peroxidase (HRP) tract tracing. Arrows in (A, B) point to the SB as it courses through the anterior limb of the IC. Arrow heads in (C, D) show the CB encapsulated by the posterior cingulate cortex.

Short- and long-association cortico-cortical bundles

DP-generated fibers based on ROIs in the ventral and dorsal banks of the principal sulcus formed a short, U-shaped association cortico-cortical bundle in the WM beneath the fundus of the principal sulcus (Fig. 4). Note that no fibers traversed the geodesic between the two ROIs despite their proximity.

FIG. 4.

Principal sulcus (PS) from ventral (v) to dorsal (d) banks (K=200). Coronal views shown in bottom row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

Long-association bundles for SLFIII and AF were generated (Fig. 5). In the latter case, it was necessary to use K=800 and to prune the inconsistent fibers that either due to high anisotropy bias or geodesic tendency were ranked higher; indeed, the highest rank of the fiber in the edited AF bundle (K=64) was 302. In addition, two bundles were generated with the SLFII-FrontalB and AF-Frontal as two separate start ROIs and the SLFII-ParietalB and AF-Temporal combined as one end ROI. The two edited bundles (K=104 and 82, respectively) are seen to merge, that is, cross or kiss as they approach the SLFII-ParietalB ROI (Fig. 6). In each case, K=400 was initially used with some of the fibers similar to those in Figure 5 edited out.

FIG. 5.

Top: SLFIII Frontal to Parietal (K=200). Middle: AF Frontal to Temporal (K=64). CS and STG are marked as reference points. Bottom: Coronal views of AF Frontal to Temporal; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

FIG. 6.

Two SLFII-FrontalB (light green, K=104) and AF-Frontal (light blue, K=82) bundles tracked to a combined SLFII-ParietalB and AF-Temporal. Contrast with the AF bundle from Figure 7. CS is marked as a reference point. Here, the color code for the connection metric is not used. Coronal views shown in bottom row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

Association bundles for the UF (Fig. 7), FOF (Fig. 8), and fornix (Fig. 9) were also generated. The latter two proved to be challenging for several reasons. Both have relatively lower FA than adjacent bundles (corpus callosum, CB for the FOF, and, in addition, the anterior commissure, IC, and external medullary lamina of the thalamus for the fornix). So, it was necessary to lower the FA threshold from 0.4 to 0.1, which reflects the fact that FOF carries multiple tracts between different cortices (see Chapter 27 in Schmahmann and Pandya, 2006). Then, waystations were needed to prevent shunting into adjacent bundles. For the FOF, three waystations at coronal sections 241, 255, and 270, which respectively corresponded to sections 93, 97, and 105 in the SP atlas, were used with the four initial values of K=1000, 400, 200, and 400 reduced to 31, 56, 95, and 349, respectively, after editing. The result is a bundle that is consistent with Figure 19-2 of Schmahmann and Pandya (2006). From the PO, the fibers begin in the SLFI, separate into FOF proper, and then, become sandwiched between the corpus callosum and corona radiata while coursing rostrally toward the frontal lobe; see pages 234–236, 268, 272–276, and chapter 27 in Schmahmann and Pandya (2006).

FIG. 7.

Uncinate fasciculus (UF) temporal to frontal (K=131). CS is marked as reference point. Here, the color code for the connection metric is not used. Coronal views shown in bottom row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

FIG. 8.

FOF-frontal to FOF-PO via four different bundles and three waystations. Here, the color code for the connection metric is not used. Coronal views shown in bottom row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

FIG. 9.

Fornix–Fimbria to mammillary bodies (MMB) via three different bundles and two waystations. Here, the color code for the connection metric is not used. Coronal views shown in bottom row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

For the fornix, which is technically a cortico-SB as the hippocampus is old cortex, two waystations were employed as suggested by Concha et al. (2005) and Saunders and Aggleton (2007). The waystations were located within the fornix in coronal sections 192 and 252 (corresponding to SP sections 69 and 97, respectively) and were used with the three initial values of K=200, 200, and 150 reduced to 8, 48, and 40, respectively, after editing. An additional difficulty was the proximity of the fornix to the mid-sagittal plane, which coincided with the boundary of the image volume.

Cortico-striatal bundles

Motor corticostriatal tracts were generated from six distinct cortical ROIs defined by SP Motor Cases 24 through 29 to the putamen (Fig. 10). Some overlapping of fiber bundles was observed, that is, kissing and crossing fibers, as well as a bifurcating and reconnecting bundle.

FIG. 10.

Tracts to the putamen from SP Motor Case 24 (light green, K=200), Case 25 (light blue, K=200), Case 26 (pink, K=161), Case 27 (green, K=159), Case 28 (blue, K=191), and Case 29 (purple, K=118). Here, the color code for the connection metric is not used. Coronal views shown in bottom row; see also movie at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

Discussion

A probabilistic method for tracking of fiber bundles that utilizes DP minimization of a quadratic function based on the Gaussian form of the full DTI tensor has been developed. DP has been shown to generate realistic tracts, including crossing and kissing fibers in DTI images without need for multiple gradients. DP is computationally less intensive than other probabilistic tracking methods that have emerged in recent years due to its focus on generating tracts between user-specified ROIs and, thus, is hypothesis driven rather than exploratory (cf. Croxson et al., 2005). DP also uses the scoring method based on a connectivity metric, such as in Sherbondy et al. (2008).

The DP approach can be applied to more advanced, high angular-resolution (HAR) techniques albeit via a different cost function based on a probabilistic model. Although the limitations of Gaussian diffusion have been acknowledged, it continues to be widely used. This is because of several practical issues of the HAR approach: the large number of required diffusion encoding directions and, more importantly, the large b-values (typically more than 3000 sec/mm²). As demonstrated by Ben-Amitay et al. (2012), such values cause much of the MR signal to decay to the noise floor, resulting in poor SNR. More importantly, the large b-value leads to higher motion sensitivity and additional eddy current distortion, which are difficult to detect and correct because of low SNR. So, it is likely that research and clinical studies will continue to use standard DTI scans.

The DP approach is attractive with regard to DTI data; deterministic approaches such as those mentioned earlier are known to generate severe false-negative results, because if there is a voxel with an incorrect fiber orientation due to noise or complex anatomy, deviation from the real path occurs. Since the deterministic approaches are usually combined with knowledge-based anatomical constraints such as at least using at least two ROIs with known connectivity, the paths that deviate from real ones are eliminated. Likewise, the DP approach is driven by anatomical knowledge; the two locations with known connections are specified and the most probable path is searched, reducing the chance of false positives. This approach is more robust against noise, partial volume effect, and problematic regions with complex tract anatomy. It is, however, important to examine whether the reconstructed paths agree with known trajectories between the two specified locations.

DP falls within a class of computationally efficient pathfinding strategies that are beginning to be useful for tracking bundles. One is the A* approach (Richter et al., 2013), which was found to be more efficient than global optimization strategies such as Gibbs tracking (Reisert et al., 2011). In fact, DP was able to completely traverse the 2D spiral phantom, which could not be done with Gibbs tracking (see Fig. 3 in Richter et al., 2013). The strong principle of optimality used in DP states that all optimal paths comprise optimal sub-paths, and it is this fact which permits the reduction in complexity of the problem (Bertsekas, 1987). The cost function assumes an additive cost that depends locally only on the current node and path direction, which is equivalent to the strong principle of optimality. As mentioned earlier, the number of paths evaluated in a brute force approach would be in the order of MV^N , where V is the connectedness, which in this case is 26. DP iteratively examines MV optimal sub-paths of ever-increasing length by assuming the strong principle of optimality, reducing the complexity to MVN and thereby making the problem tractable. Reducing V from 26 to 6 will still yield optimal paths, however at the price of a higher cost as exemplified by the transition costs for a skewed directional vector in the last row of Table 1. Since a larger neighborhood of 74 was used in the A* approach, it is possible to use a larger search space in the current software. However, the computational efficiency of the A* approach is dependent on the selection of an appropriate heuristic function. A poor choice can result in a brute force search. With DP, the reduction in search complexity is direct and independent of heuristics or other subjective conditions.

It might be argued that the feasibility of DP should be evaluated using phantom data such as PISTE (http://cubric.psych.cf.ac.uk/commondti) or FiberCup (www.lnao.fr/spip.php?rubrique79). However, it is perhaps more meaningful to generate tracts between anatomically defined ROIs using real data. So, a hemi-brain from a macaque was used as a viable means of understanding the capabilities and limitations of DP. Data used in this study for comparison or development of more sophisticated algorithms are being made available at http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/

Indeed, the macaque brain has been previously used as a model for testing methods for studying WM tracts. Parker et al. (2002a) was the first to use the rhesus monkey to evaluate tractography methods, albeit at a low resolution. More recently, Hofer and Frahm (2008) used a 2.9T scanner to study tracts in a macaque DTI scan that was upsampled to 0.75 mm³ isotropic voxels. Using an FA threshold of 0.15, they used the FACT algorithm in MRIStudio (Mori et al., 1999) to generate several tracts. While the ROI seeds they used are different from those used here, the generated SLFII and somatosensory U-fibers are similar. Unlike DP, FACT was unable to deal with kissing and crossing fibers seen with HAR (Wedeen et al., 2008, 2012, ). Adluru et al. (2012) also used a deterministic approach to track several bundles with an FA threshold of 0.1–0.15 in a template of 0.5 mm³ resolution generated as a population average from 271 macaques from a 3T scanner.

DP has several limitations. First, there is inherent bias to high anisotropy and length. As can be seen from Figure 1, it is important to visualize the tracts for anatomical consistency as well as ensure that the ROIs are defined appropriately. Perhaps not surprisingly, due to the cost function, false tracts can be generated especially in the neighborhood of bundles with higher anisotropy as in the case of the FOF and the fornix. This problem is not uncommon even with probabilistic methods (Croxson et al., 2005). This confound can be addressed in one of two ways. One is by using a scoring system based on connectivity metric to ensure anatomically consistent bundles are generated; it is quite straightforward to show that the connectivity metric is approximately inversely proportional to the fractional anisotropy, that is, the lower the cost the higher the anisotropy. The other way is to exclude adjacent bundles from the search space as was done for the deterministic approach by Zhang et al. (2010). For example, generating the FOF and fornix was resolved by excluding adjacent bundles and using waystations with a lower FA threshold of 0.1. Second, the image resolution is insufficient to deal with fiber bundles emanating from designated ROIs that come in very close proximity to each other; hence, the importance of checking the fibers for consistency with known fiber pathways. Third, the method is dependent on high resolution or smooth data. Due to the nature of the 26-voxel neighborhood, it is possible for individual paths to turn 35–45° when passing through noisy or isotropic regions, but these can be masked by generating multiple fibers. Here, it will be helpful to know Contrast-to-Noise Ratio or SNR of the image being studied (Farrell et al., 2007; Kingsley, 2006). In addition, an option for smoothing fibers via tensioned B-splines was implemented but it was not used.

Though no direct comparisons with other methods were made, it is worth noting that Figure 2 could not be replicated by FACT because of the crossing fibers in the CB. In addition, given the FOF-Frontal and FOF-PO ROIs, FACT with a threshold of 0.1 generated a bundle that was more inferior to the true FOF going through the middle longitudinal fasciculus and Extreme Capsule (see Figs. 14-2 and 15-1 in Schmahmann and Pandya, 2006) similar to those generated by deterministic approaches (Adluru et al., 2012; Hofer and Frahm, 2008). This led to the inference of a macaque homologue of the human inferior FOF that does not exist (Catani et al., 2007; Schmahmann and Pandya, 2007).

Last but not the least, it is critically important to validate tracts with regard to staining of tracts via histological sections (Annese, 2012; Dauguet et al., 2006; Seehaus et al., 2013). In this regard, there was very close correspondence between tracts generated with DP and in vivo anterograde transport of HRP in a macaque brain. Generated tracts were also consistent with those from the Schmahmann and Pandya (2006) atlas. It should be possible to reproduce the tractography for the right hemisphere, but the focus on the left has been motivated by an ongoing study of histology and MRI analyses in a macaque model of fetal exposure to radiation in which only the left hemisphere is available (Ceritoglu et al., 2010; Selemon et al., 2013). It is worth noting that the vast majority of published tracer studies have traditionally focused on one side.

Extensions and modifications of DP are possible in several ways. First, a larger neighborhood than the 26-voxel one used here could allow for acute turns (Merhof et al., 2006a, 2006b; Richter et al., 2013). The approach to modify the DP state space to control curvature and torsion (Lal, 2004) can be applied to construct a Markov model in which the transition cost depends on the previous transition along the path. Thus, a Markovian approach (Iglesias et al., 2012) could be used to incorporate bending energy into the cost function (Poupon et al., 2000) for regularization. Adjusting the state space and cost function to be based on voxel triples (Collins et al., 2011) allows for more sophisticated regularization. It is also possible to incorporate Bayesian priors based on tensor coherence or dot-product with regard to a neighboring voxel (Fout et al., 2005; Merhof et al., 2006a, 2006b). Second, for lower-resolution data such as human DTI data, it may be necessary to subparcellate the initial and end ROIs into m and n smaller ROIs to generate m×n bundles each with K fibers, which would be similar to the example of subdividing the motor cortex to generate cortico-striatal tracts. To minimize the need for editing fibers or setting K too large, one could use a WM parcellation atlas such as that developed by Huang et al. (2006) to exclude regions from the state space to generate anatomically consistent bundles as suggested earlier. Furthermore, the connectivity metric provides a single value for each fiber, which could be helpful in studying bundles disrupted by disease and may provide a better alternative to calculating the average fractional anisotropy in the bundle. These limitations and extensions will be explored more fully, particularly in human data in subsequent papers.

Footnotes

Acknowledgments

The authors thank Joseph Hennessey, Hao Huang, Can Ceritoglu, Michael Bowers, Elizabeth Postell, Carolyn Zin, and Geoffrey Gunter for technical assistance. Research was supported by NIH grants (R01-AG020012, R01-EB003543, P41-EB015909, R01-EB000975, and P50-MH71616).

Author Disclosure Statement

No competing financial interests exist.

References

Acar

, Yörük

. 2009. DT-MRI connectivity and/or tractography? Two new algorithms. Aja-Fernandez

, de Luis Garcia

, Tao

, Li

. Tensors in Image Processing and Computer Vision. London, UK: Springer-Verlag, 335–353.

Adluru

, Zhang

, Fox

, Shelton

, Ennis

, Bartosic

, Oler

, Tromp do

, Zakszewski

, Gee

, Kalin

, Alexander

. 2012. A diffusion tensor brain template for rhesus macaques. Neuroimage, 59:306–318.

Alexander

, Hasan

, Kindlmann

, Parker

, Tsurudal

. 2000. A geometric analysis of diffusion tensor measurements of the human brain. Mag Reson Med, 44:283–291.

Annese

. 2012. The importance of combining MRI and large-scale digital histology in neuroimaging studies of brain connectivity and disease. Front Neuroinform, 6:13.

Basser

, Jones

. 2002. Diffusion-tensor MRI: theory, experimental design and data analysis—a technical review. NMR Biomed, 15:456–467.

Basser

, Mattiello

, Le Bihan

. 1994a. Estimation of the effective self-diffusion tensor from the NMR spin echo. J Magn Reson B, 103:247–254.

Basser

, Mattiello

, Lebihan

. 1994b. MR diffusion tensor spectroscopy and imaging. Biophys J, 66:259–267.

Basser

, Pajevic

, Pierpaoli

, Duda

, Aldroubi

. 2000. In vitro fiber tractography using DT-MRI data. Magn Reson Med, 44:625–632.

Beaulieu

. 2002. The basis of anisotropic water diffusion in the nervous system—a technical review. NMR Biomed, 15:435–455.

10.

Behrens

, Johansen-Berg

, Woolrich

, Smith

, Wheeler-Kingshott

, Boulby

, Barker

, Sillery

, Sheehan

, Ciccarelli

, Thompson

, Brady

, Matthews

. 2003. Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging. Nat Neurosci, 6:750–757.

11.

Ben-Amitay

, Jones

, Assaf

. 2012. Motion correction and registration of high b-value diffusion weighted images. Magn Reson Med, 67:1694–1702.

12.

Bertsekas

. 1987. Dynamic Programming. Englewood Cliffs, NJ: Prentice-Hall, Inc.

13.

Catani

, Allin

, Husain

, Pugliese

, Mesulam

, Murray

, Jones

. 2007. Symmetries in human brain language pathways correlate with verbal recall. Proc Natl Acad Sci USA, 104:17163–17168.

14.

Ceritoglu

, Wang

, Selemon

, Csernansky

, Miller

, Ratnanather

. 2010. Large deformation diffeomorphic metric mapping registration of reconstructed 3D histological section images and in vivo MR images. Front Hum Neurosci, 4:43.

15.

Collins

, Singh

, Alexander

. 2011. Network connectivity via inference over curvature-regularizing line graphs. Lect Notes Comp Sci, 6492:65–78.

16.

Concha

, Gross

, Beaulieu

. 2005. Diffusion tensor tractography of the limbic system. AJNR Am J Neuroradiol, 26:2267–2274.

17.

Conturo

, Lori

, Cull

, Akbudak

, Snyder

, Shimony

, McKinstry

, Burton

, Raichle

. 1999. Tracking neuronal fiber pathways in the living human brain. Proc Nat Acad Sci USA, 96:10422–10427.

18.

Croxson

, Johansen-Berg

, Behrens

, Robson

, Pinsk

, Gross

, Richter

, Kastner

, Rushworth

. 2005. Quantitative investigation of connections of the prefrontal cortex in the human and macaque using probabilistic diffusion tractography. J Neurosci, 25:8854–8866.

19.

Dauguet

, Peled

, Berezovskii

, Delzescaux

, Warfield

, Born

, Westin

C-F

. 2006. 3D histological reconstruction of fiber tracts and direct comparison with diffusion tensor MRI tractography. Med Image Comput Comput Assist Interv, 4190:109–116.

20.

Everts

, Bekker

, Roerdink Jos

BTM

. 2009. Visualizing White Matter Structure of the Brain Using Dijkstra's Algorithm. Zinterhof

, Lončarić

, Uhl

, Carini

. Proceedings of the 6th International Symposium on Image and Signal Processing and Analysis (ISPA 2009), Salzberg, Austria, 575–580.

21.

Farrell

, Landman

, Jones

, Smith

, Prince

, Van Zijl

, Mori

. 2007. Effects of signal-to-noise ratio on the accuracy and reproducibility of diffusion tensor imaging-derived fractional anisotropy, mean diffusivity, and principal eigenvector measurements at 1.5 T. J Magn Reson Imaging, 26:756–767.

22.

Fout

, Huang

, Ding

. 2005. Visualization of Neuronal Fiber Connections from DT-MRI with Global Optimization. SAC’05: Proceedings of the 2005 ACM Symposium on Applied Computing, 2005 New York, NY: ACM Press, 1200–1206.

23.

Frandsen

, Hobolth

, Ostergaard

, Vestergaard-Poulsen

, Vedel Jensen

. 2007. Bayesian regularization of diffusion tensor images. Biostatistics, 8:784–799.

24.

Friman

, Westin

C-F

. 2005. Uncertainty in white matter fiber tractography. Lect Notes Comp Sci, 8:107–114.

25.

Goldman-Rakic

, Selemon

, Schwartz

. 1984. Dual pathways connecting the dorsolateral prefrontal cortex with the hippocampal formation and parahippocampal cortex in the rhesus monkey. Neuroscience, 12:719–743.

26.

Hageman

, Toga

, Narr

, Shattuck

. 2009. A diffusion tensor imaging tractography algorithm based on Navier-Stokes fluid mechanics. IEEE Trans Med Imaging, 28:348–360.

27.

Hofer

, Frahm

. 2008. In vivo mapping of fiber pathways in the rhesus monkey brain. Open Med Imaging J, 2:32–41.

28.

Huang

, Walker

, Ying

, Zee

, Van Zijl

, Mori

. 2006. 3D DT-MRI atlas of macaque brain. Proc Int Soc Magn Reson Med, 14:440.

29.

Iglesias

, Thompson

, Liu

, Tu

. 2012. Fast approximate stochastic tractography. Neuroinformatics, 10:5–17.

30.

Iturria-Medina

, Canales-Rodriguez

, Melie-Garcia

, Valdes-Hernandez

, Martinez-Montes

, Aleman-Gomez

, Sanchez-Bornot

. 2007. Characterizing brain anatomical connections using diffusion weighted MRI and graph theory. Neuroimage, 36:645–660.

31.

Jackowski

, Kao

, Qiu

, Constable

, Staib

. 2004. Estimation of anatomical connectivity by anisotropic front propagation and diffusion tensor imaging. Med Image Comput Comput Assist Interv, 3217:663–671.

32.

Jbabdi

, Johansen-Berg

. 2011. Tractography: where do we go from here? Brain Connect, 1:169–183.

33.

Jbabdi

, Woolrich

, Andersson

, Behrens

. 2007. A Bayesian framework for global tractography. Neuroimage, 37:116–129.

34.

Jiang

, Van Zijl

, Kim

, Pearlson

, Mori

. 2006. DTIStudio: resource program for diffusion tensor computation and fiber bundle tracking. Comput Methods Programs Biomed, 81:106–116.

35.

Johansen-Berg

, Rushworth

. 2009. Using diffusion imaging to study human connectional anatomy. Annu Rev Neurosci, 32:75–94.

36.

Jonasson

, Bresson

, Hagmann

, Cuisenaire

, Meuli

, Thiran

. 2005. White matter fiber tract segmentation in DT-MRI using geometric flows. Med Image Anal, 9:223–236.

37.

Kingsley

. 2006. Introduction to diffusion tensor imaging mathematics: part III. tensor calculations, noise, simulations, and optimization. Concepts Magn Reson Part A, 28A:155–179.

38.

Lal

. 2004. Tracking in Diffusion Tensor Imaging. MSE Thesis. http://cis.jhu.edu/data.sets/macaque_hemi-brain_fiber_tracking/ The Johns Hopkins University.

39.

Lazar

, Alexander

. 2003. An error analysis of white matter tractography methods: synthetic diffusion tensor field simulations. Neuroimage, 20:1140–1153.

40.

Le Bihan

, Mangin

, Poupon

, Clark

, Pappata

, Molko

, Chabriat

. 2001. Diffusion tensor imaging: concepts and applications. J Magn Reson Imaging, 13:534–546.

41.

Le Bihan

, Poupon

, Amadon

, Lethimonnier

. 2006. Artifacts and pitfalls in diffusion MRI. J Magn Reson Imaging, 24:478–488.

42.

Le Bihan

, Turner

, Macfall

. 1989. Effects of intravoxel incoherent motions (IVIM) in steady-state free precession (SSFP) imaging: application to molecular diffusion imaging. Magn Reson Med, 10:324–337.

43.

Lifshits

, Tamir

, Assaf

. 2009. Combinatorial fiber-tracking of the human brain. Neuroimage, 48:532–540.

44.

Marigonda

, Orlandi

. 2011. A mathematical model for neuronal fibers. Commun Appl Ind Math, 2:e-363.

45.

Merhof

, Enders

, Hastreiter

, Ganslandt

, Fahlbusch

, Nimsky

, Stamminger

. 2006a. Neuronal fiber connections based on A*-pathfinding. Manduca

, Amini

. Medical Imaging 2006: Physiology, Function, and Structure from Medical Images. San Diego, CA: SPIE, 61431S-8.

46.

Merhof

, Richter

, Enders

, Hastreiter

, Ganslandt

, Buchfelder

, Nimsky

, Greiner

. 2006b. Fast and accurate connectivity analysis between functional regions based on DT-MRI. Med Image Comput Comput Assist Interv, 4191:225–233.

47.

Moldrich

, Pannek

, Hoch

, Rubenstein

, Kurniawan

, Richards

. 2010. Comparative mouse brain tractography of diffusion magnetic resonance imaging. Neuroimage, 51:1027–1036.

48.

Mori

, Crain

, Chacko

, Van Zijl

. 1999. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann Neurol, 45:265–269.

49.

Mori

, Van Zijl

. 2002. Fiber tracking: principles and strategies–a technical review. NMR Biomed, 15:468–480.

50.

Mukherjee

, Chung

, Berman

, Hess

, Henry

. 2008. Diffusion tensor MR imaging and fiber tractography: technical considerations. AJNR Am J Neuroradiol, 29:843–852.

51.

Oishi

, Huang

, Yoshioka

, Ying

, Zee

, Zilles

, Amunts

, Woods

, Toga

, Pike

, Rosa-Neto

, Evans

, Van Zijl

, Mazziotta

, Mori

. 2011. Superficially located white matter structures commonly seen in the human and the macaque brain with diffusion tensor imaging. Brain Connect, 1:37–47.

52.

Pace

, Pierpaoli

. 1963. Photodynamic effect of diidroquinidine on isolated hearts of mammals. Nature, 199:915–916.

53.

Parker

, Stephan

, Barker

, Rowe

, Macmanus

, Wheeler-Kingshott

, Ciccarelli

, Passingham

, Spinks

, Lemon

, Turner

. 2002a. Initial demonstration of in vivo tracing of axonal projections in the macaque brain and comparison with the human brain using diffusion tensor imaging and fast marching tractography. Neuroimage, 15:797–809.

54.

Parker

, Wheeler-Kingshott

, Barker

. 2002b. Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging. IEEE Trans Med Imaging, 21:505–512.

55.

Pichon

, Westin

, Tannenbaum

. 2005. A Hamilton-Jacobi-Bellman approach to high angular resolution diffusion tractography. Med Image Comput Comput Assist Interv, 8:180–187.

56.

Pierpaoli

, Jezzard

, Basser

, Barnett

, Chiro

. 1996. Diffusion tensor MR imaging of human brain. Radiology, 201:637–648.

57.

Poupon

, Clark

, Frouin

, Regis

, Bloch

, Le Bihan

, Mangin

. 2000. Regularization of diffusion-based direction maps for the tracking of brain white matter fasicles. Neuroimage, 12:184–195.

58.

Poynton

, Lal

, Ratnanather

, Mori

, Boatman

, Miller

. 2005. Probabilistic Tracking of Fiber Pathways Using Dynamic Programming. Human Brain Mapping Conference, Toronto.

59.

Ramnani

, Behrens

, Penny

, Matthews

. 2004. New approaches for exploring anatomical and functional connectivity in the human brain. Biol Psychiatry, 56:613–619.

60.

Reisert

, Mader

, Anastasopoulos

, Weigel

, Schnell

, Kiselev

. 2011. Global fiber reconstruction becomes practical. Neuroimage, 54:955–962.

61.

Richter

, Zolal

, Ganslandt

, Buchfelder

, Nimsky

, Merhof

. 2013. Evaluation of diffusion-tensor imaging-based global search and tractography for tumor surgery close to the language system. PLoS One, 8:e50132.

62.

Saunders

, Aggleton

. 2007. Origin and topography of fibers contributing to the fornix in macaque monkeys. Hippocampus, 17:396–411.

63.

Schmahmann

, Pandya

. 2006. Fiber Pathways of the Brain. Oxford, New York: Oxford University Press.

64.

Schmahmann

, Pandya

. 2007. The complex history of the fronto-occipital fasciculus. J Hist Neurosci, 16:362–377.

65.

Seehaus

, Roebroeck

, Chiry

, Kim

, Ronen

, Bratzke

, Goebel

, Galuske

. 2013. Histological validation of DW-MRI tractography in human postmortem tissue. Cereb Cortex, 23:442–450.

66.

Selemon

, Ceritoglu

, Ratnanather

, Wang

, Harms

, Aldridge

, Begovic

, Csernansky

, Miller

, Rakic

. 2013. Distinct abnormalities of the primate prefrontal cortex caused by ionizing radiation in early or midgestation. J Comp Neurol, 521:1040–1053.

67.

Selemon

, Goldman-Rakic

. 1985. Longitudinal topography and interdigitation of corticostriatal projections in the rhesus monkey. J Neurosci, 5:776–794.

68.

Selemon

, Goldman-Rakic

. 1988. Common cortical and subcortical targets of the dorsolateral prefrontal and posterior parietal cortices in the rhesus monkey: evidence for a distributed neural network subserving spatially guided behavior. J Neurosci, 8:4049–4068.

69.

Sherbondy

, Dougherty

, Ben-Shachar

, Napel

, Wandell

. 2008. ConTrack: finding the most likely pathways between brain regions using diffusion tractography. J Vis, 8:15.1–16.

70.

Tournier

, Calamante

, King

, Gadian

, Connelly

. 2002. Limitations and requirements of diffusion tensor fiber tracking: an assessment using simulations. Magn Reson Med, 47:701–708.

71.

Tournier

, Mori

, Leemans

. 2011. Diffusion tensor imaging and beyond. Magn Reson Med, 65:1532–1556.

72.

Tuch

. 2004. Q-ball imaging. Magn Reson Med, 52:1358–1372.

73.

Wedeen

, Rosene

, Wang

, Dai

, Mortazavi

, Hagmann

, Kaas

, Tseng

. 2012. The geometric structure of the brain fiber pathways. Science, 335:1628–1634.

74.

Wedeen

, Wang

, Schmahmann

, Benner

, Tseng

, Dai

, Pandya

, Hagmann

, D'Arceuil

, De Crespigny

. 2008. Diffusion spectrum magnetic resonance imaging (DSI) tractography of crossing fibers. Neuroimage, 41:1267–1277.

75.

Xue

, Van Zijl

, Crain

, Solaiyappan

, Mori

. 1999. In vivo three-dimensional reconstruction of rat brain axonal projections by diffusion tensor imaging. Magn Reson Med, 42:1123–1127.

76.

Zalesky

. 2008. DT-MRI fiber tracking: a shortest paths approach. IEEE Trans Med Imaging, 27:1458–1471.

77.

Zalesky

, Fornito

. 2009. A DTI-derived measure of cortico-cortical connectivity. IEEE Trans Med Imaging, 28:1023–1036.

78.

Zhang

, Richards

, Yarowsky

, Huang

, Van Zijl

, Mori

. 2003. Three-dimensional anatomical characterization of the developing mouse brain by diffusion tensor microimaging. Neuroimage, 20:1639–1648.

79.

Zhang

, Zhang

, Oishi

, Faria

, Jiang

, Li

, Akhter

, Rosa-Neto

, Pike

, Evans

, Toga

, Woods

, Mazziotta

, Miller

, Van Zijl

, Mori

. 2010. Atlas-guided tract reconstruction for automated and comprehensive examination of the white matter anatomy. Neuroimage, 52:1289–1301.