图书馆订阅: Guest
Begell Digital Portal Begell 数字图书馆 电子图书 期刊 参考文献及会议录 研究收集
生物医学的图像处理,计算和显示


ISSN 在线: 2162-3511

Archives: Volume 1, 2012 to Volume 2, 2013

生物医学的图像处理,计算和显示

DOI: 10.1615/VisualizImageProcComputatBiomed.2012004964

FAST CONSTRAINED CANONICAL CORRELATION ANALYSIS FOR FMRI

Mingwu Jin,1,2 Rajesh Nandy,3 & Dietmar Cordes2

1Department of Physics, University of Texas at Arlington, Arlington, Texas 76019, USA

2Department of Radiology, School of Medicine, University of Colorado Denver, Aurora, Colorado 80045, USA

3Department of Radiology, School of Medicine, University of Colorado Denver, Aurora, Colorado 80045, USA

Abstract.

Local canonical correlation analysis (CCA) methods can increase the detection power in functional magnetic resonance imaging (fMRI) data analysis. Spatial constraints are usually needed to overcome model overfitting and loss of spatial specificity. The resulting constrained CCA (cCCA) method suffers from a marked increase of computation time, even though a computationally efficient branch-and-bound method (cCCA-BB) was developed that avoids an exhaustive search. In this study, a fast implementation of cCCA using a novel region-growing strategy (cCCA-RG) is proposed. This method can significantly shorten the computational time without compromising accuracy. Using simulated and real fMRI data, the performance (accuracy and speed) between cCCA-RG, cCCA-BB, and the conventional general linear model (GLM) methods (GLM with or without fixed Gaussian smoothing) is compared. Results demonstrate that cCCA-RG achieves similar detection power to cCCA-BB, and both cCCA methods have a higher detection power than any of the GLM methods. The cCCA-RG method is at least 12 times faster than the cCCA-BB method. The proposed fast implementation using cCCA-RG enables a broad application of cCCA for fMRI data analysis.

Keywords:

fMRI, general linear model (GLM), constrained canonical correlation analysis (cCCA), fast implementation, branch-and-bound method, region-growing method


1. Introduction

The general linear model (GLM) is the dominant functional magnetic resonance imaging (fMRI) data analysis method used to determine activations in the brain because of its efficiency in estimation and good detection sensitivity [1]. The estimated parameters and their variances are used to construct various contrast statistics, either t or F, for inference of regional effects. The mass-univariate GLM routinely uses fixed isotropic Gaussian spatial smoothing to suppress noise variations and to control the family-wise error parametrically, based on the theory of random fields [2, 3]. The fixed spatial smoothing causes blurring of the activation, and could potentially eliminate weak activations, such as areas involved in episodic memory functions, if the filter kernel is larger than the activated area. Moreover, the spatial correlation of fMRI noise can diminish the noise suppression effect from Gaussian smoothing [4]. The use of adaptive spatial smoothing may alleviate these problems of isotropic Gaussian smoothing.

Local multivariate methods such as canonical correlation analysis (CCA) [5] and its variants [6–8] have been demonstrated to increase significantly the detection power of fMRI activations due to their adaptive spatial smoothing nature that can boost the signal-to-noise ratio. In conventional CCA [5], neighborhood-specific spatial filtering has been applied on a 3 × 3 group of in-plane voxels and the maximum canonical correlation was assigned to the center voxel of the 3 × 3 region. Though this pioneering method increased the detection sensitivity, it also increased the number of false positives due to more freedom in finding favorable linear combinations with non-active voxel time series, thus leading to a decrease in specificity. To solve this problem, [6] proposed constrained CCA using a linear combination of four spatially adaptive smoothing kernels. A better assignment scheme [7] was proposed in order to improve the spatial specificity of conventional CCA. Recently, we systematically studied the effectiveness of various constraints used in cCCA and found as high as a 20% increase of detection of fMRI activations [8]. However, the detection improves at the expense of increased computation time because the optimal solution was found by either an exhaustive search [8] or a moderately accelerated branch-and-bound method (cCCA-BB) [6, 9].

Another problem with previous CCA methods [5–7] is lack of a general test statistic allowing computation of estimable linear contrasts of interest as in GLM. This problem was addressed by [10] through the connection between the multivariate multiple-regression (MVMR) model and CCA. The weights of fMRI neighboring voxels determined by cCCA were treated as a linear transformation of the dependent variables in the MVMR model. Then, the weights of temporal regressors from a complicated design matrix (i.e., the weighing parameters of the explanatory variables) could be used to test any estimable linear contrast of temporal regressors after one-time cCCA estimation similar to the MVMR model. Such a treatment avoided reparametrization and reestimation for each contrast, where the “reparametrization” means that the original design matrix can be linearly transformed so that one particular transformed regressor represents a linear contrast of several original regressors and is orthogonal to other transformed regressors ([8], Appendix B). This novel test statistic for cCCA [10] is necessary for the fast implementation developed in this study.

