Magnetic resonance diffusion-weighted imaging coupled with fiber tractography (DFT) is the only non-invasive method for measuring white matter pathways in the living human brain. DFT is often used to discover new pathways. But there are also many applications, particularly in visual neuroscience, in which we are confident that two brain regions are connected, and we wish to find the most likely pathway forming the connection. In several cases, current DFT algorithms fail to find these candidate pathways. To overcome this limitation, we have developed a probabilistic DFT algorithm (ConTrack) that identifies the most likely pathways between two regions. We introduce the algorithm in three parts: a sampler to generate a large set of potential pathways, a scoring algorithm that measures the likelihood of a pathway, and an inferential step to identify the most likely pathways connecting two regions. In a series of experiments using human data, we show that ConTrack estimates known pathways at positions that are consistent with those found using a high quality deterministic algorithm. Further we show that separating sampling and scoring enables ConTrack to identify valid pathways, known to exist, that are missed by other deterministic and probabilistic DFT algorithms.

*R*

_{1}to

*R*

_{2}should equal the score of the same pathway measured from

*R*

_{2}to

*R*

_{1}. Independence means that the score of a pathway between

*R*

_{1}and

*R*

_{2}depends only on the data along that pathway. Therefore, the score along the green pathway between

*R*

_{1}and

*R*

_{2}should be the same in the presence or absence of another pathway. For example, if the orange tensors in Figure 2 become isotropic due to a disease, thereby lowering the score of a pathway traversing them, the score along the

*R*

_{1}–

*R*

_{2}pathway should remain unchanged.

*R*

_{1}to

*R*

_{2}is unlikely to be identified by a deterministic algorithm with typical stopping rules based on a fractional anisotropy (FA) threshold; the FA of the green tensor in the fourth column is very low and would stop the growth. But if we know

*R*

_{1}is connected to

*R*

_{2}, the green pathway is the most likely route.

*s*

_{i}). The pathway score is determined by the diffusion data at each node and the angles between the linear segments.

*Q*(

*s*), as,

*s*is a pathway and

*D*are the data along the pathway. The scoring function is defined in this way to separate the data-dependent portion of the score,

*p*(

*D*∣

*s*), from the data independent portion,

*p*(

*s*). Both the

*p*(

*D*∣

*s*) and

*p*(

*s*) components of the pathway score are a product of the local scores along all the nodes of the pathway, as described in the following paragraphs.

*p*(

*D*∣

*s*) measures how well the existence of pathway

*s*is supported by the tensors

*D*derived from diffusion measurements. The pathway score aggregates local terms multiplicatively, i.e.,

*n*-nodes that define the pathway

*s*with diffusion tensors

*D*

_{i}and tangent vectors

*t*

_{i}. The local score for the diffusion data given the pathway tangent vector is computed using a Bingham distribution (Cook, Alexander, & Parker, 2004; Kaden, Knösche, & Anwander, 2007; Mardia & Jupp, 2000),

*v*

_{1},

*v*

_{2}, and

*v*

_{3}are the three ordered eigenvectors of the

*i*th diffusion tensor,

*D*

_{i};

*C*(

*σ*

_{3},

*σ*

_{2}) is a normalizing constant so that the function integrates to one over the sphere; and

*σ*

_{i}is the decay rate in the direction of the eigenvector,

*v*

_{i}. The two decay rate parameters are determined by summing two dispersion parameters that are based on (a) uncertainty in the data (

*σ*

_{m}) and (b) the ellipsoid shape (

*σ*

_{i}*). That is,

*σ*

_{i}=

*σ*

_{m}+

*σ*

_{i}*.

*σ*

_{m}) from repeated measures of the principal direction of diffusion or PDD (

*v*

_{1}). The repeatability of the data sets a lower bound on fiber direction uncertainty. We estimate

*σ*

_{m}using a bootstrap procedure (Efron, 1979): a tensor is fit to 1000 permutations of the raw DWI. The bootstrap procedure estimates the mean PDD and the dispersion about that mean. Other bootstrap procedures, such as the wild bootstrap (Liu, 1988), can be used if the DWI data lacks an adequate number of repetitions (Whitcher, Tuch, & L., 2005). We summarize the bootstrapped dispersion,

*σ*

_{m}, using the Watson distribution on the 3D sphere (Mardia & Jupp, 2000; Schwartzman, Dougherty, & Taylor, 2005). The Watson distribution is identical to a Bingham distribution with radial symmetry, i.e., when

*σ*

_{2}=

*σ*

_{3}.

*C*

_{L}> 0.4) in all of our data sets. The linearity index is a measure of anisotropy (Peled, Gudbjartsson, Westin, Kikinis, & Jolesz, 1998) and is the positive difference between the largest two eigenvalues of the diffusion tensor divided by the sum of its eigenvalues. In our algorithm, increasing the linearity criterion does not change the minimum dispersion value; decreasing the linearity criterion increases the minimum dispersion to a value larger than 4°. The precise value of this parameter is not critical to the performance of ConTrack.

*σ*

_{i}*, from the diffusion shape. This dispersion term is necessary because the PDD of an oblate or nearly isotropic diffusion ellipsoid might be highly reliable, yet such voxels clearly have high fiber direction uncertainty. Reliable measurements of oblate and spherical ellipsoids occur when voxels include many crossing fibers, such as where the tapetum interdigitates with the inferior longitudinal fasciculus.

*δ*is a dispersion factor calculated from

*C*

_{L}the tensor linearity index, and a user parameter

*η,*discussed below.

*σ*

_{m}= 4° and

*λ*

_{1}=

*λ*

_{2}=

*λ*

_{3}, would have a uniform distribution about the sphere (

*σ*

_{2}=

*σ*

_{3}= 54°).

*δ*varies between the maximum added dispersion and zero within ∼0.15 units of

*C*

_{L}. Adjusting these constants did not alter the results significantly.

*p*(

*s*) (Equation 1) measures how well the pathway represents prior knowledge about neuronal connections. This term currently encodes assumptions such as the smoothness of white matter pathways, length preference, as well as the possible locations for endpoints of pathways. We expect that over the years our knowledge will increase based on studies of the neuroanatomy and the function can evolve to incorporate this knowledge.

*θ*

_{i}, see Figure 3), the location of pathway nodes, and the number of nodes, we assess the score of a pathway (independent of the data) as

*p*

_{curve}, defined on the local intra-pathway angles, is a Watson distribution (Cook et al., 2004; Kaden et al., 2007; Mardia & Jupp, 2000)

*θ*

_{i}∣ ≤

*π*/2, i.e., only directions on the hemisphere are possible,

*σ*

_{c}is the angular dispersion parameter in degrees, and

*C*(

*σ*

_{c}) is a normalizing constant. The Watson distribution of axes on a sphere is easily transformed into a distribution of directions on a hemisphere and this was a motivation for selecting this distribution. The possible directions are distributed on the hemisphere because it is assumed that the pathways have a node spacing small enough such that an intra-pathway angle >90° would be too sharp to reconstruct the desired anatomy. When

*σ*

_{c}is less than 54° (uniform dispersion), a pathway is penalized for having many large angular deviations.

*p*

_{length}and

*p*

_{end}are defined on the pathway nodes.

*λ is*a manually chosen length scoring parameter and the white matter mask is defined as those voxels that are within the brain mask and have an FA >0.15 and either mean diffusivity (MD) less than 1.1

*μ*m

^{2}/ms or FA >0.4. This mask is dilated by one voxel to allow for partial voluming at the edges of the white matter. Using FA and MD to create the white matter mask is a common method (e.g., Basser et al., 2000; Behrens, Woolrich, et al., 2003; Mori, Crain, Chacko, & van Zijl, 1999; Parker et al., 2003). As with other algorithms, the mask prevents pathways from entering CSF or gray matter. The accuracy of the white matter mask was ensured by visual inspection.

*σ*

_{3}< 14°), then the next step direction is chosen according to

*p*(

*D*

_{i}∣

*t*

_{i}); (b) otherwise the step direction is drawn from the

*p*

_{curve}(

*θ*

_{i}) distribution. The step distance is currently 1 mm. The pathway is terminated if the current node is outside of the white matter mask or within an ROI voxel. Only pathways with endpoints in both ROIs are retained. The process repeats until a desired number (∼10

^{5}) of pathway samples is collected.

