Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2896997

Formats

Article sections

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 July 28.

Published in final edited form as:

Published online 2010 March 22. doi: 10.1109/TMI.2010.2040625

PMCID: PMC2896997

NIHMSID: NIHMS187341

Hongtu Zhu^{*}

Department of Biostatistics, and Biomedical Research Imaging Center, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA

Martin Styner

Department of Psychiatry, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA

Niansheng Tang

Department of Statistics, Yunnan University, Kunming, Yunnan, P.R.China

Zhexing Liu

Department of Psychiatry, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA

Weili Lin

Department of Radiology, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA

Department of Psychiatry, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

Diffusion tensor imaging (DTI) provides important information on the structure of white matter fiber bundles as well as detailed tissue properties along these fiber bundles *in vivo*. This paper presents a functional regression framework, called FRATS, for the analysis of multiple diffusion properties along fiber bundle as functions in an infinite dimensional space and their association with a set of covariates of interest, such as age, diagnostic status and gender, in real applications. The functional regression framework consists of four integrated components: the local polynomial kernel method for smoothing multiple diffusion properties along individual fiber bundles, a functional linear model for characterizing the association between fiber bundle diffusion properties and a set of covariates, a global test statistic for testing hypotheses of interest, and a resampling method for approximating the *p*-value of the global test statistic. The proposed methodology is applied to characterizing the development of five diffusion properties including fractional anisotropy, mean diffusivity, and the three eigenvalues of diffusion tensor along the splenium of the corpus callosum tract and the right internal capsule tract in a clinical study of neurodevelopment. Significant age and gestational age effects on the five diffusion properties were found in both tracts. The resulting analysis pipeline can be used for understanding normal brain development, the neural bases of neuropsychiatric disorders, and the joint effects of environmental and genetic factors on white matter fiber bundles.

Water molecules in the human brain diffuse preferentially along the major axis of white matter fiber bundles. Thus, Diffusion Tensor Imaging (DTI), which can track the effective diffusion of water in the human brain in vivo, is used to accurately map in vivo the structure and orientation of fiber tracts in the white matter of the brain [1], [2]. The directional dependence of water diffusion in each voxel can be characterized by a matrix, called a diffusion tensor (DT), and the degree of diffusivity can be quantified by using the three eigenvalue-eigenvector pairs of DT and its related parameters, such as fractional anistropy (FA) [3], [4], [5], [6]. There has been a wealth of neuroimaging studies using a set of water diffusion related parameters including FA, mean diffusivity (MD), and the three eigenvalues of diffusion tensor as a marker for white matter tract maturation and integrity found in [7], [8], [9], [10], and many others. Analyzing these water diffusion related parameters across subjects requires specific statistical methods for the group analysis of diffusion imaging data.

In the current literature, three major approaches to the group analysis of diffusion imaging data are region-of-interest (ROI) analysis, voxel based analysis, and fiber tract based analysis [11], [12], [13]. The region-of-interest (ROI) method used in some neuroimaging studies [14], [15] primarily averages diffusion properties in some manually drawn ROIs for each subject and then creates a single statistic per ROI [13]. Some major drawbacks of ROI analysis are the difficulty in identifying meaningful ROIs, particularly the long curved structures common in fiber tracts, the instability of statistical results obtained from ROI analysis, and the partial volume effect in relative large ROIs. The ROI method is also based on a stringent assumption that diffusion properties in all voxels of the same ROI are essentially homogeneous, which is largely false for DTI.

Voxel based analysis is used more commonly than ROI analysis in neuroimaging studies [16], [17], [18], [19]. The first step involves fitting a statistical model to the smoothed and registered diffusion properties imaging data from multiple subjects at each voxel to generate a parametric map of test statistics (or *p*-values). The second step is to correct for the multiple comparisons across the many voxels of the imaging volume [20], [21], [22]. It suffers from issues of alignment quality and the arbitrary choice of smoothing extent [23], [20], [11], [24]. As shown in Jones et al. (2005), the final results of voxel-based analysis can strongly depend on the amount of smoothing in the smoothed diffusion imaging data.

With the drawbacks mentioned of ROI and voxel based analysis, there is a growing interest in the DTI literature in developing fiber tract based analysis of diffusion properties, such as eigenvalues and fractional anisotropy (FA) values [11], [12], [25], [26]. For instance, Smith and his coauthors develop a *tract-based spatial statistics* framework for constructing local diffusion properties along the white matter skeleton followed by performing pointwise hypothesis tests on the skeleton [11]. Yushkevich and coauthors propose a model-based framework for the statistical analysis of diffusion properties on the medial manifolds of fiber tracts followed by testing pointwise hypotheses on the medial manifolds [25]. Similar to the work proposed in this paper, Goodlett and coauthors propose to use a functional data analysis method to compare a univariate diffusion property, such as fractional anisotropy, across two (or more) populations for a single hypothesis test per tract [26]. Their method is limited to a univariate diffusion property and cannot control for other covariates of interest, such as age, gender and behavioral variables. Moreover, the permutation test and the Hotelling *T*^{2} statistic used in [26] may be statistically over-sensitive, because only the data in the reduced principal component analysis space are permuted. Statistically, such a procedure, which ignores substantial noise in the original data, can lead to misleading results.

What these three methods do not account for is the comparison of fiber bundle diffusion properties across groups and the development of fiber bundle diffusion properties along time, while controlling for other covariates of interest, such as gender [16], [14], [11], [17], [18], [19]. Making these comparisons requires a regression modeling framework for the analysis of fiber bundle diffusion properties and a set of covariates of interest, such as age, diagnostic status and gender. This paper presents a functional regression analysis of DTI tract statistics, called FRATS, for modeling the relationship between fiber bundle diffusion properties and covariates of interest. Compared with [26] and other existing literature, there are four methodological contributions in this paper. First, the local polynomial kernel method is used to regularize multiple diffusion properties along individual fiber bundles. Second, a functional linear model is developed to characterize the association between fiber bundle diffusion properties and any covariate of interest. Third, a global test statistic is proposed for testing hypothesis of interest. Fourth, a resampling method is developed for estimating the *p*-value of the global test statistic. The proposed methodology is used to characterize the development of five diffusion properties along the splenium and right internal capsule tracts in a clinical study of neurodevelopment in the University of North Carolina at Chapel Hill.

This paper introduces a set of statistical tools for the group analysis of fiber bundle diffusion properties while controlling for a set of covariates of interest. Comparing diffusion properties in populations of DTIs requires imaging methods for establishing correspondence among regions of anatomy in DTIs. For this purpose, we use the DTI atlas building methods proposed in [26], which are included here for completion. The focus of this paper is to develop the functional regression model pipeline, called FRATS, for assessing the association between fiber bundle diffusion properties and covariates of interest.

Before applying FRATS, a quality control pipeline is applied to each of the subjects' diffusion weighted images (DWIs) to check if the DWIs contain large slice brightness artifacts, intra-gradient Venetian blind artifacts, and motion artifacts using DTIPrep (www.ia.unc.edu/dev). Skull stripping is then performed automatically via a brain mask computed from both the non-diffusion weighted baseline image and isotropic diffusion weighted image using a free software itkEMS (www.ia.unc.edu/dev/download/itkems). Masking out the non-brain tissues stabilizes the subsequent DTI estimation and registration step. Diffusion tensors are estimated for each subject from a series of cleaned diffusion weighted images using a weighted least square tensor estimation [27], [2].

After DTI estimation, registration occurs and its accuracy is crucial for identifying meaningful group differences. We use a nonlinear fluid deformation based high-dimensional, unbiased atlas computation method to carry out a large deformation non-linear registration [28]. The atlas building procedure is initialized by an affine registration and followed by a nonlinear registration of a feature image which is sensitive to the geometry of white matter. Further reference of the DTI atlas building procedure has been described in [26]. Then, we warp each of the tensor images into the unbiased space and reorient them using the finite strain approximation proposed by [29] to get the aligned DTIs. After averaging all the registered DTIs, a study specific unbiased DTI atlas is created. The tensor warping and averaging are performed using a novel Log-Euclidean geometry [30].