In this study, we propose a fast implementation of cCCA methods using the region-growing method (a forward search technique). Using receiver operating characteristic (ROC) techniques [11–13] applied to pseudo-real data [7, 14, 15], we quantitatively compare the sensitivity and specificity of the proposed method with the GLM without and with fixed Gaussian spatial smoothing and the branch-and-bound cCCA method (cCCA-BB). We also apply a nonparametric approach [16] to estimate the family-wise error rate for all methods using resampled resting-state data and show the activation maps from real fMRI data for a simple visual cortex activation paradigm and a more complicated memory paradigm.

2. METHODS

For a group of K local neighboring voxels in fMRI data, the multivariate multiple-regression model can be written as

(1)

where Y = (y1,y2,…,yK) is a n × K matrix containing n temporal observations of K neighboring voxels from fMRI time series, X is a predefined n × p design matrix, B = (β1, β2,…,βK) is a p × K parameter matrix to be estimated, and E = (ε1, ε2,…,εK) is the n × K error matrix. In order to increase detection power of weak activations, local spatial smoothing is usually applied to decrease the noise variance. Let α be the vector containing the spatial smoothing coefficients. Multiplication of both sides of Eq. (1) with α gives

(2)

where β and ε.

In conventional fixed Gaussian smoothing, both Y and α are fixed and treated as known, and only β has to be estimated. To determine local adaptive spatial weights in the formalism of CCA, our goal is to find a subset Y* of Y with given constraints on α (see next section, “Fast Computational Approaches for solving the cCCA problem”) and to estimate vectors and (for α and β) to have minimal fitting errors. The latter can be achieved by either the least-squares (LS) method or by CCA. Both of these solutions are equivalent but differ in scale. However, this scaling difference will cancel out in the calculation of statistics for hypothesis testing [10].

Although conventional CCA applied to a local neighborhood can increase the sensitivity, it will also increase the number of false positives if there is no constraint on the spatial coefficients α enforced. In this article, we apply two important constraints for the spatial weights and focus more on a fast implementation of cCCA. Other elaborate investigations of different constraints for cCCA can be found in [8]. First, the basic requirement is that all of the spatial weights need to be positive to meet the requirements of being a low-pass (smoothing) filter. The local adaptiveness of the cCCA method allows versatile filtering kernels with flexible shape and size within a prespecified region (for example, a 3 × 3 pixel region). Note that if both signs of and are flipped, the canonical correlation will not change. Therefore, the actual constraint used in our methods is that has to have the same sign for all components. Second, a constraint is needed to account for the following artifact: When the center voxel has a weight that is only weakly positive (and thus the voxel is very likely to be non-active) but most neighboring voxels have large spatial weights, this center voxel could be incorrectly assigned a large significance value if activations of these neighboring voxels are significant, due to the multivariate property of CCA. This scenario is very likely to impact (non-active) boundary regions of activations. To reduce the number of false positives near the boundary of activations, an additional constraint for the spatial weight of the center voxel using an empirically determined parameter δ in [0,1] is used according to

(3)

This parameter allows adjusting the tradeoff between sensitivity and specificity of the detection of activation. When δ = 0, there is no constraint on the center voxel. When δ = 1, the center voxel has the largest weight in the neighborhood. Thus, with δ increasing from 0 to 1, the probability that an inactive center voxel is labeled as “active” decreases at the expense of the sensitivity of the method.

For each subset Y* of a prescribed search area Y, CCA has to be applied to find the solutions of Eq. (2) satisfying the constraints. Although an exhaustive search could provide the optimal solution [8], this is not always feasible because the high dimensionality of the problem due to the number of configurations for each voxel and the number of voxels in whole-brain fMRI data leads to a high cost in computation time. For example, a 3 × 3 pixel region requires an exhaustive search of 28 = 256 different pixel configurations. This number increases to 226 if the search is extended to a 3D 3 × 3 × 3 voxel neighborhood. The computations will increase exponentially with the number of voxels defining the local neighborhood. Heuristic methods for significantly reducing computation time without sacrificing accuracy will be indispensible for practical applications of fMRI data analysis.

2.1. Fast Computational Approaches for Solving the cCCA Problem

We use the notation to label the time series of a local neighborhood of K voxels as Y = {y1,…,yK–1,yK} and the center voxel as y1. In other words, each column of Y denotes the time-series data of a voxel measured by fMRI. The center voxel of the neighborhood is of interest because the cCCA test statistic is assigned to it, instead of to the whole neighborhood, to increase specificity of the method. Our goal is to find a submatrix Y* of Y with y1Y* so that the spatial weight vector α of Y* satisfies the constraints and the correlation of Y*α and is maximized. Two methods of finding the submatrix Y* and corresponding vectors α and β are described as follows.

