Evaluation of Diffusion-Tensor Imaging-Based Global Search and Tractography for Tumor Surgery Close to the Language System

Pre-operative planning and intra-operative guidance in neurosurgery require detailed information about the location of functional areas and their anatomo-functional connectivity. In particular, regarding the language system, post-operative deficits such as aphasia can be avoided. By combining functional magnetic resonance imaging and diffusion tensor imaging, the connectivity between functional areas can be reconstructed by tractography techniques that need to cope with limitations such as limited resolution and low anisotropic diffusion close to functional areas. Tumors pose particular challenges because of edema, displacement effects on brain tissue and infiltration of white matter. Under these conditions, standard fiber tracking methods reconstruct pathways of insufficient quality. Therefore, robust global or probabilistic approaches are required. In this study, two commonly used standard fiber tracking algorithms, streamline propagation and tensor deflection, were compared with a previously published global search, Gibbs tracking and a connection-oriented probabilistic tractography approach. All methods were applied to reconstruct neuronal pathways of the language system of patients undergoing brain tumor surgery, and control subjects. Connections between Broca and Wernicke areas via the arcuate fasciculus (AF) and the inferior fronto-occipital fasciculus (IFOF) were validated by a clinical expert to ensure anatomical feasibility, and compared using distance- and diffusion-based similarity metrics to evaluate their agreement on pathway locations. For both patients and controls, a strong agreement between all methods was observed regarding the location of the AF. In case of the IFOF however, standard fiber tracking and Gibbs tracking predominantly identified the inferior longitudinal fasciculus that plays a secondary role in semantic language processing. In contrast, global search resolved connections in almost every case via the IFOF which could be confirmed by probabilistic fiber tracking. The results show that regarding the language system, our global search is superior to clinically applied conventional fiber tracking strategies with results similar to time-consuming global or probabilistic approaches.


