White matter fiber tracking directed by interpolating splines and a methodological framework for evaluation

Image-based tractography of white matter (WM) fiber bundles in the brain using diffusion weighted MRI (DW-MRI) has become a useful tool in basic and clinical neuroscience. However, proper tracking is challenging due to the anatomical complexity of fiber pathways, the coarse resolution of clinically applicable whole-brain in vivo imaging techniques, and the difficulties associated with verification. In this study we introduce a new tractography algorithm using splines (denoted Spline). Spline reconstructs smooth fiber trajectories iteratively, in contrast to most other tractography algorithms that create piecewise linear fiber tract segments, followed by spline fitting. Using DW-MRI recordings from eight healthy elderly people participating in a longitudinal study of cognitive aging, we compare our Spline algorithm to two state-of-the-art tracking methods from the TrackVis software suite. The comparison is done quantitatively using diffusion metrics (fractional anisotropy, FA), with both (1) tract averaging, (2) longitudinal linear mixed-effects model fitting, and (3) detailed along-tract analysis. Further validation is done on recordings from a diffusion hardware phantom, mimicking a coronal brain slice, with a known ground truth. Results from the longitudinal aging study showed high sensitivity of Spline tracking to individual aging patterns of mean FA when combined with linear mixed-effects modeling, moderately strong differences in the along-tract analysis of specific tracts, whereas the tract-averaged comparison using simple linear OLS regression revealed less differences between Spline and the two other tractography algorithms. In the brain phantom experiments with a ground truth, we demonstrated improved tracking ability of Spline compared to the two reference tractography algorithms being tested.


