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. 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 (more ...)
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 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 (more ...)
Symmetry means that the pathway score measured from R1 to R2 should equal the score of the same pathway measured from R2 to R1. Independence means that the score of a pathway between R1 and R2 depends only on the data along that pathway. Therefore, the score along the green pathway between R1 and R2 should be the same in the presence or absence of another pathway. For example, if the orange tensors in become isotropic due to a disease, thereby lowering the score of a pathway traversing them, the score along the R1–R2 pathway should remain unchanged.
Notice also that the pathway from R1 to R2 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 R1 is connected to R2, the green pathway is the most likely route.
Consider a single sampled pathway with endpoints in the two regions of interest (ROIs) of . We call the sample points along the pathway “nodes.” illustrates a short pathway composed of two linear segments connecting three nodes (si). 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(si)) is labeled for each node (si). Pathway tangent vectors at interior nodes (t2) are computed as the average of the line segments from the two connecting nodes. Tangent vectors at endpoints (t1, t3 (more ...)
We define the score, Q(s), as,
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
), from the data independent portion, p
). Both the p
) and p
) 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.,
, and v3
are the three ordered eigenvectors of the i
th diffusion tensor, Di
) 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, vi
. 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
We calculate one dispersion parameter (σm
) from repeated measures of the principal direction of diffusion or PDD (v1
). 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
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 (CL
>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.
The value δ is a dispersion factor calculated from CL 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 λ1 = λ2 = λ3, would have a uniform distribution about the sphere (σ2 = σ3 = 54°).
The constant 0.015 is chosen so that δ varies between the maximum added dispersion and zero within ~0.15 units of CL. Adjusting these constants did not alter the results significantly.
Four examples of local distribution functions are shown in . 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 CL 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 (more ...)
Scoring the pathway shape independent of the data
The function p
) (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 ), the location of pathway nodes, and the number of nodes, we assess the score of a pathway (independent of the data) as
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
) 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, plength and pend are defined on the pathway nodes.
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 μ
/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.
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 (σ3 < 14°), then the next step direction is chosen according to p(Di | ti); (b) otherwise the step direction is drawn from the pcurve(θ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 (~105) 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.
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.