Introduction
The neurological network of the human brain comprises complex functional units such as the language system, which consists of multiple functional areas connected via white matter pathways. These subcortical networks are of major interest for preoperative neurosurgical planning and intra-operative guidance close to eloquent structures to prevent neurological deficits such as post-operative aphasia.
Insight into the general architecture of the language system is provided by various anatomo-functional connectivity studies with invasive and non-invasive techniques. Intra-operative subcortical mapping (ISM) locates pathways by inducing temporal neurological dysfunctions with electrical stimulation [1]. In contrast, with image data from functional magnetic resonance imaging (fMRI), functional nodes and their causal associations can be identified by various types of correlation analyses [2][3][4].
Both types of functional connectivity studies have been supplemented by diffusion tensor imaging (DTI) [3,5], which is a magnetic resonance imaging (MRI) technique that measures the diffusion characteristics of water molecules [6]. In brain parenchyma, anisotropic diffusion is observed in aligned microstructures such as white matter pathways, whereas isotropic diffusion is dominant in gray matter and other tissues. To characterize these differences, DTI encodes information about the direction and intensity of diffusion by a 2 nd -order tensor D. Based on D, the trajectories of neuronal pathways can be reconstructed by fiber tractography.
For the detection of language related white matter tracts, ISM can be applied during awake surgery. In such a procedure, the patient undergoes repetitive language testing while the subcortical surface of the resection cavity is being stimulated. This invasive method can also be used to verify the results of DTI fiber tractography [1,5]. In contrast, with the combination of fMRIbased statistical functional mapping and DTI-based neuronal fiber reconstruction, functional and anatomical connectivity can be resolved without additional invasive steps. Recently, this combination has become an established method for pre-operative planning and stereotactic neuronavigation in brain tumor patients [7][8][9] because it provides patient-specific information about individual anatomy and displacement effects resulting from tumor growth [7,10]. Due to brain shift and tumor resection during surgical procedures, shifting of fibers may occur which requires intra-operative tractography to ensure maximum tumor resection and preservation of neurological function [11,12]. To meet clinical requirements, accurate and fast fiber reconstruction methods that can be categorized into local and global approaches are required: Several local approaches based on the principle of fiber tracking such as streamline propagation (SP) [13], tensor deflection (TD) [14] and the FACT algorithm [15] have been devised. They make restricted use of local tensor information by deterministically computing a fiber direction at each considered point and propagate fibers sequentially in independent steps. A general limitation of this concept is that minor inaccuracies in the chain of local decisions can accumulate and significantly affect the final result. Therefore, local approaches have difficulties with crossing or branching pathways or with pathways in regions of low anisotropic diffusion.
In contrast, global approaches consider a larger subspace of the diffusion tensor field at once to infer global characteristics of the diffusion process and to model uncertainty in the diffusion data. For instance, fast marching methods evolve a diffusion front according to the underlying tensor profiles and compute fibers with an assigned level of confidence by backward processing based on gradient descent [16,17]. Another approach applies Monte-Carlo simulation of fiber tracking to provide a probabilistic interpretation of tensors [18]. Multiple fibers originating from one seed point are created by randomly sampling probability distributions defined on the diffusion tensor. From probability maps computed by counting the number of fibers crossing each voxel, conclusions on the architecture of the neuronal network can be derived. In [19], two connectivity maps are computed for two seed regions of interest (ROIs). Based on the number of fibers that cross each voxel in the same or the opposing direction of the major eigenvector, areas of connecting and merging pathways can be separated. By computing a probability index for connections to both seed regions (PIBS), this approach (denoted in the following as ConProb) enforces the expressiveness of probability maps for structural connectivity as compared to probabilistic tractography using one seed region. Another approach based on particle simulation models the diffusion field by external energies and the location and direction of and connection between short fiber segments by internal energies [20]. Optimized by a Markov chain Monte Carlo method and a tailored local tracking algorithm, this approach (denoted as GibbsT) produces optimally distributed fibers along the white matter structures [21,22].
Either of these global approaches computes fibers or probability maps of high quality, however, the computational effort is beyond a clinically acceptable time frame. Other global approaches introduce the concept of target-orientation explicitly by reconstructing connections between two predefined ROIs. This reduces the amount of data while the level of quality is maintained. The guided fiber approach [23] extends a standard fiber tracking algorithm by a random sampling strategy that is constrained to a narrow region surrounding the expected location of the pathway. Since multiple alternative directions are considered, complex pathway arrangements in regions of low anisotropic diffusion can be resolved. Even more alternative directions are considered by a search for minimum cost paths that is combined with a Bayesian framework for tensor field regularization [24].
Our previously presented approach [25] applies a global search (GS) to evolve optimal connections between two functional areas according to a diffusion dependent cost function. As this function is sampled for a large number of directions at each point of the search space, complex pathway configurations can be resolved. Moreover, GS is based on the A * shortest path algorithm, whose time complexity is optimal [26], i.e. it is guaranteed that the best possible solution is found with the smallest computational effort. For this reason, GS is particularly suitable for clinical application.
In this work, GS was compared with the aforementioned fiber tracking algorithms, SP and TD, which are commonly used in clinical practice [7,8]. Two tractography algorithms, ConProb and GibbsT, which employ stronger concepts such as explicit regularization to model noise in diffusion data and global optimization with particle simulation, were additionally included in the evaluation.
To evaluate the reconstruction capabilities of the considered methods, two neuronal pathways, the arcuate fasciculus (AF) and the inferior fronto-occipital fasciculus (IFOF), were selected. These pathways are of particular interest for neurosurgical interventions as they connect the Broca and Wernicke areas by a dual route that is part of the language processing system. The AF, which is a part of the superior longitudinal fasciculus, is generally regarded as the dorsal language pathway responsible for phonological processing [3,27,28]. In contrast, the exact anatomical correlate of the ventral semantic pathway is not clearly established. Various fiber tracts have been described, including the inferior longitudinal fascicle (ILF) [3,29], middle longitudinal fascicle [3] and the IFOF to subserve this connection [29]. The function of the ILF for language processing was recently described as dispensable by a subcortical stimulation study [30], whereas the function of the IFOF was found to be crucial. Moreover, for both the existence of the IFOF in humans and its role in language processing, there is a wealth of anatomical [31,32], neurosurgical [28,30,33], and diffusion tensor imaging studies [34,35]. The relevance of this tract for semantic processing was additionally confirmed by a study on young adults with dyslexia showing significant correlation with DTI [36]. In a dissection study of brains of rhesus monkey, the existence of this pathway has been challenged [37], but as white matter tracts serving the language processing might be underdeveloped, conclusions following the transition between species might not be strong enough to counter recent research.
Despite of this controversy, the IFOF was selected to challenge the above tractography approaches as its reconstruction in the context of semantic language processing has been shown to be a difficult task. A combined fMRI and DTI study for language relevant cortical sites and white matter pathways of healthy subjects identified problems of fiber tracking to reveal the terminals of the IFOF in the middle temporal gyrus that are relevant for semantic processing [29]. This was explained by the stronger occipital terminals and the tendency of streamline tractography to follow the other major neuronal pathways in this region. Due to this difficulties, assessing the capabilities of GS, SP, TD, GibbsT and ConProb in the context of brain tissue affected by tumors is a challenging research question.
For this purpose, fiber tracts reconstructed by all methods on image data of both patients undergoing tumor surgery and healthy controls were assessed by a neuromedical expert under the qualitative aspect of anatomical feasibility. Then, the agreement of the reconstruction methods on the location of neuronal pathways was evaluated quantitatively by distance-and diffusion-based similarity metrics. In addition, all methods were validated on ground-truth data of simulated DTI software phantoms comprising fiber configurations of varying complexity.

Ethics Statement
This study was approved by the ethics committee of the University of Erlangen-Nuremberg according to the ethics committee proposal number 3578. This included all clinical investigations with pre-and intra-operative MRI and functional neuronavigation, and subsequent data analyses for all patients and controls. Written informed consent was obtained from all patients and controls. All clinical investigations were conducted according to the Declaration of Helsinki.