Major fiber bundles are tracked in the atlas space within 3D Slicer (www.slicer.org). First, ROIs are drawn with the Editor module and then fibers are tracked using a stream line method with the Tractography module. Usually, the fibers directly generated in this way contain some fibers not belonging to the bundle that we are interested in due to either the noise in DTI data or not proper ROI delineation. We use FiberViewer to visualize and clean the fibers. Tractography of atlas fiber bundles are more reliable than individual ones because of the improved signal-to-noise ratio of the atlas.

With the fiber bundles in atlas space, each subject's DTI data is transformed into the atlas space. When transforming the diffusion properties from native space to the atlas space, tri-linear interpolation is used to interpolate diffusion properties from the native space. Finally, we get a set of individual tracts with corresponding geometry but varying diffusion properties. Within FiberViewer, a cutting plane is first placed perpendicular to the fibers at a location where the fibers are organized neatly to define the origin of arc length. The diffusion properties along fiber tracts are then measured and plotted as functions of arc length. Using FiberViewer, these plots from each subject are generated. This allows us to conduct a group comparison analysis of the along-bundle tensor properties. Except for 3D Slicer, all DTI tools and itkEMS are part of the open source UNC NeuroLib software repository (https://www.ia.unc.edu/dev/).

After spatial normalization of tensor images, we propose to use a functional regression model to analyze diffusion properties along the same fiber bundle from multiple subjects. FRATS consists of four major components: a regularization of individual tracts, a functional linear model, a global test statistic for hypothesis testing, and a resampling method for estimating the *p*-value of the global test statistic. A schematic overview of FRATS is in Figure 1. We describe each of these components in detail below.

A schematic overview of FRATS: a nonparametric model for regularizing individual tracts, a functional linear model, a global test statistic for hypothesis testing, and a resampling method for estimating the *p*-value of the global test statistic.

To regularize individual tracts, we develop a *nonparametric model* for multiple diffusion properties (e.g., FA) along the same fiber bundle as a smooth function of arc length from each subject as follows. For the *i*-th subject, we consider a *m* × 1 vector of diffusion properties, denoted by **y**_{i,j} = (*y*_{ij,1}, , *y*_{ij,m})^{T}, and its associated arc length *s*_{j} for the *j*-th location grid point on the fiber bundle for *j* = 1, ,*L*_{0} and *i* = 1, , *n*, where *L*_{0} and *n* denote the numbers of grid points and subjects, respectively. The nonparametric model is given by

$${\mathbf{y}}_{i,j}={\mathbf{f}}_{i}\left({s}_{j}\right)+{\u03f5}_{i,j},$$

(1)

where **f**_{i}(*s*) = (*f*_{i,1}(*s*), , *f*_{i,m}(*s*))^{T} is an *m* × 1 vector of continuous functions with second-order continuous derivative, *E*[**y**_{i,j}|**f**_{i}(*s*_{j})] = **f**_{i}(*s*_{j}), and Cov[**y**_{i,j}|**f**_{i}(*s*_{j})] = Σ(*s*_{j}). Our aim is to reconstruct the continuous function vector **f**_{i}(*s*) based on observed fiber bundle diffusion properties.

To simultaneously construct all individual function vectors **f**_{i}(*s*), we develop an adaptive *local polynomial kernel smoothing* technique [31], [32], [33], [34], [35], [36]. Specifically, using Taylor's expansion, we can expand **f**_{i}(*s*_{j}) at *s* to obtain

$${\mathbf{f}}_{i}\left({s}_{j}\right)={\mathbf{f}}_{i}\left(s\right)+{\stackrel{.}{\mathbf{f}}}_{i}\left(s\right)({s}_{j}-s)={\mathbf{A}}_{i}{\mathbf{z}}_{j},$$

(2)

where **z**_{j} = (1, *s*_{j} − *s*)^{T} and ${\mathbf{A}}_{i}=\left[{\mathbf{f}}_{i}\left(s\right)\phantom{\rule{thickmathspace}{0ex}}{\stackrel{.}{\mathbf{f}}}_{i}\left(s\right)\right]$ is an *m* × 2 matrix with ${\stackrel{.}{\mathbf{f}}}_{i}\left(s\right)=d{\mathbf{f}}_{i}\left(s\right)\u2215ds$. We develop an algorithm to estimate **A**_{i} as follows.

Step (1.1). Step (1.1) is to construct an initial estimate of **f**_{i}(*s*) for each *i*. Let **a**_{i;k} be the *k*-th row of **A**_{i} and *K*(·) be a kernel function, such as the Gaussian, Epanechnikov, or uniform kernels [31]. Without loss of generality, we assume that the kernel function *K*(·) is always a symmetric probability density function. Throughout the paper, we take *K*(*t*) = 0.75(1 − *t*^{2})**1**(|*t*| ≤ 1), where **1**(·) is an indicator function. For each *k* and a fixed bandwidth *h*_{k}, we estimate **a**_{i;k} by minimizing an objective function given by

$$\sum _{j=1}^{{L}_{0}}{({y}_{ij,k}-{\mathbf{a}}_{i;k}^{T}{\mathbf{z}}_{j})}^{2}{K}_{{h}_{k}}({s}_{j}-s),$$

(3)

