Archives: Volume 1, 2012 to Volume 2, 2013 |
## 生物医学的图像处理，计算和显示DOI: 10.1615/VisualizImageProcComputatBiomed.2012004964
## FAST CONSTRAINED CANONICAL CORRELATION ANALYSIS FOR FMRI
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.
fMRI, general linear model (GLM), constrained canonical correlation analysis (cCCA), fast implementation, branch-and-bound method, region-growing method ## 1. IntroductionThe 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 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. METHODSFor a group of
where
where In conventional fixed Gaussian smoothing, both 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
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, 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 2^{8} = 256 different pixel
configurations. This number increases to 2^{26} 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 ProblemWe use the notation to label the time series of a local neighborhood of
Y
with y_{1} ∈ Y so that the
spatial weight vector ^{*}α of Y satisfies the constraints and the correlation
of ^{*}Y and ^{*}αXβ 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 -
Try CCA using
**Y**containing all^{*}*K*column vectors of**Y**and**X**. -
Try CCA using all possible
**Y**containing^{*}*K*–1 column vectors of**Y**(including the center voxel**y**_{1}) and**X**. -
Try CCA using all possible
**Y**containing^{*}*K*–2 column vectors of**Y**(including the center voxel**y**_{1}) and**X**. ⋮
If at any step the vector 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
where 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 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 DataFunctional 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 ParadigmFor 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 ParadigmWe 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. PreprocessingAll 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 ## 2.4. Construction of Simulated and Pseudo-Real DataTo quantitatively determine the performance of different methods, we constructed both simulated and pseudo-real data in the following way:
In this equation ## 2.5. Methods of fMRI Data AnalysisFour 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 ## 3. RESULTS AND DISCUSSION## 3.1. Coefficient Estimation for Simulated DataTo 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 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
## 3.2. Area under ROC Curves for Pseudo-Real DataTo 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 (
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.
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.
For the strong activations of the contrast V-F, cCCA-RG achieves the smallest TFF at
## 3.3. Activation Maps Using Real DataTo compare the four analysis methods using real visual and memory activation
data, it is necessary to get the proper 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.
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.
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.
## 3.4. Influence of the Constraint of the Center VoxelThe 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 (
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 TimeOf 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. CONCLUSIONIn 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. ## REFERENCES1. Worsley, K. and Friston, K., Analysis of fMRI time-series revisited-again, 2. Worsley, K. J., Evans, A. C., Marrett, S., and Neelin, P., A three-dimensional
statistical analysis for CBF activation studies in human brain, 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, 4. Jin, M., Nandy, R., and Cordes, D., Effectiveness of Gaussian smoothing on spatially
correlated noise: A 3T case analysis, 5. Friman, O., Cedefamn, J., Lundburg, P., Borga, M., and Knutsson, H.,
Detection of neural activity in functional MRI using canonical correlation
analysis, 6. Friman, O., Borga, M., Lundberg, P., and Knutsson, H., Adaptive analysis of fMRI data,
7. Nandy, R. and Cordes, D., Improving the spatial specificity of CCA in fMRI, 8. Cordes, D., Jin, M., Curran, T., and Nandy, R.,
Optimizing the performance of local canonical correlation
analysis in fMRI using spatial constraints, 9. Das, S. and Sen, P., Restricted canonical correlations, 10. Jin, M., Nandy, R., Curran, T., and Cordes, D., Extending local canonical correlation
analysis to handle general linear contrasts for fMRI data, 11. Metz, C., Basic principle of ROC analysis, 12. Sorenson, J. and Wang, X., ROC methods for evaluation of fMRI techniques, 13. Skudlarski, P., Constable, T., and Gore, J., ROC analysis of statistical methods used in functional MRI:
Individual subjects, 14. Nandy, R. and Cordes, D., Novel ROC-type method for testing the efficiency
of multivariate statistical methods in fMRI, 15. Nandy, R. and Cordes, D., A novel nonparametric approach to canonical correlation
analysis with applications to low CNR functional MRI data, 16. Nandy, R. and Cordes, D., A semi-parametric approach to estimate the family-wise error rate in
fMRI using resting-state data, 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, 18. Prichard, D. and Theiler, J., Generating surrogate data for time-series with several
simultaneously measured variables, 19. Jin, M. and Cordes, D., Validation of resampling methods for fMRI data, |

Begell Digital Portal | Begell 数字图书馆 | 电子图书 | 期刊 | 参考文献及会议录 | 研究收集 |