ConTrack: Finding the most likely pathways between brain regions using diffusion tractography

Abstract

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.

Introduction

Magnetic resonance diffusion-weighted imaging coupled with fiber tractography (DFT) provides important information about white matter in the living human brain (Basser, Mattiello, & LeBihan, 1994; Catani, Howard, Pajevic, & Jones, 2002; Conturo et al., 1999; Mori & Barker, 1999; Mori & van Zijl, 2002). There has been significant progress during the past decade in both diffusion measurement methods and the tractography algorithms designed to identify and quantify human white matter tracts (Dougherty, Ben-Shachar, Bammer, Brewer, & Wandell, 2005; Wakana, Jiang, Nagae-Poetscher, van Zijl, & Mori, 2004). These advances have generated new discoveries about human anatomy and promise clinical utility. 

While these methods do identify many known tracts, they do not identify several large tracts that are known to exist, and there is need for further algorithm development. Because of their diversity, the visual pathways offer excellent targets for algorithm development. Some pathways are easy to find (occipital–frontal) (Catani et al., 2002) and occipital–callosal (Dougherty, Ben-Shachar, Bammer, et al., 2005; Dougherty, Ben-Shachar, Deutsch, et al., 2005); other known pathways represent challenging targets for identification (e.g., LGN–MT, or callosal–MT). Further, the ability to find certain visual pathways has immediate clinical application (Ebeling & Reulen, 1988; Powell et al., 2005; Yamamoto et al., 2007). 

In this paper, we describe a probabilistic DFT algorithm (ConTrack) developed to trace known pathways. We developed and applied this algorithm to pathways in the visual system of the living human brain. We show that ConTrack identifies pathways in visual cortex that are known to exist but that are missed by other algorithms. 

Tractography algorithms can be grouped loosely into two types: deterministic and probabilistic. Deterministic tractography algorithms, such as Streamlines Tracing Technique (STT) (Basser, Pajevic, Pierpaoli, Duda, & Aldroubi, 2000; Conturo et al., 1999; Mori & Barker, 1999) and Tensorlines/TEND (Lazar et al., 2003; Weinstein, Kindlmann, & Lundberg, 1999), generate pathways by making a sequence of discrete, locally optimal decisions. These algorithms have the advantage of simplicity and speed. Further, they generate discrete tract estimates that are easily visualized and analyzed. However, they have several limitations. First, they do not account for uncertainty in the data: a pathway either exists or not. The binary assessment conceals valuable information about the reliability of the underlying data. Second, current deterministic algorithms do not explore thoroughly the entire space of possible white matter pathways because decision making is strictly local. A pathway that passes through even a small region that violates deterministic rules is not considered. Such local thresholding makes these algorithms vulnerable to small regional aberrations in the data such as those caused by measurement artifacts and noise, as well as unresolved features such as crossing fiber tracts. While there are proposals using high angular resolution diffusion imaging to resolve crossing tracts (Hagmann et al., 2004; Tournier, Calamante, Gadian, & Connelly, 2004; Tuch et al., 2002), there is no solution to the problem of measurement uncertainty. 

Several groups have developed probabilistic algorithms to address the limits of deterministic algorithms (Anwander, Tittgemeyer, von Cramon, Friederici, & Knösche, 2007; Behrens, Woolrich, et al., 2003; Bjornemo, Brun, Kikinis, & Westin, 2002; Friman & Westin, 2005; Hagmann et al., 2003; Hosey, Williams, & Ansorge, 2005; Lazar & Alexander, 2002; Parker, Haroon, & Wheeler-Kingshott, 2003; Perrin et al., 2005). The advantage of these probabilistic methods is that they (a) expand the pathway search space beyond that explored by deterministic algorithms and (b) have principled ways of ordering the likelihood among the possible pathways. However, there remain several limitations to current probabilistic tractography algorithms. 

One limitation concerns the primary output of these algorithms, referred to interchangeably as the “probability of a dominant streamline passing through any single voxel” (Behrens, Berg, Jbabdi, Rushworth, & Woolrich, 2007, p. 4), or more commonly as “connection probabilities” (Behrens, Johansen-Berg, et al., 2003; Croxson et al., 2005; Parker et al., 2003; Powell et al., 2005; Rushworth, Behrens, & Johansen-Berg, 2006). In our view, these algorithms are useful because they use probabilistic reasoning to search for anatomically valid pathways, not because they compute an accurate probability of brain connections. For example, in a normally sighted individual, we can be certain that the lateral geniculate nucleus (LGN) is connected to primary visual cortex (V1 and V2). However, the estimated connection probability in the widely used FMRIB's diffusion tractograhpy (FDT) tool is substantially less than one (Behrens, Johansen-Berg, et al., 2003; see Figure 11). 

As an extreme example of inappropriate connection probability estimation, current probabilistic algorithms have also been criticized for failing to identify connecting pathways (i.e., connection probability is 0), even when these are known to exist (Behrens et al., 2007; Hosey et al., 2005; Jbabdi, Woolrich, Andersson, & Behrens, 2007; Sherbondy et al., 2006). Authors have proposed ways to force algorithms to identify at least some connecting pathways (Jbabdi et al., 2007; Sherbondy et al., 2006). But these methods do not address the fundamental source of this problem, which is that they link pathway sampling with the estimation of anatomical validity (scoring). Specifically, current probabilistic algorithms sample pathways at a frequency that is strictly proportional to their score. Consequently, highly valid anatomical pathways with relatively low scores are sampled many fewer times and these valid pathways can be missed entirely. 

ConTrack finds pathways invisible to other tractography algorithms by separating the pathway sampling and scoring processes. Rather than sampling strictly according to the pathway score, the pathway generation algorithm identifies a broad range of potential pathways, including some that will have rather modest validity. A separate scoring stage rates the pathway anatomical validity. Separating the pathway sampling and scoring has the advantage that pathway scores do not depend on (a) the location of the seed-point used to generate the pathway (symmetry) or (b) the scores of other pathways (independence). ConTrack should be useful for applications in which the user wishes to measure the properties of pathways that are known to exist rather than to estimate the probability that these pathways exist. This arises in various clinical applications, such as identifying visual pathways prior to surgical resection (Powell et al., 2005; Yamamoto et al., 2007). 

In the next section, we describe the principles and implementation of the ConTrack algorithm. We then describe experiments that clarify its useful characteristics as well as its limitations. Finally, we compare the contributions of the ConTrack approach to those of other fiber tractography algorithms. 

ConTrack algorithm

ConTrack comprises three essential parts. First, scoring procedures encode knowledge about pathway likelihood; the scoring algorithm is a key component in finding anatomically correct pathways. Second, the vast number of possible pathways makes it impossible to score every pathway, so ConTrack includes a sampling procedure to select potential pathways for scoring. Third, the algorithm specifies a simple inference method that defines how to use sampled and scored pathways to select the most likely pathways that connect two regions. Figure 1 illustrates the ConTrack algorithm by showing example results from each of the algorithm's three stages: ROI specification, pathway sampling, and pathway scoring and selection. 

Figure 1

 

Example outputs of the three stages of the ConTrack algorithm. (A) ROI selection: One ROI was placed to cover the corpus callosum (cyan) and a second ROI covers a region of occipital cortex (green). (B) Pathway sampling: Pathways (yellow curves) are generated that have an endpoint in each ROI. (C) Pathway selection: The highest scoring 1% of the sampled pathways is retained as representing the most likely pathways connecting the two regions. Crosshair: 1 × 1 cm.

Scoring procedures

Principles

Pathway scoring encapsulates our knowledge about the likely properties of white matter pathways. The scoring used in ConTrack also incorporates two basic principles that have not been adopted by other algorithms: symmetry and independence (Figure 2). 

Figure 2

 

Independence and symmetry in pathway scoring. Glyphs represent tensor data. The shape varies from isotropic and round (low FA) to thin cylinders (high FA). In this figure, the tensors are generally isotropic and randomly oriented (gray) apart from a row of tensors between R₁ and R₂ (green) and a column of tensors rising from R₁ (orange). Symmetry asserts that the green pathway should have the same score whether measured from R₁ to R₂ or from R₂ to R₁. Independence asserts that the green pathway score should be the same whether the orange pathway is present or replaced by a random set of tensors.

Symmetry means that the pathway score measured from R₁ to R₂ should equal the score of the same pathway measured from R₂ to R₁. Independence means that the score of a pathway between R₁ and R₂ depends only on the data along that pathway. Therefore, the score along the green pathway between R₁ and R₂ 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₁–R₂ pathway should remain unchanged. 

Notice also that the pathway from R₁ to R₂ 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₁ is connected to R₂, the green pathway is the most likely route. 

Algorithm

Consider a single sampled pathway with endpoints in the two regions of interest (ROIs) of Figure 2. We call the sample points along the pathway “nodes.” Figure 3 illustrates a short pathway composed of two linear segments connecting three nodes (s_i). The pathway score is determined by the diffusion data at each node and the angles between the linear segments. 

Figure 3

 

Pathway scoring variables. The diffusion tensor (D(s_i)) is labeled for each node (s_i). Pathway tangent vectors at interior nodes (t₂) are computed as the average of the line segments from the two connecting nodes. Tangent vectors at endpoints (t₁, t₃) are in the direction of the one connecting segment. The angle between adjacent pathway segments is π − θ₂ or θ₂ if the segments are aligned tail to tail.

We define the score, Q(s), as,  

$Q\left(s\right)=p(D|s)p\left(s\right),$

(1)

where 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. 

Scoring the pathway fit to the data

The function 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.,  

$p(D|s)={\prod }_{i=1}^{n}p\left(D_i|t_i\right),$

(2)

for the 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),  

$p\left(D_i|t_i\right)=C({\sigma }_3,{\sigma }_2)\exp \{-{\left(\frac{v_3^T{t}_i}{\sin {\sigma }_3}\right)}^2-{\left(\frac{v_2^T{t}_i}{\sin {\sigma }_2}\right)}^2\},$

(3)

where v₁, v₂, and v₃ are the three ordered eigenvectors of the ith diffusion tensor, D_i; C(σ₃, σ₂) 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*. 

We calculate one dispersion parameter (σ_m) from repeated measures of the principal direction of diffusion or PDD (v₁). 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 σ₂ = σ₃. 