2.1.1. Branch-and-Bound Method (cCCA-BB)

The branch-and-bound approach was originally proposed by [9] to find non-negative spatial weights for CCA. This method uses the monotonicity of the canonical correlation as a function of the degrees of freedom of the spatial functions, i.e., number of components in Y. An implementation of this method to find Y* is as follows:

  1. Try CCA using Y* containing all K column vectors of Y and X.
  2. Try CCA using all possible Y* containing K–1 column vectors of Y (including the center voxel y1) and X.
  3. Try CCA using all possible Y* containing K–2 column vectors of Y (including the center voxel y1) and X.

If at any step the vector α satisfies the positivity constraint (or equivalently all components of α have the same sign) and also the center voxel constraint, then the search is stopped and the optimal submatrix Y* is obtained. If there are several Y* found satisfying the constraints, then the one Y* corresponding to the largest canonical correlation provides the optimal solution.

Such a branch-and-bound method was first adopted by [6] for fMRI data analysis. They composed the filter kernel space of one isotropic center basis filter and several oriented parts (so called “steerable filters”). By weighting them differently, anisotropic low-pass filters in any orientation can be constructed. In order to obtain a smoothing (low-pass filter) effect, only positive weights were allowed when combining the basis filters. The weights of these basis filters were determined in a two-step process: (1) Perform cCCA without the center basis filter to determine the best orientation of a spatial filter; (2) perform cCCA with both the center basis filter and the resulting oriented filter from step (1) to obtain the final correlation coefficient to the center voxel.

Our branch-and-bound approach is similar to what [6] proposed in essence. However, we apply this method differently in several aspects. First, more flexible filter kernels are allowed in our method. The steerable filter kernels used in their work [6] can only generate spatial filters that are symmetric to the center voxel. Second, we explicitly define the center voxel dominance condition as an optional constraint and use a one-step cCCA estimation. Third, we use a classic hemodynamic response function in a complex experimental design and use our developed novel cCCA statistic [10] instead of the correlation coefficient for the computation of contrasts of interest. Our generalized test statistic is

(4)

where c is the contrast vector that allows us to test any linear combination of experimental effects , and DF = n–p–K* specifies the degrees of freedom (DOF) given n observations (i.e., time samples), p nonconstant regressors (i.e., linear equations for β), and K* voxels selected in cCCA.

2.1.2. Novel Local Region-Growing Method (cCCA-RG)

In order to provide a more efficient solution of the cCCA spatial problem, we propose a novel local region-growing method starting from the center voxel y1 to find the optimal subset Y*. The branch-and-bound method and the proposed region-growing method essentially are based on two opposite search strategies. The former prunes the voxel space from the largest set of voxels, while the latter grows it from the smallest set, i.e., the center voxel. The region-growing method is implemented as follows: starting at the center voxel y1, add one voxel in each iteration that gives the highest canonical correlation and also meets the positivity constraint for α and also the center voxel constraint α1 δ max(αc). The iteration stops if no such voxel exists in searching the neighborhood of y1 and outputs the current Y* as the optimal subset.

Both cCCA-BB and cCCA-RG methods yield solutions that could fall in a local maximum because the search schemes are based on a necessary condition for cCCA (i.e., CCA solutions satisfying constraints). For example, a particular configuration of voxels could achieve a global maximum canonical correlation when constrained to have positive spatial weights, but when unconstrained will give both positive and negative spatial weights for a higher correlation. Such a configuration of voxels will be excluded from the search algorithms. Despite the issue of no guarantee to obtain the global optimal cCCA solution, the proposed cCCA-RG method is computationally efficient and can yield substantially better detection power than the GLM, as shown in the “Results” section. Thus, cCCA-RG can serve as a good substitute of the conventional GLM method when additional detection power is desired.

2.2. Imaging Data

Functional MRI (fMRI) was performed at the Brain Imaging Center of the University of Colorado, Denver, in a 3.0 T GE HDx MRI scanner equipped with an eight-channel head coil and parallel imaging technology. Stimulus presentation was done with a rear projection system (AVOTEC, Inc.). Two different paradigms (visual paradigm and memory paradigm) were performed on five healthy adult subjects and fMRI data were collected according to local IRB approval. The pulse sequence to collect fMRI data was echo-planar imaging (EPI) with the following parameters: ASSET = 2, repetition time (TR)/echo time (TE) = 2 s/30 ms, flip angle (FA) = 70°, field of view (FOV) = 22 cm × 22 cm, slice thickness = 4 mm, gap = 1 mm, 25 slices, in-plane matrix 96 × 96. For the visual paradigm we prescribed axial slices (slice orientation perpendicular to the AC-PC line) and collected 150 volumes, whereas for the memory paradigm we prescribed coronal oblique slices perpendicular to the long axis of the hippocampus and collected 288 volumes. The first 5 volumes were discarded to establish signal equilibrium of the imaging sequence. In addition, we acquired a co-planar standard high-resolution T2-weighted anatomical scan (FOV 22 cm, resolution 256 × 256, TR/TE = 3 s/85 ms, NEX 2, slice thickness 4 mm, gap 1 mm) that we used as an underlay of the functional data. Here we briefly describe the paradigms programed in EPRIME (Psychological Software Tools, Inc.). More details can be found in [10, 17].

2.2.1. Visual Paradigm

For each subject we acquired two fMRI data sets. The first data set was collected during resting state where the subject tried to relax and refrain from executing any overt task with eyes closed. The second data set was collected with an event-related paradigm where the subject was looking at a flashing checkerboard (10 Hz flashing frequency, duration 2 s) which alternated with a fixation period of random (jittered) duration (2–10 s, uniformly distributed). The jittered null events are commonly employed in event-related tasks to prevent the subject from guessing when the next task activation will occur (which when guessed correctly would provide unwanted spikes in attention). Furthermore, jittering leads to different linear overlaps of the hemodynamic response of the activation and allows accurate and efficient isolation of the hemodynamic response function for the visual event. During the fixation period a black screen containing a small white cross (about 1 in. in size) at the center was shown and the subject was instructed to focus on this cross.

2.2.2. Memory Paradigm

We collected two fMRI data sets for each subject. The first set contained resting-state data, and the second set was acquired while the subject performed a memory task. The memory paradigm started with a fixation period of 16 s followed by six 89-s-long cycles of “5 s instruction”, “21 s encoding condition”, “5 s instruction”, “11 s control condition”, “5 s instruction”, and “42 s recognition condition.” It ended with another fixation period of 16 s. The short “instruction period” consisted of a single sentence and reminded the subject of the following task to be performed. The “encoding condition” consisted of a series of novel outdoor pictures, where each picture was displayed for 3 s for the subject to memorize it. During the “control condition” the subject saw the letters “Y” or “N” appearing in random order every 100 ms. The subject was instructed to press as fast as possible the right button when Y appears or the left button when N appears. The purpose of the “control” condition was to keep the subject′s attention away from the previously viewed pictures without producing any activations in regions associated with the memory circuit (hippocampal complex, posterior cingulate cortex, precuneus, fusiform gyrus). The “recognition condition” is a recognition task where the subject saw a series of randomly arranged pictures where half of the pictures were novel and the other half of the pictures were identical to the pictures from the previous “encoding period.” Each picture was displayed on the screen for 3 s. The subject was instructed to press the right button if the picture was seen in the previous encoding period and to press the left button if the picture was not seen in the previous encoding session. The four conditions of “instruction,” “encoding,” “recognition,” and “control” are denoted as “I”, “E”, “R”, and “C”, respectively. Due to the complexity of the memory paradigm, all subjects were trained on a computer before fMRI scanning in a quiet room with the memory paradigm using a different set of images.

2.3. Preprocessing

All data were preprocessed in SPM5 using realignment to correct for motion artifacts, slice timing correction to correct for differences in image acquisition time between slices, and high-pass filtering using T = 150 s to remove low-frequency components and signal drifts. The classic two-gamma hemodynamic response function was used to construct the design matrix. In the next section, we give examples for the contrasts “Visual minus Fixation” (denoted as “V-F”) for visual data and “Encoding minus Control” (denoted as “E-C”) for memory data for clarity and omit other contrasts of interest.

2.4. Construction of Simulated and Pseudo-Real Data

To quantitatively determine the performance of different methods, we constructed both simulated and pseudo-real data in the following way:

(5)

In this equation x is the vector representing the time series of a voxel with activation contribution xact and noise contribution xnull. The noise fraction parameter f is a fixed number to adjust the noise level in the data vector x given that xact and xnull have the same or similar power. Now, to construct simulated data we used for xact an ideal hemodynamic response function and for xnull white Gaussian noise. To use more realistic data, we constructed another pseudo-real data from fMRI activation data and resting-state data. To obtain pseudo-real data according to Eq. (5) we used for xact activation fMRI data where the active set was obtained by setting a very high significance threshold so that the activations obtained are almost certainly true activations. The vector xnull was chosen to be Fourier resampled resting-state fMRI data [18, 19]. The advantage of constructing pseudo-real data from real active data and resampled resting-state data is that the spatial and temporal correlations of both the activations and the noise are similar to real data and the location of active and non-active voxels are known, which allows ROC techniques to be applied on the data for quantitative comparison of difference methods. For pseudo-real data, we used two values for f: 0.55 representing the low noise case, and 0.65 representing the high noise case as proposed in our previous publication [10].

2.5. Methods of fMRI Data Analysis

Four methods were investigated in this study: (i) GLM without smoothing, denoted as “GLM-NS,” (ii) Gaussian smoothing followed by GLM, denoted as “GLM-GS,” (iii) cCCA with the branch-and-bound method, denoted as “cCCA-BB,” and (iv) cCCA with the region-growing method, denoted as “cCCA-RG. The search neighborhood for the cCCA methods is an in-plane 3 × 3 pixel area for each voxel. The full width at half maximum (FWHM) of Gaussian smoothing in the GLM was chosen as 2.24 pixels. This number is not only falling in the generally recommended smoothing size (two to three times of the spatial resolution) in fMRI data analysis, but is also equal to the average size of all possible 256 configurations within a 3 × 3 pixel area that includes the center pixel. For the contrasts of V-F of visual data and E-C of memory data, the t statistic was used in GLM-NS and GLM-GS, while the novel cCCA statistic developed in [10] was used in cCCA-BB and cCCA-RG.

3. RESULTS AND DISCUSSION

3.1. Coefficient Estimation for Simulated Data

To define different patterns of spatial activations, 100,000 randomly shaped activations within a 3 × 3 grid of pixels having a size of 2–9 pixels were generated. The center pixel was always assigned to be active. The corresponding time courses for the activated pixels xact were simulated to be linear combinations of the four memory regressors (i.e., I, E, R, and C) used in the memory paradigm with random amplitudes βI, βE, βR, βC, uniformly distributed in [0,1]. Different levels of white Gaussian noise were used for xnull and added to the 3 × 3 grid of pixels. Both xnull and xact were normalized to have unit variance before the mixture. Since the center pixel is always active, we used δ = 0 to enforce maximum adaptive smoothing.

We computed the mean square errors (MSE) between the estimated coefficients of the linear combination and the original ones generating the data (Table 1) for a random noise fraction parameter f uniformly distributed in [0,1]. Both proposed cCCA-BB and cCCA-RG outperform GLM-NS and GLM-GS. The GLM-GS method is inferior to GLM-NS due to the small and irregularly defined activations. The cCCA-BB method performs best and has an improvement of about 20% in estimating the coefficients based on MSE demonstrating the superior estimation power of the proposed cCCA method.

Table 1. Mean square errors (MSE) of estimated coefficients for different methods. The mean square errors (MSE) between the original simulated amplitudes of activation and estimated ones are shown for different methods. Both cCCA methods achieve about 20% less MSE compared to the GLM-NS. The GLM-GS method is worse than GLM-NS due to the small and irregularly defined activation patterns.

  ΔβI2 ΔβE2 ΔβR2 ΔβC2
GLM-NS 0.3273 0.5488 0.7549 0.4102
GLM-GS 0.3454 0.5793 0.7933 0.4335
cCCA-BB 0.2638 0.4460 0.6170 0.3317
cCCA-RG 0.2596 0.2596 0.6053 0.3255

3.2. Area under ROC Curves for Pseudo-Real Data

To quantitatively test both sensitivity and specificity of different methods in a more realistic setting, we computed the area under the ROC curve (called “AUR”) for a false positive fraction (FPF) less than 0.1 as an index of method performance. This AUR quantity provides a measure of detection power for FPF < 0.1, which is relevant to the low FPF requirement in detection of fMRI activations in order to control the family-wise error rate.

The AUR quantities for the contrast V-F of the visual data are shown in Fig. 1, where curves of δ [0,1] for both cCCA methods along with the values of GLM-NS and GLM-GS are plotted. Since the induced activations at the visual cortex are strong and spatially extended, GLM-GS achieves nearly the best performance. Our proposed methods, cCCA-BB and cCCA-RG, perform similarly to each other and only outperform GLM-GS marginally for δ < 0.1. Note that the GLM-NS method works very poorly in the high noise case (f = 0.65) and increasing δ leads to a decreasing performance of cCCA methods.

(a) (b)

Figure 1. Performance of different data analysis methods showing the area under ROC curves (AUR), integrated over FPF ∈[0,0.1], as a function of the parameter δ in cCCA-BB and cCCA-RG using the contrast V-F of the visual paradigm for pseudo-real data: (a) Low noise case when f = 0.55 and (b) high noise case when f = 0.65. The horizontal lines correspond to GLM-NS (solid lines) and GLM-GS (dashed lines). Since the induced activations at the visual cortex are strong and spatially extended, GLM-GS has nearly optimum performance. The cCCA-BB and cCCA-RG methods perform similarly to each other and only outperform GLM-GS marginally for δ < 0.1. Note that the GLM-NS method works very poorly in the high noise case (f = 0.65) whereas both cCCA methods are very powerful for small δ. Note also, that an increase of δ leads to a decrease in performance of both cCCA methods.