*k*-space acquisition). We acquired 48–54 axial, 2-mm thick slices (no skip) for two

*b*-values,

*b*= 0 and

*b*= 800 s/mm

^{2}. The high

*b*-value was obtained by applying gradients along 12 different diffusion directions. Two gradient axes were energized simultaneously to minimize TE and the polarity of the effective diffusion-weighting gradients was reversed for odd repetitions to reduce cross-terms between diffusion gradients and both imaging and background gradients.

*b*-spline algorithm based on code from SPM5 (Friston & Ashburner, 2004). Also, the eddy-current intensity correction (Rohde et al., 2004) was applied to the diffusion-weighted images at this stage. We note that the 7th-order

*b*-spline interpolation does not require image variance correction (Rohde, Barnett, Basser, & Pierpaoli, 2005) due to the large support kernel. Preserving the signal variance structure in the interpolated data is crucial for an accurate bootstrap variance estimate. After rotating the diffusion-weighting gradient directions to preserve their orientation with respect to the rotated diffusion images, the tensors were fit using a least-squares algorithm. We confirmed that this registration technique aligns the DTI and T1 images to within a few millimeters (except in regions prone to susceptibility artifacts, such as orbito-frontal and anterior temporal regions). All the custom image processing software is available as part of our open-source mrDiffusion package available for download from our Web site (http://sirl.stanford.edu/software/).

*η, λ,*and

*σ*

_{c}.

*σ*

_{c}parameter was empirically set by taking a collection of STT pathways seeded throughout the entire brain with seed spacing 1 mm, step distance 1 mm, FA threshold of 0.15, and angle threshold of 30° and fitting a Watson distribution on the hemisphere to the set of adjacent step directions within all the STT pathways. The fitted value was 14°.

*η*parameter was set based on an examination of the

*C*

_{L}value of the voxels in regions of suspected fasciculi crossings. This value was set to 0.175.

*λ*parameter was chosen using a search procedure. We generated a database of pathways with the sampling algorithm on the same set of sixteen anatomical problems as described in the Results section. We selected the thousand pathways with the highest scores and created a binary volume indicating voxels that contain at least one pathway. We measured pathway similarity by computing the correlation coefficient between two binary volumes estimated using two sets of parameters. We repeated these pathway estimates multiple times varying

*λ*between values

*e*

^{−4}and

*e*

^{1}. As

*λ*increased to the positive edge of this range the estimated pathways looped and were inconsistent with known anatomy. At

*λ*=

*e*

^{−2}, there was no looping. Moreover, small changes near this value produced very high correlation coefficients (Figure 5). Hence, we used

*λ*=

*e*

^{−2}.

*η*and

*σ*

_{c}parameters as there is little change in the resulting pathway sets across the range of attempted values (see Figure 5). The user may want to adjust

*λ*depending on the distance between the anatomical regions of interest, although we never observed a large effect (Figure 5).

STT L–R (DOC) | ConTrack L–R (DOC) | ConTrack and STT (DOC) | ConTrack L–R (MT+) | |
---|---|---|---|---|

S1 | 0 | 6 | 0 | 30 |

S2 | 30 | 35 | 31 | 15 |

S3 | 35 | 27 | 27 | 20 |

S4 | 17 | 25 | 40 | 33 |

*R*

_{1}and

*R*

_{2}, that is,

*p*(

*R*

_{1}→

*R*

_{2}) ≠

*p*(

*R*

_{2}→

*R*

_{1}).

*R*

_{1}to

*R*

_{2}divided by the total number of pathways from

*R*

_{1}(Behrens, Woolrich, et al., 2003; Friman & Westin, 2005; Parker & Alexander, 2003). In principle, one should not evaluate the probability of a connection between

*R*

_{1}and

*R*

_{2}based on the score (likelihood) of unrelated pathways. ConTrack avoids this failure of independence entirely.

^{3}and large enough to cover the entire LGN (Andrews, Halpern, & Purves, 1997). The estimated optic radiation forms a thin sheet that includes short direct connections and longer connections that match the expectations of Meyer's loop (Behrens, Johansen-Berg et al., 2003; Parker et al., 2003).

*λ*value. The other two parameters did not noticeably alter the top scoring pathways.