The bootstrap sometimes produces an unreasonably low dispersion parameter. Hence, we set a minimum dispersion to 4°—the mean dispersion of all brain voxels with a high linearity index (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. 

We derive a second dispersion term, σ_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. 

A user-specified function relates ordered tensor eigenvalues (λ_i) to the shape dispersion term.  

$\begin{array}{c}{\sigma }_2^{*}=\frac{{\lambda }_2}{{\lambda }_2+{\lambda }_3}\delta \\ \\ {\sigma }_3^{*}=\frac{{\lambda }_3}{{\lambda }_2+{\lambda }_3}\delta \\ \\ \delta ={100}^{°}/(1+\exp \{-(\eta -C_{\mathrm{L}})/0.015\}).\end{array}$

(4)

 

The value δ is a dispersion factor calculated from C_L the tensor linearity index, and a user parameter η, discussed below. 

The constant, 100°, is maximum total dispersion. A perfect spherical tensor, where σ_m = 4° and λ₁ = λ₂ = λ₃, would have a uniform distribution about the sphere (σ₂ = σ₃ = 54°). 

The constant 0.015 is chosen so that δ varies between the maximum added dispersion and zero within ∼0.15 units of C_L. Adjusting these constants did not alter the results significantly. 

Four examples of local distribution functions are shown in Figure 4. When the linearity index is small, the shape of the tensor adds significant uncertainty to the local distribution function. This uncertainty may be concentrated along one axis for more oblate-shaped ellipsoids or distributed uniformly across the sphere for more spherical ellipsoids. When the ellipsoid has a very high linearity index (i.e., a prolate ellipsoid), the total fiber direction uncertainty equals the bootstrap dispersion estimate. 

Figure 4

 

Scoring function parameters. (A) The sigmoidal dependence of δ on C_L when η = 0.175 (Equation 4). (B) These four spheres plot four extremes of the local distribution on the pathway tangent to the diffusion data. The conditions are prolate (λ₂ = λ₃) and oblate (λ₂ > λ₃) tensors with either low or high linearity values.

Scoring the pathway shape independent of the data

The function 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. 

Using the angles of the pathway as it enters and exits each node (θ_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\left(s\right)=p_{\mathrm{e}\mathrm{n}\mathrm{d}}\left(s_1\right)p_{\mathrm{e}\mathrm{n}\mathrm{d}}\left(s_n\right){\prod }_{i=2}^{n-1}{p}_{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{v}\mathrm{e}}\left({\theta }_i\right)p_{\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}}\left(s_i\right).$

(5)

 

The function p_curve, defined on the local intra-pathway angles, is a Watson distribution (Cook et al., 2004; Kaden et al., 2007; Mardia & Jupp, 2000)  

$p_{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{v}\mathrm{e}}\left({\theta }_i\right)=C\left({\sigma }_c\right)\exp \left(\frac{{\cos }^2{\theta }_i}{{\sin }^2{\sigma }_c}\right),$

(6)

where we set 0 ≤ ∣θ_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. 

The two functions, p_length and p_end are defined on the pathway nodes.  

$p_{\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}}\left(s_i\right)=\{\begin{array}{c}\lambda ,\mathrm{i}\mathrm{f}{s}_i\mathrm{i}\mathrm{n}\mathrm{w}\mathrm{h}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\\ \\ 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array},$

(7)

 

$p_{\mathrm{e}\mathrm{n}\mathrm{d}}\left(s_i\right)=\{\begin{array}{c}1,\mathrm{i}\mathrm{f}{s}_i\mathrm{i}\mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{R}\mathrm{O}\mathrm{I}\\ \\ 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array},$

(8)

where λ 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²/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. 

Pathway sampling

The pathway sampling algorithm identifies a large set of potential pathways that have endpoints in the two ROIs. The pathway sampling procedure differs from procedures used in other probabilistic tractography algorithms (Behrens, Woolrich, et al., 2003; Friman & Westin, 2005; Parker & Alexander, 2005). Rather than sampling strictly proportionally to estimated anatomical validity (likelihood), the ConTrack sampler includes pathways with a broad range of scores. 

The algorithm begins by placing a seed point randomly within one of the voxels in an ROI. A pathway is grown from this seed position by iteratively choosing step directions and moving a unit step distance along that direction. A step direction is chosen in one of two ways; (a) if this is the first step or if the fiber direction uncertainty is small (σ₃ < 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⁵) of pathway samples is collected. 

While the generation of a single pathway must start in one ROI, pathways are seeded equally from each of the two ROIs. Hence, the population of pathways is symmetric and not biased to either ROI. 

Inference methodology

ConTrack scores the likelihood of each sampled pathway. We infer that the highest scoring pathways are the most likely to exist; the selection criterion for anatomical validity is based on a combination of visual inspection of the pathways and comparison of the anatomical validity scores for pathways identified by STT. We derive various quantitative measures from the highest scoring pathways, including the pathway destination (projection zone), mean fractional anisotropy, pathway length, and so forth. 

This inference—that the highest scoring pathways are the most likely to exist—is a simple and perhaps obvious design choice for our application; yet, the reader should note that the inference methodology used in ConTrack differs significantly from that used by probabilistic tractography algorithms, which measure the connection probability between two regions. In the Discussion section, we describe why the combination of separating sampling from scoring and the difference in inference strategy enables ConTrack to find pathways that other algorithms miss. 

Experimental methods

Subjects

Data were obtained from four healthy subjects with no history of neurological disease, head injury, attention deficit/hyperactivity disorder, or psychiatric disorder. The Stanford Panel on Human Subjects in Medical and Non-Medical Research approved all procedures. The four subjects (S1: female aged 11, S2: male aged 11, S3: female aged 8, and S4: female aged 9) were recruited from the San Francisco Bay region of California and were paid for their participation. The data were collected as part of a longitudinal study of reading (Ben-Shachar, Dougherty, Deutsch, & Wandell, 2007; Dougherty, Ben-Shachar, Bammer, et al., 2005; Dougherty, Ben-Shachar, Deutsch, et al., 2005). Written informed consent/assent was obtained from all parents and children. 

Image acquisition

The DTI protocol involved averaging 8–10 repetitions of a 90-second whole-brain scan. The pulse sequence was a diffusion-weighted single-shot spin-echo, echo planar imaging sequence (63 ms TE; 6 s TR; 260 mm FOV; 128 × 128 matrix size; ±110 kHz bandwidth; partial 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². 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. 

We also collected high-resolution T1-weighted anatomical images using an eight-minute sagittal 3D-SPGR sequence (approximately 1 × 1 × 1 mm voxel size). The following anatomical landmarks were defined in the T1 images: the anterior commissure (AC), the posterior commissure (PC), and the mid-sagittal plane. With these landmarks, we utilized a rigid-body transform to convert the T1 data to the conventional AC–PC aligned space. 

Eddy current distortions and subject motion in the diffusion-weighted images were removed by a 14-parameter non-linear co-registration that is based on the expected pattern of eddy-current distortions given the phase-encode direction of the acquired data (Rohde, Barnett, Basser, Marenco, & Pierpaoli, 2004). Each diffusion-weighted image was registered to the mean of the (motion-corrected) non-diffusion-weighted images using a two-stage coarse-to-fine approach that maximized the normalized mutual information. The mean of the non-diffusion-weighted images was also automatically aligned to the T1 image using a rigid-body mutual information algorithm. All raw images from the diffusion sequence were then resampled to 2 mm isotropic voxels by combining the motion correction transform, eddy-current correction transform, and anatomical alignment transform into one transform and resampling using a 7th-order 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/). 

For each subject, we also acquired functional MR data to define the ROIs used in the Results section. BOLD responses were measured using a spiral pulse sequence with 21–30 obliquely oriented slices acquired every 2.4 seconds (30 ms TE; 1.2 s TR; 2 interleaves; 70° flip angle; 2 × 2 × 3 mm effective voxel size). Functional scan sessions lasted about 45 minutes, broken into 8 four-minute scans, allowing the subject to rest between scans. Visual field maps were measured in the subject using rotating wedge and expanding ring stimuli. For all subjects, we parametrically varied the contrast of moving luminance gratings while subjects were engaged in a two-interval, forced-choice speed discrimination task. We measured BOLD responses as a function of stimulus contrast in motion sensitive regions (MT+, V1, and V3A). 

Algorithm parameters

There are three user-defined parameters in the scoring procedure: η, λ, and σ_c. 

The σ_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°. 

The η 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. 

The λ 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⁻⁴ and e¹. As λ increased to the positive edge of this range the estimated pathways looped and were inconsistent with known anatomy. At λ = e⁻², there was no looping. Moreover, small changes near this value produced very high correlation coefficients (Figure 5). Hence, we used λ = e⁻². 

Figure 5

 

Scoring algorithm sensitivity to parameters. The dotted light gray line summarizes pathways connecting dorsal occipital cortex with the corpus callosum (CC) and the solid black line summarizes (MT+)–CC pathways. The correlation coefficient is described in the text. (A) Varying the length penalty λ between e⁻⁴ and e¹ per node of pathway. (B) Varying the smoothness dispersion parameter, σ_c, from 10° to 45°. (C) Varying the position of the sigmoid mediating the effect of the diffusion tensor on the score, η, between 0.05 and 0.35. Both η and σ_c have little effect on the relative pathway scoring, while altering λ produces a smooth difference across a large range of values. (D) The effect of the number of initial samples on the similarity of the output pathway set for η = 0.175, σ_c = 14°, and λ = e⁻². All pathway sets were compared to the set retained after collecting 100,000 pathways.

Finally, using these three parameter values, we examined the scoring function's sensitivity to changes in these values. Each parameter was searched independently while the others were held at a fixed value. There were sixteen anatomical problems presented to ConTrack—one for each subject, hemisphere, and section of anatomy, as described in the Results section. The analyses confirmed that there is little value in modifying the η 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). 

Algorithm time requirements

The tensor fitting and bootstrapping of the tensor fits to our data required approximately 30 minutes on a computer with a 2.2-GHz, 64-bit CPU, and 8 GB of RAM. Pathway generation required an additional 0.5–1.5 hours for each of the cases presented in the Results section on the same machine. Pathway scoring and thresholding could be done in approximately 10 minutes. 

ConTrack software, documentation as well as supporting example data will be available for download at https://simtk.org/home/contrack. 

Results

Algorithm validation

The first experiment was designed to confirm that highly reliable white matter pathways produced by STT are also recovered by ConTrack. In previous work, we found that pathways from dorsal occipital visual areas that pass through the corpus callosum (CC) were reliably estimated by STT (Dougherty, Ben-Shachar, Bammer, et al., 2005). We determined the extent to which ConTrack reliably identifies these pathways. 

We tested STT by comparing how homologous fiber groups from the two hemispheres converge in the callosum (Dougherty, Ben-Shachar, Bammer, et al., 2005). We use this method again to evaluate ConTrack pathway estimates. Specifically, we identified two regions of interest in homologous left and right dorsal occipital cortex (DOC; see Figure 6A). We also identified the entire corpus callosum (CC; see Figure 6A). ConTrack was run once between the left DOC and CC and then again between the right DOC and CC. 

Figure 6

 