where *K*_{hk}(·) = *K*(·/*h*_{k})*h*_{k} is a rescaled kernel function. With some calculation, it can be shown that

$${\widehat{\mathbf{a}}}_{i;k}={\left[\sum _{j=1}^{{L}_{0}}{\mathbf{z}}_{j}{K}_{{h}_{k}}({s}_{j}-s){\mathbf{z}}_{j}^{T}\right]}^{-1}\sum _{j=1}^{{L}_{0}}{K}_{{h}_{k}}({s}_{j}-s){\mathbf{z}}_{j}{y}_{ij,k}.$$

(4)

Let **e**_{1,2} = (1, 0)^{T}. Then, since ${\mathbf{a}}_{i;k}={({f}_{i,k}\left(s\right),{\stackrel{.}{f}}_{i,k}\left(s\right))}^{T}$, *f*_{i,k}(*s*) can be estimated by

$${\widehat{f}}_{i,k}\left(s\right)={\mathbf{e}}_{1,2}^{T}{\widehat{\mathbf{a}}}_{i;k}=\sum _{j=1}^{{L}_{0}}{\stackrel{~}{K}}_{{h}_{k}}^{0}({s}_{j}-s,s){y}_{ij,k},$$

(5)

where ${\stackrel{~}{K}}_{{h}_{k}}^{0}(\cdot ,\cdot )$ are the empirical equivalent kernels [31]. Thus, ${\widehat{\mathbf{f}}}_{i,k}={({\widehat{f}}_{i,k}\left({s}_{1}\right),\cdots ,{\widehat{f}}_{i,k}\left({s}_{{L}_{0}}\right))}^{T}={S}_{i,k}{\mathbf{y}}_{i:,k}$, where **y**_{i:,k} = (*y*_{i1,k}, , *y*_{iL0,k})^{T} and *S*_{i,k} is the smoother matrix for the *k*-th measurement of the *i*-th subject. For each *k*, we pool the data from all *n* subjects and select the optimal bandwidth *h*_{k}, denoted by ${\widehat{h}}_{k,\mathit{opt}}^{\left(1\right)}$, by minimizing the generalized cross-validation score given by

$${\mathrm{GCV}}_{k}\left({h}_{k}\right)={\left(n\right)}^{-1}\sum _{i=1}^{n}\sum _{j=1}^{{L}_{0}}\frac{{[{y}_{ij,k}-{\widehat{f}}_{i,k}\left({s}_{j}\right)]}^{2}}{1-{L}_{0}^{-1}\mathrm{tr}\left({S}_{i,k}\right)}.$$

(6)

Based on the optimal ${\widehat{h}}_{k,\mathit{opt}}^{\left(1\right)}$, we can estimate ${\widehat{f}}_{i,k}\left(s\right)$ for all *i*.

Step (1.2). Step (1.2) is to construct an estimator of the covariance matrix Σ(*s*_{j}) at *s*_{j}. Specifically, we consider the unbiased sample covariance matrix at *s*_{j} given by

$$\widehat{\Sigma}\left({s}_{j}\right)={(n-m)}^{-1}\sum _{i=1}^{n}{[{\mathbf{y}}_{i,j}-{\widehat{\mathbf{f}}}_{i}\left({s}_{j}\right)]}^{\otimes 2},$$

(7)

where ${\widehat{\mathbf{f}}}_{i}\left(s\right)={({\widehat{f}}_{i,1}\left(S\right),\cdots ,{\widehat{f}}_{i,m}\left(s\right))}^{T}$ and **a**^{2} = **aa**^{T} for any vector **a**. It can be shown that $\widehat{\Sigma}\left({s}_{j}\right)$ converges to the true Σ(*s*_{j}) in probability as both *n* and *L*_{0} go to infinity.

Step (1.3). Step (1.3) is to compute an adaptive estimator of **f**_{i}(*s*) for all *i* using the initial results from Steps (1.1) and (1.2). For all *k* and a fixed bandwidth *h*, we estimate **A**_{i} by minimizing an objective function given by

$$\sum _{j=1}^{{L}_{0}}{({\mathbf{y}}_{ij}-{\mathbf{A}}_{i}{\mathbf{z}}_{j})}^{T}\widehat{\Sigma}{\left({s}_{j}\right)}^{-1}({\mathbf{y}}_{ij}-{\mathbf{A}}_{i}{\mathbf{z}}_{j}){K}_{h}({s}_{j}-s).$$

(8)

Let ${\mathbf{Z}}_{j}=\text{block diagonal}({\mathbf{z}}_{j}^{T},{\mathbf{z}}_{j}^{T},\cdots ,{\mathbf{z}}_{j}^{T})$ be an *m* × 2*m* matrix and ${\mathbf{B}}_{i}={({\mathbf{a}}_{i;1}^{T},\cdots ,{\mathbf{a}}_{i;m}^{T})}^{T}$. It can be shown that

$${\widehat{\mathbf{B}}}_{i}={\left[\sum _{j=1}^{{L}_{0}}{K}_{h}({s}_{j}-s){\mathbf{Z}}_{j}^{T}\widehat{\Sigma}{\left({s}_{j}\right)}^{-1}{\mathbf{Z}}_{j}\right]}^{-1}\sum _{j=1}^{{L}_{0}}{K}_{h}({s}_{j}-s){\mathbf{Z}}_{j}^{T}\widehat{\Sigma}{\left({s}_{j}\right)}^{-1}{\mathbf{y}}_{ij},$$

(9)

which leads to a new estimator of *f*_{i,k}(*s*), denoted by ${\widehat{f}}_{i,k}{\left(s\right)}^{\mathit{sec}}$ for each *i* and *k*. Let ${\stackrel{~}{S}}_{i,k}$ be the smoother matrix for the *k*-th measurement of the *i*-th subject such that ${\widehat{\mathbf{f}}}_{i,k}^{\mathit{sec}}={({\widehat{f}}_{i,k}{\left({s}_{1}\right)}^{\mathit{sec}},\cdots ,{\widehat{f}}_{i,k}{\left({s}_{{L}_{0}}\right)}^{\mathit{sec}})}^{T}={\stackrel{~}{S}}_{i,k}{\mathbf{y}}_{i:,k}$. We pool the data from all *n* subjects and *m* measurements and select the optimal bandwidth *h*, denoted by ${\widehat{h}}_{\mathit{opt}}$ by minimizing the generalized cross-validation score, denoted by GCV(*h*), given by

$$\sum _{i=1}^{n}\sum _{j=1}^{{L}_{0}}\frac{{[{\mathbf{y}}_{ij}-{\widehat{\mathbf{f}}}_{i}{\left({s}_{j}\right)}^{\mathit{sec}}]}^{T}\widehat{\Sigma}{\left({s}_{j}\right)}^{-1}[{\mathbf{y}}_{ij}-{\widehat{\mathbf{f}}}_{i}{\left({s}_{j}\right)}^{\mathit{sec}}]}{n-n{\left(m{L}_{0}\right)}^{-1}\sum _{k=1}^{m}\mathrm{tr}\left({\stackrel{~}{S}}_{i,k}\right)},$$

(10)

where ${\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}={({\widehat{f}}_{i,1}{\left(s\right)}^{\mathit{sec}},\cdots ,{\widehat{f}}_{i,m}{\left(s\right)}^{\mathit{sec}})}^{T}$. Based on the optimal ${\widehat{h}}_{\mathit{opt}}$, we can estimate ${\widehat{f}}_{i,k}{\left(s\right)}^{\mathit{sec}}$ for all *i* and *k*. Similar to the arguments in [35], it can be shown that when Σ(*s*) varies across *s*, ${\widehat{f}}_{i,k}{\left(s\right)}^{\mathit{sec}}$ based on the optimal ${\widehat{h}}_{\mathit{opt}}$ is more accurate than ${\widehat{f}}_{i,k}\left(s\right)$ obtained from Step (1.1).

Step (1.4). Step (1.4) is to estimate the mean function **f**(*s*) and the covariance function Γ(*s*, *t*) of **f**_{i}(*s*). Specifically, following [36], [34], we can estimate **f**(*s*) and Γ(*s*, *t*) by using their empirical counterparts of the estimated ${\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}$ as

$$\widehat{\mathbf{f}}\left(s\right)={n}^{-1}\sum _{i=1}^{n}{\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}\phantom{\rule{1em}{0ex}}\text{and}$$

(11)

$$\widehat{\Gamma}(s,t)={n}^{-1}\sum _{i=1}^{n}[{\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}-\widehat{\mathbf{f}}\left(s\right)]{[{\widehat{\mathbf{f}}}_{i}{\left(t\right)}^{\mathit{sec}}-\widehat{\mathbf{f}}\left(t\right)]}^{T}$$

(12)

The diagonal of Γ(*s*, *s*) reflects the variance of **f**_{i}(*s*) at the location *s*.

We develop a *functional linear model* to characterize the relationship between all diffusion properties along fiber tracts and a set of covariates of interest, such as age, group, and gender. We assume that

$${\mathbf{f}}_{i}\left(s\right)=B\left(s\right){\mathbf{x}}_{i}+{\eta}_{i}\left(s\right),\phantom{\rule{1em}{0ex}}i=1,\cdots ,n,$$

(13)

where *B*(*s*) is a *m* × *p* matrix of functions of *s*, **x**_{i} is a *p* × 1 vector of covariates of interest, and η_{i}(*s*) satisfies *E*[η_{i}(*s*)|**x**_{i}] = 0 and Cov[η_{i}(*s*), η_{i}(*t*)|**x**_{i}] = Γ_{η}(*s*, *t*). *B*(*s*) characterizes the association between fiber bundle diffusion properties and the covariates of interest **x**_{i}.

As an illustration, in our clinical study on early brain development, we are interested in studying the evolution of the three eigenvalues λ_{i} of diffusion tensor (λ_{1} ≥ λ_{2} ≥ λ_{3}) along two selected fiber tracts in 128 healthy pediatric subjects (Figs. 2 (a) and 5 (a)). See clinical data and results section for details. We consider a functional linear model of three eigenvalues along a specific tract as follows:

$$\begin{array}{cc}\hfill {f}_{i,1}\left(s\right)& =({\beta}_{11}\left(s\right),{\beta}_{12}\left(s\right),{\beta}_{13}\left(s\right),{\beta}_{14}\left(s\right)){\mathbf{x}}_{i}+{\eta}_{i1}\left(s\right),\hfill \\ \hfill {f}_{i,2}\left(s\right)& =({\beta}_{21}\left(s\right),{\beta}_{22}\left(s\right),{\beta}_{23}\left(s\right),{\beta}_{24}\left(s\right)){\mathbf{x}}_{i}+{\eta}_{i2}\left(s\right),\hfill \\ \hfill {f}_{i,3}\left(s\right)& =({\beta}_{31}\left(s\right),{\beta}_{32}\left(s\right),{\beta}_{33}\left(s\right),{\beta}_{34}\left(s\right)){\mathbf{x}}_{i}+{\eta}_{i3}\left(s\right),\hfill \end{array}$$

(14)

where *f*_{i,k}(·) equals the smoothed λ_{k} curve from the *i*-th subject for *k* = 1, 2, 3, and **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T}, in which *g*_{i} and Gage_{i} denote gender and gestational age, respectively. Moreover, gestational age is the age of a newborn infant, both gestational age and age are continuous covariates, and gender is a discrete covariate. In this case, *m* = 3, *B*(*s*) = (β_{jk}(*s*)) is a 3 × 4 matrix, and **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T}.

Splenium tract and diffusion properties along the splenium tract: (a) the splenium tract extracted from the tensor atlas with color representing mean FA value; (b) FA; (c) MD; (d) λ_{1}; (e) λ_{2}; (f) λ_{3}. The diffusion properties in **...**

Results from the analysis of FA and MD on the splenium tract: the −log_{10}(*p*) values of test statistics *S*_{n}(*s*_{j}) for testing gender effect in panel (a), gestational age effect in panel (b), and age effect in panel (c) on FA; the −log_{10}(*p*) **...**

We develop an estimation algorithm to estimate *B*(*s*) and Γ_{η}(*s*, *t*) as follows.

Step (2.1). Step (2.1) is to estimate *B*(*s*). Let *B*_{k}(*s*) be the *k*th row of *B*(*s*). Then, based on the estimated function vectors ${\widehat{\mathbf{f}}}_{i}\left(s\right)$, we calculate the least-squares estimator of *B*(*s*), denoted by $\widehat{B}\left(s\right)$, by minimizing an objective function given by

$$\sum _{i=1}^{n}{[{\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}-B\left(s\right){\mathbf{x}}_{i}]}^{T}[{\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}-B\left(s\right){\mathbf{x}}_{i}].$$

(15)

Specifically, the least-squares estimator of *B*_{k}(*s*), denoted by ${\widehat{B}}_{k}\left(s\right)$, is given by

$${\widehat{B}}_{k}{\left(s\right)}^{T}={\left(\sum _{i=1}^{n}{\mathbf{x}}_{i}^{\otimes 2}\right)}^{-1}\sum _{i=1}^{n}{\mathbf{x}}_{i}{\widehat{f}}_{i,k}{\left(s\right)}^{\mathit{sec}}$$

(16)

for *k* = 1, , *m*.

Step (2.2). Step (2.2) is to estimate Γ_{η}(*s*, *t*). Let ${\widehat{\eta}}_{i}\left(s\right)={\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}-\widehat{B}\left(s\right){\mathbf{x}}_{i}$. Then, the covariance matrix Γ_{η}(*s*, *t*) can be estimated by

$${\widehat{\Gamma}}_{\eta}(s,t)={(n-m)}^{-1}\sum _{i=1}^{n}{\widehat{\eta}}_{i}\left(s\right){\widehat{\eta}}_{i}{\left(t\right)}^{T}.$$

(17)

The covariance matrix Γ_{η}(*s*, *t*) characterizes the variation of η_{i}(*s*), which is different from Γ(*s*, *t*).

To carry out statistical inference on *B*(*s*), we need some asymptotic results for $\widehat{B}\left(s\right)$. Let vec(*B*(*s*)) = (*B*_{1}(*s*)^{T}, ,*B*_{m}(*s*)^{T})^{T}. Theoretically, following the arguments in [36], it can be shown that under some regularity conditions,

$$\sqrt{n}[\mathrm{vec}\left(\widehat{B}\left(s\right)\right)-\mathrm{vec}\left(B\left(s\right)\right)+O\left({\widehat{h}}_{\mathit{opt}}^{2}\right)]{\to}^{L}\mathrm{GP}(0,{\Gamma}_{B}),$$

(18)

where →^{L} denotes the convergence in distribution and GP(**0**, Γ_{B}) is a Gaussian process with mean zero and covariance function ${\Gamma}_{B}(s,t)={\Gamma}_{\eta}(s,t)\otimes {\Omega}_{X}^{-1}$, in which denotes the Kronecker product of two matrices and ${\Omega}_{X}={\mathrm{lim}}_{n\to \infty}{n}^{-1}{\sum}_{i=1}^{n}{\mathbf{x}}_{i}^{\otimes 2}$. The asymptotic result in (18) ensures the validity of constructing the global test statistic detailed below.

We develop a *global test statistic* to test linear hypotheses of *B*(*s*) in order to answer various scientific questions involving a comparison of fiber bundle diffusion properties along fiber bundles across two (or more) diagnostic groups and the development of fiber bundle diffusion properties along time. We can formulate these questions as linear hypotheses of *B*(*s*) as follows:

$${H}_{0}:\mathbf{C}\mathrm{vec}\left(B\left(s\right)\right)={\mathbf{b}}_{0}\left(s\right)\phantom{\rule{1em}{0ex}}\text{for all}\phantom{\rule{1em}{0ex}}s\phantom{\rule{1em}{0ex}}\mathrm{vs}.$$

(19)

$${H}_{1}:\mathbf{C}\mathrm{vec}\left(B\left(s\right)\right)\ne {\mathbf{b}}_{0}\left(s\right),$$

(20)

where **C** is a *r* × *mp* matrix of full row rank and **b**_{0}(*s*) is a given *r* × 1 vector of functions.

As an illustration, we are interested in testing the age effect on the evolution of the three eigenvalues of the diffusion tenors along the two selected fiber tracts in 128 healthy pediatric subjects in our clinical study on early brain development. Statistically, for model (1), testing the age effect can be formularized as follows:

$${H}_{0}:{\beta}_{14}\left(s\right)={\beta}_{24}\left(s\right)={\beta}_{34}\left(s\right)=0\phantom{\rule{1em}{0ex}}\text{for all}\phantom{\rule{1em}{0ex}}s\phantom{\rule{1em}{0ex}}\mathrm{vs}.$$

$${H}_{1}:{\beta}_{14}\left(s\right)\ne 0\phantom{\rule{1em}{0ex}}\text{or}\phantom{\rule{1em}{0ex}}{\beta}_{24}\left(s\right)\ne 0\phantom{\rule{1em}{0ex}}\text{or}\phantom{\rule{1em}{0ex}}{\beta}_{34}\left(s\right)\ne 0.$$

In this case, we have

$$\mathbf{C}=\left(\begin{array}{cccccccccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right)$$

and **b**_{0}(*s*) (0, 0, 0)^{T} for all *s*. The use of multiple diffusion properties in model (13) allows us to compare different functions in *B*(*s*) associated with different diffusion properties. For instance, suppose that we want to test

$${H}_{0}:{\beta}_{14}\left(s\right)={\beta}_{24}\left(s\right)\phantom{\rule{1em}{0ex}}\text{for all}\phantom{\rule{1em}{0ex}}s\phantom{\rule{1em}{0ex}}\text{vs.}\phantom{\rule{1em}{0ex}}{H}_{1}{\beta}_{14}\left(s\right)\ne {\beta}_{24}\left(s\right).$$

In this case, we have

$$\mathbf{C}=\left(\begin{array}{cccccccccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right)$$

and **b**_{0}(*s*) (0) for all *s*.

We test the null hypothesis *H*_{0} : **C**vec(*B*(*s*)) = **b**_{0}(*s*) using a global test statistic *S _{n}* defined by

$${S}_{n}=n{\int}_{0}^{{F}_{0}}\mathbf{d}{\left(s\right)}^{T}{\left[\mathbf{C}({\widehat{\Gamma}}_{\eta}(s,s)\otimes {\widehat{\Omega}}_{X}^{-1}){\mathbf{C}}^{T}\right]}^{-1}\mathbf{d}\left(s\right)ds,$$

(21)

where ${\widehat{\Omega}}_{X}={n}^{-1}{\sum}_{i=1}^{n}{\mathbf{x}}_{i}^{\otimes 2},\mathbf{d}\left(s\right)=\mathbf{C}\mathrm{vec}\left(\widehat{B}\left(s\right)\right)-{\mathbf{b}}_{0}\left(s\right)$, and *F*_{0} is the whole arc length of a specific fiber bundle. In order to use *S _{n}* as a test statistic, we need an asymptotic result. Specifically, similar to the arguments in [36], we can show that under some conditions and

In addition, at a given grid point *s _{j}* on a specific tract, we can also test the local null hypothesis

$${S}_{n}\left({s}_{j}\right)=n\mathbf{d}{\left({s}_{j}\right)}^{T}{\left[\mathbf{C}({\widehat{\Gamma}}_{\eta}({s}_{j},{s}_{j})\otimes {\widehat{\Omega}}_{X}^{-1}){\mathbf{C}}^{T}\right]}^{-1}\mathbf{d}\left({s}_{j}\right).$$

(22)

Under some conditions and *H*_{0}(*s _{j}*), $\sqrt{n}\mathbf{d}\left({s}_{j}\right)$ and

We develop a *resampling method* (or wild bootstrap method) to approximate the *p*-value of *S _{n}* [37], [38]. The resampling method has four key steps as follows.

Step (3.1): Fit the functional linear model ${\widehat{\mathbf{f}}}_{i}\left(s\right)=B\left(s\right){\mathbf{x}}_{i}+{\eta}_{i}\left(s\right)$, *i* = 1, , *n*, under the null hypothesis *H*_{0}, which yields ${\widehat{B}}^{\ast}\left(s\right)$ and ${\widehat{\eta}}_{i}^{\ast}\left(s\right)={\widehat{\mathbf{f}}}_{i}\left(s\right)-{\widehat{B}}^{\ast}\left(s\right){\mathbf{x}}_{i}$.

Step (3.2): Generate a random sample ${\tau}_{i}^{\left(g\right)}$ from a *N*(0, 1) generator for *i* = 1, , *n* and then construct ${\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\left(g\right)}={\widehat{B}}^{\ast}\left(s\right){\mathbf{x}}_{i}+{\tau}_{i}^{\left(g\right)}{\widehat{\eta}}_{i}^{\ast}\left(s\right)$. Then, based on ${\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\left(g\right)}$, we calculate

$${\widehat{B}}_{k}{\left(s\right)}^{\left(g\right)T}={\left(\sum _{i=1}^{n}{\mathbf{x}}_{i}^{\otimes 2}\right)}^{-1}\sum _{i=1}^{n}{\mathbf{x}}_{i}{\widehat{\mathbf{f}}}_{i,k}{\left(s\right)}^{\left(g\right)}$$

(23)

for *k* = 1, , *m*, where ${\widehat{B}}_{k}{\left(s\right)}^{\left(g\right)T}$ and ${\widehat{\mathbf{f}}}_{i,k}{\left(s\right)}^{\left(g\right)}$ are, respectively, the *k*th row of $\widehat{B}{\left(s\right)}^{\left(g\right)}$ and ${\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\left(g\right)}$. Finally, let $\mathbf{d}{\left(s\right)}^{\left(g\right)}=\mathbf{C}\mathrm{vec}\left(\widehat{B}{\left(s\right)}^{\left(g\right)}\right)-{\mathbf{b}}_{0}\left(s\right)$, we compute ${S}_{n}^{\left(g\right)}=n{\int}_{0}^{{F}_{0}}\mathbf{d}{\left(s\right)}^{\left(g\right)T}{\left[\mathbf{C}({\widehat{\Gamma}}_{\eta}(s,s)\otimes {\widehat{\Omega}}_{X}^{-1}){\mathbf{C}}^{T}\right]}^{-1}\mathbf{d}{\left(s\right)}^{\left(g\right)}ds$ and ${S}_{n}{\left({s}_{j}\right)}^{\left(g\right)}=n\mathbf{d}{\left({s}_{j}\right)}^{\left(g\right)T}{\left[\mathbf{C}({\widehat{\Gamma}}_{\eta}({s}_{j},{s}_{j})\otimes {\widehat{\Omega}}_{X}^{-1}){\mathbf{C}}^{T}\right]}^{-1}\mathbf{d}{\left({s}_{j}\right)}^{\left(g\right)}$ for *j* = 1, , *L*_{0}.

Step (3.3): Repeat Step (3.2) *G* times to obtain $\{{S}_{n,\mathrm{max}}^{\left(g\right)}={\mathrm{max}}_{1\le j\le {L}_{o}}{S}_{n}{\left({s}_{j}\right)}^{\left(g\right)}:g=1,\cdots ,G\}$ and calculate

$$p\left({s}_{j}\right)={G}^{-1}\sum _{g=1}^{G}1({S}_{n,\mathrm{max}}^{\left(g\right)}\ge {S}_{n}\left({s}_{j}\right))$$

for each *s _{j}*. The

Step (3.4): Repeat Step (3.2) *G* times to obtain $\{{S}_{n}^{\left(g\right)}:g=1,\cdots ,G\}$ and calculate

$$p={G}^{-1}\sum _{g=1}^{G}1({S}_{n}^{\left(g\right)}\ge {S}_{n}).$$

If *p* is smaller than a pre-specified value α, say 0.05, then we reject the null hypothesis *H*_{0}.

We note several advantages of using the resampling method in the above test procedure. Computationally, the above procedure only requires the repeated calculation of ${\widehat{B}}_{k}^{\left(g\right)}$ and ${S}_{n}^{\left(g\right)}$. Thus, because computing ${\widehat{B}}_{k}^{\left(g\right)}$ is very fast, the proposed test procedure is computationally efficient. The proposed resampling method also performs better than parametric bootstrap. Specifically, the parametric bootstrap requires parametric assumptions for η_{i}(*s*). Moreover, for functional linear model (13), permutation method is not directly applicable without very strong assumption such as complete exchangeability.

We conducted a set of Monte Carlo simulations to evaluate the Type I and II error rates of the global test statistic *S*_{n} based on the resampling method. We simulated FA and MD measures along the right internal capsule tract (Fig. 7(a)) obtained from our clinical data according to

$$\begin{array}{c}\hfill {({\mathrm{FA}}_{i}\left({s}_{j}\right),{\mathrm{MD}}_{i}\left({s}_{j}\right))}^{T}={({f}_{i,1}\left({s}_{j}\right),{f}_{i,2}\left({s}_{j}\right))}^{T}+{\u03f5}_{i,j},\hfill \\ \hfill {f}_{i,1}\left(s\right)=({\beta}_{11}\left(s\right),{\beta}_{12}\left(s\right),{\beta}_{13}\left(s\right),{\beta}_{14}\left(s\right)){\mathbf{x}}_{i}+{\eta}_{i1}\left(s\right),\hfill \\ \hfill {f}_{i,2}\left(s\right)=({\beta}_{21}\left(s\right),{\beta}_{22}\left(s\right),{\beta}_{23}\left(s\right),{\beta}_{24}\left(s\right)){\mathbf{x}}_{i}+{\eta}_{i2}\left(s\right),\hfill \end{array}$$

where **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T}, η_{i}(*s*) = (η_{i1}(*s*), η_{i2}(*s*))^{T} is a 2 × 1 vector of Gaussian process with zero mean and covariance matrix Γ_{η}(*s*, *t*) and ϵ_{i,j} is a 2×1 vector of Gaussian random variables with zero mean and covariance matrix Σ(*s _{j}*). To mimic imaging data in practice, we used FA and MD measures along the right internal capsule tract from the 128 subjects in our clinical data to estimate $\widehat{\Sigma}\left({s}_{j}\right)$ of Σ(

Right internal capsule tract and diffusion properties along the tract: (a) the right internal capsule tract extracted from the tensor atlas with color representing mean FA value; (b) FA; (c) MD; (d) λ_{1}; (e) λ_{2}; (f) λ_{3}. The diffusion **...**

We tested the hypotheses *H*_{0} : β_{12}(*s*) = β_{22}(*s*) = 0 for all *s* along the right internal capsule tract against *H*_{1} : β_{12}(*s*) ≠ 0 or β_{22}(*s*) ≠ 0 for at least a *s* on the tract. We assumed *c* = 0 to assess the Type I error rates for the global test statistic, and then we assumed *c* = 0.1, 0.2, 0.3, and 0.4 to examine the Type II error rates for *S _{n}*. In both cases,

$$\mathbf{C}=\left(\begin{array}{cccccccc}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right)\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{\mathbf{b}}_{0}\left(s\right)\equiv \left(\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)$$

for all *s*. We set *n* = 128 and 64. For *n* = 128, we used the same values of age, gestational age, and gender from all 128 subjects in our clinical data. For *n* = 64, we randomly chose 32 males and 32 females from our clinical data and used their values of age, gestational age, and gender to simulate the values of FA and MD along the right internal capsule tract.

For each simulation, the significance levels were set at α = 5% and 1%, and 1,000 replications were used to estimate the rejection rates. For a fixed α, if the Type I rejection rate is smaller than α, then the test is conservative, whereas if the Type I rejection rate is greater than α, then the test is anticonservative, or liberal.

This study was approved by the Institute Review Board of the University of North Carolina at Chapel Hill. A total of 128 healthy full-term subjects (75M and 53F) during the first year of age were recruited and written informed consent was obtained from their parents before imaging acquisition. The subjects were taken from a larger study designated to investigating early brain development at our institution. A total of 128 neonates with a mean postnatal age of 25±17.9 days (range: 9 to 144 days) were included in this study. Efforts were made to ensure the subjects slept comfortably inside the MR scanner, and thus none of the subjects was sedated during imaging session. All subjects were fed and calmed to sleep on a warm blanket with proper ear protection.

All images were acquired on a 3T Allegra head only MR system (Siemens Medical Inc., Erlangen, Germany) with a maximal gradient strength of 40 mT/m and a maximal slew rate of 400 mT/(m·msec). A single shot EPI DTI sequence (TR/TE=5400/73 msec) with eddy current compensation was used to obtain DTI images. Diffusion gradients with a b-value of 1000 s/mm^{2} were applied in six non-collinear directions, (1,0,1), (−1,0,1), (0,1,1), (0,1,−1), (1,1,0), and (−1,1,0). The reference scan (b=0) was also acquired for the construction of diffusion tensor matrices. Contiguous slices with slice thickness of 2mm were scanned to cover the whole brain. The voxel resolution was isotropic 2mm, and the in-plane field of view was 256mm in both directions. A total of five scans were acquired and averaged to improve the signal-to-noise ratio of the images.

A weighted least square estimation method was used to construct the diffusion tensors [27], [2]. Then the image processing steps in the DTI atlas building were used to process all 128 DTI data sets and compute diffusion properties along all fiber tracts of interest. For the sake of space, we chose two tracts of interest including the splenium of the corpus callosum tract and the right internal capsule tract (Figs. 2(a) and 7(a)) and then computed fractional anistropy (FA), mean diffusivity (MD), and the three eigenvalues of the diffusion tensors, denoted by λ_{1} ≥ λ_{2} ≥ λ_{3}, at each grid point on both tracts for each of the 128 subjects. FA denotes the inhomogeneous extent of local barriers to water diffusion, while MD measures the averaged magnitude of local water diffusion. The three eigenvalues of diffusion tensor may, respectively, reflect the magnitude of water diffusivity along and perpendicular to the long axis of white matter fibers [39].

We applied FRATS to the joint analysis of FA and MD values along the splenium tract as follows. We fitted the functional linear model (13) to these smoothed FA and MD functions from all 128 subjects, in which **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T} and *m* = 2, and then we used equation (16) to estimate the function of regression coefficient vector $\widehat{\beta}\left(s\right)$. Secondly, we constructed the global test statistic *S*_{n} to test the effects of all the three covariates for FA alone, MD alone, and joint FA and MD, respectively, and performed hypothesis testing on the whole splenium tract. The *p*-value of *S _{n}* was approximated using the resampling method with

Similar to the analysis of the splenium tract, we also applied FRAST to the analysis of FA and MD values along the right internal capsule tract. Then, we fitted the functional linear model (13) to these smoothed FA and MD functions from all 128 subjects, in which **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T} and *m* = 2. We performed formal hypothesis testings in order to examine the effects of gender, gestational age, and age on FA and MD values along the right internal capsule tract. We further analyzed the three eigenvalues of diffusion tensors along the right internal capsule tract.

Overall, the rejection rates for *S*_{n} based on the resampling method were accurate for all sample sizes (*n* = 64, or 128) at both significance levels (α = 0.01 or 0.05) (Figs. 3(a) and (b)). Consistent with our expectations, the statistical power for rejecting the null hypothesis increased with the sample size.

We used the adaptive local polynomial kernel smoothing technique to simultaneously smooth FA and MD function vectors for each of 128 subjects, which leads to the smoothed curves (Figs. 4(a) and (b)). We computed the sample covariance matrix and obtained correlation coefficient between FA and MD at each grid point. We observed negative correlations at most grid points (Fig. 4 (c)). Based on the smoothed FA and MD functions, we estimated their mean functions and their covariance functions using equations (11) and (12) (Figs. 4 (d) and (e)). Then, we fitted the functional linear model (13) to the smoothed FA and MD functions from all 128 subjects, in which **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T} and *m* = 2. The elements of $\widehat{\beta}\left(s\right)$ corresponding to FA are shown in Fig. 4(f).

Results from the analysis of FA and MD on the splenium tract: reconstructed curves ${\widehat{\mathbf{f}}}_{i}{\left(s\right)}^{\mathit{sec}}$ for FA in panel (a) and MD in panel (b); (c) estimated correlation between FA and MD along the tract; estimated covariance matrices $\widehat{\Gamma}(s,t)$ for FA in **...**

From our model, we performed formal hypothesis testings in order to examine the effects of gender, gestational age, and age on FA and MD values along the splenium tract. We first considered tests based on FA and MD, respectively, and performed hypothesis testing at each grid point along the splenium tract (Figs. 5). For FA alone, no significant effect of any covariate was found, even though the −log_{10}(*p*) value of *S _{n}*(

For FA alone, no significant effects of all covariates were found in Table 1, which agrees with our previous findings based on the local test statistics *S _{n}*(

The *p*-values of *S*_{n} for testing the effects of gender, gestational age (Gage), and age on the splenium and right internal capsule tracts

Along the splenium tract, we observed positive correlations among the three eigenvalues (Fig. 6(d)). Particularly, large correlations between λ_{2} and λ_{3} along the splenium tract indicate the small differences between λ_{2} and λ_{3} on the splenium tract (Figs. 2 (e) and (f)), whereas relatively weak correlations between λ_{1} and λ_{3} indicate large differences between λ_{1} and λ_{3} along the splenium tract (Fig. 2(d)). Then, we fitted the functional linear model (13) to the smoothed three eigenvalues functions from all 128 subjects, in which **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T} and *m* = 3.

Results from the analysis of the three eigenvalues of diffusion tensor on the splenium tract: the −log_{10}(*p*) values of test statistics *S*_{n}(*s*_{j}) for testing gender effect in panel (a), gestational age effect in panel (b), and age effect in panel (c) **...**

We examined association between all covariates of interest and all eigenvalues along the splenium tract by performing formal hypothesis testing at each grid point. No significant gender effect was found for all the three eigenvalues (Fig. 6 (a)), whereas significant gestational age effect was observed at the head and tail of the splenium tract (Fig. 6 (b)). We observed a similar association pattern between gestational age and the three eigenvalues of DT (Fig. 6 (b)). We picked a grid point with a significant *p* value of *S _{n}*(

We simultaneously smoothed the FA and MD function vectors of all 128 subjects and then estimated their mean functions (Figs. 8 (a) and (d)) and their covariance functions (Figs. 8 (b) and (e)). Then, we fitted the functional linear model (13) to the smoothed FA and MD functions from all 128 subjects, in which **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T} and *m* = 2. We considered tests based on FA alone, MD alone, and (FA, MD), respectively, and performed hypothesis testing at each grid point along the right internal capsule tract (Figs. 8(c) and (f)). No significant gender effect was found for all tests (Table 1). For FA alone, significant gestational age and age effects were observed for the right internal capsule tract (Figs. 8 (c) and (f)). We observed negative correlations between FA and MD at most grid points (Fig. 8 (g)). We picked a grid point with the significant *p* value of *S _{n}*(

Results obtained from the analysis of FA and MD on the right internal capsule tract: the estimated mean function and covariance function for FA in panels (a) and (b), respectively; the estimated mean function and covariance function for MD in panels (d) **...**

We also constructed the global test statistic *S _{n}* to assess the covariates of interest based on FA alone, MD alone, and (FA, MD), respectively, and performed hypothesis testing on the right internal capsule tract. The

We also observed positive correlations among the three eigenvalues (Fig. 9(b)). We observed large correlation between λ_{2} and λ_{3} and relatively weak correlation between λ_{1} and λ_{3}. This agrees with the small differences between λ_{2} and λ_{3} along the right internal capsule tract (Figs. 7 (d)–(f)) and the relatively large differences between λ_{1} and λ_{3} in the middle of the right internal capsule tract (Figs. 7(d)–(f)). Then, we fitted the functional linear model (13) to the smoothed three eigenvalues functions from all 128 subjects, in which **x**_{i} = (1, g_{i}, Gage_{i}, age_{i})^{T} and *m* = 3. We performed formal hypothesis testings by examining the association between all covariates of interest and the three eigenvalues along the right internal capsule tract at each grid point. No significant gender effect was found for all the three eigenvalues (Fig. 9 (d)). Gestational age was strongly associated with all the three eigenvalues (Fig. 9 (e)). The scatter plot at a selected grid point shows the decreasing pattern for all eigenvalues with the gestational age. We observed a moderate age effect on all the three eigenvalues (Fig. 9 (f)).

The contributions of our work are two folds. Technically, we have presented a novel functional regression framework, called FRATS, for analyzing multiple fiber bundle diffusion properties with a set of covariates of interest, such as age, diagnostic status and gender, in real applications. FRATS can be used to understand normal brain development, the neural bases of neuropsychiatric disorders, and the joint effects of environmental and genetic factors on white matter fiber bundles. From the application end, FRATS is applied to characterizing the development of diffusion properties on fiber bundles in a clinical study of neurodevelopment. Our approach is able to reveal the complex inhomogeneous spatiotemporal maturation patterns as the apparent changes in FA, MD, and the three eigenvalues of DT. Specifically, our results suggest that white matter maturation patterns are different in different white matter regions.

There are still limitations that need to be addressed. First, we have focused on the analysis of a set of water diffusion related parameters based on diffusion tensor image rather than those based on high angular resolution diffusion image (HARDI), because diffusion tensor image is commonly used for characterizing major white matter fiber bundles. In the future, we will apply FRATS to water diffusion related parameters obtained from HARDI [42], [43], [44]. Second, we have focused on tensor derived univariate measures. It is also interesting to extend FRATS to principal directions and full diffusion tensors on fiber bundles [45], [46], [47], [48], [49]. Third, we have focused on analyzing fiber bundle diffusion properties with a set of covariates. The proposed methodology can be readily extended to more complex fiber structures, such as the medial manifolds of fiber tracts [25]. Fourth, we have focused on cross-sectional studies. Extending FRATS to longitudinal studies and family studies requires further research [50].

We have developed FRATS for statistically analyzing multiple diffusion properties along fiber bundle and assessing their association with a set of covariates of interest in real applications. FRATS integrates four statistical tools for formally testing hypothesis of interest and carrying out statistical inference in real applications. The proposed methodology is demonstrated in a clinical study of neurodevelopment. In this study, significant age and gestational age effects on multiple diffusion properties were examined and localized in two representative tracts. We expect that this novel statistical tool will lead to new findings in our clinical applications.

The authors would like to thank the editor, an associate editor and three referees for valuable suggestions.

This work was supported in part by NSF grant BCS-08-26844 and NIH grants UL1-RR025747-01, MH086633, and AG033387 to Dr. Zhu, NIH grants MH064065, HD053000, and MH070890 to Dr. Gilmore, NIH grant R01NS055754 to Dr. Lin, Lilly Research Laboratories, the UNC NDRC HD 03110, and NIH Roadmap Grant U54 EB005149-01, NAMIC to Dr. Styner.

[1] Basser PJ, Mattiello J, LeBihan D. Mr diffusion tensor spectroscopy and imaging. Biophysical Journal. 1994;66:259–267. [PubMed]

[2] Basser PJ, Mattiello J, LeBihan D. Estimation of the effective self- diffusion tensor from the nmr spin echo. Journal of Magnetic Resonance Ser. B. 1994;103:247–254. [PubMed]

[3] Pierpaoli C, Basser PJ. Toward a quantitative assessment of diffusion anisotropy. Magn Reson Med. 1996;36:893–906. [PubMed]

[4] Hasan KM, Basser PJ, Parker DL, Alexander AL. Analytical computation of the eigenvalues and eigenvectors in dt-mri. J Magn Reson. 2001;152:41–47. [PubMed]

[5] Hasan KM, Narayana PA. Computation of the fractional anisotropy and mean diffusivity maps without tensor decoding and diagonalization: theoretical analysis and validation. Magn Reson Med. 2003;50:589–598. [PubMed]

[6] Zhu HT, Xu D, Amir R, Hao X, Zhang H, Alayar K, Ravi B, Peterson B. A statistical framework for the classification of tensor morphologies in diffusion tensor images. Magnetic Resonance Imaging. 2006;24:569–582. [PMC free article] [PubMed]

[7] Moseley M. Diffusion tensor imaging and aging-a review. NMR Biomed. 2002;15:553–560. [PubMed]

[8] Mukherjee P, McKinstry RC. Diffusion tensor imaging and tractography of human brain development. Neuroimaging Clin N Am. 2006;16:19–43. [PubMed]

[9] Cascio CJ, Gerig G, Piven J. Diffusion tensor imaging: application to the study of the developing brain. J Am Acad Child Adolesc Psychiatry. 2007;46:213–223. [PubMed]

[10] Rollins NK. Clinical applications of diffusion tensor imaging and tractography in children. Pediatr Radiol. 2007;37:769–780. [PubMed]

[11] Smith SM, Jenkinson M, Johansen-Berg H, Rueckert D, Nichols TE, Mackay CE, Watkins KE, Ciccarelli O, Cader MZ, Matthews PM, Behrens TE. Tractbased spatial statistics: voxelwise analysis of multi-subject diffusion data. NeuroImage. 2006;31:1487–1505. [PubMed]

[12] O'Donell LJ, Westin CF, Golby AJ. Tract-based morphometry for white matter group analysis. NeuroImage. 2009;45:832–844. [PMC free article] [PubMed]

[13] Snook L, Plewes C, Beaulieu C. Voxel based versus region of interest analysis in diffusion tensor imaging of neurodevelopment. NeuroImage. 2007;34:243–252. [PubMed]

[14] Bonekamp D, Nagae LM, Degaonkar M, Matson M, Abdalla WMA, Barker PB, Mori S, Horsk A. Diffusion tensor imaging in children and adolescents: reproducibility, hemispheric, and age-related differences. NeuroImage. 2008;34:733–742. [PMC free article] [PubMed]

[15] Gilmore JH, Smith LC, Wolfe HM, Hertzberg BS, Smith JK, Chescheir NC, Evans DD, Kang C, Hamer RM, Lin W, Gerig G. Prenatal mild ventriculomegaly predicts abnormal development of the neonatal brain. Biol Psychiatry. 2008;64:1069–1076. [PMC free article] [PubMed]

[16] Chen YS, An HY, Zhu HT, Stone T, Smith JK, Hall C, Bullitt E, Shen DG, Lin WL. White matter abnormalities revealed by diffusion tensor imaging in non-demented and demented hiv+ patients. NeuroImage. 2009;47:1154–1162. [PMC free article] [PubMed]

[17] Focke NK, Yogarajah M, Bonelli SB, Bartlett PA, Symms MR, Duncan JS. Voxel-based diffusion tensor imaging in patients with mesial temporal lobe epilepsy and hippocampal sclerosis. NeuroImage. 2008;40:728–737. [PubMed]

[18] Cmara E, Bodammer N, Rodrguez-Fornells A, Tempelmann C. Age-related water diffusion changes in human brain: A voxel-based approach. NeuroImage. 2007;34:1588–1599. [PubMed]

[19] Snook L, Paulson LA, Roy D, Phillips L, Beaulieu C. Diffusion tensor imaging of neurodevelopment in children and young adults. NeuroImage. 2005;26:1164–1173. [PubMed]

[20] Ashburner J, Friston KJ. Voxel-based morphometry: the methods. NeuroImage. 2000;11:805–821. [PubMed]

[21] Wager TD, Keller MC, Lacey SC, Jonides J. Increased sensitivity in neuroimaging analyses using robust regression. Neuroimage. 2005;26:99–113. [PubMed]

[22] Worsley KJ, Taylor JE, Tomaiuolo F, Lerch J. Unified univariate and multivariate random field theory. NeuroImage. 2004;23:189–195. [PubMed]

[23] Hecke WV, Sijbers J, Backer SD, Poot D, Parizel PM, Leemans A. On the construction of a ground truth framework for evaluating voxel-based diffusion tensor mri analysis methods. NeuroImage. 2009;46:692–707. [PubMed]

[24] Jones DK, Symms MR, Cercignani M, Howard RJ. The effect of filter size on vbm analyses of dt-mri data. NeuroImage. 2005;26:546–554. [PubMed]

[25] Yushkevich PA, Zhang H, Simon TJ, Gee JC. Structure-specific statistical mapping of white matter tracts. Neuroimage. 2008;41:448–461. [PMC free article] [PubMed]

[26] Goodlett CB, Fletcher PT, Gilmore JH, Gerig G. Group analysis of dti fiber tract statistics with application to neurodevelopment. NeuroImage. 2009;45:S133–S142. [PMC free article] [PubMed]

[27] Zhu HT, Zhang HP, Ibrahim JG, Peterson BG. Statistical analysis of diffusion tensors in diffusion-weighted magnetic resonance image data (with discussion) Journal of the American Statistical Association. 2007;102:1085–1102.

[28] Joshi S, Davis B, Jomier M, Gerig G. Unbiased diffeomorphic atlas construction for computational anatomy. Neuroimage. 2004;23:S151–S160. [PubMed]

[29] Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Transactions on Medical Imaging. 2001;20:1131–1139. [PubMed]

[30] Arsigny V, Fillard P, Pennec X, Ayache N. Log-euclidean metrics for fast and simple calculus on diffusion tensors. Magnetic Resonance in Medicine. 2006;56:411–421. [PubMed]

[31] Fan J, Gijbels I. Local Polynomial Modelling and Its Applications. Chapman and Hall; London: 1996.

[32] Wand MP, Jones MC. Kernel Smoothing. Chapman and Hall; London: 1995.

[33] Wu HL, Zhang JT. Nonparametric Regression Methods for Longitudinal Data Analysis. John Wiley & Sons, Inc.; Hoboken, New Jersey: 2006.

[34] Ramsay JO, Silverman BW. Functional Data Analysis. Springer-Verlag; New York: 2005.

[35] Welsh AH, Yee TW. Local regression for vector responses. Journal of Statistical Planning and Inference. 2006;136:3007–3031.

[36] Zhang J, Chen J. Statistical inference for functional data. The Annals of Statistics. 2007;35:1052–1079.

[37] Zhu HT, Ibrahim JG, Tang N, Rowe DB, Hao X, Bansal R, Peterson BS. A statistical analysis of brain morphology using wild bootstrapping. IEEE Trans Med Imaging. 2007;26:954–966. [PMC free article] [PubMed]

[38] Lin DY. An efficient monte carlo approach to assessing statistical significance in genomic studies. Bioinformatics. 2005;6:781–787. [PubMed]

[39] Song SK, Sun SW, Ju WK, Lin SJ, Cross AH, Neufeld AH. Diffusion tensor imaging detects and differentiates axon and myelin degeneration in mouse optic nerve after retinal ischemia. Neuroimage. 2003;20:1714–1722. [PubMed]

[40] Neil JJ, Shiran SI, McKinstry RC, Schefft GL, Snyder AZ, Almli CR, Akbudak E, Aronovitz JA, Miller JP, Lee BC, Conturo TE. Normal brain in human newborns: apparent diffusion coefficient and diffusion anisotropy measured by using diffusion tensor mr imaging. Radiology. 1998;209:57–66. [PubMed]

[41] Zhai G, Lin W, Wilber KP, Gerig G, Gilmore JH. Comparisons of regional white matter diffusion in healthy neonates and adults performed with a 3.0-t head-only mr imaging unit. Radiology. 2003;229:673–681. [PubMed]

[42] Lenglet C, Campbell JSW, Descoteaux M, Haro G, Savadjiev P, Wassermann D, Anwander A, Deriche R, Pike GB, Sapiro G, Siddiqi K, Thompson PM. Mathematical methods for diffusion mri processing. NeuroImage. 2009;45:S111–S122. [PMC free article] [PubMed]

[43] Leow AD, Zhu S, Zhan L, McMahon K, de Zubicaray GI, Meredith M, Wright MJ, Toga AW, Thompson PM. The tensor distribution function. Magnetic Resonance in Medicine. 2009;61:205–214. [PMC free article] [PubMed]

[44] Tuch DS, Reese TG, Wiegell MR, Makris N, Belliveau JW, Wedeen VJ. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magnetic Resonance in Medicine. 2002;48:577–582. [PubMed]

[45] Schwartzman A. PhD thesis. Stanford University; California: Jul, 2006. Random Ellipsoids and False Discovery Rates: Statistics for Diffusion Tensor Imaging Data.

[46] Lepore N, Brun CA, Chou Y, Chiang M, Dutton RA, Hayashi KM, Luders E, Lopez OL, Aizenstein HJ, Toga AW, Becker JT, Thompson PM. Generalized tensor-based morphometry of hiv/aids using multivariate statistics on deformation tensors. IEEE Transactions in Medical Imaging. 2008;27:129–141. [PMC free article] [PubMed]

[47] Schwartzman A, Dougherty RF, Taylor JE. Cross-subject comparison of principal diffusion direction maps. Magn Reson Med. 2005;53:1423–1431. [PubMed]

[48] Zhu HT, Cheng YS, Ibrahim JG, Li YM, Hall C, Lin WL. Intrinsic regression models for positive definitive matrices with applications in diffusion tensor images. Journal of the American Statistical Association. 2009;104:1203–1212. [PMC free article] [PubMed]

[49] Whitcher B, Wisco JJ, Hadjikhani N, Tuch DS. Statistical group comparison of diffusion tensors via multivariate hypothesis testing. Magnetic Resonance in Medicine. 2007;57:1065–1074. [PMC free article] [PubMed]

[50] Fang Y, Wang Y. Testing for familial aggregation of functional traits. Statistics in Medicine. 2010 [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |