Data-driven analysis methods, in particular independent com- ponent analysis (ICA) ... estimation of the task-related time courses and in the spatial activation.
FLEXIBLE COMPLEX ICA OF FMRI DATA Hualiang Li, T¨ulay AdalÕ, Nicolle Correa, Pedro A. Rodriguez and Vince D. Calhoun University of Maryland Baltimore County, Baltimore, Maryland 21250 {lihua1, adali}@umbc.edu ABSTRACT
or the link in such areas can be stronger when the phase and magnitude of the data are jointly processed. Using the phase Data-driven analysis methods, in particular independent cominformation together with magnitude is also likely to lead to ponent analysis (ICA) has proven quite useful for the analysis better localization of the origin of the signal by helping us to of functional magnetic imaging (fMRI) data. In addition, by distinguish larger vessels from smaller ones. enabling one to work in its native, complex form, complexIndependent component analysis (ICA) has been shown valued ICA algorithms provide better estimation performance to be quite effective for the analysis of fMRI data. By uscompared to the traditional approach that uses only the maging a simple generative model based on linear mixing, ICA nitude data. In the complex domain, circularity has been a can minimize the constraints imposed on the temporal—or common assumption even though most data acquisition meththe spatial—dimension of the fMRI data, and hence provides ods collect fMRI data that end up being noncircular when valuable new insights, especially when studying paradigms saved in complex form. In this paper, we show that a comfor which reliable models of brain activity are not available. plex ICA approach that does not assume circularity and also Thus, ICA emerges as a desirable alternate tool which promises adapts to the source density is the more desirable one for to be quite valuable for analysis of the natively complex fMRI performing ICA of complex fMRI data. We show that by data as well. adaptively matching the underlying fMRI density model, the Advantages of complex ICA for analyzing fMRI data in analysis performance can be improved in terms of both the its native form has been noted recently [3, 5, 11]. However, estimation of the task-related time courses and in the spatial the complex ICA algorithms implemented in these references activation. all use a ¿xed nonlinearity, which is used to implicitly genIndex Terms— Independent component analysis, fMRI, erate the higher-order statistics and assumes a ¿xed density complex analysis, ICA model for the underlying fMRI sources. Some other wellknown complex ICA algorithms, such as the complex Fas1. INTRODUCTION tICA algorithm [12] and circular Infomax algorithm [6], also assume a ¿xed nonlinearity function to match the underlying Functional magnetic resonance imaging (FMRI) has contributed density model. For all the complex ICA algorithms presented greatly to our understanding of the human brain by enabling here, the optimization criteria of them are based on maxiresearchers to study temporal and spatial changes in both the mum likelihood (ML)—which is equivalent to information healthy and the diseased brain as a function of various stimuli. maximization—and maximization of negentropy (MN) [1]. The MRI signal is acquired as a quadrature signal using two It can be shown that the two measures are equivalent when orthogonal detectors. The signal that is thus acquired in the the demixing matrix is constrained to be unitary. For both complex frequency space (k-space) is inverse Fourier transML and MN, the algorithms are optimal when the form of the formed into the complex image space. After this stage, most nonlinear function used in the cost function matches the form studies discard the phase and primarily use the magnitude imof the probability density functions of the sources. Another ages. However, processing of the fMRI data in its native comimportant consideration when performing ICA in the complex form is attractive for a number of reasons [2]. In [3], we plex domain is the circular/noncircular nature of the source show on the average 20% increased sensitivity for detection distribution. Since very little is known about the nature of the of task-related components when compared to the magnitude fMRI sources, and in particular of the complex fMRI data, it only approach. Also, by using both the magnitude and phase would be much more desirable to use an adaptive algorithm of the time courses in fMRI data, one is likely to elucidate such that the density model is estimated for each source adapdifferent networks while studying the brain connectivity as tively. One such algorithm is introduced in [14] where a genthe magnitude and phase potentially activate different areas, eralized Gaussian density model is used for the sources and the parameters of the model are estimated during the adapThis work is supported by the NSF grants NSF-CCF 0635129 and NSFtation, estimation of the demixing matrix W. In this paper, IIS 0612076.
978-1-4244-4296-6/10/$25.00 ©2010 IEEE
2050
ICASSP 2010
we show the application of the algorithm introduced in [14], the adaptive complex maximization of non-Gaussianity (ACMN) algorithm to complex fMRI data and show that indeed, it leads to improved performance. In addition, we address a number of issues in the preprocessing and visualization of ICA results of complex fMRI data and demonstrate the results using actual fMRI data for multiple subjects performing a motor task. In the next section, we ¿rst give a brief review of complex statistics and ICA, and then introduce the A-CMN algorithm in section 3. In section 4, we explain the pre- and postprocessing for the complex fMRI data and provide comparison of the performances of three competitive complex ICA algorithms when applied to fMRI analysis. We present our discussions in the last section. 2. COMPLEX ICA OF FMRI DATA The probability density function (pdf) of a complex random variable X = Xre + iXim is de¿ned through its joint density function pX (x) pXre Xim (xre , xim ). A complex random variable X is called circular in the strict-sense if X and Xejθ have the same pdf. A zero-mean random variable is called second-order circular (or proper) when its pseudo-covariance is zero, i.e., E{X 2 } = 0. In the ICA analysis of fMRI data, we assume independence of spatial brain activations (for Nspatial ICA), and write the complex ICA model as X = i=1 ai sTi where the matrix X ∈ CT ×V is formed using the fMRI data such that the ith row is formed by Àattening the volume image data of V voxels into a row, and the rows are indexed as a function of time. The assumed mixing column vector a represents the time course for the ith source and the independent source vector si is the ith spatial map. The task of the ICA algorithm is to determine a weight matrix W such that y = Wx = PΛs, where P, a permutation matrix, represents the permutation ambiguity and Λ, a diagonal matrix, represents the complexvalued scaling ambiguity of ICA. When performing complex ICA, most complex ICA algorithms assume that the sources are circular, such as the complex FastICA algorithm [12] and the Infomax algorithm that uses a circular nonlinearity [6]. It is desirable to relax the assumption of circularity so that the estimated source distributions are not constrained but also preprocessing inside the scanner is likely to lead to noncircular data [2, 8]. In addition, it is desirable to match the source distributions during the ICA estimation process. In the next section, we introduce such an adaptive ICA algorithm.
distribution (CGGD) model for the underlying sources. The form of the CGGD density model is given as β(c) −[η(c)¯sH C ¯ −1 ¯ s] c p(¯s) = e ¯ |C| ¯ is the where ¯s = [s s∗ ]T is the augmented source vector, C Γ(2/c) cΓ(2/c) covariance matrix of ¯s, η(c) = 2Γ(1/c) , β(c) = πΓ(1/c) 2 and c is the shape parameter. Using this CGGD model, we can model both sub-Gaussian (c > 1) and super-Gaussian (0 < c < 1)sources can be modeled through the Àexible shape parameter and we can model bth circular and noncircular sources through the estimation of the covariance matrix. The A-CMN cost function is de¿ned as J(y) = E{(yy ∗ )c } where y is the estimated source in a deÀationary mode. The ¯ are estimated onshape parameter c and covariance matrix C line using a maximum likelihood estimator. This on-line adaptation, modi¿es the cost function to match each source distribution and hence improves the overall performance. 4. RESULTS 4.1. fMRI data The dataset used in this paper is from 16 subjects performing a ¿nger-tapping motor task while receiving auditory instructions. The paradigm has a block design with alternating periods of 30s ON (¿nger tapping) and 30s OFF (rest). The experiments were performed on a 3T Siemens TRIO TIM system with a 12-channel radio frequency (RF) coil. The fMRI experiment used a standard Siemens gradient-echo EPI sequence modi¿ed to store real and imaginary data separately. The data was pre-processed for motion correction and spatial normalization into standard Montreal Neurological Institute space using the MATLAB Toolbox for Statistical parametric Mapping (SPM) [17]. 4.2. ICA algorithms Besides the A-CMN algorithm, for the comparison, we include the complex FastICA [12] algorithm and the circular Infomax algorithm [6] in our experiments. The cost function used in the complex FastICA algorithm is given by J(y) = log(|y|2 ). In the circular Infomax algorithm, the score func−|y| . Both algorithms assume tion is chosen as [sign(y)] 1−e 1+e−|y| that the sources are circular but have been shown to be robust choices for a number of applications including ICA of fMRI data [7]. 4.3. Preprocessing
3. ADAPTIVE COMPLEX ICA ALGORITHM The Adaptable complex maximization of non-Gaussianity (ACMN) algorithm [14] uses a complex Generalized Gaussian
2051
First, we apply quality map phase de-noising (QPMD) method presented in [15] to generate the mask to minimize the effects of noise in the phase of the fMRI data. The binary mask generated by QPMD identi¿es the good quality voxels in each
perform the clustering. We set the cluster number equal to the number of estimated components. If the ICA algorithms give consistent results over different runs, the output of ICASSO will have 26 compact clusters for each subject. Since the stability of the algorithms is not the focus of study in this paper, we use ICASSO only for determining the centrotype of 10 runs that has highest correlation with other components within one cluster. Therefore, for each subject, there are 26 centrotype estimates for each algorithm. 4.4. Performance evaluation The two components we considered for this fMRI data set are the motor and the temporal components. The motor component is de¿ned as the right motor activation area responding to the task, while the temporal component is the area activated in the bilateral temporal lobe. We observed that these two components have been consistently estimated by all three algorithms. To compare their performance, we consider three measures: 1) the number of activated voxels for the two components; 2) the maximum activation value for the two components; 3) for the motor (task-related) component, the correlation between the estimated time course and a model time course, which is generated by convolving a temporal model Fig. 1. Right motor component estimated by (a) A-CMN (b) of the on-off task with a canonical hemodynamic response complex FastICA algorithms function. To count the number of activated voxels for an activation map, we use WFU Pickatlas [16] to create the mask by selectfMRI slice for every fMRI data volume obtained at every ing different areas of the brain. We create two masks for the time point. Then the fMRI data are multiplied by the mask right motor component and the bilateral temporal component. and smoothed. The original dimension of the fMRI data is These masks include neuronal areas that are expected to be 16×153594×165 where 16 is the number of subjects, 153594 activated during the task as well as those in the temporal lobe. is the number of voxels in all the slices, 165 is the number The areas included in each of the two masks are as follows: of time points. After masking and smoothing, the data size 1) right motor task-related: Brodmann areas 1, 2, 3, 4 and 6. decreases to approximately 16 × 50000 × 165. For each sub2) bilateral temporal lobe: Brodmann areas 20, 21, 22, 37, 38 ject, we use complex principal component analysis (PCA) to and 42. further reduce the dimension of time points from 165 to 26, After the identi¿cation of the right motor and the temporal where the PCA weighting matrix is determined by selecting component, we count the number of voxels in the estimated the ¿rst 26 signi¿cant eigenvectors of the covariance matrix components that overlap with the corresponding masks for of the data. Then the three complex ICA algorithms, A-CMN, each algorithm and each subject. The average number of acticomplex FastICA and circular Infomax, are performed on the vated voxels for all the subjects are displayed in Table 1 where ¿nal data set. The resulting mixing matrix is then used to RM denotes the right motor components and TMP denotes back-reconstruct spatial maps and to obtain the time courses the temporal components. As shown in the table, the numfor each subject. ber of activated voxels estimated by A-CMN algorithm is the highest and the performance of circular Infomax is closely For each subject, we run the three algorithms 10 times behind, and both perform better than complex FastICA. with different initial conditions. Since the ICA estimates for iterative algorithms change for different runs, we use the ICASSO The maximum activation values, thresholded by a Z-score [13] scheme to calculate the centrotype of the 10 results for value 2, are displayed in Table 2. Again, the A-CMN algoeach algorithm and each subject. ICASSO collects the comrithm achieves the best performance with the circular Infomax ponents estimated from all runs and then matches components algorithm closely behind. The complex FastICA algorithm across runs by clustering components based on the absolute obtains the lowest activation value for both of the two comvalue of the correlation between squared source estimates. ponents. The correlations between the estimated time course Since we have both magnitude and phase results, here we use and the model time course are displayed in Table 3. From this the correlations between the magnitude of the estimates to Table, the improvement on the estimation of the phase time
2052
course by A-CMN algorithm is clearly observed. An example spatial map estimated by A-CMN algorithm is shown in Fig.??, where both the magnitude time course and phase time course have correlation values as high as 0.86 with the model time course. As a comparison, one sample estimate by complex FastICA algorithm is shown in Fig.1. The two ¿gures demonstrate that A-CMN algorithm can obtain much better results on the estimation of phase time course. Hence, based on the three performance indices, we note the importance of relaxing the circularity assumption of the sources and adaptively matching the underlying fMRI source model for improving the performance of the complex ICA of fMRI data. RM No of voxels
A-CMN 728 ± 275
CfastICA 622 ± 310
C-Infomax 667 ± 238
TMP No of voxels
A-CMN 1855 ± 352
CFastICA 1782 ± 369
C-Infomax 1829 ± 440
the ICA of fMRI data in its native, complex form, such as denoising and thresholding of the results for visualization. 6. REFERENCES [1] T. AdalÕ, H. Li, M. Novey, and J. Cardoso, “Complex ICA using nonlinear functions,” IEEE Trans. Signal Proc., vol. 56, no. 9, pp. 4536–4544, 2008. [2] T. AdalÕ and V. D. Calhoun, “Complex ICA of brain imaging data,” IEEE Signal Processing Magazine, vol. 24, no. 5, pp. 136–139, Sep. 2007. [3] V. D. Calhoun, T. AdalÕ, G. D. Pearlson, P. C. van Zijl, and J. J. Pekar, “Independent component analysis of fMRI data in the complex domain,” Magn Reson. Med., vol. 48, pp. 180–192, 2002. [4] V. D. Calhoun and T. AdalÕ, “Unmixing fMRI with independent component analysis,” IEEE Eng. Med. Bio. Mag., vol. 25, no. 2, pp. 79–90, 2006.
Table 1. Number of voxels in activation area of components that overlap with the corresponding mask
[5] V. D. Calhoun, T. AdalÕ, and Y. O. Li, “Independent component analysis of complex-valued functional MRI data by complex nonlinearities,” in Proc. IEEE Int. Symp. Bio. Imag. (ISBI), Arlington, VA, 2004. [6] J. Anem¨uller, T. J. Sejnowski and S. Makeig, “Complex independent component analysis of frequency-domain electroencephalographic data,” Neural Networks, vol. 16, no. 9, pp. 1311–1323, 2003.
RM Maximum Z-score
A-CMN 12 ± 2.4
CfastICA 11 ± 4
C-Infomax 12 ± 2.4
[7] N. Correa, T. AdalÕ, and V. D. Calhoun, “Performance of blind source separation algorithms for fMRI analysis using a group ICA method,” Magn. Reson. Imag., vol. 25, no. 5, pp. 684–694, June 2007.
TMP Maximum Z-score
A-CMN 11 ± 2.1
CFastICA 9 ± 2.6
C-Infomax 10 ± 2.3
[8] F. G. Hoogenraad, et al., “Quantitative differentiation between BOLD models in fMRI” Magn. Reson. Med., vol. 45, pp. 233–246, 2001.
Table 2. Maximum Z-scores of the estimated components
[9] M. J. McKeown, et al., “Analysis of fMRI data by blind separation into independent spatial components,” Human Brain Mapping, vol. 6, no. 3, pp. 160–188, 1998.
RM Magnitude
A-CMN 0.68 ± 0.24
CfastICA 0.64 ± 0.25
C-Infomax 0.68 ± 0.26
RM Phase
A-CMN 0.69 ± 0.19
CFastICA 0.56 ± 0.24
C-Infomax 0.6 ± 0.27
[10] F. Zhao, et al., “Sources of phase changes in BOLD and CBV-weighted fMRI,” Magn Reson. Med., vol. 57, pp. 520–527, 2007. [11] W. Xiong, Y. O. Li, H. Li and T. AdalÕ, “On ICA of complexvalued fMRI: Advantages and order selection,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing (ICASSP), Las Vegas, Nevada, April 2008.
Table 3. Correlation of the estimated time courses with model time courses
[12] E. Bingham an A. Hyvarinen, “A fast ¿xed-point algorithm for independent component analysis of complex valued signals,” Int. J. Neural Systems., vol. 10, pp. 1–8, Feb. 2000. [13] J. Himberg and A. Hyvarinen, “ICASSO: software for investigating the reliability of ICA estimates by clustering and visualization,” Proc. NNSP, Toulouse, France, pp. 259–268, 2003.
5. DISCUSSION ICA is one of the most fruitful methods for the study of fMRI data and the use of phase information in addition to the typically used magnitude promises to provide new insights for the analysis. In this paper, we showed that by relaxing the circularity assumption and adaptively matching the density model, we can improve the performance of the analysis results for complex fMRI data using ICA. In particular, we demonstrated that the adaptive A-CMN algorithm provides a more Àexible way to extract the intrinsic features of the noisy fMRI data. In addition, we addressed important issues when performing
2053
[14] M. Novey and T. AdalÕ, “Adaptable nonlinearity for complex maximization of nongaussianity and a ¿xed-point algorithm,” in Proc. MLSP, Maynooth, Ireland, pp. 79–84, Sep. 2006. [15] P. Rodriguez, N. Correa, T. AdalÕ and V. Calhoun “Quality map thresholding for de-noising of complex-valued fMRI data and its application to ICA of fMRI,” in Proc. MLSP, Sep. 2009. [16] J. Maldjian, P. Laurienti, R. Kraft, and J. Burdette “An automated method for neuroanatomic and cytoarchitectonic atlas-based interrogation of fMRI data set,” NeuroImage,vol. 19, pp. 1233–1239, July 2003. [17] SPM (http://www.¿l.ion.ucl.ac.uk/spm).