Callosal projection zones from dorsal occipital cortex (DOC) measured with ConTrack and STT. (A) Locations of the three regions of interest: Left and right DOC (green) and the corpus callosum (cyan). (B) The four images in this column show the callosal projection zones from left and right DOC in four subjects. Colors distinguish regions containing projections from left DOC only (dark yellow), right DOC only (medium yellow), and both (bright yellow). (C) The callosal projection zones from DOC estimated with STT are shown in the same four subjects. Projections from left DOC only (dark green), right DOC only (medium green), and both (bright green) are shown. Scale bar: 1 cm. Crosshair: 1 × 1 cm.

Many of the dorsal occipital fibers from the left and right hemispheres connect homologous regions; hence, there should be substantial overlap in the CC between the left and right tract estimates. The overlap from the independent ConTrack pathway estimates in the mid-sagittal plane through the CC is substantial and confined to the splenium (Figure 6B). The estimates for these same subjects using independent STT pathway estimates are also shown (Figure 6C; Table 1). In order to compare STT estimates with those of ConTrack, the STT pathways were clipped to the ROIs that were specified for ConTrack and only pathways that reached both ROIs were retained. STT generates fewer pathways than ConTrack. But the STT center-of-mass of the intersection agrees within 2 mm of the location estimated by ConTrack, except for S1 where STT failed to find overlap between left and right hemispheres (Figure 6C). 

Table 1

 

Percentage of the callosal projection zone overlap from left (L) and right (R) for dorsal occipital cortex (DOC) and MT as estimated by ConTrack and STT. The first two columns are the percentage overlap between left and right DOC–CC pathways as estimated by ConTrack and STT, respectively. The third column is the percentage of overlap between algorithms for the union of left and right DOC–CC pathways. The fourth column is the percentage overlap between left and right (MT+)–CC pathways as estimated by ConTrack. Measurements are in percentages.

	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

There are both similarities and differences in the pathways identified by the two algorithms. First, we compared the scores of the identified pathways (Figure 7). ConTrack identifies many more pathways than STT, including high-scoring pathways. STT pathways, however, typically contain the very highest scoring pathways. STT also produces many low-scoring pathways. 

Figure 7

 

Score histograms for pathways connecting left DOC to CC (left column) and right DOC to CC (right column) for four subjects (rows) in a log-log plot. There are far fewer pathways returned by the STT algorithm (open circles) as compared to the ConTrack algorithm (open triangles).

We also measured the distances between pathways to quantify the difference between the ConTrack and STT pathway estimates. We computed the distance for a specific STT-ConTrack pathway by finding, for each node on the STT pathway, the closest point along the ConTrack pathway. We summarize the pathway separation by using either the maximum or the mean of these distances. The STT-ConTrack pair with the smallest distance is considered the corresponding pair. For 7 out of the 8 hemispheres, the maximum distance between corresponding pairs was no greater than 5 mm; in one hemisphere the maximum distance was within 10 mm. The mean distances between corresponding pairs in all 8 hemispheres was less than 4 mm for the high scoring pathways. 

The few cases where the pathway distance exceeded 5 mm were due to differences in the projection location of the pathways within the cortex. As the pathways left core white matter and approach the ROI in the DOC, the destination voxel between the two algorithms in the ROI diverged. Hence, the results from STT and ConTrack, particularly in the white matter core, are quite similar. 

Using ConTrack to identify connections previously missed

As is noted in the literature, STT fails to identify certain types of connections. Missed connections of specific interest to us are the callosal fibers that project to area MT+, a motion-selective region located on the lateral surface at the border between the occipital and temporal lobes. Dense callosal projections to area MT+ have been demonstrated by anatomical studies in monkey (Krubitzer & Kaas, 1993; Maunsell & Van Essen, 1987) and post-mortem human brains (Clarke & Miklossy, 1990). These callosal fibers are not detected with STT (Dougherty, Ben-Shachar, Bammer, et al., 2005; Huang et al., 2005; Wakana et al., 2004), probably due to the problem of crossing fibers, as these callosal pathways likely interdigitate with the much larger inferior longitudinal fasciculus. Evidence from monkey (Pandya, Karol, & Heilbronn, 1971; Schmahmann & Pandya, 2006) suggests that MT+ callosal connections intersect the callosum at the dorsal part of the splenium. 

In this second experiment, we identified the most likely pathways between right and left MT+ (Figure 8) in the four subjects. The location of MT+ was identified using BOLD fMRI and a standard motion localizer stimulus. In all subjects, most pathways travel near the dorsal surface of the lateral ventricle (i.e., the tapetum), as expected (Clarke & Miklossy, 1990). 

Figure 8

 

ConTrack pathways between left and right MT+ (blue) and between left and right DOC (green). The pathways shown are the top scoring 1% pathways connecting the specific ROIs. The red object is the lateral ventricle, shown for reference. The four panels show estimates from four subjects. Crosshair: 1 × 1 cm.

In a third experiment, we compared the intersection points within a medial slice of the corpus callosum for left (MT+)–CC and right (MT+)–CC pathways (Figure 9A). This experiment is designed to validate the (MT+)–CC connections as in the previous work (Figure 6; see also Dougherty, Ben-Shachar, Bammer, et al., 2005). We identified three regions of interest around left and right MT+ and the corpus callosum. If the left and right MT+ communicate as we expect, the callosal intersection points from left and right will overlap (Zeki, 1980). Indeed, the intersection points from pathways of the independent data sets from left and right MT+ overlap substantially in the callosum (Figure 9B; Table 1). The position of these crossing fibers is slightly dorsal to the position of the DOC crossing fibers. 