INTRODUCTION
Diffusion-weighted MRI (DW-MRI) of the brain is a technique used for measuring white matter (WM) fiber architecture. Tractography (Mori et al., 1998;Basser et al., 2000) aims at reconstructing the fiber trajectories, and thus gives in vivo information about myelinated fiber bundle connections in the brain. Various diffusion parameters, e.g., fractional anisotropy (FA), mean diffusivity (MD), and radial diffusivity (RD), are used to quantify developmental or disease related white matter structural changes in the brain (Westlye et al., 2010;Lebel et al., 2012;Madden et al., 2012).
WM fiber tracking was originally based on the use of a voxelwise diffusion tensor (DT) representation of the recorded data. The earlier form of DW-MRI was first introduced by Basser et al. . Later, Westin et al. (1999) showed various useful measures for tract analysis, e.g., the scalar FA measure. Although highly valuable, the use of a single tensor to represent the diffusion in a voxel has limitations. The coarse image resolution will give voxels that contain fibers of possibly different directions. The diffusion ellipsoid, reconstructed from the eigenvalues and eigenvectors of the tensor, is applicable to tractography in voxels with a single fiber direction, but contain deficient information in voxels with more than one fiber direction (Jones, 2011). A recent study (Jeurissen et al., 2012) estimated that 63-90% of WM voxels contain crossing fibers, meaning that the diffusion tensors provide insufficient information about fiber directions in the majority of WM voxels, and could therefore mislead a tractography method.
To solve the challenge of complex fiber configurations, more advanced (and time-consuming) acquisition schemes have been developed, such as diffusion spectrum imaging (DSI) (Wedeen et al., 2000;Tuch et al., 2001) and Q-ball imaging (Tuch, 2004). They are referred to as "high angular resolution diffusion imaging" (HARDI). In these techniques, the diffusion is measured at higher spatial density of the diffusion sensitizing gradients, and at higher b-values. These studies also introduced an alternative way to represent the diffusion signal, different from the classical diffusion tensor. The measured diffusion signal is now mapped to a function on a sphere, denoted the orientation distribution function (ODF), using a transformation proposed by Stejskal and Tanner (1965). The ODF either reflects the directional distribution of diffusion and is referred to as a diffusion ODF (dODF) (Tuch, 2004), or estimates the fiber directions, and is referred to as fiber ODF (fODF) (Descoteaux et al., 2009). Subsequent studies have shown that these transformation techniques improve tracking in DSI (Wedeen et al., 2008;Fernandez-Miranda et al., 2012), and recent investigations have shown that the use of ODFs may also improve fiber tracking in standard clinical acquisitions (Prckovska et al., 2010). However, most HARDI reconstruction techniques are resource demanding and time-consuming and seem not to be readily available for the clinics.
In this context, we are developing a new tractography algorithm inspired by the smooth nature of WM fiber bundles (Gössl et al., 2002;Jones, 2011;Catani and de Schotten, 2012). Standard tracking methods reconstruct piecewise linear trajectories that are subsequently smoothed. Instead of subsequent smoothing, our method smooths (by spline construction) the fibers iteratively, and we investigate the effect this approach has on the reconstructed trajectories. We first evaluate the method on DW-MRI data acquired at three different time points from eight healthy elderly people participating in a longitudinal study of cognitive aging, using a 1.5 T scanner and a clinically applicable acquisition protocol with 25 diffusion sensitizing directions. In the evaluation, we focus on the corticospinal tract (CST) and the trans-callosal forceps minor (Fminor) fiber bundles. The CST has been studied previously (Westlye et al., 2010;Lebel et al., 2012), reporting only a 0.02 reduction in tract-average FA for the age span from 40 to 80 years. In contrast, the FA decline in Fminor for the same age span was found to be 0.08 in Westlye et al. (2010) and 0.04 in Lebel et al. (2012), indicating less affected FA values during aging for the CST compared to the Fminor tract. Moreover, along-tract FA analysis in e.g., (Colby et al., 2012;Hodneland et al., 2012;Yeatman et al., 2012) revealed local FA variations that were not captured by simple grand mean FA across the tract. The two well-studied fiber bundles CST and Fminor were therefore chosen as a real-data tested for comparing our Spline method to other tracking methods. For further testing the performance of our algorithm we employed DW-MRI data from a brain phantom, designed for the MICCAI 2009 Fiber Cup with known ground truth (Poupon et al., 2008(Poupon et al., , 2010Fillard et al., 2011;Côte et al., 2012), where a 3T MRI-scanner and 64 diffusion sensitizing directions had been used.
The rest of the paper is organized as follows. First we describe the orientation distribution function (dODF) in its probabilistic setting, and then we describe our Spline tracking algorithm using the voxel-wise dODF estimates. The streamline tracking algorithms which are compared to our Spline approach is then presented before more detailed anatomical descriptions of the selected fiber bundles CST and Fminor are provided. In the Materials and Methods section we also give a description of the MR imaging data and acquisition parameters being used for the experiments, and the linear mixed-effects model used in the assessment of longitudinal recordings. In the Results section we describe the outcome from processing the longitudinal healthy aging MRI-data, comparing Spline with the two other tracking methods, before processing results from the Fiber Cup test object recordings are given. In the last section we discuss our approach in a broader context and point to improvements and follow-up experiments for future work.

ORIENTATION DISTRIBUTION FUNCTION
Orientation distribution functions (ODF) (Wedeen et al., 2000;Tuch et al., 2001;Descoteaux et al., 2007;Yeh et al., 2010) are computed from the original diffusion data and can be used to obtain voxelwise information about preferred diffusion directions. The dODF reconstruction technique used in this study is implemented in TrackVis (Wang et al., 2007;Wedeen et al., 2008). From the diffusion signal S, acquired at the diffusion sensitizing directions q, a probability density function (PDF) of spin displacements p(r) is computed: where F is the discrete Fourier transform and r is the 3D displacement vector. The dODF w(v) is then computed by projection of the data in the radial direction: where v is a unit 3D vector, ρ = |r| and ρ 2 dρ is the 3D volume element. The integral is computed as a discrete sum over the range ρ ∈ [0, 5] voxels in diffusion r-space. Then, w(v) is evaluated for a set of vectors at the vertices of a regular polyhedron. Local maxima of the dODFs are estimated subsequently and used as voxelwise fiber orientations. Figure 1 displays the fiber configuration for a single voxel (Figures 1A,B), the diffusion signal (Figure 1C), the corresponding diffusion ellipsoid (Figure 1D), and the dODF ( Figure 1E). Clearly, the dODF representation contains significantly more detailed information about the local water diffusion than the diffusion tensor ellipsoid, and is able to approximate two fiber directions instead of only one (see also Figure 7).

SPLINE TRACKING
Our spline tracking algorithm is inspired by the smooth nature of WM fibers (Gössl et al., 2002;Jones, 2011;Catani and de Schotten, 2012). To achieve a similar smoothness iteratively, we make the use of cubic b-splines (Gerald and Wheatley, 1999). This has two main advantages. Firstly, the propagation is directed by the dODFs along the splines, instead of along the piecewise linear segments. Secondly, noisy voxels may terminate fibers too early using streamline tracking. The use of splines iteratively during the tracking can reduce the impact of such noisy voxels. The cubic b-splines implemented here require at least four control points, and returns a smooth curve between the middle two points. Given four control points p (i) 3 , the first element of the position vector S (i) where i indicates the iteration number and u ∈ {0ν, 1ν, 2ν, . . . , 1}. We choose ν = 0.01. The two other components of s  The initial points of a fiber (the control points), p (1) 0 , . . ., p (1) 3 , and the resulting spline (blue curve). p (1) 1 is the seed point, p (1) 2 is found by going a step-length in the direction of a local maxima in p (1) 1 , p (1) 0 is found by going in the opposite direction, and p (1) 3 is found by going a step-length in the direction of the local maxima most closely to the incoming direction in p 3 in blue are the points from the first iteration. Note that p (2) 1 from this iteration and p (1) 2 from the previous iteration are in this case not the same point. The same is the case for p (2) 2 from this iteration and p (1) 3 from the previous iteration. The resulting splines are shown for iteration 1 (blue curve) and iteration 2 (green curve).
the incoming direction, in the sense of the largest dot product. In general, we findQ where M = {q k } = arg local maxû w(p k ,û) is the set of local maxima directions of w at p (i) k . Using Equation (3) we obtain a spline that is valid between p (1) 1 and p (1) 2 . In the next iteration, we define p 2 is found a step-length from p (2) 1 in the direction of the local maximum at p (2) 1 closest to the incoming direction. In a similar way, p (2) 3 is found by considering the dODF in p (2) 2 , and going a step in the local maximum closest to the incoming direction. Both local maxima are found using Equation (4). This procedure is repeated iteratively for new points.
The maximum number of iterations n is not pre-defined. However, when the fiber reaches a voxel with FA ≤ 0.15 or the fiber takes a sharp turn, the iterations stop. A sharp turn is defined Frontiers in Neuroinformatics www.frontiersin.org July 2013 | Volume 7 | Article 13 | 3 by a large dot product between the fiber direction of consecutive iterations, when 1 be a seed point, andQ k andQ k + 1 satisfy (4). h is the step-length.

then Break End
Since we update the points as we move along the fiber, the splines are not continuous, as we see in Figure 2. When the full trajectory is extracted, we select points along the tract and compute a smooth tract using cubic b-splines similarly to Equation (3), now only using every 20 fiber coordinate.

COMPARATIVE STREAMLINE TRACKING
We performed the comparative streamline tracking using two available state-of-the-art reconstruction methods in TrackVis (Wang et al., 2007;Wedeen et al., 2008), using dODFs and DTs (denoted DTKODF and DTKTensor, respectively). Standard streamline tractography is based on the Frenet equation: Here, r(s) is the curve, s the arc length and t(s) the unit tangent vector of the curve at s. The equation is solved iteratively using Euler integration. In DTKODF, t(s) equals the local maximum of the dODF most closely aligned with the incoming direction (as in Equation 4), whereas in DTKTensor, t(s) is the principal eigenvector of the DT.
To further validate our method, we used the MICCAI 2009 Fiber Cup phantom dataset (Poupon et al., 2008(Poupon et al., , 2010Fillard et al., 2011). It resembles three coronal brain slices and was made of hydrophobic acrylic fibers with proportions similar to myelinated fibers. DW-MRI of the phantom was acquired with 3 × 3 × 3 mm 3 isotropic voxels and 64 diffusion sensitizing directions (b = 1500 s/mm 2 ) in addition to one b 0 volume on a 3 T Siemens Tim Trio MRI system, equipped with a whole body gradient coil (40 mT/m, 200 T/m/s) and a 12-channel receive only head coil, used in combination with the whole body transmit coil of the MRI system. The acquisitions were performed using a single-shot diffusion-weighted twice refocused spin echo echoplanar pulse sequence, with volume size 64 × 64 × 3 voxels TR/TE/FA/NAC = 5000 ms/94 ms/65 • /2.

SELECTION OF FIBER TRACTS
For the experiments on real data, we have chosen to study the corticospinal tract (CST), carrying information from the cerebral cortex to the spinal cord, and the forceps minor (Fminor), also named forceps anterior, which are trans-callosal fibers curving forward from the genu of the corpus callosum into the frontal lobes, carrying inter-hemispheric information. Both tracts were extracted for each subject at each of the three study waves (see above) by defining ROIs manually according to Wakana et al. (2007). For the CST, one ROI was drawn on the cerebral peduncle at the level of the decussation of the superior cerebellar peduncle. After displaying the fibers intersecting this ROI, the central sulcus and the projection to the motor cortex were identified and a ROI was delineated in the axial slice right after the bifurcation to the motor cortex and sensory cortex. The CST was defined in our experiments as the fibers intersecting these two ROIs (cf. Figure 3A). As mentioned in Wakana et al. (2007), the challenge with CST fibers continuing into the opposite brain hemisphere is solved by cutting the fibers at the inferior ROI. Fminor was extracted by drawing two ROIs in a coronal slice at the midpoint between the anterior tip of the frontal lobe and the anterior edge of the genu of the corpus callosum, one in each hemisphere, delineating the entire frontal lobe (cf. Figure 3B).

USING LINEAR MIXED-EFFECTS MODELS TO ASSESS DIFFERENCES BETWEEN TRACKING METHODS ON LONGITUDINAL DATA
For each of the tracts CST and Fminor, and each tracking method DTKODF, DTKTensor, and Spline, a linear mixed-effects (LME) model was fitted to the data, i.e., Here, FA ij , the continuous response variable in the model, is mean FA in given tract and for given tracking method in subject i (i = 1, . . . , N) at wave j (j = 1, . . . , n i ). In our case, we have N = 8 and three waves with complete data, i.e., n i = 3 for all i. Age ij is age (in years) of subject i at wave j, and a predictor variable in the model. The model parameters β 0 and β 1 are the fixed effects parameters. The variables b 0i and b 1i (i = 1, . . . , N) are the random effects parameters, assumed to be normally distributed with zero mean. They denote individual deviations in intercept and slope, respectively, from the group-level fixed effect. Finally, the random residual errors ij are assumed to be independent and normally distributed with zero mean and constant variance σ 2 . In the LME parameter estimation we used the restricted maximum likelihood (REML) criterion. Model Equation (7) allows for correlation among measurements obtained from the same subject at different waves after the fixed (population) effects have been taken into account.
We also compared the LME-model fits with simple linear OLS regression, assuming (naively) independence of between and within-subject measurements over time. The LME-model incorporates two sources of variation: a biological variation due to subjects, and a within-subject variability. In the comparison of tracking methods, we were interested in their influence on the fixed effect parameter estimates, on the variation between subjects, and on variation within subjects (residuals). For the longitudinal analysis we used R 3.0.0 (www.r-project.org) with the lmer() function of the lme4 package (http://lme4.r-forge. r-project.org) for LME model fitting, the ggplot2 package for making graphics (cf. Figure 4), and the stargazer package for generating tables (cf . Tables 2, 3).

TRACT-AVERAGED AND ALONG-TRACT ANALYSIS IN HUMAN DW-MRI
In our evaluation study we have examined longitudinal changes between wave1, wave2, and wave3 in the well-defined CST and Fminor fiber bundles. For the along-tract FA analyses, we employed the Colby MATLAB toolbox (Colby et al., 2012). Tracking was performed with step-length h = 1 mm (cf. Algorithm Spline), maximum angle θ = 60 • , and minimum FA threshold at 0.15. The results for CST and Fminor are presented in Table 1, listing mean and standard deviations across all subjects. Table 1 shows that FA values of CST extracted by the different tracking methods in general has stable values over time. DTKTensor resulted in the overall highest mean FA values across the waves, followed by DTKODF and Spline. Averaging the standard deviations, DTKODF had the lowest standard deviations, followed by Spline and DTKTensor (both 0.029). The Fminor tract showed slightly decreasing FA values during time. For this tract, DTKODF resulted in the overall highest mean FA values, followed by Spline and DTKTensor. Also in this case, DTKODF had the lowest standard deviation on average, followed by Spline and DTKTensor. To quantify the changes in FA values across the three study waves, the absolute yearly changes are listed in the two last columns. Spline and DTKODF CST FA values were slightly increasing, while DTKTensor tracking gave slightly decreasing FA values. The Fminor FA values were decreasing the most for Spline, followed by DTKODF and DTKTensor tracking, where the Spline tracking method resulted in a mean reduction of 0.006 FA units per year, or about 1.4% yearly reduction with respect to a mean age of 58 years at baseline.
The wave-by-wave summary statistics (cf. Table 1) was supplied by more detailed longitudinal data analysis in terms of linear mixed-effects model fitting as explained in the Materials and Methods section. Initially, an exploratory visual analysis was performed by producing spaghetti plots of the tract-averaged FA values versus subject age at each examination time point, separately for each tract and for each tracking method. In Figure 4 the tract-averaged FA values for wave1, wave2, and wave3, are connected with color-coded line segments for a given subject to indicate subject-specific change over time.
From Figure 4 an immediate observation is the age-resistant and slightly higher FA values in CST compared to Fminor. For the CST, the three tracking methods performed qualitatively very similar with a rather stable FA value at about 0.49, and where subject 27 was an outlier with very low FA values in all tracking methods and subject 20 had consistently high FA values. Also for the Fminor, the methods performed qualitatively similar showing (on average) a clear reduction in FA over time, and where subject 29 seemed to be a consistent outlier with generally higher FA values for age in all tracking methods. However, more differentiating results could be obtained using the LMEmodeling approach, enabling detailed quantitative assessment  linear regression gave a negative slope of −0.0013 (SE 0.0010). Moreover, mean FA in CST using Spline tracking assessed with LME model-fitting gave comparable slope as with simple OLS fitting, −0.00036 and −0.00047, respectively, while the intercepts were 0.52122 and 0.52213, respectively. Thus, the LMEmodel predicts slightly higher mean FA in the age-range 51-73 years compared to simple OLS regression. Finally, the Spline tracking method together with LME model-fitting yielded the highest mean FA of the CST in elderly subjects (>65 years) compared to the other tracking methods. Regarding the Fminor tract, both DTKODF tracking and Spline tracking gave the largest decrease in mean FA with age (fixed-effects slopes −0.00525 and −0.00474, respectively), while DTKTensor tracking yielded a more moderate change in mean FA over time with a fixed-effects slope of −0.00127. The individual linear random-effects estimates had a larger variation in slope and intercept for the Spline tracking method (Figures 4E,F) compared to the two other methods. In the order of DTKODF, DTKTensor, and Spline tracking, the random effects summary statistics are given in Table 4.
Visually and numerically, for the given LME model-fitting (Equation 7), this indicates that Spline tracking exhibits a slightly larger "sensitivity" to within-individual variation in time trajectories of mean FA for a tract, compared to DTKODF and DTKTensor tracking.
To further explore performance differences between tracking methods, we used along-tract analyses of FA. Comparing the along-tract FA curves in Figure 5, there are at least two interesting observations. Regarding the CST (Figures 5A-C), all tracking methods gave two local maxima, at about 30% and about 80%, respectively, of the normalized tract length. Furthermore, our Spline method had the lowest second local maxima. Figure 6 shows specifically the part of the fiber that is represented in the region of the second local FA maximum (encircled, one subject). Apparently, there is no striking difference between Spline and DTKODF. They both track fibers which are spreading toward the cortex. This is in contrast to the fiber bundle obtained from DTKTensor tracking, which is narrower and is much more spatially restricted. In Figure 7 we have depicted Frontiers in Neuroinformatics www.frontiersin.org July 2013 | Volume 7 | Article 13 | 7

Tract Random effect parameter
Tracking method the voxelwise reconstructed fiber directions in a comparable slice, using both tensor (upper panel) and dODF representations (lower panel). For CST proper there seems to be only minor differences, but outside the CST the dODF helps to reconstruct more (anatomically plausible) fiber directions (circumscribed areas in Figure 7). Regarding Fminor, results from all tracking methods demonstrate very similar curve shape patterns, with one maximum FA between two FA plateaus. The largest difference between the methods is the lower FA values obtained with the DTKTensor tracking.

DTKODF
The tract-averaged fiber analysis was computed by averaging the FA across all fiber coordinates. Additionally, an alternative voxel-averaged fiber analysis was conducted, where each voxel visited by at least one tract was counted only once, and the FA was averaged over these voxels. This approach showed similar results as the tract-averaged fiber analysis.

PHANTOM DATASET
Tracking in the brain phantom dataset was performed in a similar manner as for the human DW-MRI data, except from the steplength h, which was set to 1.5 mm. The evaluation of the tracking was done twofold for the phantom. First, we used a dedicated program (fibercup) for comparing the tracking results, available at the web-site of the MICCAI Fiber Cup 2009 (http://www.lnao.fr/ spip.php?rubrique79). The coordinates for 16 ground truth fibers are provided by the phantom dataset, in addition to 16 ROIs used to define fibers corresponding to the ground truth. In the fiber tracking method comparison, interpolating splines were used and the fibers were sampled using s = 1000 points. A distance measure (L 2 ) between two corresponding fibers were then computed for these points. The 16 fibers and the ROIs are depicted in Figure 8A. dODFs were computed in TrackVis, and tracking was accomplished by Spline, DTKODF and DTKTensor, respectively. The results from the comparison are presented in Figure 8B. The mean distance (i.e., average value of the L 2 metric estimated over the number of sampling points) between ground truth fibers and the fibers obtained by each of the three tracking methods, is slightly lower for our Spline method (followed by DTKTensor)-in terms of grand mean across the 16 ROIs (last row in Figure 8B). Among the fibers where our Spline method performed better, are the Fminor-mimicking tract (ROIs 8,9,and 10), and also the long-range CST-mimicking tract, crossing the midline (ROIs 11,15,and 16). Inferior results with Spline, however, occurred for tracts defined by ROIs 3 and 4.
The effect of varying the step-length h in our Spline method was assessed by performing tracking in the same data set. Figure 8C shows the distance measure (L 2 ) between the ground truth and our results using h = 1 mm, h = 1.5 mm, and h = 2 mm. The best performance using these three different step-lengths was obtained using h = 1.5 mm. A different evaluation framework for the phantom data was later proposed by Côte et al. (2012). The ROIs were then redefined and larger. The number of fibers was reduced to seven, where each fiber was defined by two ROIs (cf. Figures 9A,B),

FIGURE 7 | The differences in the reconstruction techniques applied to one of our datasets, showing the principal eigenvectors of the diffusion tensors (top), and the dODF local maxima (bottom).
Details of the rectangular inserts are shown to the right. Along the CST, there are a few neighboring voxels where the dODF reveals more than one fiber direction. These are indicated by white circles and show that the voxel-wise dODF provides additional information about fiber directions, compared to the diffusion tensor representation.
mimicking in vivo WM fiber bundles as indicated in Figure 9C. Figures 9D-F, displays the fibers resulting from tracking with Spline (D), DTKTensor (E), and DTKODF (F), connecting the pairs of ROIs shown in Figure 9B, where the underlying dODFs were computed in TrackVis as in the previous evaluation. Results obtained comparing Spline with varying step-length h (1, 1.5, and 2 mm). The phantom dataset contains 16 ROIs, and the fibers intersecting each of these ROIs constitute fiber bundle 1-16. For the comparison, the average spatial metric using 1000 points were employed (computed with interpolating splines). The lower the metric value, the closer the tracking is to ground truth. Our Spline tracking method ( Figure 9D) was able to find fibers connecting all seven pairs of ROIs, where each of them are visually comparable to the ground truth ( Figure 9A) collectively in Figure 10H, whereas the ground truth fibers are shown in Figure 10I. Two of the fiber bundles deviate somewhat from the ground truth, i.e., bundle 3 ( Figure 10C) and 6 ( Figure 10F). Our fiber 3 consists of two bundles, one following the correct path of fiber 3 in Figure 10I, the other running down along fiber 2 and up again along fiber 6. A few tracts in fiber 6 bifurcate in the middle segment and continue along fiber 7, before turning laterally, joining the end point ROI of fiber 6. Also, some of the tracts in fiber 6 had a sharp turn close to this ROI, and followed fiber 3 for a short distance. Tracking with DTKODF ( Figure 9F) resulted also in fibers connecting 7 out of the 7 pairs of ROIs, all of which are visually comparable to the ground truth fibers. Moreover, this tracking method had similar deficiency as the Spline method in that fiber 6 had a bifurcation with a branch along each of fibers 5 and 7. Finally, the DTKTensor tracking method (Figure 9E) resulted in fiber connections only between 4 out of the 7 pairs of ROIs, as fibers 2, 3, and 7 were partially missing. Thus, for the phantom dataset, the DTKTensor tracking method was clearly inferior to DTKODF and Spline tracking.

DISCUSSION
The development of HARDI acquisitions and introduction of the orientation distribution function has improved the performance and spatial resolution of in vivo investigations of "WM integrity" and the assessment of fiber bundle connections in the brain. However, the new advanced acquisition protocols are time-consuming and most state-of-the-art tracking methods are still solving the Frenet equation using Euler or Runge-Kutta schemes (Basser et al., 2000)-although the tangent tracking vector is frequently taken from the ODFs, rather than the tensor representation when recordings have many diffusion sensitizing directions. Despite these efforts, WM tractography, when applied to clinically feasible DW-MRI acquisitions, is still considered inaccurate with false negative and false positive tracts, and the underlying data seems not fully exploited (Li et al., 2012;Feigl et al., 2013). This opens up for further investigations. We have therefore designed and tested a new tracking method for WM fibers, denoted Spline.
Our tracking method is based on the assumption that gross WM consists of smooth tubular structures (Gössl et al., 2002;Jones, 2011;Catani and de Schotten, 2012). This assumption is partly implicit and partly different from working principles of methods merely solving the Frenet equation. With the traditional approach, one can experience sharp turns arising from noise and inaccuracy in the recorded data and in the derived vector fields. This might lead to termination of the tracking, or to erroneous "jumping" and spurious tracts. Normally, these methods incorporate a smoothing spline applied post-tracking. However, these approximating splines will not influence significantly the path taken by a tract across a region of noisy voxels because the major path has already been decided by the tracking algorithm. Our method was designed to counteract these deficiencies. The tracking is iteratively directed by a smoothing spline, which likely decreases the influence of noisy voxel data on the path of the tract. In our work we have investigated how this strategy affects tracking results in data obtained by clinically applicable recording schemes, and in HARDI data. In both cases, the local water diffusion was chosen to be represented by the previously established dODF as calculated by TrackVis (Wang et al., 2007).
To evaluate the performance of our approach, we compared our results with those obtained with two widely used tracking methods, using both real data and data from a phantom. The phantom data, obtained from the MICCAI 2009 Fiber Cup contest, has a known ground truth regarding its fiber connections and been processed by many competing tracking algorithms (Li et al., 2012;Pontabry et al., 2013). It could therefore be argued that the phantom data alone suffices for the evaluation of our Spline tracking method. Although mimicking a brain coronal slice, the water diffusion properties of the test object may not, however, be representative for water diffusion in white matter tissue of the living human brain. We therefore decided to also include real data in our evaluation, despite the fact that comparative tractography on real data is not straight-forward since the ground truth is not known. Also, post mortem brain and anatomical preparations can be useful in the evaluation of tracking methods (Dyrby et al., 2007;Li et al., 2012;Prats-Galino et al., 2012), but a ground truth is still lacking. On the other side, several WM fiber tracts are typically identified in wholebrain tractography examinations (Bassett et al., 2011;Catani and de Schotten, 2012;Côte et al., 2012;Guevara et al., 2012) and have been subject to extensive studies by many neuroimaging researchers. To better assess the applicability and performance of our Spline tracking method in real life situations, we selected two well-described fiber bundles for the evaluation: the cortocospinal tract, and the frontal lobe forceps minor. Moreover, these tract evaluations were done on longitudinal data, employing clinical data acquisition schemes. This enabled us to compare tracking algorithms with respect to tract detection/segmentation, and quantification of biological change (in terms of tract-wise FA values) related to aging using linear mixed-effects model fitting.
Large cross-sectional studies involving segmentation of the CST have shown a minor reduction in FA values with age (Westlye et al., 2010;Lebel et al., 2012). Westlye et al. (2010) reported a slight FA decline in CST from approximately 0.58 at the age of 40 years, to approximately 0.56 for subjects at the age of 80. The study by Lebel et al. (2012) reported a reduction in FA from ∼0.52 to ∼0.50 for the same age difference. In our experiments with longitudinal data, the tract-averaged FA obtained from eight participants with age range 52-66 years, using Spline, DTKODF and DTKTensor tracking, showed relatively constant CST FA during the 6 years of aging (mean FA in the range 0.49-0.50) and therefore supports the conclusions from these studies of near constant mean FA of the corticospinal tract into old age. The same investigators (Westlye et al., 2010;Lebel et al., 2012) reported a clear reduction in mean FA of the Fminor fiber bundle from 0.55 to 0.47 (Westlye et al., 2010) and from 0.54 to 0.50  between the age span from 40 to 80 years. Our results, using Spline and DTKODF tracking are in accordance with these findings, where we observed a reduction in mean FA of Fminor from 0.44 to 0.41 during the 6 years of aging in our sample. The unlikely finding of a slight increase in average FA of Fminor from wave1 to wave2 using DTKTensor tracking (Table 1), followed by a slight decrease from wave2 to wave3, although not statistically significant, might be due to noise sensitivity of the local tensor representation and the FACT algorithm implemented in DTKTensor, being highly dependent on a strong principal diffusion direction. Our Spline tracking method and DTKODF tracking, however, demonstrated a pattern of monotonic decrease in mean FA of the Fminor tract from wave1 to wave3.
The along-tract analysis revealed some interesting patterns. Firstly, using Spline, DTKODF, and DTKTensor tracking, respectively, the distinct local FA minimum located at approximately 60% along the length-normalized CST tract (Figures 5A-C) seems to correspond to the intersection of the CST and fibers from the corpus callosum, which could explain the reduced FA due to crossing fibers. This might advocate the use of more reliable and sensitive diffusion parameters, like generalized fractional anisotropy (GFA) (Tuch, 2004;Fritzsche et al., 2010;Jones et al., 2012). We used, however, FA as the diffusion metric since it has been shown to be satisfactory for b = 1000 (Fritzsche et al., 2010). All tracking methods resulted in similar FA curves for the first local maximum (at about 30% of trajectory length), but the second local FA maximum (at about 80% trajectory length), was clearly lower for Spline tracking compared to DTKODF and DTKTensor. This local discrepancy between tracking methods was further investigated and Figure 6 shows the 60-100% segment of the length-normalized CST for a single subject. By construction, we expect Spline tracking to allow fibers to spread and cover a larger volume than the fibers obtained with DTKTensor tracking. Since FA in general Frontiers in Neuroinformatics www.frontiersin.org July 2013 | Volume 7 | Article 13 | 11 is lower in the outskirts of the CST, closer to the cerebral cortex, this would explain the lower local mean FA values along this part of the tract. The local deviation in FA between Spline and DTKODF tracking is harder to explain since the fibers in this segment of CST seem to have a similar spreading pattern (Figures 6B,C). However, ground truth is not attainable for the detailed fiber pathways of CST (and the corresponding FA values) in a specific brain and region, and no firm conclusion can be made about local accuracy of the tracking methods being compared. Visually, Spline and DTKODF seem to give better results than DTKTensor, since tracking by the latter method produces a rather concentrated and spatially restricted bundle of CST fibers toward the cortical sheet which is anatomically less plausible. Secondly, the along-tract analysis of Fminor using Spline, DTKODF and DTKTensor tracking, respectively, gave comparable FA curves with closely corresponding peaks and valleys (Figures 5D-F). However, the DTKTensor tracking produced consistently lower along-tract FA values across the three waves compared to that obtained with the two other methods. To explore whether this was due to difference in tract density or overall volume positioning, we conducted the alternative, voxel-averaged fiber analysis. This confirmed the finding of lower FA along Fminor using DTKTensor and indicate that the DTKTensor tracking method applied to Fminor will segment a partly different white matter volume than that obtained with the two other tracking methods. Further investigations are thus needed regarding the reduced FA along Fminor using DTKTensor tracking.
Along-tract analysis have been used in several recent DWI-MRI studies (Jones et al., 2005;Colby et al., 2012;Hodneland et al., 2012;Yeatman et al., 2012;Hill et al., 2013). In our alongtract analysis of CST and Fminor, the peaks and valleys of the FA curves were very similar to those reported by Colby et al. (2012, Figure 5). A major difference, however, is that FA in our study was 0.1-0.2 units lower along the CST than that reported in Colby et al. (2012). For Fminor, the along-tract FA discrepancy was much less (<0.1 units in lower values), where the Spline and the DTKODF tracking methods produced FA curves closest to those in Colby and coworkers. It should be noted that the along-tract curves reported in Colby et al. (2012, Figure 5) are from a single subject, which can explain the discrepancies in FA. Furthermore, the DWI-MRI acquisitions in Colby et al. (2012) were performed on a 3T scanner, with different pulse sequences, and with voxel-wise tensor representations, also giving possible rise to the observed bias in FA.
For evaluation of tracking performance on the brain phantom data, two different strategies were adopted. Applying the first strategy (Fillard et al., 2011), we showed an improved fiber tracking ability of Spline, compared to the DTKODF and DTKTensor methods. Comparing our results with those from the MICCAI 2009 Fiber Cup contest, our Spline tracking method could not entirely match the best results. However, it should be noted that some of the methods being used in this contest included a fiber bundle selection (e.g., Reisert et al., 2009) where only a selected subset of the resulting fibers were used in the evaluation. Following the second strategy (Côte et al., 2012), we achieved good results with Spline and DTKODF. Both methods were able to find 7 out of 7 fibers, whereas DTKTensor was able to reconstruct only 5 out of 7 fibers. We suspect the poorer results with DTKTensor to be due to its voxel-wise diffusion tensor representation, instead of dODFs as in DTKODF and Spline, misleading the tracking in regions with crossing fibers. This observation supports the use of dODFs rather than DTs for practical tracking applications.
There have been only few reports (e.g., Prckovska et al., 2010) on the use of ODF to DWI-MRI recordings with a low to moderate angular resolution, i.e., 20-30 diffusion sensitizing directions.
We have in this work demonstrated that the orientation distribution function can be better suited to represent regions of crossing fibers than using a local tensor representation-both for the brain phantom data and for (clinical) data with relatively few (<30) diffusion sensitizing directions. In Figure 7 we have depicted the voxelwise reconstructed fiber directions in a slice from one of the subjects in the longitudinal study, using both tensor (upper panel) and dODF representations (lower panel). Most studies with DWI-MRI data similar to ours have been using diffusion tensor fiber tracking (cf. DTKtensor). This study indicate that dODF reconstruction can be employed even at low angular resolution of 25 diffusion sensitizing directions recorded within clinically feasible acquisition times.

SUMMARY AND CONCLUSION
We have proposed a new tracking method (Spline) as an alternative to the traditional streamline tracking method based on the Frenet equations. This is done by imposing a contextual constraint via an iterative smoothing spline, using voxel-wise estimation of the diffusion orientation distribution function, rather than a tensor representation. The motivation and design principle of Spline was to improve the resistance against erroneous tracking in noisy voxels. The other major contribution is the incorporation of a comprehensive evaluation framework for comparing tracking methods, spanning both real DWI-MRI recordings (longitudinal data with linear mixed-effects model fitting), detailed along-tract analysis of diffusion parameters (FA), and the use of DWI-MRI acquisitions from a test object, mimicking a brain, with ground truth for its fiber connections. Using this evaluation framework, we compared Spline with two other tracking methods (DTKODF and DTKtensor) provided in widely used software for DWI-MRI analysis. We found some major differences between Spline and two other tracking methods, that might have practical implications. This regards results from along-tract analyses of fractional anisotropy, in particular for the corticospinal tract and also for the forceps minor. Testing our method on data from the phantom with a ground truth, the performance of our method was generally better than the two widely used tracking methods being compared.
We find these design considerations and results encouraging, and further algorithmic improvements and evaluation procedures are planned. Additional metrics and evaluation strategies are explained in Fillard et al. (2011), Côte et al. (2012, and the use of these would provide a more thorough evaluation of our results. To our disposal we also have DW-MRI Frontiers in Neuroinformatics www.frontiersin.org July 2013 | Volume 7 | Article 13 | 12 data from nearly 100 participants across three study waves obtained in the aforementioned study of aging, genetics and cognition (Westlye et al., 2011;Ystad et al., 2011;Hodneland et al., 2012). Moreover, it would be important to study and reconstruct a wider collection of fiber tracts, other than the CST and Fminor, in order to quantify regional differences that can be attributed to the aging process, and to better clarify the influence of tracking method being used. For such longitudinal investigations, where also results from cognitive testing will be included, proper statistical methods are essential, and linear and non-linear mixed-effects models should be methods of choice.