- Email: [email protected]
Intersubject variability in the analysis of diffusion tensor images at the group level: fractional anisotropy mapping and fiber tracking techniques
Available online at www.sciencedirect.com
Magnetic Resonance Imaging 27 (2009) 324 – 334
Intersubject variability in the analysis of diffusion tensor images at the group level: fractional anisotropy mapping and fiber tracking techniques Hans-Peter Müller⁎, Alexander Unrath, Axel Riecker, Elmar H. Pinkhardt, Albert C. Ludolph, Jan Kassubek Department of Neurology, University of Ulm, D – 89081 Ulm, Germany Received 29 November 2007; revised 22 June 2008; accepted 1 July 2008
Abstract Introduction: Diffusion tensor imaging (DTI) provides comprehensive information about quantitative diffusion and connectivity in the human brain. Transformation into stereotactic standard space is a prerequisite for group studies and requires thorough data processing to preserve directional inter-dependencies. The objective of the present study was to optimize technical approaches for this preservation of quantitative and directional information during spatial normalization in data analyses at the group level. Methods: Different averaging methods for mean diffusion-weighted images containing DTI information were compared, i.e., region of interest-based fractional anisotropy (FA) mapping, fiber tracking (FT) and corresponding tractwise FA statistics (TFAS). The novel technique of intersubject FT that takes into account directional information of single data sets during the FT process was compared to standard FT techniques. Application of the methods was shown in the comparison of normal subjects and subjects with defined white matter pathology (alterations of the corpus callosum). Results: Fiber tracking was applied to averaged data sets and showed similar results compared with FT on single subject data. The application of TFAS to averaged data showed averaged FA values around 0.4 for normal controls. The values were in the range of the standard deviation for averaged FA values for TFAS applied to single subject data. These results were independent of the applied averaging technique. A significant reduction of the averaged FA values was found in comparison to TFAS applied to data from subjects with defined white matter pathology (FA around 0.2). Conclusion: The applicability of FT techniques in the analysis of different subjects at the group level was demonstrated. Group comparisons as well as FT on group averaged data were shown to be feasible. The objective of this work was to identify the most appropriate method for intersubject averaging and group comparison which incorporates intersubject variability of the directional information. © 2009 Elsevier Inc. All rights reserved. Keywords: Magnetic resonance imaging; Diffusion tensor imaging; Fractional anisotropy mapping; Fiber tracking
1. Introduction Diffusion in human brain white matter (WM) can be noninvasively mapped by diffusion tensor magnetic resonance imaging (DTI). By DTI, the directional dependence of diffusion in each voxel can be characterized by the so-called ⇉ Y Y Y diffusion tensor D. The eigenvectors ( v 1; ; v 2 ; v 3 ) and ⇉ eigenvalues (λ1, λ2, λ3) of the 3×3 matrix D reflect the diffusivity of water in each direction. This can be used to ⁎ Corresponding author. Tel.: +49 731 177 1206; fax: +49 731 177 1202. E-mail addresses: [email protected] (H.-P. Müller), [email protected] (J. Kassubek). 0730–725X/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2008.07.003
quantify the diffusivity by so-called fractional anisotropy (FA) maps on a voxelwise basis [1–3]. The directional information can be used as the basis for reconstruction of the interconnectivity of brain regions by following the fiber pathways by fiber tracking (FT) techniques. The basis of FT is the consecutive connection of neighbored tensors along their principal directions. Basically, the FT techniques can be divided into two groups: streamline or deterministic FT [4–7] and probabilistic FT [8–14]. The objective of the present study was to optimize the technical approach for the preservation of quantitative and directional information contained in diffusion weighted images/DTI data in order to provide a framework for the comparison of FT of subject groups with defined anatomical
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
WM alterations and groups of normal subjects. To find a difference between two subject groups may eventually assist the diagnosis in finding ways to categorize a subject to one of the group based on some discrimination metric. The prerequisite for this comparison of subject groups was the transformation of the basic DTI data into a stereotactic standard space, e.g., the Montreal Neurological Institute (MNI) space [15]. After this normalization, a comparison at the group level might be performed in several ways: (i) FA maps of different subject groups could be compared, while the directional information within these maps is not considered. (ii) FT results of group-averaged data sets could be qualitatively compared by neuroimaging experts on the basis of standardized human brain atlases [16,17]. The result images were judged in terms of the three factors directionality, cortical projection areas and homogeneity of fiber bundles. (iii) Skeletons that consist of fiber tracts could be used for quantitative comparison on the basis of the directional dependence. In order to show the validity and power of the different techniques, the following questions were to be addressed: • How could intersubject averaging and group comparison be performed and which were the most appropriate methods to quantify differences between subject groups? • How to consider intersubject variability of the directional information? Like in other advanced magnetic resonance imaging (MRI) methods, DTI- and FT-based studies pursue the ultimate goal to categorize individual patients' brain morphology in order to facilitate the diagnostic process based on some discrimination metric. In the present study, the standardization of technical approaches to group-based analyses of certain brain pathologies can be considered the prerequisite for the identification of pathological patterns of altered brain anatomy. As an example for a groupwise comparison, differences between patients with atrophy of the corpus callosum (CC) and age-matched healthy controls were mapped. As a model of CC alteration, subjects with the neurodegenerative disease of complicated hereditary spastic paraparesis were investigated as a prototype of morphological alterations (thinning) along the whole structure of the CC [18,19]. The CC was chosen as the most appropriate structure in the brain to be analyzed since it is one of the white matter structures with the most strongly directed fibers [20]. 2. Data recording and data preprocessing 2.1. Data recording and standard data preprocessing All DTI data were acquired on the same 1.5T scanner (Symphony, Siemens Medical, Erlangen, Germany). Six
325
healthy controls (three men, three women, average age 32.7±4.5 years), and six patients with thinned CC (tCC) (three men, three women, 32.5±12.1 years old) underwent the MRI protocol. All DTI acquisitions consisted of 13 volumes (45 slices, 128×128 voxels, slice thickness 2.2 mm, in-plane voxel size 1.5×1.5 mm), representing 12 gradient directions and one scan with gradient 0 (b0). Echo time (TE) and repetition time (TR) were 93 ms and 8000 ms, respectively. b was 800 s/mm2, five scans were k-space-averaged online by the Siemens SYNGO operating software. As an anatomical background, a high-resolution T1-weighted, magnetization-prepared, rapid-acquisition gradient echo (MPRAGE) sequence was used (TR=9.7 ms, TE=3.93 ms, flip angle 15°, matrix size 256×256 mm2, voxel size 1.0×0.96×0.96 mm 3 ), consisting of 160–200 sagittal partitions depending on the head size. All analyses were performed by the software package (Tensor Imaging and Fiber Tracking) [21]. Standard image processing procedures such as eddy current correction and transformation to iso-voxels and smoothing have already been described previously [21,22]. Also, the spatial normalization protocol to MNI standard space has also been described in [21]. All data sets were preprocessed and normalized accordingly. 2.2. MNI normalization and storage of major eigenvectors Basically, a complete nonlinear MNI normalization consisted of three deformation components (DC): • DC 1: a rigid brain transformation to align the basic coordinate frames. The rotation angles had to be stored ⇉ in a rotation matrix R. • DC 2: an affine deformation according to landmarks. The six affine deformation parameters for the different brain regions had to be stored in an array with 6 Y components S. • DC 3: a non-affine normalization equalizing non-linear brain shape differences. The 3D vector shifts were different for each voxel leading to a separate transformation for each voxel of the 3D voxel array Y Y Y (a 5D transformation array T with indices image row, image column, image slice, dilation matrix row, dilation matrix column) (details see below, Eq. 3). Consequently, for normalized data the resulting diffusion ⇉
tensor D i of each voxel i had to be rotated according to all the rotations listed above. • A rotation resulting from the aligning to the basic coordinate frame (corresponding to DC 1) had to be applied ⇉
⇉
⇉
D i V¼ R D i
ð1Þ Y
Y
Y
• The components of the eigenvectors (V 1 ;V 2 ;V 3 ) had to be adapted according to the six affine
326
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334 Y
normalization parameters of S (dependent on the brain region sa, a=1...6) of the affine deformation (corresponding to DC 2).
and the standard deviation of averaged FA
vw;j j ¼ sa vw;j j
Fstd
ð2Þ
with w=1...3 and j=x,y,z. After adapting, the eigenvectors could have lost their unity length and had to be normalized. • Standard trigonometry gives a rotation matrix (for each voxel independently), resulting from the 3D vector shifts following the concepts of Alexander et al. [23] in order to preserve the directional relations between eigenvectors of neighbored voxels, different shifts of two neighbored voxels result in rotations of the corresponding eigenvectors. The dilation⇉matrices were used for the alignment of the tensor D of each voxel to the surrounding voxels (corresponding to DC 3). ⇉
⇉
⇉
DiW ¼ t i Di V ⇉
ð3Þ Y Y Y
where t i are the components of T . These⇉fine-corrections (according to Eqs.2 and 3) of the tensors D i were essential for a correct FT [22], and the corresponding parameters had to be stored for each subject data set separately. With the standard normalization procedure, the corresponding MPRAGE data were also normalized to MNI space. From the average of these data sets, a display background was created. 3. Analysis techniques 3.1. Region of interest analysis of FA maps The diffusion anisotropy could be quantified by FA: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u3 k1 k 2 þ k 2 k 2 þ k 3 k 2 F¼t 2 k21 þ k22 þ k23
ð4Þ
1 XN Y Y F i;j f s v i jbrs i N
ð6Þ
where N was the ROI size. Extensive studies via ROI analyses about the preservation of the DTI-specific parameters during the process of MNI normalization had been performed previously [21]. The FA maps were calculated separately for each subject data in MNI space, making statistical comparison of the subject groups feasible. Whole brain-based statistical analysis (WBSA) was performed by t test, i.e., voxelwise statistical comparison of normal subject FA maps and FA maps from subjects with tCC. 3.2. Fiber tracking methods For the calculation of the diffusion spheroid, the eigenvector corresponding to the major eigenvalue indicated the fiber direction in WM regions. Based on this directional information, different methods and algorithms had been proposed to estimate WM connectivity. In the present study, the deterministic streamline tracking (DST) technique [24,25] and a technique based on probabilistic fiber tracking (PFT) [8,9], which models the propagation in the major eigenvector field of the brain were applied. Generally, the FT positions resulted from floating point numbers. The corresponding eigenvector direction (pointing towards the consecutive FT position) was the interpolation of the directions of the neighboring voxels weighted by the proportionate position (linear nearest neighbor interpolation) Y
v new ði; j; k Þ ¼
X8 w¼i
Y
aw v ðlw ; mw ; nw Þ:
ð7Þ
Y
where k was the arithmetic average of all eigenvalues. The intensity was related to the FA value and the color-coding was as follows: red for major eigenvector mainly in left-right direction, blue for major eigenvector mainly in inferiorsuperior direction, and green for major eigenvector mainly in posterior-anterior direction. A region of interest (ROI)-based approach of defined brain regions was the method of choice for the analysis of separate areas. An ROI was set by a sphere with an observer-defined radius rs and its center at the user-defined Y focus f s . All FA values Fi of the N voxels inside this sphere Y (position v i ) were included in the following parameterization: the average FA, i.e., average diffusion strength Favg ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 XN Y Y ¼ F F avg i i;j f s v i jbrs N 1
ð5Þ
Here, v new (i,j,k) was the resulting vector at the new position i,j,k (floating point numbers), and l,m,n were the voxel coordinates (integer numbers) of the eight neighboring voxels. The factors aw were the respective eight weighting factors for the interpolation. The threshold for the scalar product of the major eigenvectors (proportional to the angle between the major eigenvectors of two consecutive FT positions) was set to 0.9. The distance between two FT positions, i.e., the stepwidth, was set to 0.5 mm, corresponding to 0.5 voxels (this value of 0.5 mm had also been suggested by [4]). 3.2.1. Modified DST technique The DST approach (Fig. 1A) was modified in the following way: Together with the selected FT starting point, additional 26 fiber tracts started at the 26 directly neighbored voxels i. At each FT step, the mean and
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
327
Fig. 1. Methods of FT for single subject data. (A) DST; starting points were located in a ROI of the radius 2 voxels (2 mm). (B) MDST; compared to (A) runaway fiber tracts were avoided. (C) MPFT. (D) MPFT with display of the end points of fiber tracts — the color bar encodes the probability that a voxel is crossed by a fiber tract (P) in (C) and (D). The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the FA maps. Y
the standard deviation of these 27 vectors v new;j were calculated. 1 X27 ðk Þ kÞ vðavg ¼ v ð8Þ i¼1 i 27 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 X27 ðk Þ ðk Þ ðk Þ 2 vstd ¼ vi vavg ð9Þ i¼1 26 Y
where k indicates the components of the vector v i. Vectors which deviated by more than the standard deviation from the average vector were replaced by the average vector in order to avoid runaway fiber tracts (Fig. 1B). 3.2.2. Modified PFT technique PFT was modified in the following way: a volume with a radius of three voxels around each FT starting point was defined. Then, based on Monte Carlo simulations, additional 10,000 randomly chosen starting points (Gaussian-like distribution) were created. Runaway fiber tracts were avoided by the same technique as described above (Section 3.2.1, Eqs. 8 and 9). A probability map could be set up, which was weighted by the number of fiber tracts intersecting each voxel (Fig. 1C). If only connected brain regions were of interest, rather than the discrete fiber tracts, then the only regions displayed were those in which fiber tracts ended (Fig. 1D). 3.2.3. Intersubject fiber tracking Based upon the directional information for all subject data, a new FT technique could be developed, called “intersubject
FT” (ISFT). For each of the M group subject data sets (M=6 in this study), the identical starting point was defined. The average and the standard deviation of the M new vectors were calculated at each FT step. Vectors deviating more than the standard deviation from the average vector were replaced by the average vector. Note that the new FT position could be different for each subject data. In this way, runaway fiber tracts were avoided and an FT was set up that took the inter-subject variability of the eigenvectors at each FT step into account. 3.3. Tractwise fractional anisotropy statistics The fiber tracts could be used as a skeleton mask for selective statistics [9,26–30]. For that purpose, tractwise fractional anisotropy statistics (TFAS) was used [22]. Each voxel that was crossed by a fiber tract was defined as “active” for statistics and thus contributed to the skeleton mask so that the underlying FA values of the voxels contributed to the statistical t test. The average and standard deviation were calculated by 1 XN TFAS ¼ F ; iaskeleton ð10Þ Favg i¼1 i N TFAS Fstd
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 XN ¼ Favg Fi ; iaskeleton i¼1 N 1
ð11Þ
Separate skeletons were used for volunteers and for tCC subjects, i.e., FT was performed both on group averaged normal subject data and on group averaged tCC data. In
328
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
this way, two skeletons were derived, and the comparison was performed by applying skeleton 1 to group averaged healthy subject FA maps and skeleton 2 to group averaged tCC FA maps. 3.4. Averaging of DTI data 3.4.1. Storage of fine corrections Each DTI data was normalized, and all DTI data were averaged before FA mapping. The fine corrections were also arithmetically averaged (Section 2.2). This approach resulted in an averaged fine correction to be applied to the resulting averaged DTI data.
3.4.2. Storage of fine-corrected eigenvectors Major eigenvectors were stored after fine-correction, and a new major eigenvector array was created for the averaged DTI data. For this purpose, a representative major eigenvector was selected by the criterion XM Y Y v d vj ð12Þ Si ¼ j¼1; j p i i The major eigenvector which had the maximum Si was selected at a respective voxel position out of all major eigenvectors of the M different subjects. This new major eigenvector replaced the major eigenvector calculated from
Fig. 2. WBSA on a voxelwise basis. Comparison of FA maps of six healthy subjects (A) and six subjects with tCC (B). The histograms indicate the distribution of the absolute FA values within the ROI. (C) t Test analysis, display threshold Pb.05. The results are overlaid on an MPRAGE template (averaged from all 12 subjects). The ROI used for analysis (Table 1) is indicated in light blue.
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
the averaged DTI data at the corresponding voxel position. This method required increased memory space and approximately doubled computation time.
For the group comparison, averaging of results for different subjects was necessary. After all subjects' data had been transformed into MNI space, averaging was feasible. 4.1. FA map averaging without directional information WBSA was applied to FA maps of the 12 subjects contributing to the study. Fig. 2 shows the results for WBSA at Pb.05. For tCC, major differences could be detected in the CC. Therefore, a ROI with a radius of 45 mm (45 voxels) was defined including the CC, i.e., the ROI sphere included approximately 380,000 voxels. Average FA values were 0.42±0.17 for normal controls and 0.34±0.14 for tCC patients, respectively. The reduction of FA values was also observed in the histograms of the FA values within the ROI (Fig. 2A and B). The main contribution resulted from the left-right direction (red) (Table 1), compatible with the normal fiber distribution in the CC. By use of t test, approximately 51,000 voxels within the ROI were observed to pass the threshold of Pb.05 (Fig. 2C). 4.2. Single subject fiber tracking The modified PFT (MPFT) results for single subject data using three starting points in the CC according to Table 2 are displayed in Fig. 3. TFAS analysis showed averaged FA values between 0.33 and 0.43 for the different subjects (Table 3). 4.3. Averaging of DTI data from different subjects Averaged DTI data were created by the storage of fine corrections method (Section 3.4.1). The corresponding FT (starting points, cf. Table 2) is displayed in Fig. 4A and B for modified DST (MDST) and MPFT, respectively. Visual inspection showed minor differences. Differences mainly occurred in juxtacortical regions with low FA values. The method of storage of fine-corrected eigenvectors (Section
Table 1 ROI analysis and comparison for subject groups (see also Fig. 2) ROI position MNI: 0–11 mm, 18 mm; ROI radius: 45 mm, ROI size: 380,000 voxels Healthy controls
tCC patients
Favg±Fstd
Total Left-right Superior-inferior Anterior-posterior
Table 2 MNI coordinates for starting points for the fiber tracking (FT)
Starting point 1 Starting point 2 Starting point 3
4. Results
0.42±0.17 0.35±0.14 0.37±0.14 0.56±0.16
t test 51,000 voxels (13%) passing threshold of Pb.05
0.34±0.14 0.26±0.08 0.31±0.09 0.51±0.11
329
x/mm
y/mm
z/mm
0 0 0
−40 −10 21
17 29 8
3.4.2) showed the FT (starting points, cf. Table 2) displayed in Fig. 4C and D for MDST and MPFT, respectively. MDST and MPFT showed corresponding results independent of the averaging method in regions near the CC (identification on the basis of standardized human brain atlases [16,17]). Here, only descriptive comparison is feasible. The comparison was performed by four neuroimaging experts (J.K., A.U., A.R., E.P.) on the basis of standardized human brain atlases [16,17]. The resulting images were judged in terms of the three factors directionality, corresponding cortical projection areas, and homogeneity of fiber bundles (three-item scale: good–moderate–bad concordance of the results). With respect to directionality and anatomical localization of the cortical projection areas, the concordance was always judged as “good.” The homogeneity of the fiber bundles was judged as “good” in the white matter regions adjacent to the connection line between anterior commissure and posterior commissure, whereas it was judged as “moderate” in juxtacortical regions, i.e., differences between the two averaging methods mainly occurred in juxtacortical regions where the FA values were b0.2. 4.4. Intersubject fiber tracking ISFT was based on the results of DTI analysis of single subject data. By taking the intersubject variability of the eigenvectors at each FT step into account, the result displayed in Fig. 4E was achieved. 4.5. Tractwise fractional anisotropy statistics The qualitative comparison by visual analysis (see Section 4.3) was extended by a quantification of the differences of the FT results as performed by TFAS. The fiber tracts obtained by the different FT techniques could be used as skeletons for TFAS, which will be described in the following. The averaged FA values for skeletons obtained by MPFT (starting points in the CC) in each subject (Eqs. 10 and 11) are summarized in Table 3. Alternatively, skeletons based on FT on group averaged data were used for statistical comparison between the normal subject group and the tCC subject group (Fig. 5). For the healthy control group (fine-correction method, Section 4.3), the average FA values were 0.35 for MPFT, 0.39 for MDST and 0.38 for ISFT, respectively. These values did not differ significantly from the single subject results (Table 3) and were significantly different from averaged FA values for the tCC patient group (fine-correction method). The results are summarized in Table 4.
330
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
Fig. 3. Fiber tracking results for different normal subjects by MPFT. FA threshold was 0.1, display background was MPRAGE template (averaged from all 12 subjects who participated in the study). The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set.
4.6. Interconnectivity of brain regions
5. Discussion
Fiber tracking could be used to investigate interconnectivity of brain regions (Fig. 6). The end points of all fiber tracts were displayed. An averaged fiber tract (from all fiber tracts contributing to the end point distribution) was created to exhibit the interconnectivity. The course of these averaged fiber tracts as well as the regions identified to be connected in frontal, parietal and occipital brain areas correlate well with neuroanatomical knowledge of commissural fibers. Here, the averaged fiber tract connecting the frontal brain areas of both hemispheres correlated with parts of the minor forceps (frontal part of CC), while interconnectivity within occipital brain areas corresponded to the major forceps (occipital part of CC).
MNI normalization and comparison of FA maps at the group level allowed for several possibilities to quantify differences between subject groups. 5.1. FA map analysis WBSA could be applied to normalized FA maps of controls and subjects with tCC to obtain a 3D comparison visualized using maps with color coded significant regions. The main contribution to differences from a discrete region could be analyzed via histograms. Nevertheless, the directional dependence was not taken into account by these analysis methods. 5.2. FT analysis
Table 3 TFAS for normal subjects (MPFT) TFAS FTFAS avg ±Fstd
Subject 1 Subject 2 Subject 3 Subject 4 Subject 5 Subject 6
0.43±0.20 0.37±0.23 0.39±0.19 0.33±0.22 0.41±0.21 0.35±0.23
In order to include directional information, averaging techniques for directional information were applied. The results of visual comparison between single-subject FT and the FT on group averaged data showed for the two averaging techniques (Sections 3.4.1 and 3.4.2) differences mainly in regions with low FA values (FAb0.2; this threshold had also been suggested by [31]). Thus, both averaging methods were considered appropriate. Never-
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
331
Fig. 4. Different analyses of FT in the averaged DTI data from 6 normal controls, depending if fine-corrections were averaged (A, B) or storage of major eigenvector was used (C, D). (A) MDST-averaged fine-corrections. (B) MPFT-averaged fine-corrections. (C) MDST storage of major eigenvector. (D) MPFT storage of major eigenvector. (E) ISFT. The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set.
theless, the technique based on fine-corrections was considered to be superior because of its advantages in memory usage and computation time.
The application of FT to tCC subject data caused problems due to low FA values in the CC. Variations in FT between the healthy group and the tCC group were mainly
Fig. 5. ISFT: (A) calculated from averaged DTI data of 6 normal controls. (B) Calculated from averaged DTI data of 6 subjects with tCC. The FT is color coded according to FA values. The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set.
332
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
Table 4 TFAS with different skeletons, depending on the FT mode Healthy subjects
tCC patients
TFAS FTFAS avg ±Fstd
MDST MPFT ISFT
0.39±0.18 0.35±0.16 0.38±0.18
t test P
0.19±0.08 0.18±0.07 0.19±0.10
b.001 b.001 b.001
not related to changes in the diffusion direction; moreover, the low values in mean diffusion weighted images (as detected by low FA values) led to errors in the diffusion tensors, and thus, information on the directional dependency was lost. 5.3. TFAS There have been various efforts to establish diffusion anisotropy as a marker for white matter tract integrity [26–28] by use of the underlying FA maps for selective statistics. Smith et al. [29,30] developed an algorithm for an
alignment-invariant tract representation to overcome normalization problems; this approach was referred to as tract-based spatial statistics [9,29,30]. In the present study, an alternative approach (TFAS) was applied, as previously described [22]. Here, bundles of fiber tracts were used in the sense of a skeleton which represented the basis of statistical analysis of the underlying FA maps. The novel character of TFAS was its use of averaged DTI data, i.e., the processing steps of normalization and averaging were not performed on FA maps but on one newly created group averaged DTI data (for a detailed description of the TFAS method, see Ref. [22]). 5.4. Interconnectivity of brain regions Interconnectivity of brain regions has repeatedly been a topic of research with various MRI techniques such as functional MRI [32] and other DTI-based approaches. In the present study, the investigation of interconnectivity was restricted to normal subjects due to directional loss in eigenvector calculation from DTI data of patients with tCC, e.g., the FA signal in the CC was sometimes low (FAb0.2)
Fig. 6. Interconnectivity display of brain regions — ISFT with starting points in the CC. Only the end points of the FT process are displayed. Coronal (A), sagittal (B), axial (C) and 3D view (D). The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set. Interconnectivity (derived from an averaged FT) is indicated as yellow lines following the fiber tracks. The color bar encodes the probability that a voxel is an end point of a fiber tract (P).
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
for the tCC data sets and, thus, sometimes led to an erroneous FT. Kunimatsu et al. [31] also suggested an FA value of 0.2 as the optimal trackability threshold. Generally, this technique has to be applied very thoroughly because low FA values (as usually occur in juxtacortical positions) might lead to incorrect FT and, thus, to an incorrect localization of fiber tract end points. 6. Conclusions Intersubject averaging of DTI data could be performed with respect to diffusion amplitude (by FA information) and diffusion direction (by FT). Averaging of FA maps allowed for the statistical comparison of subject groups with the WBSA method. Comparison of methods for maintaining directional information during FT data processing was performed. This work gives an instruction how intersubject averaging and group comparison could be performed and how intersubject variability of the directional information could be considered. In the following, the two questions from the introduction are to be addressed in detail. (i) How could intersubject averaging and group comparison be performed and which were the most appropriate methods to quantify differences between subject groups? The averaging methods showed similar results; therefore, systematic errors during the spatial normalization process could be excluded. Differences between the method of averaged fine-corrections and the method of stored major eigenvectors occurred in juxtacortical brain positions where the FA values were known to be b0.2 so that correct FT was impossible. The major eigenvector method has the drawback of huge memory consumption, i.e., in case of large patient number studies. (ii) How to consider intersubject variability of the directional information? ISFT is a new technique which takes into account inter-subject variability of directional information during FT. All FT techniques used (i.e., MDST, MPFT, ISFT) showed similar results for TFAS analysis. The results were in agreement with those of the TFAS analysis in single subjects. Therefore, DTI data averaging and FT calculated from these averaged DTI data was considered to be applicable. In summary, the averaging methods, FT techniques and FA quantification approaches (ROI analyses as well as TFAS) were demonstrated to be applicable to the comparison between different subjects. In this study, the use of the technical approach in alterations of the CC was specifically addressed and future studies will be necessary to demonstrate the applicability to other WM brain alterations. References [1] Basser PJ, Mattiello J, LeBihan D. Estimation of the effective selfdiffusion tensor from the NMR spin echo. J Magn Reson B 1994;103:247–54. [2] Basser PJ, Mattiello J, LeBihan D. MR Diffusion Tensor Spectroscopy and Imaging. Biophys J 1994;66:259–67.
333
[3] Pierpaoli C, Jezzard P, Basser PJ, Barnett A, DiChiro G. Diffusion tensor MR imaging of the human brain. Radiology 1996;201:637–48. [4] Conturo, Lori NF, Cull TS, Akbudak E, Snyder AZ, Shimony JS, et al. Tracking neuronal fibre pathways in the living human brain. Proc Nat Acad Sci U S A 1999;96:10422–7. [5] Lori NF, Akbudak E, Shimony JS, Cull TS, Snyder AZ, Guillory RK, et al. Diffusion tensor fibre tracking of human brain connectivity: acquisition methods, reliability analysis and biological results. NMR Biomed 2002;15:494–515. [6] Mori S, Kaufmann WE, Davatzikos C, Stieltjes B, Amodei L, Fredericksen K, et al. Imaging cortical association tracts in the human brain using diffusion-tensor-based axonal tracking. Magn Reson Med 2002;47:215–23. [7] Lazar M, Weinstein DM, Tsuruda JS, Hasan KM, Arfanakis K, Meyerand ME, et al. White matter tractography using diffusion tensor deflection. Hum Brain Mapp 2003;18:306–21. [8] Behrens TE, Johansen-Berg HJ, Woolrich MW, Smith SM, WheelerKingshott CA, Boulby PA, et al. Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging. Nat Neurosci 2003;6:750–7. [9] Behrens TE, Berg HJ, Jbabdi S, Rushworth MF, Woolrich MW. Probabilistic diffusion tractography with multiple fibre orientations: what can we gain? Neuroimage 2007;34:144–55. [10] Hagmann P, Thiran JP, Jonasson L, Vandergheynst P, Clarke S, Maeder P, et al. DTI mapping of human brain connectivity: statistical fibre tracking and virtual dissection. Neuroimage 2003;19:545–54. [11] Parker GJ, Haroon HA, Wheeler-Kingshott CA. A framework for a streamline-based probabilistic index of connectivity (PICo) using a structural interpretation of MRI diffusion measurements. J Magn Reson Imaging 2003;18:242–54. [12] Parker GJ, Alexander. Probabilistic anatomical connectivity derived from the microscopic persistent angular structure of cerebral tissue. Philos Trans R Soc Lond B Biol Sci 2005;360:893–902. [13] Jones DK, Pierpaoli C. Confidence mapping in diffusion tensor magnetic resonance imaging tractography using a bootstrap approach. Magn Reson Med 2005;53:1143–9. [14] Mori S, van Zijl PC. Fibre tracking: principles and strategies — a technical review. NMR Biomed 2002;15:468–80. [15] Brett M, Johnsrude IS, Owen AM. The problem of functional localization in the human brain. Nat Rev Neurosci 2002;3:243–9. [16] Talairach J. Tournoux. Coplanar stereotaxic atlas of the human brain. New York: Thieme Medical; 1988. [17] Mori S, van Zijl PCM, Wakana S, Nagae-Poetscher LM. MRI Atlas of Human White Matter. Amsterdam: Elsevier; 2005. [18] Kassubek J, Juengling FD, Baumgartner A, Ludolph AC, Sperfeld AD. Different regional brain volume loss in pure and complicated hereditary spastic paraparesis: a voxel-based morphometry study. Amyotrophic Lateral Sclerosis 2007;8:328–36. [19] Fink JK. Hereditary spastic paraplegia. Curr Neurol Neurosci Rep 2006;6:65–76. [20] Sullivan EV, Adalsteinsson E, Pfefferbaum A. Selective age-related degradation of anterior callosal fibre bundles quantified in vivo with fibre tracking. Cerebral Cortex 2006;16:1030–9. [21] Müller HP, Unrath A, Ludolph AC, Kassubek J. Preservation of diffusion tensor properties during spatial normalization by use of tensor imaging and fibre tracking on a normal brain database. Phys Med Biol 2007;52:N99–N109. [22] Müller HP, Unrath A, Sperfeld AD, Ludolph AC, Riecker A, Kassubek J. Diffusion tensor imaging and tractwise fractional anisotropy statistics: quantitative analysis in white matter pathology. Biomed Eng OnLine 2007;6:42. [23] Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans Med Imaging 2001;20:1131–9. [24] Basser PJ, Jones DK. Diffusion-tensor MRI: theory, experimental design and data analysis – a technical review. NMR Biomed 2002;15:456–67.
334
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
[25] Mori S, Cain B, Chacko VP, van Zijl PCM. Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann Neurol 1999;45:265–9. [26] Pagani E, Filippi M, Rocca MA, Horsfield MA. A method for obtaining tract-specific diffusion tensor MRI measurements in the presence of disease: Application to patients with clinically isolated syndromes suggestive of multiple sclerosis. Neuroimage 2005;26:258–65. [27] Kanaan RA, Shergill SS, Barker GJ, Catani M, NG VW, Howard R, et al. Tract-specific anisotropy measurements in diffusion tensor imaging. Psychiatry Res 2006;146:73–82. [28] Wilson M, Tench CR, Morgan PS, Blumhardt LD. Pyramidal tract mapping by diffusion tensor magnetic resonance imaging in multiple sclerosis: improving correlations with disability. J Neurol Neurosurg Psychiatry 2003;74:203–7.
[29] Smith SM, Jenkinson M, Johansen-Berg H, Rueckert D, Nichols, Mackay CE, et al. Tract-based spatial statistics: voxelwise analysis of multi-subject diffusion data. Neuroimage 2006;31:1487–505. [30] Smith SM, Johansen-Berg H, Jenkinson M, Rueckert D, Nichols TE, Miller KL, et al. Acquisition and voxelwise analysis of multi-subject diffusion data with tract-based spatial statistics. Nat Protoc 2007;2:499–503. [31] Kunimatsu A, Aoki S, Masutani Y, Abe O, Hayashi N, Mori H, et al. The optimal trackability threshold of fractional anisotropy for diffusion tensor tractography of the corticospinal tract. Magn Reson Med Sci 2004;3:11–7. [32] Rogers BP, Morgan VL, Newton AT, Gore JC. Assessing functional connectivity in the human brain by fMRI. Magn Reson Imaging 2007;25:1347–57.
Magnetic Resonance Imaging 27 (2009) 324 – 334
Intersubject variability in the analysis of diffusion tensor images at the group level: fractional anisotropy mapping and fiber tracking techniques Hans-Peter Müller⁎, Alexander Unrath, Axel Riecker, Elmar H. Pinkhardt, Albert C. Ludolph, Jan Kassubek Department of Neurology, University of Ulm, D – 89081 Ulm, Germany Received 29 November 2007; revised 22 June 2008; accepted 1 July 2008
Abstract Introduction: Diffusion tensor imaging (DTI) provides comprehensive information about quantitative diffusion and connectivity in the human brain. Transformation into stereotactic standard space is a prerequisite for group studies and requires thorough data processing to preserve directional inter-dependencies. The objective of the present study was to optimize technical approaches for this preservation of quantitative and directional information during spatial normalization in data analyses at the group level. Methods: Different averaging methods for mean diffusion-weighted images containing DTI information were compared, i.e., region of interest-based fractional anisotropy (FA) mapping, fiber tracking (FT) and corresponding tractwise FA statistics (TFAS). The novel technique of intersubject FT that takes into account directional information of single data sets during the FT process was compared to standard FT techniques. Application of the methods was shown in the comparison of normal subjects and subjects with defined white matter pathology (alterations of the corpus callosum). Results: Fiber tracking was applied to averaged data sets and showed similar results compared with FT on single subject data. The application of TFAS to averaged data showed averaged FA values around 0.4 for normal controls. The values were in the range of the standard deviation for averaged FA values for TFAS applied to single subject data. These results were independent of the applied averaging technique. A significant reduction of the averaged FA values was found in comparison to TFAS applied to data from subjects with defined white matter pathology (FA around 0.2). Conclusion: The applicability of FT techniques in the analysis of different subjects at the group level was demonstrated. Group comparisons as well as FT on group averaged data were shown to be feasible. The objective of this work was to identify the most appropriate method for intersubject averaging and group comparison which incorporates intersubject variability of the directional information. © 2009 Elsevier Inc. All rights reserved. Keywords: Magnetic resonance imaging; Diffusion tensor imaging; Fractional anisotropy mapping; Fiber tracking
1. Introduction Diffusion in human brain white matter (WM) can be noninvasively mapped by diffusion tensor magnetic resonance imaging (DTI). By DTI, the directional dependence of diffusion in each voxel can be characterized by the so-called ⇉ Y Y Y diffusion tensor D. The eigenvectors ( v 1; ; v 2 ; v 3 ) and ⇉ eigenvalues (λ1, λ2, λ3) of the 3×3 matrix D reflect the diffusivity of water in each direction. This can be used to ⁎ Corresponding author. Tel.: +49 731 177 1206; fax: +49 731 177 1202. E-mail addresses: [email protected] (H.-P. Müller), [email protected] (J. Kassubek). 0730–725X/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2008.07.003
quantify the diffusivity by so-called fractional anisotropy (FA) maps on a voxelwise basis [1–3]. The directional information can be used as the basis for reconstruction of the interconnectivity of brain regions by following the fiber pathways by fiber tracking (FT) techniques. The basis of FT is the consecutive connection of neighbored tensors along their principal directions. Basically, the FT techniques can be divided into two groups: streamline or deterministic FT [4–7] and probabilistic FT [8–14]. The objective of the present study was to optimize the technical approach for the preservation of quantitative and directional information contained in diffusion weighted images/DTI data in order to provide a framework for the comparison of FT of subject groups with defined anatomical
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
WM alterations and groups of normal subjects. To find a difference between two subject groups may eventually assist the diagnosis in finding ways to categorize a subject to one of the group based on some discrimination metric. The prerequisite for this comparison of subject groups was the transformation of the basic DTI data into a stereotactic standard space, e.g., the Montreal Neurological Institute (MNI) space [15]. After this normalization, a comparison at the group level might be performed in several ways: (i) FA maps of different subject groups could be compared, while the directional information within these maps is not considered. (ii) FT results of group-averaged data sets could be qualitatively compared by neuroimaging experts on the basis of standardized human brain atlases [16,17]. The result images were judged in terms of the three factors directionality, cortical projection areas and homogeneity of fiber bundles. (iii) Skeletons that consist of fiber tracts could be used for quantitative comparison on the basis of the directional dependence. In order to show the validity and power of the different techniques, the following questions were to be addressed: • How could intersubject averaging and group comparison be performed and which were the most appropriate methods to quantify differences between subject groups? • How to consider intersubject variability of the directional information? Like in other advanced magnetic resonance imaging (MRI) methods, DTI- and FT-based studies pursue the ultimate goal to categorize individual patients' brain morphology in order to facilitate the diagnostic process based on some discrimination metric. In the present study, the standardization of technical approaches to group-based analyses of certain brain pathologies can be considered the prerequisite for the identification of pathological patterns of altered brain anatomy. As an example for a groupwise comparison, differences between patients with atrophy of the corpus callosum (CC) and age-matched healthy controls were mapped. As a model of CC alteration, subjects with the neurodegenerative disease of complicated hereditary spastic paraparesis were investigated as a prototype of morphological alterations (thinning) along the whole structure of the CC [18,19]. The CC was chosen as the most appropriate structure in the brain to be analyzed since it is one of the white matter structures with the most strongly directed fibers [20]. 2. Data recording and data preprocessing 2.1. Data recording and standard data preprocessing All DTI data were acquired on the same 1.5T scanner (Symphony, Siemens Medical, Erlangen, Germany). Six
325
healthy controls (three men, three women, average age 32.7±4.5 years), and six patients with thinned CC (tCC) (three men, three women, 32.5±12.1 years old) underwent the MRI protocol. All DTI acquisitions consisted of 13 volumes (45 slices, 128×128 voxels, slice thickness 2.2 mm, in-plane voxel size 1.5×1.5 mm), representing 12 gradient directions and one scan with gradient 0 (b0). Echo time (TE) and repetition time (TR) were 93 ms and 8000 ms, respectively. b was 800 s/mm2, five scans were k-space-averaged online by the Siemens SYNGO operating software. As an anatomical background, a high-resolution T1-weighted, magnetization-prepared, rapid-acquisition gradient echo (MPRAGE) sequence was used (TR=9.7 ms, TE=3.93 ms, flip angle 15°, matrix size 256×256 mm2, voxel size 1.0×0.96×0.96 mm 3 ), consisting of 160–200 sagittal partitions depending on the head size. All analyses were performed by the software package (Tensor Imaging and Fiber Tracking) [21]. Standard image processing procedures such as eddy current correction and transformation to iso-voxels and smoothing have already been described previously [21,22]. Also, the spatial normalization protocol to MNI standard space has also been described in [21]. All data sets were preprocessed and normalized accordingly. 2.2. MNI normalization and storage of major eigenvectors Basically, a complete nonlinear MNI normalization consisted of three deformation components (DC): • DC 1: a rigid brain transformation to align the basic coordinate frames. The rotation angles had to be stored ⇉ in a rotation matrix R. • DC 2: an affine deformation according to landmarks. The six affine deformation parameters for the different brain regions had to be stored in an array with 6 Y components S. • DC 3: a non-affine normalization equalizing non-linear brain shape differences. The 3D vector shifts were different for each voxel leading to a separate transformation for each voxel of the 3D voxel array Y Y Y (a 5D transformation array T with indices image row, image column, image slice, dilation matrix row, dilation matrix column) (details see below, Eq. 3). Consequently, for normalized data the resulting diffusion ⇉
tensor D i of each voxel i had to be rotated according to all the rotations listed above. • A rotation resulting from the aligning to the basic coordinate frame (corresponding to DC 1) had to be applied ⇉
⇉
⇉
D i V¼ R D i
ð1Þ Y
Y
Y
• The components of the eigenvectors (V 1 ;V 2 ;V 3 ) had to be adapted according to the six affine
326
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334 Y
normalization parameters of S (dependent on the brain region sa, a=1...6) of the affine deformation (corresponding to DC 2).
and the standard deviation of averaged FA
vw;j j ¼ sa vw;j j
Fstd
ð2Þ
with w=1...3 and j=x,y,z. After adapting, the eigenvectors could have lost their unity length and had to be normalized. • Standard trigonometry gives a rotation matrix (for each voxel independently), resulting from the 3D vector shifts following the concepts of Alexander et al. [23] in order to preserve the directional relations between eigenvectors of neighbored voxels, different shifts of two neighbored voxels result in rotations of the corresponding eigenvectors. The dilation⇉matrices were used for the alignment of the tensor D of each voxel to the surrounding voxels (corresponding to DC 3). ⇉
⇉
⇉
DiW ¼ t i Di V ⇉
ð3Þ Y Y Y
where t i are the components of T . These⇉fine-corrections (according to Eqs.2 and 3) of the tensors D i were essential for a correct FT [22], and the corresponding parameters had to be stored for each subject data set separately. With the standard normalization procedure, the corresponding MPRAGE data were also normalized to MNI space. From the average of these data sets, a display background was created. 3. Analysis techniques 3.1. Region of interest analysis of FA maps The diffusion anisotropy could be quantified by FA: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u3 k1 k 2 þ k 2 k 2 þ k 3 k 2 F¼t 2 k21 þ k22 þ k23
ð4Þ
1 XN Y Y F i;j f s v i jbrs i N
ð6Þ
where N was the ROI size. Extensive studies via ROI analyses about the preservation of the DTI-specific parameters during the process of MNI normalization had been performed previously [21]. The FA maps were calculated separately for each subject data in MNI space, making statistical comparison of the subject groups feasible. Whole brain-based statistical analysis (WBSA) was performed by t test, i.e., voxelwise statistical comparison of normal subject FA maps and FA maps from subjects with tCC. 3.2. Fiber tracking methods For the calculation of the diffusion spheroid, the eigenvector corresponding to the major eigenvalue indicated the fiber direction in WM regions. Based on this directional information, different methods and algorithms had been proposed to estimate WM connectivity. In the present study, the deterministic streamline tracking (DST) technique [24,25] and a technique based on probabilistic fiber tracking (PFT) [8,9], which models the propagation in the major eigenvector field of the brain were applied. Generally, the FT positions resulted from floating point numbers. The corresponding eigenvector direction (pointing towards the consecutive FT position) was the interpolation of the directions of the neighboring voxels weighted by the proportionate position (linear nearest neighbor interpolation) Y
v new ði; j; k Þ ¼
X8 w¼i
Y
aw v ðlw ; mw ; nw Þ:
ð7Þ
Y
where k was the arithmetic average of all eigenvalues. The intensity was related to the FA value and the color-coding was as follows: red for major eigenvector mainly in left-right direction, blue for major eigenvector mainly in inferiorsuperior direction, and green for major eigenvector mainly in posterior-anterior direction. A region of interest (ROI)-based approach of defined brain regions was the method of choice for the analysis of separate areas. An ROI was set by a sphere with an observer-defined radius rs and its center at the user-defined Y focus f s . All FA values Fi of the N voxels inside this sphere Y (position v i ) were included in the following parameterization: the average FA, i.e., average diffusion strength Favg ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 XN Y Y ¼ F F avg i i;j f s v i jbrs N 1
ð5Þ
Here, v new (i,j,k) was the resulting vector at the new position i,j,k (floating point numbers), and l,m,n were the voxel coordinates (integer numbers) of the eight neighboring voxels. The factors aw were the respective eight weighting factors for the interpolation. The threshold for the scalar product of the major eigenvectors (proportional to the angle between the major eigenvectors of two consecutive FT positions) was set to 0.9. The distance between two FT positions, i.e., the stepwidth, was set to 0.5 mm, corresponding to 0.5 voxels (this value of 0.5 mm had also been suggested by [4]). 3.2.1. Modified DST technique The DST approach (Fig. 1A) was modified in the following way: Together with the selected FT starting point, additional 26 fiber tracts started at the 26 directly neighbored voxels i. At each FT step, the mean and
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
327
Fig. 1. Methods of FT for single subject data. (A) DST; starting points were located in a ROI of the radius 2 voxels (2 mm). (B) MDST; compared to (A) runaway fiber tracts were avoided. (C) MPFT. (D) MPFT with display of the end points of fiber tracts — the color bar encodes the probability that a voxel is crossed by a fiber tract (P) in (C) and (D). The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the FA maps. Y
the standard deviation of these 27 vectors v new;j were calculated. 1 X27 ðk Þ kÞ vðavg ¼ v ð8Þ i¼1 i 27 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 X27 ðk Þ ðk Þ ðk Þ 2 vstd ¼ vi vavg ð9Þ i¼1 26 Y
where k indicates the components of the vector v i. Vectors which deviated by more than the standard deviation from the average vector were replaced by the average vector in order to avoid runaway fiber tracts (Fig. 1B). 3.2.2. Modified PFT technique PFT was modified in the following way: a volume with a radius of three voxels around each FT starting point was defined. Then, based on Monte Carlo simulations, additional 10,000 randomly chosen starting points (Gaussian-like distribution) were created. Runaway fiber tracts were avoided by the same technique as described above (Section 3.2.1, Eqs. 8 and 9). A probability map could be set up, which was weighted by the number of fiber tracts intersecting each voxel (Fig. 1C). If only connected brain regions were of interest, rather than the discrete fiber tracts, then the only regions displayed were those in which fiber tracts ended (Fig. 1D). 3.2.3. Intersubject fiber tracking Based upon the directional information for all subject data, a new FT technique could be developed, called “intersubject
FT” (ISFT). For each of the M group subject data sets (M=6 in this study), the identical starting point was defined. The average and the standard deviation of the M new vectors were calculated at each FT step. Vectors deviating more than the standard deviation from the average vector were replaced by the average vector. Note that the new FT position could be different for each subject data. In this way, runaway fiber tracts were avoided and an FT was set up that took the inter-subject variability of the eigenvectors at each FT step into account. 3.3. Tractwise fractional anisotropy statistics The fiber tracts could be used as a skeleton mask for selective statistics [9,26–30]. For that purpose, tractwise fractional anisotropy statistics (TFAS) was used [22]. Each voxel that was crossed by a fiber tract was defined as “active” for statistics and thus contributed to the skeleton mask so that the underlying FA values of the voxels contributed to the statistical t test. The average and standard deviation were calculated by 1 XN TFAS ¼ F ; iaskeleton ð10Þ Favg i¼1 i N TFAS Fstd
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 XN ¼ Favg Fi ; iaskeleton i¼1 N 1
ð11Þ
Separate skeletons were used for volunteers and for tCC subjects, i.e., FT was performed both on group averaged normal subject data and on group averaged tCC data. In
328
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
this way, two skeletons were derived, and the comparison was performed by applying skeleton 1 to group averaged healthy subject FA maps and skeleton 2 to group averaged tCC FA maps. 3.4. Averaging of DTI data 3.4.1. Storage of fine corrections Each DTI data was normalized, and all DTI data were averaged before FA mapping. The fine corrections were also arithmetically averaged (Section 2.2). This approach resulted in an averaged fine correction to be applied to the resulting averaged DTI data.
3.4.2. Storage of fine-corrected eigenvectors Major eigenvectors were stored after fine-correction, and a new major eigenvector array was created for the averaged DTI data. For this purpose, a representative major eigenvector was selected by the criterion XM Y Y v d vj ð12Þ Si ¼ j¼1; j p i i The major eigenvector which had the maximum Si was selected at a respective voxel position out of all major eigenvectors of the M different subjects. This new major eigenvector replaced the major eigenvector calculated from
Fig. 2. WBSA on a voxelwise basis. Comparison of FA maps of six healthy subjects (A) and six subjects with tCC (B). The histograms indicate the distribution of the absolute FA values within the ROI. (C) t Test analysis, display threshold Pb.05. The results are overlaid on an MPRAGE template (averaged from all 12 subjects). The ROI used for analysis (Table 1) is indicated in light blue.
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
the averaged DTI data at the corresponding voxel position. This method required increased memory space and approximately doubled computation time.
For the group comparison, averaging of results for different subjects was necessary. After all subjects' data had been transformed into MNI space, averaging was feasible. 4.1. FA map averaging without directional information WBSA was applied to FA maps of the 12 subjects contributing to the study. Fig. 2 shows the results for WBSA at Pb.05. For tCC, major differences could be detected in the CC. Therefore, a ROI with a radius of 45 mm (45 voxels) was defined including the CC, i.e., the ROI sphere included approximately 380,000 voxels. Average FA values were 0.42±0.17 for normal controls and 0.34±0.14 for tCC patients, respectively. The reduction of FA values was also observed in the histograms of the FA values within the ROI (Fig. 2A and B). The main contribution resulted from the left-right direction (red) (Table 1), compatible with the normal fiber distribution in the CC. By use of t test, approximately 51,000 voxels within the ROI were observed to pass the threshold of Pb.05 (Fig. 2C). 4.2. Single subject fiber tracking The modified PFT (MPFT) results for single subject data using three starting points in the CC according to Table 2 are displayed in Fig. 3. TFAS analysis showed averaged FA values between 0.33 and 0.43 for the different subjects (Table 3). 4.3. Averaging of DTI data from different subjects Averaged DTI data were created by the storage of fine corrections method (Section 3.4.1). The corresponding FT (starting points, cf. Table 2) is displayed in Fig. 4A and B for modified DST (MDST) and MPFT, respectively. Visual inspection showed minor differences. Differences mainly occurred in juxtacortical regions with low FA values. The method of storage of fine-corrected eigenvectors (Section
Table 1 ROI analysis and comparison for subject groups (see also Fig. 2) ROI position MNI: 0–11 mm, 18 mm; ROI radius: 45 mm, ROI size: 380,000 voxels Healthy controls
tCC patients
Favg±Fstd
Total Left-right Superior-inferior Anterior-posterior
Table 2 MNI coordinates for starting points for the fiber tracking (FT)
Starting point 1 Starting point 2 Starting point 3
4. Results
0.42±0.17 0.35±0.14 0.37±0.14 0.56±0.16
t test 51,000 voxels (13%) passing threshold of Pb.05
0.34±0.14 0.26±0.08 0.31±0.09 0.51±0.11
329
x/mm
y/mm
z/mm
0 0 0
−40 −10 21
17 29 8
3.4.2) showed the FT (starting points, cf. Table 2) displayed in Fig. 4C and D for MDST and MPFT, respectively. MDST and MPFT showed corresponding results independent of the averaging method in regions near the CC (identification on the basis of standardized human brain atlases [16,17]). Here, only descriptive comparison is feasible. The comparison was performed by four neuroimaging experts (J.K., A.U., A.R., E.P.) on the basis of standardized human brain atlases [16,17]. The resulting images were judged in terms of the three factors directionality, corresponding cortical projection areas, and homogeneity of fiber bundles (three-item scale: good–moderate–bad concordance of the results). With respect to directionality and anatomical localization of the cortical projection areas, the concordance was always judged as “good.” The homogeneity of the fiber bundles was judged as “good” in the white matter regions adjacent to the connection line between anterior commissure and posterior commissure, whereas it was judged as “moderate” in juxtacortical regions, i.e., differences between the two averaging methods mainly occurred in juxtacortical regions where the FA values were b0.2. 4.4. Intersubject fiber tracking ISFT was based on the results of DTI analysis of single subject data. By taking the intersubject variability of the eigenvectors at each FT step into account, the result displayed in Fig. 4E was achieved. 4.5. Tractwise fractional anisotropy statistics The qualitative comparison by visual analysis (see Section 4.3) was extended by a quantification of the differences of the FT results as performed by TFAS. The fiber tracts obtained by the different FT techniques could be used as skeletons for TFAS, which will be described in the following. The averaged FA values for skeletons obtained by MPFT (starting points in the CC) in each subject (Eqs. 10 and 11) are summarized in Table 3. Alternatively, skeletons based on FT on group averaged data were used for statistical comparison between the normal subject group and the tCC subject group (Fig. 5). For the healthy control group (fine-correction method, Section 4.3), the average FA values were 0.35 for MPFT, 0.39 for MDST and 0.38 for ISFT, respectively. These values did not differ significantly from the single subject results (Table 3) and were significantly different from averaged FA values for the tCC patient group (fine-correction method). The results are summarized in Table 4.
330
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
Fig. 3. Fiber tracking results for different normal subjects by MPFT. FA threshold was 0.1, display background was MPRAGE template (averaged from all 12 subjects who participated in the study). The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set.
4.6. Interconnectivity of brain regions
5. Discussion
Fiber tracking could be used to investigate interconnectivity of brain regions (Fig. 6). The end points of all fiber tracts were displayed. An averaged fiber tract (from all fiber tracts contributing to the end point distribution) was created to exhibit the interconnectivity. The course of these averaged fiber tracts as well as the regions identified to be connected in frontal, parietal and occipital brain areas correlate well with neuroanatomical knowledge of commissural fibers. Here, the averaged fiber tract connecting the frontal brain areas of both hemispheres correlated with parts of the minor forceps (frontal part of CC), while interconnectivity within occipital brain areas corresponded to the major forceps (occipital part of CC).
MNI normalization and comparison of FA maps at the group level allowed for several possibilities to quantify differences between subject groups. 5.1. FA map analysis WBSA could be applied to normalized FA maps of controls and subjects with tCC to obtain a 3D comparison visualized using maps with color coded significant regions. The main contribution to differences from a discrete region could be analyzed via histograms. Nevertheless, the directional dependence was not taken into account by these analysis methods. 5.2. FT analysis
Table 3 TFAS for normal subjects (MPFT) TFAS FTFAS avg ±Fstd
Subject 1 Subject 2 Subject 3 Subject 4 Subject 5 Subject 6
0.43±0.20 0.37±0.23 0.39±0.19 0.33±0.22 0.41±0.21 0.35±0.23
In order to include directional information, averaging techniques for directional information were applied. The results of visual comparison between single-subject FT and the FT on group averaged data showed for the two averaging techniques (Sections 3.4.1 and 3.4.2) differences mainly in regions with low FA values (FAb0.2; this threshold had also been suggested by [31]). Thus, both averaging methods were considered appropriate. Never-
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
331
Fig. 4. Different analyses of FT in the averaged DTI data from 6 normal controls, depending if fine-corrections were averaged (A, B) or storage of major eigenvector was used (C, D). (A) MDST-averaged fine-corrections. (B) MPFT-averaged fine-corrections. (C) MDST storage of major eigenvector. (D) MPFT storage of major eigenvector. (E) ISFT. The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set.
theless, the technique based on fine-corrections was considered to be superior because of its advantages in memory usage and computation time.
The application of FT to tCC subject data caused problems due to low FA values in the CC. Variations in FT between the healthy group and the tCC group were mainly
Fig. 5. ISFT: (A) calculated from averaged DTI data of 6 normal controls. (B) Calculated from averaged DTI data of 6 subjects with tCC. The FT is color coded according to FA values. The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set.
332
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
Table 4 TFAS with different skeletons, depending on the FT mode Healthy subjects
tCC patients
TFAS FTFAS avg ±Fstd
MDST MPFT ISFT
0.39±0.18 0.35±0.16 0.38±0.18
t test P
0.19±0.08 0.18±0.07 0.19±0.10
b.001 b.001 b.001
not related to changes in the diffusion direction; moreover, the low values in mean diffusion weighted images (as detected by low FA values) led to errors in the diffusion tensors, and thus, information on the directional dependency was lost. 5.3. TFAS There have been various efforts to establish diffusion anisotropy as a marker for white matter tract integrity [26–28] by use of the underlying FA maps for selective statistics. Smith et al. [29,30] developed an algorithm for an
alignment-invariant tract representation to overcome normalization problems; this approach was referred to as tract-based spatial statistics [9,29,30]. In the present study, an alternative approach (TFAS) was applied, as previously described [22]. Here, bundles of fiber tracts were used in the sense of a skeleton which represented the basis of statistical analysis of the underlying FA maps. The novel character of TFAS was its use of averaged DTI data, i.e., the processing steps of normalization and averaging were not performed on FA maps but on one newly created group averaged DTI data (for a detailed description of the TFAS method, see Ref. [22]). 5.4. Interconnectivity of brain regions Interconnectivity of brain regions has repeatedly been a topic of research with various MRI techniques such as functional MRI [32] and other DTI-based approaches. In the present study, the investigation of interconnectivity was restricted to normal subjects due to directional loss in eigenvector calculation from DTI data of patients with tCC, e.g., the FA signal in the CC was sometimes low (FAb0.2)
Fig. 6. Interconnectivity display of brain regions — ISFT with starting points in the CC. Only the end points of the FT process are displayed. Coronal (A), sagittal (B), axial (C) and 3D view (D). The images are glass-brain plots of fiber tracts overlaid on glass-brain plots of the averaged MPRAGE data set. Interconnectivity (derived from an averaged FT) is indicated as yellow lines following the fiber tracks. The color bar encodes the probability that a voxel is an end point of a fiber tract (P).
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
for the tCC data sets and, thus, sometimes led to an erroneous FT. Kunimatsu et al. [31] also suggested an FA value of 0.2 as the optimal trackability threshold. Generally, this technique has to be applied very thoroughly because low FA values (as usually occur in juxtacortical positions) might lead to incorrect FT and, thus, to an incorrect localization of fiber tract end points. 6. Conclusions Intersubject averaging of DTI data could be performed with respect to diffusion amplitude (by FA information) and diffusion direction (by FT). Averaging of FA maps allowed for the statistical comparison of subject groups with the WBSA method. Comparison of methods for maintaining directional information during FT data processing was performed. This work gives an instruction how intersubject averaging and group comparison could be performed and how intersubject variability of the directional information could be considered. In the following, the two questions from the introduction are to be addressed in detail. (i) How could intersubject averaging and group comparison be performed and which were the most appropriate methods to quantify differences between subject groups? The averaging methods showed similar results; therefore, systematic errors during the spatial normalization process could be excluded. Differences between the method of averaged fine-corrections and the method of stored major eigenvectors occurred in juxtacortical brain positions where the FA values were known to be b0.2 so that correct FT was impossible. The major eigenvector method has the drawback of huge memory consumption, i.e., in case of large patient number studies. (ii) How to consider intersubject variability of the directional information? ISFT is a new technique which takes into account inter-subject variability of directional information during FT. All FT techniques used (i.e., MDST, MPFT, ISFT) showed similar results for TFAS analysis. The results were in agreement with those of the TFAS analysis in single subjects. Therefore, DTI data averaging and FT calculated from these averaged DTI data was considered to be applicable. In summary, the averaging methods, FT techniques and FA quantification approaches (ROI analyses as well as TFAS) were demonstrated to be applicable to the comparison between different subjects. In this study, the use of the technical approach in alterations of the CC was specifically addressed and future studies will be necessary to demonstrate the applicability to other WM brain alterations. References [1] Basser PJ, Mattiello J, LeBihan D. Estimation of the effective selfdiffusion tensor from the NMR spin echo. J Magn Reson B 1994;103:247–54. [2] Basser PJ, Mattiello J, LeBihan D. MR Diffusion Tensor Spectroscopy and Imaging. Biophys J 1994;66:259–67.
333
[3] Pierpaoli C, Jezzard P, Basser PJ, Barnett A, DiChiro G. Diffusion tensor MR imaging of the human brain. Radiology 1996;201:637–48. [4] Conturo, Lori NF, Cull TS, Akbudak E, Snyder AZ, Shimony JS, et al. Tracking neuronal fibre pathways in the living human brain. Proc Nat Acad Sci U S A 1999;96:10422–7. [5] Lori NF, Akbudak E, Shimony JS, Cull TS, Snyder AZ, Guillory RK, et al. Diffusion tensor fibre tracking of human brain connectivity: acquisition methods, reliability analysis and biological results. NMR Biomed 2002;15:494–515. [6] Mori S, Kaufmann WE, Davatzikos C, Stieltjes B, Amodei L, Fredericksen K, et al. Imaging cortical association tracts in the human brain using diffusion-tensor-based axonal tracking. Magn Reson Med 2002;47:215–23. [7] Lazar M, Weinstein DM, Tsuruda JS, Hasan KM, Arfanakis K, Meyerand ME, et al. White matter tractography using diffusion tensor deflection. Hum Brain Mapp 2003;18:306–21. [8] Behrens TE, Johansen-Berg HJ, Woolrich MW, Smith SM, WheelerKingshott CA, Boulby PA, et al. Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging. Nat Neurosci 2003;6:750–7. [9] Behrens TE, Berg HJ, Jbabdi S, Rushworth MF, Woolrich MW. Probabilistic diffusion tractography with multiple fibre orientations: what can we gain? Neuroimage 2007;34:144–55. [10] Hagmann P, Thiran JP, Jonasson L, Vandergheynst P, Clarke S, Maeder P, et al. DTI mapping of human brain connectivity: statistical fibre tracking and virtual dissection. Neuroimage 2003;19:545–54. [11] Parker GJ, Haroon HA, Wheeler-Kingshott CA. A framework for a streamline-based probabilistic index of connectivity (PICo) using a structural interpretation of MRI diffusion measurements. J Magn Reson Imaging 2003;18:242–54. [12] Parker GJ, Alexander. Probabilistic anatomical connectivity derived from the microscopic persistent angular structure of cerebral tissue. Philos Trans R Soc Lond B Biol Sci 2005;360:893–902. [13] Jones DK, Pierpaoli C. Confidence mapping in diffusion tensor magnetic resonance imaging tractography using a bootstrap approach. Magn Reson Med 2005;53:1143–9. [14] Mori S, van Zijl PC. Fibre tracking: principles and strategies — a technical review. NMR Biomed 2002;15:468–80. [15] Brett M, Johnsrude IS, Owen AM. The problem of functional localization in the human brain. Nat Rev Neurosci 2002;3:243–9. [16] Talairach J. Tournoux. Coplanar stereotaxic atlas of the human brain. New York: Thieme Medical; 1988. [17] Mori S, van Zijl PCM, Wakana S, Nagae-Poetscher LM. MRI Atlas of Human White Matter. Amsterdam: Elsevier; 2005. [18] Kassubek J, Juengling FD, Baumgartner A, Ludolph AC, Sperfeld AD. Different regional brain volume loss in pure and complicated hereditary spastic paraparesis: a voxel-based morphometry study. Amyotrophic Lateral Sclerosis 2007;8:328–36. [19] Fink JK. Hereditary spastic paraplegia. Curr Neurol Neurosci Rep 2006;6:65–76. [20] Sullivan EV, Adalsteinsson E, Pfefferbaum A. Selective age-related degradation of anterior callosal fibre bundles quantified in vivo with fibre tracking. Cerebral Cortex 2006;16:1030–9. [21] Müller HP, Unrath A, Ludolph AC, Kassubek J. Preservation of diffusion tensor properties during spatial normalization by use of tensor imaging and fibre tracking on a normal brain database. Phys Med Biol 2007;52:N99–N109. [22] Müller HP, Unrath A, Sperfeld AD, Ludolph AC, Riecker A, Kassubek J. Diffusion tensor imaging and tractwise fractional anisotropy statistics: quantitative analysis in white matter pathology. Biomed Eng OnLine 2007;6:42. [23] Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans Med Imaging 2001;20:1131–9. [24] Basser PJ, Jones DK. Diffusion-tensor MRI: theory, experimental design and data analysis – a technical review. NMR Biomed 2002;15:456–67.
334
H.-P. Müller et al. / Magnetic Resonance Imaging 27 (2009) 324–334
[25] Mori S, Cain B, Chacko VP, van Zijl PCM. Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann Neurol 1999;45:265–9. [26] Pagani E, Filippi M, Rocca MA, Horsfield MA. A method for obtaining tract-specific diffusion tensor MRI measurements in the presence of disease: Application to patients with clinically isolated syndromes suggestive of multiple sclerosis. Neuroimage 2005;26:258–65. [27] Kanaan RA, Shergill SS, Barker GJ, Catani M, NG VW, Howard R, et al. Tract-specific anisotropy measurements in diffusion tensor imaging. Psychiatry Res 2006;146:73–82. [28] Wilson M, Tench CR, Morgan PS, Blumhardt LD. Pyramidal tract mapping by diffusion tensor magnetic resonance imaging in multiple sclerosis: improving correlations with disability. J Neurol Neurosurg Psychiatry 2003;74:203–7.
[29] Smith SM, Jenkinson M, Johansen-Berg H, Rueckert D, Nichols, Mackay CE, et al. Tract-based spatial statistics: voxelwise analysis of multi-subject diffusion data. Neuroimage 2006;31:1487–505. [30] Smith SM, Johansen-Berg H, Jenkinson M, Rueckert D, Nichols TE, Miller KL, et al. Acquisition and voxelwise analysis of multi-subject diffusion data with tract-based spatial statistics. Nat Protoc 2007;2:499–503. [31] Kunimatsu A, Aoki S, Masutani Y, Abe O, Hayashi N, Mori H, et al. The optimal trackability threshold of fractional anisotropy for diffusion tensor tractography of the corticospinal tract. Magn Reson Med Sci 2004;3:11–7. [32] Rogers BP, Morgan VL, Newton AT, Gore JC. Assessing functional connectivity in the human brain by fMRI. Magn Reson Imaging 2007;25:1347–57.