Image Data
Image data originated from 25 right-handed patients (11 female, 14 male, with mean age of 41.6+12.6 years) and 6 right-handed controls (1 female, 5 male with mean age of 34.7+9.9 years, 4 with tumors in right hemisphere, 2 healthy probands), all with lefthemisphere dominance. The two groups did not significantly differ in age (p~0:22 by non-parametric tests of Mann-Whitney [38] with significance level a~0:05). The patient group comprised subjects with primary brain tumors (WHO grade I-IV). To support the grading, the volume of tumors was assessed by manual or semi-automatic intensity-based segmentation. For enhancing tumors, post-contrast T1-MPRAGE scans, acquired as mentioned below, provided images with appropriate contrast. For nonenhancing tumors, images were acquired with T2-weighted scans.
All tumors were located in the white matter of the left dominant hemisphere along the course of the language pathways. Even with larger distances between tumors and neuronal pathways, the influence of tumors can be strong in terms of both fiber displacement and impairment of associated functions. Since functional deficits were observed also in these cases, patients were included with distances from 1 up to 30 mm. The controls were considered to assess the similarity of fiber tracts reconstructed with different methods in unaffected brain tissue.
For all subjects, the image data was acquired between 2005 and 2009 with a Magnetom Sonata Maestro Class 1.5 Tesla scanner (Siemens Medical Solutions, Erlangen, Germany) equipped with a gradient system with a field strength of up to 40 mT/m (effective 69 mT/m) and a slew rate of up to 200 T/m/s (effective 346 T/ m/s).
Data Acquisition -fMRI. During the acquisition of functional data using fMRI, visual stimulation was performed because of the reduced signal-to-noise ratio as compared to auditory stimulation. For each subject, 30 rest and activation measurements were obtained according to a block paradigm (''boxcar'') with a total duration of 15 minutes. Activities included reading, naming of objects, answering questions and inferring verbs and sentences from given nouns [9].
The functional data was acquired using an echo-planar imaging (EPI) sequence (parameters: TR = 2470 ms, TE = 60 ms, image matrix 64|64, slice thickness 3 mm, resolution 3|3 mm, 16-25 slices). These measurements were interpolated to 75 slices with 1 mm in plane resolution. An image based prospective acquisition correction for head motion was used [39]. Hereby, interpolation was done in the k-space.
Statistical functional mapping on a voxel-by-voxel basis of activated functional areas was performed using BrainVoyager [2]. Activation maps were individually determined by analyzing the correlation for each voxel between signal intensity and a square wave reference function according the paradigm. Voxels exceeding a threshold (typically 0.40) were displayed if at least six contiguous voxels of the measured slices built a cluster to eliminate isolated voxels. The functional slices were aligned to the MPRAGE images with 160 slices of 1 mm slice thickness and resolution obtained in the same patient position.
According to this mapping, areas associated with language functions were typically located in the frontal and temporoparietal lobes of the left hemisphere (right-handed population) that correspond to the Broca and Wernicke areas. These areas could be determined for all subjects successfully and were selected as ROIs for all tractography approaches.
Data Acquisition -DTI. For the DTI measurements, a single-shot spin-echo diffusion weighted EPI sequence was used (parameters: TR = 9200 ms, TE = 86 ms, b high = 1000 s/mm 2 , b low = 0 s/mm 2 , field of view 240 mm, 1500 Hz/Px bandwidth, measurement time: 5 min and 31 sec). For each slice, five measurements were acquired and averaged by the MRI scanner. The image resolution was 128|128|60 voxels with a respective voxel size of 1:875|1:875|1:9 mm 3 . One reference image without diffusion sensitization as well as six images with diffusion sensitization in non-collinear directions were obtained according to the oblique double gradient encoding scheme ((+1,1,0), (+1,0,1), and (0,1,+1)) [40]. Compared to numerically optimized gradient directions, this scheme provides relatively low error rates for the computation of tensors and FA values, but a slightly higher sensitivity to measurement noise [41]. Moreover, compared to acquisitions using 30 diffusion gradient directions, the influences on diffusion values and tractography results are within an acceptable range [42].
Additionally, if the tumor exhibited contrast enhancement on the preoperative images, the T1-weighted axial spin echo sequence was repeated after intravenous application of contrast media (0.2 ml/kg gadolinium diethylenetriamine penta-acetic acid). In nonenhancing tumors, an additional T2-weighted turbo spin echo sequence (parameters: slice thickness 4 mm, field of view 230 mm, TR = 6490 ms, TE = 98 ms, scan time: 5 min 59 sec at three acquisitions) was acquired.

Simulated Data
To validate the fiber reconstruction methods on ground-truth data, three DTI software phantoms were selected from a database with fiber configurations of different complexity including branching and kissing pathways as well as a spiral with more than three turns [43]. While the character of branching and kissing pathways is similar to major pathways such as the pyramidal and optical tracts, strong curvatures are typical for pathways with cortical terminations such as the language pathways along the AF and IFOF.
According to a study on the influence of signal-to-noise ratio (SNR) on the pathway reconstruction [44] that used a similar scanner type and DTI protocol, SNR varies between 15 and 25. Therefore, phantom images were selected with a simulated SNR of 15 and 30. To allow for a direct comparison of similarities between reconstructed tracts, the voxel size was adjusted to comply with the image data.