In Fig. 2, we show the same plot of AUR results for the contrast E-C of memory data. Because memory activations are more subtle and localized, GLM-GS performs worse than GLM-NS whereas our proposed cCCA methods perform similarly to GLM-NS in the low noise case but gain much more detection power in the high noise case. Please also note that smaller values of δ result in more smoothing and thus lead to larger AUR quantities in the high noise case and thus better performance.

(a) (b)

Figure 2. Performance of different data analysis methods showing the area under ROC curves (AUR), integrated over FPF ∈[0,0.1], as a function of the parameter δ in cCCA-BB and cCCA-RG for the contrast E-C of the memory paradigm using pseudo-real data: (a) low noise case when f = 0.55 and (b) high noise case when f = 0.65. Because memory activations are more subtle and spatially localized, GLM-GS performs worse than GLM-NS whereas both cCCA-BB and cCCA-RG perform similarly and are optimal. Note that in the low noise case GLM-NS performs very close to the cCCA methods whereas in the high noise case both cCCA methods outperform GLM-NS. Please also note that smaller values of δ results in more local smoothing and thus leads to a better performance of the cCCA methods in the high noise case.

In addition, we plotted the curves of total false fraction (TFF) including both false positives and false negatives versus false positive fraction (FPF) in Fig. 3 (for the contrast V-F of the visual data) and Fig. 4 (for the contrast E-C of the memory data). This measurement provides another aspect of view of the detection performance of different methods. The δ in each case was chosen to have the smallest TFF.

(a) (b)

Figure 3. The total false fraction (TFF) including false positives and false negatives versus the false positive fraction (FPF) for the contrast V-F of the visual paradigm for pseudo-real data: (a) low noise case when f = 0.55 and (b) high noise case when f = 0.65. Note that all TFF curves have minima in the interval [0.001,0.01]. The GLM-GS performs poorly in the low noise case (a) but is very effective in the high noise case (b). Both cCCA methods especially cCCA-RG perform nearly optimal in both cases by achieving the minimum value of TFF.

(a) (b)

Figure 4. The total false fraction (TFF) including false positives and false negatives versus the false positive fraction (FPF) for the contrast E-C of the memory paradigm for pseudo-real data: (a) low noise case when f = 0.55 and (b) high noise case when f = 0.65. Both cCCA methods are optimal in both the low and high noise cases. Note that GLM-GS does not achieve less TFF than GLM-NS. This demonstrates that fixed Gaussian spatial smoothing does not have advantages when used on data with weak and spatially localized activations. Spatially adaptive smoothing as implemented in both cCCA methods is far more effective than any fixed spatial smoothing method.

For the strong activations of the contrast V-F, cCCA-RG achieves the smallest TFF at f = 0.55 [Fig. 3(a)], followed by cCCA-BB, GLM-NS and GLM-GS. The GLM-GS method is effective in the high noise case [Fig. 3(b), f = 0.65] and performs similarly to both cCCA methods. In Fig. 4, for the weak and small activations of the contrast E-C, the proposed cCCA methods become the optimum in both the low and high noise cases. The GLM-GS method works poorly even in the high noise case. This demonstrates that it is destructive to apply fixed Gaussian spatial smoothing on the data with weak and small activations. Adaptive spatial smoothing implemented by the proposed cCCA methods is more reliable and thus a better alternative to detect activations.

3.3. Activation Maps Using Real Data

To compare the four analysis methods using real visual and memory activation data, it is necessary to get the proper p values for the corresponding t- or novel cCCA statistics that are adjusted for multiple comparisons. The family-wise error rate (FWE) can be computed by using a nonparametric technique that does not rely on random field theory [16]. A nonparametric technique is necessary for a reliable comparison between different analysis methods because the parametric distribution of both cCCA methods is intractable due to the data-adaptive spatial filtering kernel. The FWE was calculated using Fourier resampled-resting-state data using bootstrapping of the order statistics [16]. All images in the figures are in radiological convention (left is right and vice versa) with the corrected p < 0.05 in activation maps.

In Fig. 5, we show the activation maps of the contrast V-F of visual data for different methods from one representative subject. It can be seen that GLM-GS yields the smoothest activation map at the expense of loss of the visual cortex structures, and GLM-NS preserves these folded structures much better but looks more unappealing because of broken links. Activation maps of cCCA-BB and cCCA-RG provide a good compromise between the smoothness of activations and preservation of fine cortical structure.

Figure 5. Activation maps for the contrast V-F of the visual paradigm using corrected p < 0.05. The GLM-GS method yields the smoothest activation map at the expense of showing activations reaching into white matter. The GLM-NS method preserves activations in gray matter much better but the map looks more unappealing because of apparent broken links between activated voxels in gray matter. Activation maps using cCCA-BB and cCCA-RG provide a compromise between a smooth appearance of activations and preservation of fine cortical structure localized to gray matter.