Figure 9

 

Callosal projection zones from MT+ measured with ConTrack. (A) Locations of the three regions of interest: Left and right MT+ (blue) and the corpus callosum (cyan). (B) Callosal projection zones from left (dark yellow), right (medium yellow), and both (bright yellow) for four subjects are shown. The intersection of the DOC crossing fiber positions in Figure 6 is shown for each subject (red). Scale bar: 1 cm. Cross hair: 1 × 1 cm.

Discussion

Related work

The goal of ConTrack—to estimate pathways connecting two regions—differs from the goal of most other probabilistic tractography algorithms, whose goal is to compute the probability of a connection between a seed region and other voxels. Thus, rather than producing a map of connection probabilities, ConTrack produces a set of pathways and their anatomical validity scores. We use ConTrack to trace a known pathway so that we will be able to measure the diffusion properties along those pathways with a high anatomical validity score. 

Principles

There are several differences in the scoring principles used by ConTrack and other probabilistic methods. First, scoring symmetry is not embedded into many tractography algorithms. Friman and Westin (2005) specifically point out that their method and most other tractography algorithms do not satisfy symmetry for the probability of connection between R₁ and R₂, that is, p(R₁ → R₂) ≠ p(R₂ → R₁). 

Second, ConTrack has a much greater degree of independence than other algorithms. The failure of independence in other probabilistic algorithms arises from two sources. The first failure is due to the way connection probability is calculated. Specifically, the connection probability is defined as the number of pathways from region R₁ to R₂ divided by the total number of pathways from R₁ (Behrens, Woolrich, et al., 2003; Friman & Westin, 2005; Parker & Alexander, 2003). In principle, one should not evaluate the probability of a connection between R₁ and R₂ based on the score (likelihood) of unrelated pathways. ConTrack avoids this failure of independence entirely. 

The second failure is due to the way pathways are sampled. In most algorithms, the probability of sampling a pathway strongly depends on the properties of other pathways. A relatively low likelihood pathway is sampled much less frequently or even not at all when there exists a higher likelihood pathway between these two regions. In principle, one should not sample a pathway less often simply because other anatomically valid pathways exist. The ConTrack sampling parameters are adjusted so that the decline in sampling probability as a function of score is much more gradual than in other samplers. Hence, by separating sampling from scoring, ConTrack can achieve a higher degree of independence. 

Two examples illustrating the advantage of increasing independence are described below. 

Probabilistic tracking using way points and enforced connectivity

It is possible to use “way point” methods with existing probabilistic tractography algorithms (Behrens, Woolrich, et al., 2003; Parker & Alexander, 2005) and new methods that enforce connectivity (Jbabdi et al., 2007) to find pathways between two regions known to be connected. We compare the results of both of these methods with ConTrack in two applications to the visual pathways. First, we consider pathways between the corpus callosum and a motion-responsive region on the lateral occipital lobe (MT+). Then we consider the optic radiation that connects the lateral geniculate nucleus (LGN) and primary visual cortex (V1). 

Callosal–MT+ pathways

We applied the FDT (Behrens, Woolrich, et al., 2003) and PICo (Parker et al., 2003) algorithms, using the way points method, to the same four data sets described in Experimental methods with the same two regions of interest: MT+ and the CC. 

Specifically, we downloaded the FDT algorithm included in the FSL version 3.3 from the FMRIB Web site, www.fmrib.ox.ac.uk/fsl/. For each of the four data sets, the tractography algorithm was run with the default parameters (number of samples per seed = 5000, curvature threshold = 80°, step length = 0.5 mm, loop check = true, use fractional anisotropy = true, maximum number of steps per pathway = 2000). The inverse of the white matter mask used by the ConTrack algorithm was provided as an exclusion-with-truncation mask for the FDT algorithm. The algorithm was run in “Two-Masks—Symmetric” mode with the two masks defined as the MT+ and CC ROIs. 

The PICo algorithm, as included in version 2 revision 340 of CAMINO, was downloaded from www.cs.ucl.ac.uk/research/medic/camino/. The PICo tractography algorithm was run with the following parameters that were chosen to be default when possible (SNR = 14, iterations = 5000, curvature threshold = 80°, distribution = Bingham, discard loops = true, step length = 0.2 mm). The inverse of the white matter mask used by the ConTrack algorithm was provided as an exclusion-with-truncation mask for the PICo algorithm. The algorithm was run such that the MT+ and CC ROIs were used as seed regions and way point regions, forcing the algorithm to only show paths intersecting both and starting from either. 

For both algorithms, voxels with a connection probability of at least 0.01 were included as voxels containing anatomically valid pathways. This means that at least 1% of the identified fibers intersected the voxels in the valid pathways. No connections between left MT+ and the posterior portion of corpus callosum were found in any of the four subjects for either the FDT or PICo algorithms. Right (MT+)–CC connections were found in three subjects; an example for left and right is shown for the subject in which the FDT algorithm performed the best in Figure 10. The single stray pathway connecting left MT+ (Figure 10A) is far from the expected crossing region in the splenium, which is based on the extensive monkey literature and the post-mortem human work (de Lacoste, Kirkpatrick, & Ross, 1985; Krubitzer & Kaas, 1990). No meaningful overlap analysis could be performed on the projection of the pathways onto the CC. 

Figure 10

 

Comparison of methods for identifying pathways between MT+ and the CC within subject S1. (A) Maximum intensity projections of FDT pathway estimates shown superimposed on T1 axial image. The MT+ ROIs are shown as blue and the CC ROI as cyan. Purple pathway connects left MT+ and the CC (only one found) and truncated just right of CC; light blue represent right (MT+)–CC pathways (any voxel with FDT probability of connectivity >0.01) and truncated just left of CC. (B) Three-dimensional visualization of the FDT pathways intersecting right MT+ and the CC only. These pathways are not truncated after the CC as was done in A. (C) ConTrack estimates of the top 1% scoring pathways intersecting right MT+ and the CC. Crosshairs: 1 × 1 cm.

Notice that the FDT analysis mistakenly identifies connections between the medial surface of the occipital cortex and the CC as being right (MT+)–CC connections (Figure 10C). The mistake occurs because the right MT+ ROI extends into the white matter through this pathway. It is impossible for the user to fully appreciate this error when using the inference methods of FDT, which discard the pathway information and examine voxelwise connection probability maps. The same way point analysis was carried out using PICo; the results were quite similar to FDT (see Supplementary Figure 1). 

In addition to STT, FDT, and PICo, we tested Tensorlines, a deterministic algorithm (Lazar et al., 2003; Weinstein et al., 1999). The parameters were similar to those used for STT and the puncture coefficient was set to 0.2 (Weinstein et al., 1999). Tensorlines found (MT+)–(MT+) pathways similar to those found by ConTrack in all of our subjects. The core principles in Tensorlines are similar to those of ConTrack. First, in regions of fiber direction uncertainty (as measured by tensor shapes that are more oblate or spherical rather than prolate), Tensorlines relies on a smoothness constraint. Thus, at these locations Tensorlines weights the local pathway history more than the local PDD. This is similar to ConTrack's fiber direction uncertainty function (Equation 3), which is also based on tensor shape. 

While Tensorlines was successful in finding (MT+)–(MT+) fiber tracts, the strategy of deterministically tracing a single pathway in regions of uncertainty is guaranteed to fail in situations where the correct pathway makes a relatively sharp turn. For such fiber tracts, probabilistic methods that produce many paths of differing smoothness (like ConTrack) are more likely to find the correct pathways. The optic radiation is an example and will be discussed below. 

Optic radiation

An important application of fiber tracking is to identify major pathways that should be avoided during neurosurgical resection. The optic radiation is an important example in temporal lobe resection for epilepsy (Powell et al., 2005; Yamamoto et al., 2007). A portion of the optic radiation extends into the temporal lobe white matter before curving back in a posterior direction toward V1 (calcarine sulcus). This part of the optic radiation is called Meyer's loop (Ebeling & Reulen, 1988; Krolak-Salmon et al., 2000) and surgeons must avoid cutting it during temporal lobe resection to preserve vision. 

Behrens, Johansen-Berg, et al. (2003) describe the difficulty in finding Meyer's loop with FDT. They report finding only direct connections when using conventional procedures and seeding the entire LGN; although these connections do not appear to arrive at the calcarine sulcus (Figure 11A). To identify the loop, they report that it is necessary to move the FDT seed region away from the LGN (Behrens, Johansen-Berg et al., 2003). Moving the seed region reduces the influence of the shorter direct pathways and permits FDT to sample the real, but lower likelihood, Meyer's loop. Similar seed placement methods, using knowledge about the course of the pathway, are required by other probabilistic algorithms, i.e., PICo, for identifying Meyer's loop (Powell et al., 2005). We confirmed these observations using our data (not shown). 

Figure 11

 

Comparison of methods for identifying the optic radiation. (A) The estimated course of the optic radiation (OR) using FDT (from Behrens, Johansen-Berg, et al., 2003). These fibers do not connect to the medial wall of occipital lobe or the occipital pole, where V1 is located in this slice. (B) Maximum intensity projection from a more recent algorithm based on FDT (from Jbabdi et al., 2007). These pathways connect to the expected location of the far-peripheral representation of V1 and do not appear to reach the occipital pole. (C) Three-dimensional visualization of the top scoring 1% of ConTrack pathways connecting the LGN (red) and V1 in subject S3. Pathways viewed from the bottom of the brain. Crosshair: 1 × 1 cm.

Jbabdi et al. (2007) proposed a probabilistic algorithm designed to constrain pathways ensuring a connection between regions is found. Much as we proposed (Sherbondy et al., 2006), they argue that this constraint increases the robustness of the pathway estimates and identifies connections that were previously invisible. Their algorithm builds on the methods in FDT, but it does not include our emphasis on independence and symmetry. They illustrate LGN–V1 connections for FDT and their new algorithm (Figure 11B). As shown for the original FDT algorithm (Figure 11A), these estimates also fail to reach the main portion of V1, and they fail to identify Meyer's loop. 

The ConTrack results, within subject S3, with ROIs surrounding the LGN and calcarine are shown in Figure 11C, on a different subject than Figure 11A or 11B. Our LGN ROI was approximately 270 mm³ 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). 

To be certain that the differences shown in Figure 11 are not due to differences in the data, we ran the FDT algorithm in “Two-Masks—Symmetric” mode (see above), using LGN and V1 as way points, on subject S3 Figure 11C. The FDT results on our data agree with the results reported by Behrens et al. (2007) and Jbabdi et al. (2007; Figures 11A and 11B). At the time of this writing, the technique proposed by Jbabdi et al. was not available. 