Global search (GS)
The approach of GS to anatomical connectivity [25], that is compared in this study with SP, TD, GibbsT and ConProb, searches for minimum cost paths by evolving a search tree on a structured grid between two functional areas identified as ROIs by fMRI. The cost function originally minimized in [25] is based on the propagation probability function p(v,D)~(DDvDD{l 3 )=l 1 with v denoting the current propagation direction and l i being the eigenvalues of D in descending order. To prevent a bias towards isotropic tensors, the isotropic fraction l 3 is subtracted, while the shape of the probability profile is still maintained. The cost for each path considered by the A * shortest path algorithm is then defined by accumulating the negative proportional function c(v,D)~1{p(v,D) over all path segments.
In this work, this cost function is extended in two ways: At first, the divergence div(v,D) from the principal diffusion direction (Eq. 1) is incorporated to improve guidance through regions of low anisotropic diffusion close to functional areas and tumor infiltrated tissue. Since the principal diffusion direction is not well-defined in isotropic tissue, the geometric shape of tensor ellipsoids [45] has to be considered. For linear ellipsoids, divergence is defined as the cosine of the angle between v and the first eigenvector e 1 , whereas for planar and spherical ellipsoids the cosine of the angle between v and the plane spanned by the first two eigenvectors is considered. The second extension includes fractional anisotropy FA(D) (Eq. 2) [6] to support guidance in anisotropic regions. All metrics are then combined into a new cost function c'(v,D) (Eq. 3).

div(v,D)~c
os(%(v,e 1 )) : lin: cos(p=2{%(v,e 3 )) : pl:=sph: ð1Þ The cost function is additionally supported by thresholds for the maximum bending angle between adjacent fiber segments and the minimum FA value at sample points. Thereby, the reconstruction of anatomically unlikely fibers of strong curvature or of fibers intersecting isotropic regions can be avoided. Moreover, anatomical knowledge can be incorporated to reconstruct specific pathways by restricting the search space to an adjustable bounding box. This is especially useful for the dual route of the language system, since anisotropic diffusion is different along the AF and IFOF [46] and may also depend on the location of the tumor [10].

Deterministic fiber tracking
The algorithms SP and TD evaluated in this study are clinically approved and implemented in various DTI applications [7,8]. Fiber tracking is initialized at seed points with sufficiently high FA value selected from the whole volume. From the resulting set of fibers, only those that intersect both ROIs are kept. In this way, propagation becomes more robust, as it is not initialized in functional areas with unclear diffusion direction, but in regions of anisotropic diffusion referring to white matter.
Fibers are propagated in both ortho-and retrograde directions and are terminated as soon as the FA value drops below a predefined threshold. To propagate fibers from the current location, different schemes are applied by SP and TD: In case of TD, the diffusion tensor is multiplied with the fiber direction of the previous fiber segment to obtain the next fiber point. In regions of high anisotropy, fibers are deflected towards the major eigenvector, whereas in isotropic regions, the incoming direction is maintained or only slightly affected by the tensor [14]. In contrast, fiber tracking with SP applies the 4 th -order Runge-Kutta streamline integration scheme that linearly combines the major eigenvectors of the interpolated diffusion tensors at four sample points to determine the next fiber point [13].

Similarity metrics
Several similarity metrics for reconstructed fiber tracts were proposed in the literature, which are either applied to cluster fibers or to evaluate approaches for pathway reconstruction by comparing the divergence of fibers from ground-truth data or from the results of other methods. In general, the similarity of two fibers can be quantified by their spatial displacement [21,[47][48][49][50], by the angular deviation of a direction [21], or by geometric shape properties such as curvature [21,50]. Since SP, TD and GibbsT inherently generate smooth fibers, whereas fibers of GS are gridbound, assessing the curvature or angular deviation would always provide favorable results for the former methods. Therefore, the similarity metric S avg (F ,G) presented in this section is based on the distance between fibers.
Geometric similarity. The spatial similarity of the resulting fiber tracts F and G of two methods, defined as S avg (F ,G) (Eq. 4), measures their deviation by averaging distances among a set of closest fiber pairs defined by relation R f (F ,G) (Eq. 6). The set contains all fibers in F , each one paired with the closest fiber in G, and vice versa, whereas closeness is measured in terms of the similarity metric for fiber pairs S p (f i ,g j ) and the trimming method tr(f ,g) presented later on. Since duplicate fiber pairs are avoided by R f (F ,G), the effect of outliers is emphasized as compared to using all fiber pairs. In addition to the averaging metric S avg (F ,G), S min (F ,G) (Eq. 5) determines the closest fiber pair.
S min (F ,G)~min To measure the distance between two fibers f and g, the relation of closest points R p (f ,g) (Eq. 7) assigns to each point a i on fiber f the closest point b j on fiber g and vice versa. The mean of closest point distances S p (f ,g) (Eq. 8) is then defined as the average distance between these closest points with d(a i ,b j ) denoting the Euclidean distance DDa i {b j DD. In contrast to the mean of closest point distances defined in [50] that creates symmetry by averaging two similarity values of a bi-directionally applied nonsymmetric metric, R p (f ,g) excludes duplicates from the resulting set of closest points to emphasize outliers.
Normalization of fiber pairs. Prior to applying the similarity metric S p (f ,g) to every fiber pair considered by S avg (F ,G), two normalization steps are carried out to account for two geometric differences in the results obtained by GS, SP, TD and GibbsT. Thereby, the trimming method tr(f ,g) reduces each fiber pair (f ,g) to (f t ,g t ).
While both deterministic fiber tracking approaches propagate in arbitrary directions with constant step length, fibers of GS are grid-bound with varying step length and fibers of GibbsT are constructed from points with arbitrary distances. Therefore, fibers of GS and GibbsT are resampled first to reduce the deviation of point sets introduced by the different propagation methods. This is achieved by replacing the original points of the fibers of GS by points of equal distance corresponding to the constant step length of SP and TD. In contrast to the smoothing of fibers with B-splines in [21], this normalization approach does not spatially modify the results.
Secondly, overlapping fiber segments have to be identified to avoid biasing of the similarity metric, as fibers of GS terminate in ROIs, whereas fibers of SP, TD and GibbsT are not restricted by them. For two fibers f and g selected from two fiber tracts F and G, the overlapping segments f t and g t are determined similarly to the trimming approach in [47]. At each end of the fibers, the closest point pairs are considered as illustrated in Figure 1  This approach is suited to cut protruding endings of otherwise coherent fibers, as the overall structure of the fibers is not altered. However, in case of diverging fibers as illustrated in Figure 1 (bottom), more fiber segments than necessary would be excluded. Therefore, the trimming approach was extended as follows: Again, the two closest point pairs at each end are considered, but trimming is applied to at most one of the fibers at each end under the condition that only one of the closest point pairs refers to an intermediate point, while the other refers to an end point.
Diffusion similarity. In addition to geometric properties, similarity in terms of diffusion is evaluated for each pair of fiber tracts as FA avg (F ,G) (Eq. 9) to reveal differences between groups, methods and pathways. For both fiber tracts F and G, the average FA values are computed for their trimmed and correlated fibers according to

Experimental Design
For detailed assessment of the capabilities to reconstruct anatomically feasible fiber tracts in the presence of nearby tumors, the clinically approved methods SP and TD, and the global approaches GibbsT and GS were compared with each other using the similarity metrics defined above. These methods were additionally compared with ConProb by visually assessing the overlap of the fibers with the probability maps.
All methods were applied in native space for direct comparison. Fiber tracking in standard space, provided by a registration method capable to deal with tumor data, was not applied because this study is intended to provide results of direct relevance for state-of-the-art neurosurgery.
After an evaluation of SP, TD, GibbsT and GS using DTI software phantoms, all methods were applied on two pathways of the language system, the AF and the IFOF. Results for the patient and control group were considered separately to investigate the influence of nearby tumors. Each pair of methods was evaluated in terms of the presented similarity metric S avg and the average FA value FA avg . To highlight the challenges that each of these reconstruction problems pose, both, the two pathways and the two groups, were considered separately. By applying non-parametric tests of Mann-Whitney [38] with a significance level of av0:05, statistical significance was tested for the defined groups.
Details about the parameters and adjustments to the connectivity problems for either method are presented in the following: Generally, for both the simulated and the image data, the step length of SP and TD was set to 0.5 mm. The search space of GS was defined by a cubic grid with 74 neighboring points and an average step length of GS to 2 mm, resulting in a smallest and largest step length of 1.3 and 2.3 mm, respectively. The threshold on the bending angle of GS was set to 75 0 . All methods sampled For the software phantoms, ROIs covering the full diameter of the fiber tracts were manually positioned at both ends of the simulated pathways. As the diffusion characteristics were different for each DTI phantom, FA thresholds were defined individually and set to 20% below the average FA value in the main part of the pathways for all methods. Restriction of the search space using a bounding box was not applied for the simulated data.
For the experiments based on image data of patients and healthy controls, ROIs corresponding to the speech areas obtained by fMRI were defined. The fibers created by SP and TD were separated into the dorsal and the ventral pathway by two manually defined ROIs placed at the middle of each tract.
Due to the aforementioned drawbacks of deterministic fiber tracking methods in regions of low FA, they may not succeed in locating a connection between these ROIs. Therefore, two modifications to the parameters and ROIs were applied: At first, tracking was performed with FA thresholds of 0.25 and on failure with 0.2. Secondly, ROIs were incrementally expanded by at most 4 voxels or 7.6 mm in each dimension until a connection was found with either FA threshold. In the following, this process is denoted as ROI padding. The decision rule about which fiber tract to select for the similarity metric was defined as follows: If fiber tracking with an FA threshold of 0.25 required no more than two additional ROI paddings as compared to fiber tracking with an FA threshold of 0.2, the fiber tracts with FA threshold of 0.25 were used. Otherwise, the fiber tracts with FA threshold of 0.2 contributed to the similarity value.
GibbsT was applied using the MITK implementation as presented in [51]. Given only six gradient diffusion data sets, the approximation of the orientation distribution function (ODF) with 252 directions was accomplished using spherical harmonics of order 4 according to [52]. As the ODFs are based on far less diffusion directions as usually required for high angular resolution diffusion imaging (at least 30), resolving complex fiber arrangements as in [20] might not be possible. Nevertheless, GibbsT based on the major diffusion direction is expected to produce fibers of at least the same quality as deterministic fiber tracking.
The fiber reconstruction itself was applied using the following parameters: 5:0|10 8 iterations, start and end temperatures of T~0:1 and T~0:001, 40 mm minimum fiber length, equally weighted internal and external energies, a curvature threshold of 45 0 and subtraction of the ODF mean. The particle parameters were defined as length~1:5 mm, width~0:5 mm and weight~0:0017. Fibers of GibbsT were imported into MedAlyVis (Medical Analysis and Visualization, Computer Graphics, University of Erlangen, Germany), the platform for generation of fiber tracts of SP, TD, and GS, and their visualization and comparison. Similarly to SP and TD, GibbsT applies fiber tracking to the whole brain and tends to generate fibers not reaching both ROIs. Therefore, ROI padding was applied as well.
Regarding ConProb, connection probability maps for both ROIs (without padding) were generated in the Matlab-based DTI & FiberTools [53] using N~10 4 randomly computed fibers with a maximum length of 150 voxels and an exponent of 4 on eigenvalues to sharpen diffusion tensors. The resulting two probability maps were combined to the ''probability index of forming part of the bundle of interest'' (PIBI) as described in [19], and imported into MedAlyVis for visual comparison.
In order to search for anatomical connections with GS, the bounding box was manually placed for both the AF and IFOF by considering anatomical features and activation centers. In case of the AF, the bounding box was placed to include the white matter of the left hemisphere without the rostral part of the temporal lobe. Regarding the IFOF, the upper boundary of the bounding box was placed above the frontal horn of the lateral ventricle, and the inner boundary directly medial to the external capsule. In all cases, the FA threshold was initially set to 0.3 and reduced to 0.15 in case of failure. In contrast to SP, TD and GibbsT, ROI padding was not applied for GS.
Prior to statistical evaluation, each pair of reconstructed fiber tracts was inspected visually to ensure that the expected anatomical space of the language site is covered. For this purpose, the neuromedical expert (coauthor A. Z.) investigated slices of the T1 MR images with an overlaid display of fiber tracts, fMRI activation clusters, and tumor segmentation (if applicable) as illustrated in Figure 2. Pairs of fiber tracts were excluded from the similarity evaluation if at least one of the fiber tracts did not cover main parts of the AF or the IFOF.

Results
The evaluation of the methods on ground truth data using DTI phantoms is presented in the next section. The results of the evaluations carried out on the controls with unaffected brain anatomy are explicated in the second section, on the brain tumor patients in the third section.

Evaluation on simulated data
For a qualitative evaluation on simulated data, success and failure to connect the ROIs was assessed in a first step. All methods successfully reconstructed branching and kissing pathways for images with both SNR 15 and 30, leading to similarity results. In case of the simulated spiral, only about two thirds of the first cycle was covered by SP and TD before the fibers terminated at the border of the pathway (Figure 3). Fiber segments generated by GibbsT covered the whole spiral, however connected fibers achieved not much more than one cycle. In contrast, GS successfully reconstructed a trajectory along the full extent of the spiral.
Based on the previously defined similarity metric S min (F ,G) and images with SNR 30, the deviation between fiber tracts of SP and TD amounted to 0:03 mm in the branching case and 0:04 mm in the kissing case. The deviations between SP and GS on the same phantoms were 2:68 mm and 1:38 mm, respectively. The alignment between SP and GS could be improved by increasing the angular resolution of the grid. Using a cubic grid with 154 neighboring points and an average step length of 2 mm, these similarity values amounted to 1:05 mm and 0:82 mm. As a result, the course of the closest fibers was almost identical for both methods. Regarding GibbsT compared with GS, deviations of 2:90 mm and 1:73 mm were observed in the branching and the kissing case, respectively.

Evaluation of control subjects
Reconstruction of pathways. On image data of control subjects, both SP and TD were able to reconstruct neuronal pathways, although relaxation of the FA threshold and incremental ROI padding had to be applied in some cases to obtain fibers. In one case, the FA threshold had to be reduced to 0. Visual inspection. After visual inspection, the fiber tracts generated by all methods along the dorsal pathway were confirmed to follow the AF anatomically in all cases of the control group. Regarding the match between the ventral fibers and the anatomy of the IFOF, GS was successful in all cases, while SP and TD established connections along the correct pathway in five and three cases, respectively. In the remaining cases, at least one of the methods reconstructed fibers corresponding to the ILF. Similar results were created by GibbsT for the ventral fibers. Here, fibers did not match the IFOF, but the ILF for all cases which is exemplified for the control patient in Figure 4. From visual comparison of the PIBI maps by ConProb, a strong agreement between ConProb and GS on the anatomy of the dorsal and the ventral pathway was found which is highlighted in Figure 4. The AF and the IFOF were found in all cases, whereas pathways corresponding to the ILF could be established in two cases only.  Similarity metrics. Results of pairwise comparison between all the methods excluding ConProb are summarized in the boxplots in Figure 5 with similarity metrics S avg averaged for corresponding fiber tracts. N represents the number of pairs of fiber tracts that cover the same anatomical space visually identified as either AF or IFOF.
For the control group, complete agreement of all methods was found for the AF, whereas less matching pairs between GS and SP and between GS and TD were found for the IFOF, and none for GS and GibbsT. Regarding the geometric similarities along the AF, fiber tracts of GS and SP were on average closer than of GS and TD and of GS and GibbsT. Along the IFOF, fiber tracts of GS and TD were closer than those of GS and SP. For GibbsT, no fiber was found for the IFOF.
To elucidate the effect of methods on FA, values of FA avg averaged across subjects are summarized in the barplots in Figure 6 for the AF and the IFOF reconstructed with different methods for both groups. Regarding the different pathways in the control group, larger average FA values were observed for the AF as compared to the IFOF (p~0:003). No clear differences could be observed for the methods.  GS reconstructed the dorsal phonological pathway in all but four patients with an FA threshold of 0.3. The remaining cases had either tumors of WHO grade III or IV, or the distance between pathways and the tumor was smaller than 5 mm. By reducing the FA threshold to 0.15, connections could be found in these cases as well. For the ventral semantic pathway, fiber tracts were found by GS in all but two patients with an FA threshold of 0.3. In the remaining cases, patients with a WHO grade IV tumor in the left temporal lobe and a WHO grade II tumor in the temporal stem, GS applied with an FA threshold of 0.15 found no fiber tracts.
Visual inspection. After visual inspection, the fiber tracts along the dorsal pathway were confirmed as anatomically feasible in the majority of all patients for SP, TD and GibbsT: 15, 17 and 14 out of 25 cases with ROI padding averages of 1.0, 0.81, and 1.64 respectively. The results of GS were confirmed to follow the AF for all subjects in the patient group.
In contrast to the positive results obtained for all methods dorsally, fiber tracts reconstructed by SP and TD along the ventral pathway correlated in some cases with the ILF rather than the IFOF. This problem becomes apparent in Figures 7,8,9 and 10 where either SP, TD or both fail to reconstruct the IFOF. In Figures 11 and 12, SP agreed with GS on the IFOF.
By replacing the frontal ROI with the rostral part of the external capsule, the focus of SP and TD could be changed to the IFOF. However, even after this modification, only some fiber tracts with anatomical correspondence to the IFOF could be reconstructed. Regarding the anatomic feasibility of the results of GS, all successfully reconstructed fiber tracts followed the trajectories of the IFOF.
For only one case, ventral fibers of GibbsT followed the IFOF (ROI padding of 3), and the ILF in eight cases (average ROI padding of 2.7) which is shown for the patients in Figures 4, 7, 10, and 11. Regarding the PIBI maps shown in Figures 4,5,6,7,8,9,10,11 and 12, the dorsal pathway was assigned to the AF in all cases except for two subjects both with a left glioblastoma multiforme (WHO grade IV) located between Broca and Wernicke areas. The ventral pathway corresponded to the IFOF in all but two cases, whereas the ILF was marked in all but nine cases.
Similarity metrics. Apparently, when observing the results in Figure 5, there is a higher correspondence and increased geometric similarity for fiber tracts along the AF than along the IFOF when GS is compared with either deterministic fiber tracking method. Across all fiber reconstruction methods, significantly lower values of S avg were observed for AF than for IFOF (p~0:005).
Regarding the comparison of diffusion values, higher average FA values were observed along the AF than along the IFOF in all compared fiber tract pairs except for the pair GS and GibbsT (compare barplots in Figure 6. Across all methods, a p-value of 0.005 was observed for the patient group, whereas significant differences belonged especially to the pairs GS-TD (p~0:031) and SP-TD (p~0:011).
In case of the AF, visual assessment revealed that the increased diameter of this pathway compared to the IFOF causes a stronger deviation between the methods despite their agreement on the anatomical space of the AF.

Group differences
Regarding differences of S avg between groups, no statistically significant effects could be established. In the controls, slightly lower values were observed in case of the AF for all compared methods, and in case of the IFOF for GS compared with TD. However, this trend is too weak to show statistical significance. In terms of group effects for diffusion values, significantly reduced average FA values were detected in the patient group in comparison to the controls (pv0:001). In addition, the average FA values were higher for GS than for SP and TD (p~0:058), whereas no preference between SP and TD (p~0:218) could be observed.

Computational aspects
All methods were executed on a Intel Pentium Core i7 (3.0 GHz, 12 GByte RAM) and timings were recorded for each subject including loading of data and ROIs, preprocessing such as computation of the diffusion tensors, execution of the methods and automatic filtering with ROIs. SP and TD required between 5 and 10 minutes, and GS took 20 minutes on average. The computation of the connection probability maps amounted to between 15 and 45 minutes. In contrast, the optimization of the Gibbs tracking approach took about two hours.

Discussion
In this study, five different methods for reconstructing neuronal pathways of the language system were compared: GS and GibbsT as representatives for global approaches, two clinically applied fiber tracking techniques, namely SP and TD, and ConProb, a connection-oriented probabilistic tractography approach.
For the simulated DTI data, GS was able to obtain valid fiber tracts for all configurations, whereas SP, TD, and GibbsT failed to determine connecting fiber tracts in case of the spiral. This clearly indicates the methodological advantages of GS for connectivity problems with known start and end regions. The relatively large spatial deviation between GS and SP or TD can be explained by the grid dependency of GS that becomes obvious as the alignment is improved with grids of higher angular resolution.
Regarding the pathways in the image data, smaller numbers of corresponding pairs of fiber tracts and generally higher distance values were observed for the IFOF as compared to the AF. In addition, average FA values along the IFOF were significantly lower in comparison to the AF. This explains the weaker agreement between the considered methods on the location of the IFOF, and implies that this pathway is more difficult to reconstruct. An additional explanation for these difficulties might be a lower diameter of this pathway and the associated risk to lose track more easily.  In general, a stronger agreement between methods was observed in controls compared to the brain tumor patients as shown by fewer excluded fiber tract pairs. These results clearly indicate the influence of tumors on white matter tracts, which is also reflected by smaller FA values in the presence of nearby tumors. This is a possible reason for the more complicated reconstruction of feasible fibers in brain tumor patients, which is also in line with the detailed analysis of the influence of tumors on FA values in [10]: Accordingly, in white matter infiltrated by tumors or affected by the peritumoral edema, decreased FA values can be observed, whereas displacement and compression of pathways may lead to increasing FA values along pathways. In general, neurological diseases are accompanied by a reduction of FA values. Additionally, the slightly higher mean age of the patient group compared to the control group may lead to decreased FA values [46].
With respect to methods, the application of probabilistic fiber tracking [19] to both groups provided additional information about the structural connectivity problem of the language system. While the overlap of the covered anatomical space by ConProb and GS was generally high for both AF and IFOF, connections along the ILF found in the majority by SP and TD were rated as less likely by ConProb.
Regarding the strong concept of Gibbs tracking, the results for the dorsal pathway strongly correlated with the results of the other methods. This is in line with a successful application of this method on different tract systems of the brain [22]. The reason for only one successful reconstruction via the IFOF and only some reconstructions via the ILF could be the limited number of gradient directions. For a higher number, the method has shown excellent performance previously [22].
Regarding the global search, SP and TD, slightly better average FA values were observed for GS than for SP or TD, which might be due to the inclusion of FA into the cost function. As a result, the search space is preferably expanded to regions with higher anisotropic diffusion, which may result in fibers located closer towards the center of pathways.
While GS was able to reconstruct the AF and in particular the IFOF with initial settings in most cases, both SP and TD required a reduction of the FA threshold and padding of the ROI volume in order to establish a connection at all. Moreover, visual inspection of fiber tracts resulting from reconstruction with SP and TD  showed a tendency to erroneously reconstruct fibers along the ILF rather than along the IFOF. The difficulties of both deterministic approaches to reconstruct the AF, and the IFOF in particular, can be attributed to the design of the algorithms, which only take the local tensor information into account for fiber reconstruction. Especially for pathways with multiple terminations such as the IFOF [54], this local interpretation of the tensor data is not sufficient, which makes deterministic approaches susceptible to tracking errors.
On the other hand, global approaches like GS make use of the entire tensor information and establish connections which are optimal in respect of the available diffusion data between two ROIs. The underlying concept of cost minimization is justified since it is anticipated that functional areas are connected by preferably short and direct pathways as a result of the morphological evolution of white matter pathways in the brain [55]. Additionally, since multiple potential fiber directions are considered simultaneously, GS appropriately handles uncertainty in case of splitting and merging fibers. In this sense, GS [25] is similar to global approaches like probabilistic fiber tracking [19] or the shortest path approach using a Bayesian framework [24], which all implicitly model the uncertainty of DTI data.
However, since GS avoids time-consuming computations, it can be considered as a compromise between these approaches and deterministic fiber tracking. As anatomical connectivity can be resolved with sufficient accuracy and results are provided in a fast way, GS is well suited for clinical application. In particular, in combination with DTI protocols that acquire smaller numbers of gradient images and with clearly defined start and end ROIs located by fMRI [10], GS can provide information about the location of neuronal pathways for both pre-operative planning and intra-operative neurosurgery. Moreover, since the approach itself is independent of the underlying diffusion model, its extension to more complex diffusion models based on high angular resolution diffusion imaging such as q-ball imaging [56] or higher order tensors [57] is straight forward.
Despite of its considerable advantages, one limitation of the GS approach remains. In case of functional areas that are connected via multiple alternative pathways, GS requires anatomical knowledge to appropriately define the search space with the bounding box. However, as this mechanism can be used to draw the focus to either pathway, GS is the superioralternative to investigate anatomical connectivity, since SP and TD are not able to take advantage of such an approach.

Conclusion
In this study, two deterministic approaches to fiber tracking, streamline propagation (SP) and tensor deflection (TD), a dedicated approach for anatomical connectivity analysis (GS), Gibbs tracking (GibbsT) and a connection-oriented probabilistic tractography approach (ConProb) were compared with respect to their capability to reconstruct neuronal pathways of the language system. Based on DTI data of a patient group with supratentorial tumors close to pathways and functional areas, and of a control group, two pathways of the semantic system were considered. In addition to visual validation based on anatomical knowledge, the evaluation of diffusion values and distances between fiber tracts showed that GS produces superior results as compared to SP, TD, and GibbsT. In particular, the ventral language connections via the IFOF were revealed by GS in most cases which was in line with the probability maps of ConProb. In contrast, SP, TD, and GibbsT often failed with standard settings for ROIs and thresholds on fractional anisotropy (FA). Despite relaxation of FA thresholds and expansion of ROIs, SP, TD, and GibbsT still failed in some cases. Overall, GS proved to be more reliable with less failures and less adjustments required.