The activation maps of the contrast E-C of memory data from another representative subject are shown in Fig. 6. The slice shown contains part of the hippocampal complex. Symmetrical hippocampal and parahippocampal activations were detected by GLM-NS, cCCA-BB, and cCCA-RG. The missing activation at the left hippocampus of GLM-GS demonstrates the undesirable effects of a fixed isotropic Gaussian spatial smoothing on localized and weak activation patterns. Compared to GLM-NS, both cCCA methods detected more active voxels.

Figure 6. Activation maps for the contrast E-C of the memory paradigm using corrected p < 0.05. The slice chosen shows very localized activations in the anterior portion of the hippocampal complex. Note that GLM-GS, cCCA-BB and cCCA-RG lead to symmetric (left and right) activation patterns in the hippocampus and parahippocampal gyrus. It appears that GLM-GS is not able to detect activation at the left hippocampus (compare the GLM-GS map with all other maps at the location of the white arrow). This missing activation is likely to be caused by the isotropic Gaussian spatial smoothing technique that will eliminate localized and weak activations. Compared to GLM-NS, both cCCA methods detect more active voxels.

In Fig. 7 we show a more posterior slice for the E-C contrast. Please note memory encoding activation is obtained in the posterior cingulate cortex and precuneus. Using the GLM-GS method activations appear bulgy and have some unlikely connections through white matter (shown by the black arrow). Also, weak activations in the posterior cingulate cortex (see white arrows) could not be detected by GLM-GS. Both proposed cCCA methods yielded more connected activations in gray matter than GLM-NS. The GLM-GS method not only leads to missing activations but also to artifactual activations where a large fraction of false activation shows up in white matter and CSF regions due to the spherical (nondirectional) smoothing kernel.

Figure 7. Activation maps for the contrast E-C of the memory paradigm using corrected p < 0.05. The slice chosen shows very localized activations in the anterior portion of the hippocampal complex. Note that GLM-GS, cCCA-BB and cCCA-RG lead to symmetric (left and right) activation patterns in the hippocampus and parahippocampal gyrus. It appears that GLM-GS is not able to detect activation at the left hippocampus (compare the GLM-GS map with all other maps at the location of the white arrow). This missing activation is likely to be caused by the isotropic Gaussian spatial smoothing technique that will eliminate localized and weak activations. Compared to GLM-NS, both cCCA methods detect more active voxels.

3.4. Influence of the Constraint of the Center Voxel

The parameter δ in our cCCA methods is used to prevent the fitting of non-active voxels in the vicinity of active voxels from being erroneously boosted and thus prevents a loss in specificity. In our experiment, based on the AUR results, δ = 0 is good for localized and weak activations. However, the trade-off between specificity and sensitivity may require an adaptive assignment of δ. As shown in Fig. 8, the ROC curves for the contrast E-C at different δ values (0 and 0.5) in the high noise case (f = 0.65) display different behaviors. A value δ = 0 boosts the detection power for FPF > 10—3, while δ = 0.5 achieves better performance for smaller FPF (< 10—3). How to select the optimal δ remains a challenge. For this particular pseudo-real data, the selection of δ depends on what FPF is to be achieved. For FPF > 10—3, it is safe to use δ = 0 and more detection power is obtained. For FPF < 10—3, the adaptive smoothing has to be more conservative and a δ that is not too large can be used (0 < δ ≤ 0.7). Once δ > 0.7, the performance of cCCA converges to the performance of GLM-NS and thus cCCA lost its benefits. However, the selection of δ may be data dependent and as such δ is difficult to quantitatively optimize in real applications. We found that δ = 0 worked very well for our real visual and memory data (see Figs. 7–9), and there were no significant differences for small δ values (δ < 0.3).

(a) (b)

Figure 8. ROC curves for the contrast E-C of the memory paradigm for pseudo-real data using different δ for the high noise case (f = 0.65): (a) δ = 0 and (b) δ = 0.5. A small value of δ (i.e., δ = 0) increases the detection power for FPF > 10—3, while a larger δ (i.e., δ = 0.5) achieves better performance for smaller FPF (< 10—3). Insertions use a logarithmic scale of FPF to show more details for small values of FPF.

Although the current assignment of δ to a fixed value can achieve a notable improvement over traditional GLM methods, an adaptive assignment scheme of δ is worth exploring in the future, so that smaller δ can be used for stronger signals (where FPF is very small) and larger δ for weaker signals (where FPF is large).

3.5. Computation Time

Of all four methods, GLM-NS requires the least computational time, followed by GLM-GS, cCCA-RG, and cCCA-BB. For the estimation of a 2D slice with 6317 in-brain pixels, using MATLAB on a computer equipped with Intel Core 2 2.4 GHz CPU and 4 GB memory, it takes 2.0 s for GLM-NS, 2.2 s for GLM-GS, 24 s for cCCA-RG, and 308 s for cCCA-BB (averaged for 10 runs), both with δ = 0. When δ increases, the constraints become more stringent and cause a smaller size of neighborhood that satisfies the constraints. Consequently, the computation time decreases for cCCA-RG and increases for cCCA-BB, because the former grows from the center voxel and the latter trims voxels from the whole 3 × 3 area. Therefore, cCCA-RG is at least 12 times faster than cCCA-BB with similar detection performance. Although it takes much more time than the conventional GLM methods, cCCA-RG can be finished in a reasonable amount of time (less than 10 min for a 3D volume processing in our study) and is thus feasible for real applications.

4. CONCLUSION

In this article, we proposed a fast algorithm, cCCA-RG, for cCCA for fMRI data. This technique solves the problem of slow computation intrinsic to cCCA methods and can produce more accurate activation maps by detecting weak and subtle activations in fMRI data over traditional univariate GLM methods with or without isotropic Gaussian spatial smoothing. The cCCA-RG method is much faster than the previously reported branch-and-bound method and does not compromise accuracy, as determined by simulations and real fMRI data analysis. The proposed new cCCA-RG technique will become a preferred fMRI data analysis method since fMRI technology is moving toward high-resolution imaging where the number of voxels is large, the voxel size is small, and the MR signal is usually weak.

REFERENCES

1. Worsley, K. and Friston, K., Analysis of fMRI time-series revisited-again, NeuroImage, 2:173–181, 1995.

2. Worsley, K. J., Evans, A. C., Marrett, S., and Neelin, P., A three-dimensional statistical analysis for CBF activation studies in human brain, J. Cereb. Blood Flow Metab., 12:900–918, 1992.

3. Worsley, K. J., Marrett, S., Neelin, P., Vandal, A. C., Friston, K. J., and Evans, A. C., A unified statistical approach of determining significant signals in images of cerebral activation, Human Brain Mapping, 4:58–73, 1996.

4. Jin, M., Nandy, R., and Cordes, D., Effectiveness of Gaussian smoothing on spatially correlated noise: A 3T case analysis, In Proc. of International Society for Magnetic Resonance in Medicine 17th Scientific Meeting (ISMRM), p. 1709, 2009.

5. Friman, O., Cedefamn, J., Lundburg, P., Borga, M., and Knutsson, H., Detection of neural activity in functional MRI using canonical correlation analysis, Magn. Reson. Med., 45:323–330, 2001.

6. Friman, O., Borga, M., Lundberg, P., and Knutsson, H., Adaptive analysis of fMRI data, NeuroImage, 19:837–845, 2003.

7. Nandy, R. and Cordes, D., Improving the spatial specificity of CCA in fMRI, Magn. Reson. Med., 52:947–952, 2004.

8. Cordes, D., Jin, M., Curran, T., and Nandy, R., Optimizing the performance of local canonical correlation analysis in fMRI using spatial constraints, Human Brain Mapping, 2012.

9. Das, S. and Sen, P., Restricted canonical correlations, Linear Algebra Appl., 210:29–47, 1994.

10. Jin, M., Nandy, R., Curran, T., and Cordes, D., Extending local canonical correlation analysis to handle general linear contrasts for fMRI data, Int. J. Biomed. Imaging, p. 574971, 2012a.

11. Metz, C., Basic principle of ROC analysis, Semin. Nucl. Med., 8:283-98, 1978.

12. Sorenson, J. and Wang, X., ROC methods for evaluation of fMRI techniques, Magn. Reson. Med., 36:737–744, 1996.

13. Skudlarski, P., Constable, T., and Gore, J., ROC analysis of statistical methods used in functional MRI: Individual subjects, NeuroImage, 9:311–329, 1999.

14. Nandy, R. and Cordes, D., Novel ROC-type method for testing the efficiency of multivariate statistical methods in fMRI, Magn. Reson. Med., 49:1152–1162, 2003a.

15. Nandy, R. and Cordes, D., A novel nonparametric approach to canonical correlation analysis with applications to low CNR functional MRI data, Magn. Reson. Med., 49:354–365, 2003b.

16. Nandy, R. and Cordes, D., A semi-parametric approach to estimate the family-wise error rate in fMRI using resting-state data, NeuroImage, 34:1562–1576, 2007.

17. Jin, M., Pelak, V., Nandy, R., Curran, T., and Cordes, D., A preliminary study of functional abnormalities in aMCI subjects during different episodic memory tasks, Magn. Reson. Imaging, 30:459-70, 2012b.

18. Prichard, D. and Theiler, J., Generating surrogate data for time-series with several simultaneously measured variables, Phys. Rev. Lett., 73:951–954, 1994.

19. Jin, M. and Cordes, D., Validation of resampling methods for fMRI data, Neuroimage, 41:S525, 2008.