This example illustrates a core problem with current probabilistic algorithms: the probability of connection is relative and sampling is proportional to this probability. In FDT and PICo, moving the seed reduces the likelihood of the direct pathway. This increases the probability that the sampler finds Meyer's loop. Further, the optic radiation exposes a limit of deterministic methods, such as Tensorlines. These methods do not explore enough pathways with varying levels of smoothness. Consequently, both Tensorlines and STT fail to find the Meyer's loop portion of the optic radiation (Supplemental Figure 2). 

ConTrack limitations

While developing ConTrack, we identified several aspects of probabilistic algorithms that could be improved. Also, during the development process, we routinely compared our methods with other published algorithms and identified several differences that may not be obvious upon first examination. In the next sections, we discuss areas for improvement and comment on published algorithms that are closely related to our work. 

Parameters

Parameters encapsulate neuroscientific beliefs about likely pathways; these parameters influence the scoring algorithm. ConTrack contains three parameters set by the user. These are 

•  
  η, influences the orientation uncertainty due to tensor shape;
•  
  λ, influences the value of adding or removing pathway nodes; and
•  
  σ_c, the angular dispersion on the hemisphere pertaining to the local pathway curvature.

A limitation of ConTrack is the absence of a simple way to indicate the complex relationship between these parameters and pathway estimates. An analysis of the influence of these parameters (Figures 5A–5C) suggests that the user should examine the consequences of varying the λ value. The other two parameters did not noticeably alter the top scoring pathways. 

Convergence

We do not have a theorem that proves how many pathway samples are required during the initial generation phase, such that the same top scoring pathways would be recovered by ConTrack. Our results show that for the anatomical regions tested in our four subjects, the set of top scoring pathways are not altered much beyond generating the first twenty thousand sample pathways (Figure 5D). Also, other probabilistic algorithms such as FDT and PICo suffer from this convergence limitation. 

Connections known not to exist

The ConTrack algorithm was designed to estimate connections within the white matter that are known to exist. However, it is clear that investigators also use DFT to explore potential connectivity within the white matter. The inference methodology suggested in this manuscript relies on a user defined score threshold for retaining likely pathways between two regions. At this time, we do not rely on a global significance of these scores and instead only rely on the ability of the score to rank order the pathways that connect two particular regions. The ability to compare pathway scores to determine the likelihood that a pathway exists may come with more knowledge of the full distribution of scores for all potential pathways within the brain. We suggest that using ConTrack to first find and score the pathways that are known to exist is a good start to this aim as well. 

User interaction

An important strength of FDT and PICo is that they have a modest number of user-settable parameters. For example, in using FDT the user normally establishes only seed regions, exclusion regions, step-size, number of pathway samples, anisotropy threshold, loop checks, and curvature threshold. Once these are set, the algorithm produces a spatial distribution of connection probabilities, which the user often uses to then specify a voxelwise threshold on connection probability. ConTrack has a similar number of user-settable parameters, but it requires more user involvement. Specifically, the user visualizes the pathways and their scores. This increased user involvement may make some applications more challenging. However, at this point in algorithm development, there is some consensus that human oversight should be encouraged. For the near future, a viable strategy might be to facilitate user involvement by providing pathway visualization tools (Sherbondy, Akers, Mackenzie, Dougherty, & Wandell, 2005). 

Summary

We define and implement an algorithm that uses diffusion-weighted imaging data to identify the most likely pathways connecting two regions. The ConTrack algorithm includes scoring, sampling, and inferencing components. In a series of experiments estimating major visual pathway fiber tracts, we verify that in conditions where STT performs well, ConTrack identifies the same pathways. In certain conditions where other deterministic and probabilistic algorithms fail to find pathways, ConTrack identifies excellent candidates. ConTrack is well suited to applications in which other DFT algorithms fail to identify known pathways, yet one would like to identify the most likely location of these pathways and study their properties. 

Supplementary Materials

Supplementary Figure 1 - Supplementary Figure 1 

Supplementary Figure 1. PICo pathways identified between MT+ (blue) and the CC (cyan) within Subject S1. The maximum intensity projections of PICo pathway estimates are shown superimposed on a T1 axial image. Two purple pathways connect left MT+ and the CC and truncated just right of CC (not shown); light blue represent right (MT+) − CC pathways (voxels with PICo connectivity probability >0.01 are colored) and truncated just left of CC. Crosshairs: 1cm × 1cm. 

Supplementary Figure 2 - Supplementary Figure 2 

Supplementary Figure 2. Comparison of optic radiation estimates using TEND and ConTrack. Three-dimensional visualizations of pathways connecting the LGN (red) and V1 are shown for both left and right hemispheres. The background image is a T1 anatomical, viewed from the bottom of the brain. Subject S3. (A) TEND pathways (blue) do not include the Meyer's loop portion of the OR. (B) Top scoring 1% of ConTrack pathways (yellow) discovered using the same data (see also Figure 11). These pathways show Meyer's loop in addition to the pathways estimated by TEND. Crosshairs: 1cm × 1cm. 

Acknowledgments

The authors would like to thank the Richard M. Lucas Center and especially Roland Bammer for assistance with data acquisition. We also are grateful for the support from NIH Roadmap for Medical Research U54 GM072970 and NIH EY015000 grants. 

Commercial relationships: none. 

Corresponding author: Anthony J. Sherbondy. 

Email: [e-mail]. 

Address: 450 Serra Mall Room 490, Jordan Hall Stanford, CA 94305, USA.