Group statistics of DTI fiber bundles using spatial functions of tensor measures

We present a framework for hypothesis testing of differences between groups of DTI ber tracts. An anatomical, tract-oriented coordinate system provides a basis for estimating the distribution of diffusion properties. The parametrization of sampled, smooth functions is normalized across a population using DTI atlas building. Functional data analysis, an extension of multivariate statistics to continuous functions is applied to the problem of hypothesis testing and discrimination. B-spline models of fractional anisotropy (FA) and Frobenius norm measures are analyzed jointly. Plots of the discrimination direction provide a clinical interpretation of the group differences. The methodology is tested on a pediatric study of subjects aged one and two years.

Group Statistics of DTI Fiber Bundles Using Spatial Functions of Tensor Measures Casey B. Goodlett1;2, P. Thomas Fletcher1;2, John H. Gilmore3, and Guido Gerig1;2 1 School of Computing, University of Utah 2 Scientic Computing and Imaging Institute, University of Utah 3 Department of Psychiatry, University of North Carolina? Abstract. We present a framework for hypothesis testing of dierences between groups of DTI ber tracts. An anatomical, tract-oriented coordi- nate system provides a basis for estimating the distribution of diusion properties. The parametrization of sampled, smooth functions is nor- malized across a population using DTI atlas building. Functional data analysis, an extension of multivariate statistics to continuous functions is applied to the problem of hypothesis testing and discrimination. B-spline models of fractional anisotropy (FA) and Frobenius norm measures are analyzed jointly. Plots of the discrimination direction provide a clinical interpretation of the group dierences. The methodology is tested on a pediatric study of subjects aged one and two years. 1 Introduction The diusion properties of white matter tracts measured by DTI provide a novel and important source of information for group comparison and regression in clin- ical neuroimaging studies. Signicant challenges remain in the development of an automatic framework for testing signicance of group dierences in a manner which provides clinically relevant results. Previous work has shown the impor- tance of modeling the diusion properties of a ber tract as functions sampled by arc length along the axis of the bundle [1, 2]. The major challenge in applying this type of analysis is the need for a consistent parametrization of ber bundles across a large population of images. Deformable registration has been proposed as a method of mapping a population to a reference atlas coordinate system [3{5]. Most of the analysis using atlas building has focused on independent voxelwise tests, which can be challenging to interpret and require sophisticated multiple comparison correction. Most studies have also analyzed fractional anisotropy (FA) or mean diusivity (MD) values independently. We propose to combine the anatomically relevant coordinate system of tract statistics with the popu- lation coordinate system provided by atlas building. The combination of the ? This work is part of the National Alliance for Medical Image Computing (NA-MIC), funded by the National Institutes of Health through Grant U54 EB005149 as well as support from the NIMH Silvio O. Conte Center Grant MH064065. Dr. Fletcher is partially funded through an Autism Speaks Mentor-Based Postdoctoral Fellowship. UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript tract coordinate system with atlas building provides an automated, clinically in- terpretable framework for understanding group dierences. The closest related work has been done using nonlinear registration and projection onto a skele- ton representation of FA [6]. Another proposed approach uses ber clustering to compute correspondence across a population [7]. We use deformable registration to estimate and remove shape variability in a population of images. Analysis of shape normalized ber bundles is performed in an anatomically relevant coordinate system based on ber tractography. The at- las normalized diusion measures are treated as a continuous smooth function of arc length, and statistical tests are applied for the joint analysis of the orthogonal FA and Frobenius norm measures. The framework provides a single multivariate hypothesis test between groups eliminating the need for multiple comparison correction and incorporating the joint information of tensor anisotropy and size. Visualization of the linear discriminant provides a clinically meaningful interpre- tation of the group dierences as shown in an example study of pediatric images. Fig. 1 shows an overview of the analysis procedure. Images Atlas Atlas Tract Mapped Tracts Sampled Functions Functional Statistics Fig. 1. Schematic overview of the tract analysis procedure. 2 Atlas Parametrization and Fiber Extraction Atlas building based on dieomorphic registration estimates a set of transfor- mations such that each image in the population can be mapped into the atlas coordinate system. In the DTI atlas building framework, each tensor image Ii is mapped into a common atlas space by a transformation Ti with inverse T􀀀1 i using appropriate measures to account for reorientation and interpolation of ten- sors. For the study presented here, the procedure of Goodlett et al. based on the atlas building procedure of Joshi et al. is applied [4, 8]. In our framework, images from two groups are combined to produce a single pooled atlas. Our assumption UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript is that the overall topology of the images in the two groups are similar enough to allow all images to be combined into one atlas, but dierences may occur in the diusion properties of ber tracts. Thus, we use registration to normalize the image shapes and perform statistics on the diusion properties of the nor- malized ber bundles. The set of tensor images are averaged in atlas space to produce an atlas tensor image with improved signal-to-noise ratio (SNR). The average tensor volume allows reliable extraction of tracts even in populations of images with low SNR such as pediatric images. The dieomorphic transforms guaranteed by the atlas building procedure allow atlas tracts to be mapped back into each individual subject. Fiber tracts are extracted in the average tensor image using a standard Runge-Kutta streamline integration technique based on the principal eigenvector eld. Source and target regions are manually developed to extract each bundle of interest. For each subject, the data within the ber bundle is modeled as a sampled function of arc length using a method similar to that described in Corouge et al. [1]. The result of the procedure is a set of sampled functions parametrized by arc length tj 2 [􀀀a; b] from the atlas ber tract. The atlas- normalized parametrization of the curves is possible because of the smooth, invertible nature of the transformations Ti, T􀀀1 i . That is the samples from each subject are obtained by FAi(Ti(tj)) for each sample point tj in the atlas tract. A sampled function is created for each tensor scalar measure such as fractional anisotropy, mean diusivity (MD), Frobenius norm kDk, etc. For the purpose of this work we chose FA and the Frobenius norm as orthogonal anisotropy and size measures respectively [9]. Fig. 2 shows the sampled curves extracted for the genu ber bundle for our example study. (a) Genu atlas tract −40 −30 −20 −10 0 10 20 30 40 50 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 one year two year (b) All FA curves −40 −30 −20 −10 0 10 20 30 40 50 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 x 10−3 one year two year (c) All norm curves Fig. 2. (a) Genu tract extracted from the tensor atlas colored by FA value. The diusion values are sampled along the atlas-normalized arc length for each individual in the study for FA and Frobenius norm values. The sampled FA and Frobenius norm functions for the two groups are shown in (b) and (c) respectively. The one year old subjects are the dashed red lines and the two year old subjects are the solid blue lines. The spikes in the center of the Frobenius norm functions are likely partial voluming with the ventricles. UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript 3 Functional Data Analysis Image sampling as well as the ber tract extraction process create a sampled representation of the ber bundle diusion properties. However, there exists a continuous underlying biology which generates these samples. Therefore, statis- tical analysis of the sampled diusion functions must account for the underly- ing continuity and spatial correlation of the samples. We compute statistics of the diusion curves as an innite-dimensional extension to multivariate statis- tics known as functional data analysis [10]. The simplest extensions of ordinary statistics to the functional setting are the sample mean function f(t) and the sample variance-covariance function v(s; t), which is the bivariate function given by f(t) = 1 N XN i=1 fi(t), and v(s; t) = 1 N 􀀀 1 XN i=1 (fi(s) 􀀀 fi(s))(fi(t) 􀀀 fi(t)) (1) The diagonal of the function, v(t; t), is the pointwise variance. Hypothesis testing and discriminant analysis of the space of functions has an inherent high dimension low sample size problem, because of the innite dimensional space of continuous functions. Regularization methods are, therefore, essential in the computation of functional statistics. To enforce regularity, B-spline tting and functional principal components analysis (PCA) is used for data driven smooth- ing where the number of retained PCA modes acts as a smoothing parameter. In order to make computations tractable smooth basis functions are t to the sampled diusion curves. B-splines were selected as basis functions due to the nonperiodic nature of the data, the compact support of the B-spline basis, and the ability to enforce derivative continuity. A large number of B-spline bases are rst t to the sampled functions using a least squares approach. The number of basis functions is chosen subjectively to maintain local features while providing some smoothing. Computation of the mean function is computed by the sam- ple mean of the B-spline coecients. Computation of the variance-covariance function requires accounting for mapping between basis coecients and func- tion values. Let fi(t) be the B-spline function t to the samples from subject i. In matrix notation we express all functions fi(t) as a matrix of coecients C times the basis functions f (t) = C(t): (2) Similarly, the variance-covariance function of f (t) can be written as v(s; t) = 1 N 􀀀 1(s)TCTC(t): (3) Principal component analysis (PCA) of the functions fi(t) decomposes v(s; t) into the orthogonal unit eigenfunctions (t) which satisfy Z v(s; t)(t)dt = (s): (4) UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript The B-Spline basis is not orthonormal resulting in a non-symmetric eigenvalue problem to solve (4). As shown in Ramsay and Silverman [10], this minimization can solved by the symmetric eigenvalue problem for the basis coecients b, with the change of variable W1=2u = b as W1=2CTCW1=2u = u; (5) where W is the matrix of basis function inner products with entries Wij = Z i(t)j(j): (6) In our analysis we consider jointly the analysis of FA and tensor norm functions with basis coecients C1 and C2 respectively. We therefore compute PCA from the eigenanalysis of , where ij =W1=2CTi CjW1=2; and = 11 12 21 22 : (7) Hypothesis testing and discriminant analysis is performed on the projection into the rst K PCA modes, where K serves as a smoothing parameter. Let xi and yi be the projection of the curves from the two population of functions fi(t) and gi(t) into the PCA space. In this space the basis mapping has already been incorporated and standard multivariate analysis can be applied. The normal parametric hypothesis test for mean dierences is the Hotelling T2 statistic, T2 = nxny nx + ny (x 􀀀 y)S􀀀1(x 􀀀 y)T (8) where S is the pooled covariance matrix. In order to relax the normality assump- tions associated with the parametric test, we apply a permutation test based on the T2 statistic to compute the p-value. The T2 statistic is proportional to the group mean dierences projected on the Fisher linear discriminant (FLD), ! = S􀀀1(x 􀀀 y)T : (9) The linear discriminant, therefore, provides a direction for interpreting the group dierences. The coecients of the discriminant can be expanded into the original function basis so that FLD(t) = (t)! is a function whose inner product with the original data provides maximal separation between the groups. 4 Pediatric Data Application and Validation We have tested the methodology on a study of pediatric DTI images. A popula- tion of 22 one year old subjects and 30 two year old subjects were chosen from UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript a database of pediatric DTI. In this example we expected to nd large dier- ences between the two groups, and the purpose of this study is to illustrate the methodology rather than to determine clinical results. Each image was acquired with 2x2x2mm3 isotropic voxels, 10 repetitions of a six direction protocol, and a b-value of 1000s=mm2. We selected as representative ber bundles the genu of the corpus callosum and the left motor tract. An atlas was computed from the combined set of 52 images, and tractography was performed to extract the two tracts. Sampled functions of FA and tensor norm parametrized by atlas-normalized arc length were computed in the genu and left motor tracts. For the genu curves, a B-spline basis with 60 basis functions was used to provide preliminary smoothing and smooth curve estimation. For the motor tract, 80 basis functions were used. Functional joint PCA of FA and Frobenius norm was then estimated for the whole population. The number of PCA modes was selected to retain 90% of the total variance. For this study 7 and 11 PCA modes were retained for the genu and motor tracts respectively. The mean function plus the rst two principal modes for the genu tract are shown below in Fig. 3. −40 −30 −20 −10 0 10 20 30 40 50 0.3 0.4 0.5 0.6 arc length fractional anisotropy −40 −30 −20 −10 0 10 20 30 40 50 0.8 0.9 1 1.1 1.2 x 10−3 arc length Frobenius norm (a) Genu mean functions −30 −20 −10 0 10 20 30 40 0.3 0.4 0.5 0.6 fractional anisotropy arc length −30 −20 −10 0 10 20 30 40 0.9 1 1.1 1.2 x 10−3 Frobenius norm arc length (b) Genu PCA mode 1 −30 −20 −10 0 10 20 30 40 0.3 0.4 0.5 0.6 fractional anisotropy arc length −30 −20 −10 0 10 20 30 40 0.9 1 1.1 x 10−3 Frobenius norm arc length (c) Genu PCA mode 2 Fig. 3. Visualization of the PCA modes for the joint analysis of FA and Frobenius norm in the genu tract. The (a) mean functions for the combined population are shown with (b) the rst and (c) second PCA modes. The rst PCA mode accounts for a large percentage of the variability and shows an overall constant change in FA and an anti- correlated constant change in norm. The Hotelling T2 statistic was then computed in PCA space. The genu tract test was extremely signicant with a T2 statistic of 133.1 and parametric p-value of 3.3e-8. The motor tract was also extremely signicant with T2 statistic of 93.8 and a parametric p-value of 2.7e-6. In this case there was such a large dierence between groups that the permutation test did not result in any permutations with a statistic greater than the original. The p-values are uncommonly low because of the strong dierences in the test data and the relatively large sample size. Visualization of the discriminant direction provides an interpretation of the detected dierences and is shown in Fig. 4. The discriminant direction for the UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript −40 −30 −20 −10 0 10 20 30 40 50 0 2 4 6 8 x 10−3 arc−length Fractional Anisotropy −40 −30 −20 −10 0 10 20 30 40 50 −20 −15 −10 −5 0 5 x 10−6 arc−length Frobenius Norm (a) Genu discriminant functions −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Projection on discimination direction Group frequency 1 year 2 year (b) Data functions projected on FLD −60 −50 −40 −30 −20 −10 0 10 20 30 40 50 −2 0 2 4 6 x 10−3 arc−length Fractional Anisotropy −60 −50 −40 −30 −20 −10 0 10 20 30 40 50 −2 −1 0 1 2 3 x 10−5 arc−length Frobenius Norm (c) Motor tract discriminant functions −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Projection on discimination direction Group frequency 1 year 2 year (d) Data functions projection on FLD Fig. 4. Linear discriminants from one to two years for the (a) genu and (c) left motor tracts expanded into original functional basis. These are the functions integrated with the original data that maximally separate the groups. In the genu tract the FA values increase from one to two years, and the Frobenius norm values decrease. For the motor tract, the results are similar for FA, but the norm increases at the base of the tract and decreases towards the top. The projection of the (b) genu and (d) motor tract curves onto the discrimination direction shows the strong separation between the groups. genu tract shows the dierence from one to two year old groups is caused by an overall increase in FA and correlated decrease in Frobenius norm. Furthermore, the increased value of FA in the center of the tract indicates the central region of the tract provides more discriminative power between the two groups. These results are similar to dierences which have been found between neonates and one year old subjects in the same tract [11]. The results in the motor tract indicate a similar constant increase in FA across the whole tract, and the Frobenius norm increases towards the inferior region of the tract, and decreases at a specic location in the superior region of the tract. 5 Conclusions and Discussion Computing ber tract statistics as a function of arc length provides a sensitive mechanism for detecting and understanding changes in ber tract properties be- UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript tween populations. Our framework avoids the problems of multiple comparison correction by providing a single nonparametric hypothesis test for each ber bundle. Furthermore, the discrimination information contained within the hy- pothesis test can be visualized to provide a clinically relevant interpretation of the group dierences. The framework presented here is closely related to previ- ous work on shape analysis using PCA, and we intend to explore in more detail how tools from shape analysis can be applied to this problem. We are currently applying the methodology to a study of Schizophrenia in adults. References 1. Corouge, I., Fletcher, P.T., Joshi, S., Gouttard, S., Gerig, G.: Fiber tract-oriented statistics for quantitative diusion tensor MRI analysis. Medical Image Analysis 10(5) (2006) 786{798 2. Fletcher, P.T., Tao, R., Joeng, W.K., Whitaker, R.: A Volumetric Approach to Quantifying Region-to-Region White Matter Connectivity in Diusion Tensor MRI. In: Information Processing in Medical Imaging. Volume 4584 of LNCS. (2007) 346{358 3. Peyrat, J.M., Sermesant, M., Pennec, X., Delingette, H., Chenyang Xu, McVeigh, E., Ayache, N.: A Computational Framework for the Statistical Analysis of Cardiac Diusion Tensors: Application to a Small Database of Canine Hearts. Medical Imaging, IEEE Transactions on 26(11) (2007) 1500{1514 4. Goodlett, C., Davis, B., Jean, R., Gilmore, J., Gerig, G.: Improved Correspondance for DTI population studies via unbiased atlas building. In: Medical Image Com- puting and Computer Assisted Intervention (MICCAI). Volume 4191 of LNCS., Springer-Verlag (2006) 260{267 5. Zhang, H., Yushkevich, P.A., Rueckert, D., Gee, J.C.: Unbiased white matter atlas construction using diusion tensor images. In: Medical Image Computing and Computer Assisted Intervention (MICCAI). Volume 4791 of LNCS., Springer- Verlag (2007) 211{218 6. Smith, S.M., Jenkinson, M., Johansen-Berg, H., Rueckert, D., Nichols, T.E., Mackay, C.E., Watkins, K.E., Ciccarelli, O., Cader, M.Z., Matthews, P.M., Behrens, T.E.: Tract-based spatial statistics: Voxelwise analysis of multi-subject diusion data. NeuroImage 31 (2006) 1487{1505 7. O'Donnell, L., Westin, C.F., Golby, A.: Tract-Based Morphometry. In: Medical Image Computing and Computer Assisted Intervention (MICCAI). Volume 4792 of LNCS., Springer-Verlag (2007) 161{168 8. Joshi, S., Davis, B., Jomier, M., Gerig, G.: Unbiased dieomorphic atlas construc- tion for computational anatomy. NeuroImage 23(Supplement1) (2004) S151{S160 9. Ennis, D.B., Kindlmann, G.: Orthogonal tensor invariants and the analysis of diusion tensor magnetic resonance images. Magnetic Resonance in Medicine 55(1) (2006) 136{146 10. Ramsay, J., Silverman, B.: Functional Data Analysis. Second edn. Springer (2005) 11. Gilmore, J.H., Lin, W., Corouge, I., Vetsa, Y.S., Smith, J.K., Kang, C., Gu, H., Hamer, R.A., Lieberman, J.A., Gerig, G.: Early Postnatal Development of Corpus Callosum and Corticospinal White Matter Assessed with Quantitative Tractogra- phy. American Journal of Neuroradiology 28(9) (2007) 1789{1795 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript