- Email: [email protected]
Region growing method for the analysis of functional MRI data
NeuroImage 20 (2003) 455– 465
www.elsevier.com/locate/ynimg
Region growing method for the analysis of functional MRI data Yingli Lu,a Tianzi Jiang,a,b,* and Yufeng Zanga a
National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, Beijing 100080, People’s Republic of China b Department of Computer Science, University of Houston, Houston, TX 77204-3010, USA Received 14 April 2003; revised 3 June 2003; accepted 5 June 2003
Abstract Existing analytical techniques for functional magnetic resonance imaging (fMRI) data always need some specific assumptions on the time series. In this article, we present a new approach for fMRI activation detection, which can be implemented without any assumptions on the time series. Our method is based on a region growing method, which is very popular for image segmentation. A comparison of performance on fMRI activation detection is made between the proposed method and the deconvolution method and the fuzzy clustering method with receiver operating characteristic (ROC) methodology. In addition, we examine the effectiveness and usefulness of our method on real experimental data. Experimental results show that our method outperforms over the deconvolution method and the fuzzy clustering method on a number of aspects. These results suggest that our region growing method can serve as a reliable analysis of fMRI data. © 2003 Elsevier Inc. All rights reserved.
Introduction Functional magnetic resonance imaging (fMRI) has become one powerful tool for research into human brain functions due to its noninvasiveness and high spatial-temporal resolution property. The signals detected with fMRI are mostly based on blood oxygenation level dependence (BOLD) (Ogawa et al., 1990). Therefore, fMRI can capture the spatial-temporal neural activity related to sensory stimulus. However, the functional signals are rather small and are contaminated by systematic and physiological noise. The complicated structure of the mixed signals presents a formidable challenge to analysis of fMRI data. A number of methods have been proposed for activation detection. Among them, we can classify two categories: model-driven methods and data-driven methods. The typical model-based method is the general linear model (GLM) (Friston et al., 1995). GLM is a framework that includes simple t-test, ANOVA, ANCOVA, and multiple regressions. In GLM, once the design matrix is specified, multiple regressions can be used to estimate the parameters. Statistics for the obtained parameters can be
* Corresponding author. Fax: ⫹86-10-62551993. E-mail address: [email protected] (T. Jiang). 1053-8119/03/$ – see front matter © 2003 Elsevier Inc. All rights reserved. doi:10.1016/S1053-8119(03)00352-5
considered to constitute a statistical parameter map. Then a hypothesis test is made under the null hypothesis that the voxel is nonactivated. The power of GLM is that the number of false positives can be limited through correctly assessing the P-value by utilizing the random field theory. On the other hand, limitations of the GLM are that it depends on a number of assumptions, as follows: (1) linearity is assumed (2) residuals are normally distributed, (3) the estimated parameters are normally distributed, and (4) the number of variables are chosen, i.e., how to specify the design matrix. Typical data-driven methods are clustering analysis (Baumgartner et al., 1997; Filzmoser et al., 1999; Goutte et al., 1999; Scarth et al., 1995), principal component analysis (PCA) (Hansen et al., 1999; Lai and Fang, 1999), and independent component analysis (ICA) (Calhoun et al., 2001; Jung et al., 2001; McKeown et al., 1998). These techniques can separate different types of hemodynamic responses without the hypothesis of paradigm or hemodynamic response function. The power of these methods is that they can identify nonanticipated or transient task-related components. These methods are beneficially complementary to model-based approaches. For clustering analysis, the underlying assumption is that the pattern of activation actually has a structure and can be classified into a few types of similar activations. Clustering methods have several
456
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
drawbacks; as follows: (1) a priori definition of the number of clusters, i.e., the “cluster validity” problem; (2) ill-balanced data problem, i.e., the activated regions represent a small proportion of the brain and can be embedded in the large amount of non-activated voxels (Fadili et al., 2000); (3) with larger cluster numbers, the computational complexity is high. PCA finds the orthogonal patterns capturing the greatest variance in the data. However, if task-related fMRI changes are only a small part of the total signal variance, orthogonal patterns with the greatest variance may lose the information of task-related activations or other interest signals. In addition, PCA is based only on second-order statistics. It cannot extract sources that are higher than secondorder dependencies. ICA is conceptually similar to PCA, and the only difference is that ICA assumes the produced patterns are mutually statistically independent components. In practice, for real experimental data, the orthogonality or statistical independent assumptions will not be fully satisfied. For fMRI data analysis, each model has its advantages and disadvantages. Therefore, it is important to realize that no single model can reflect all features of the data, and the assumptions underlying different models are different. Therefore, model selection for fMRI activation detection is an explorative procedure under the property of the actual data. In this article, we present a novel data-driven approach for fMRI data analysis based on the idea of region growing method. Region growing methods are applied on image segmentation literature (Morel and Solimini, 1995). These methods postulate that neighboring pixels within the same region have similar intensity values. For fMRI data, the inherent character is its tendency toward clustered activation. That is, fMRI activation detection is more likely to occur in clusters of several spatially contiguous voxels than in a single voxel (Katanota et al., 2002; Tononi et al., 1998). Therefore, it is natural to use the region growing method to capture the cluster feature of fMRI activation. The proposed method is suitable for the actual spatial and temporal feature of fMRI data and can realize the reliable and sensitive detection of clustered activation without any specific assumptions.
Materials and methods Data sets Null experiment Imaging was acquired on a 1.5 T scanner Siemens Sonata equipped with high-speed gradients. For functional imaging acquisitions, technical parameters were as follows: 2000/60 ms (TR/TE), 20 slices, 64 ⫻ 64 matrix, 90° flip angle, 22-cm FOV, 5-mm section thickness, and 2-mm gap. The null data set consisted of 180 volumes. During the acquisition, the subject was instructed to rest, be as motionless as
possible, and perform no specific cognitive task. We used actual baseline data from real fMRI scans instead of computer-simulated data, since the former more closely matches the noise structure of the real fMRI data. Simulated activations were added to the null datasets, as mentioned below. Auditory experiment This data set is from the Welcome Department of Cognitive Neurology and is used with the kind permission of Functional Imaging Laboratory (FIL) methods group. This experiment is a study about auditory bisyllabic stimulation. A modified 2 T Siemens scanner system was used to acquire the whole brain BOLD/EPI images with TR ⫽ 7 s. The volume size is 64 ⫻ 64 ⫻ 64, with the voxel size 3 ⫻ 3 ⫻ 3 mm. The paradigm design starts with rest and alternates between rest and stimulation. Each block is 42 s, and 84 scans are selected, i.e., there are 84/12 ⫽ 7 cycles. Structural images’ volume size is 256 ⫻ 256 ⫻ 54, with the voxel size 1 ⫻ 1 ⫻ 1 mm. Preprocessing Preprocessing was performed on both the null and the auditory experimental datasets. Images are first realigned with the least-squares approach and a six-parameter rigid transformation to reduce the effect of head motion. Then, a spatial smoothing is made to increase the signal-to-noise ratio with an isotropic Gaussian filter, in which the full width half-maximum (FWHM) was set to 6 mm. Last, time series were detrended using the first-order polynomial detrending method. The above-mentioned three procedures were all implemented within the SPM99 software. Simulated activation For the ROC analysis described below, simulated activations were added to the null data sets. The time series of the “activated” area was defined as a boxcar (starts with rest and alternates between rest and simulation, 10 sample points for both rest and stimulation, and totals 8 cycles) convolved with the hemodynamic response function (HRF). Here, the HRF was a combination of two ␥-functions (Friston et al., 1998a). Parameters of the two ␥-functions are the defaults of SPM99 software and the convolving of HRF with boxcar was implemented in the SPM99 software too. Then, the simulated activated time series were mixed with Gaussian white noise. Five simulated data sets were generated with the contrast-to-noise ratio (CNR) (Langer, 1996) varying among 0.2, 0.4, 0.6, 0.8, and 1.0, respectively. The corresponding amplitude of the activations was slightly smaller than a typical signal change observed in the BOLD contrast (Bandettini et al., 1992; Fadili et al. 2000; Filzmoser et al., 1999; Kwong et al., 1992). The volume size-of-the-presumed four activated regions corresponding to a 3 ⫻ 3 ⫻ 4, 4 ⫻ 4 ⫻ 4, 5 ⫻ 5 ⫻ 4, 6 ⫻ 6 ⫻ 4 cuboid, respectively (refer to Fig. 1).
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
457
the seeds by adding neighboring voxels that satisfy the homogeneity criterion. During the region selection process, some regions are selected from the obtained regions based on a predefined selection criterion. The rest of the regions are not considered.
Fig. 1. One slice of the simulated fMRI data. The white box represents the sites of artificial activated region. The top right-hand graph depicts the time series of an “activated” voxel, and the bottom right-hand graph displays the time series of a “nonactivated” voxel.
Region growing method The region growing method is a well-developed technique for image segmentation. It postulates that neighboring pixels within the same region have similar intensity values. The general idea of the region growing method is to group pixels with the same or similar intensities to one region according to a given homogeneity criterion. More precisely, the region growing method starts with a set of prespecified seed pixel(s) and grows from these seeds by merging neighboring pixels whose properties are most similar to the premerged region. Typically, the homogeneity criterion is defined as the difference between the intensity of the candidate pixel and the average intensity of the premerged region. If the homogeneity criterion is satisfied, the candidate pixel will be merged to the premerged region. The procedure is iterative: at each step, a pixel is merged according to the homogeneity criterion. This process is repeated until no more voxels are assigned to the region. In an in vivo fMRI experiment, the true neural activation typically tends to occur in a functional cluster (Katanota et al., 2002; Tononi et al., 1998). In other words, true fMRI activation is more likely to occur in clusters of some spatial connected voxels than a single voxel. Based on this feature, intuitively, we can utilize the idea of region growing method for fMRI activation detection. However, the following key questions remain: (1) How to set the homogeneity criterion? Which should be based on the spatial time series? (2) How to set the initial seed voxel(s)? Different from three-dimensional image, the fMRI data is typically a spatial-temporal signal, which is represented as a four-dimention array F ⫽ {fxyzt}MNOT, where (x, y, z) is the spatial 3D coordinates of brain voxels, and t is the time index. M ⫻ N ⫻ O denotes the total number of brain voxels in an image scan and T represents the total number of image scans. fxyzt indicates the image intensity of the brain voxel (x, y, z) at the tth instance of time. The proposed region growing method solves the abovementioned two questions successfully. It is divided into two procedures: the region growing and the region selection. In the region growing procedure, each voxel in the image was set as the initial seed and regions are separately grown from
Region growing During initialization, each voxel vi, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} is specified as the initial seed of region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. For each region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}, we can define a set E as the unassigned voxels which are boundary points of region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}, E ⫽ 兵q ⰻ A i 兩 N共q兲 艚 A i ⫽ 0其 where N(q) is the set of immediate neighbors of voxel q. For each voxel p in E, if the predefined homogeneity criterion is satisfied, voxel p is added to region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. This grown procedure is repeated until no more voxels are being assigned to region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. The final result is that for each voxel vi, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}, there is a corresponding grown region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. Region selection Set A ⫽ {A1, A2, . . . , AM⫻N⫻O}, S ⫽ null, i.e., S is defined as a null family of sets. Step 1: Find the maximal region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} in A, where maximal region indicates the region that has the largest voxel number. Reserve region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} in S, i.e., set S ⫽ S 艛 {Ai}, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. Then, regions derived from the seeds that belong to region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} are removed from A. Mathematically, A is updated by set A ⫽ A ⫺ {Aj}J, J 傺 {1, 2, . . . , M ⫻ N ⫻ O}. Here, index J is defined as J ⫽ {j 兩 vj 僆 Ai}. Step 2: Repeat step 1 until the volume size of the remaining regions in A is smaller than a predefined threshold T2. This means that only clusters larger than T2 are significant. The rest are not considered. The final result of region selection procedure is saved in S. The number of regions in S is greatly smaller than M ⫻ N ⫻ O. During region selection procedure, why do we need to record the maximal region on each repeat step? It is known that in an in vivo fMRI experiment, the activation typically tends to constitute a cluster with a peak amplitude and a gradually descending scope around the peak. We can define the peak as the activation center, i.e., the point that has the strongest coherence with the simulation task. During the region growing procedure, if the initial seed is the activation center of activation cluster, the grown region will be the
458
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
maximal one. On the other hand, if the initial seed voxel deviates from the activation center, the grown region will be a smaller one. In each repeat step, it is obvious that only the maximal region can capture the feature of activated clusters.
Therefore, in each repeat step, only the maximal region was reserved. The implementation of region growing method for fMRI time series is summarized in the following pseudocode.
Region growing Specify each voxel as the initial seed of region A1, A2, . . . , AM⫻N⫻O For i ⫽ 1: M ⫻ N ⫻ O Repeat For each voxel p in the E set If the homogeneity criterion is satisfied Add p to region Ai End If End For Until no more voxels are being assigned to region Ai End For
Region selection Set A ⫽ {A1, A2, . . . , AM⫻N⫻O}, find the maximal region Ai in A, set S ⫽ S 艛 {Ai}, A ⫽ A ⫺ {Aj}J. Repeat. Find the maximal region Ai in A, set S ⫽ S 艛 {Ai}, A ⫽ A ⫺ {Aj}J. Until the volume size of remained regions in A is smaller than T2. Practical issues In the region growing procedure, it is critical to know how many and which neighboring voxels should be considered. For a 3D volume, 6-connected, 18-connected, or 26 connected neighborhood in three orthogonal directions can be considered. In this article, 26-connected neighborhood is selected. In the region growing procedure, the definition of homogeneity criterion is another critical problem. In this article it is defined as corr 共 f xyz:, M i兲 ⬎ T 1 where fxyz: is the time series of voxel (x, y, z), Mi is the averaged time series of voxels in the premerged region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. corr(,) is the Pearson correlation coefficient. T1 is the predefined threshold. The final result of the region growing method results in some clusters. Based on the clustered feature of fMRI activation, the task-related clusters are assuredly survived from the region selection procedure. That is to say, the result of region selection will not lose any task-related clusters. The problem is that besides task-related clusters, there are some other homogenous clusters reserved. To extract the task-related clusters, we have to distinguish among responses of the different anatomical structures. The procedure mentioned below will help us to this end. As a data-driven method, it is clear that region growing method cannot stand for a fMRI activation detection step. Some
additional treatments are necessary to extract the task-related clusters. If paradigm information can be utilized, one could test for the significant periodicity via Fourier transform methods (for periodic paradigm) or for the significant correlation to the predefined reference (Fadili et al., 2000). If the paradigm information cannot be utilized, some a priori neurophysiological and anatomical knowledge must be employed. In this article, a test of significant correlation was exploited. First, a reference is specified. Then, Pearson correlation coefficient was calculated between the reference and each cluster’s averaged time series. Clusters with a correlation coefficient greater than the predefined threshold TRG is considered task-related. In addition, the task-related clusters may be intersectant. Such intersectant clusters are merged into a single cluster. The final result of this procedure results in unconnected clusters (components), which are what we expect. ROC analysis on simulated data We conducted a receiver operating characteristic (ROC) analysis on the simulated data sets to evaluate the effectiveness of the region growing method and to compare its performance with that of other methods. The following three statistical models are compared: (1) the deconvolution method, (2) fuzzy c-mean clustering analysis (FCA), and (3) the region growing method. The reason the deconvolution and the FCA method were chosen for comparison is twofold. On one hand, the deconvolution method is a popular hypothesis-driven method for fMRI activation detection (Ward, 2002). On the other hand, FCA has been proven to be a successful data-driven analysis method in analyzing fMRI data (Baumgartner et al., 1997, 2000; Moser et al., 1997). The deconvolution method is an optimal program of AFNI software (Cox, 1996; Ward, 2002). It first estimates
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
the impulse response function(s). The estimated impulse response function(s) is then convolved with the stimulus paradigm to yield the estimated response. Then, various statistics are calculated to indicate the “goodness” of the fit. Finally, according to a predefined level of statistics, those activated voxels can be identified. FCA partitions the volumetric data into a number of fuzzy clusters. To obtain the task-related clusters, Pearson correlation coefficients were calculated between a predefined reference and the centroid of each cluster. Clusters with a correlation coefficient greater than the predefined threshold TFCA are considered task-related. For each cluster, FCA produces a cluster membership map for the volumetric data. By thresholding a membership map at a predefined level of threshold, those voxels that belong to the cluster can be identified. This is analogous to- but not equivalent to, selecting a significant level of a statistical test. To give an ROC analysis, four steps were processed: First, the three methods were applied to the simulated datasets. Second, according to the known number of simulated activation, the true positive ratio in the activated region and the false positive ratio in the nonactivated region were calculated. Third, an ROC curve was calculated by plotting the true positive ratio on the false positive ratio. The curve corresponding to a certain method closest to the upper left corner should be the best. Fourth, as the conventional ROC analysis, the area under the ROC curve was taken as the detectability measure of the different methods. In the implementation of the deconvolution method, the maximum time lag is set according to the known HRF, and the degree of the polynomial in the baseline model is set to 1. For FCA, the fuzziness exponent weight is set to 2 according to previous work of Fadili (2000). There is no theory to determine the number of clusters prior to FCA. In this article, for each of the five simulated datasets, the number ranging between 2 and 30 was set as the candidate number of clusters. The number that obtained the best results, i.e., the maximal area under the ROC curve, was recorded as the cluster number. In addition, the distance measure is the hyperbolic correlation measure (Golay et al., 1998),
459
Analysis of in vivo data sets The auditory dataset was analyzed with the three methods referred to above. For the decovolution method, the maximum time lag is set to 5 and the degree of the polynomial in the baseline model is set to 1. For the FCA method, the fuzziness exponent weight was set to 2, and the number of clusters was set to 30 (this is a typically used value for FCA) (Baumgartner et al., 2000). The hyperbolic correlation measure is utilized as the distance measure. To overcome the ill-balanced data problem, prior to FCA, noisy pixels were removed by preprocessing (for all voxel time series, we calculated their 1-lag-shifted autocorrelation function, and P ⫽ 0.01 was selected for the significance level) (Fadili et al., 2000). For the region growing method, T2 is set to 4. TRG and TFCA are all set to 0.5. To extract the task-related clusters of FCA and the region growing method, the reference time series is the paradigm convolved with the HRF (combinations of two ␥-functions, SPM99). For all of the three methods, only the brain voxels are considered, which are obtained from a low-level threshold estimated from histogram. In addition, for deconvolution and FCA method, only the clusters with a cluster size larger than four voxels were considered significant. The rest of the clusters were regarded as noise. The aim of this is to stay consistent with the T2 threshold of the region growing method. To make a comparison among the three methods, a quantitative method was implemented. The focus of quantitative comparison is to select the threshold of different activation maps, as the threshold determines whether a voxel is activated or not. Moreover, the threshold determines the fairness of comparison. The scheme of the proposed quantitative comparison method is as follows: (1) some or other activated area obtained with the deconvolution method is selected as the reference area; (2) thresholds of the deconvolution (P value), the FCA (membership value), and the region growing (T1) method are adjusted by keeping the same cluster size of the reference area; (3) the volume size and the averaged time series of the nonreference area are taken as measures of effectiveness and sensitiveness of different methods.
d共 x k,v i兲 ⫽ 共1 ⫺ c ik兲/共1 ⫹ c ik兲, Results where cik denotes the Pearson correlation coefficient between the time series xk and cluster centroid vi. The reference time series is the simulated activated time series without mixing Gaussian noise. TFCA is set to 0.5. For the region growing method, the reference time series is same with that utilized in FCA. TRG is set to 0.5 (this corresponding to TFCA ⫽ 0.5). In addition, T2 is set to 2. Accordingly, for the other two methods, only clusters with a cluster size larger than two voxels were considered significant. The rest clusters were regarded as noise and hence will not be considered.
Simulated data Fig. 2 shows the area under the ROC curve for each of the three statistical models plotted as a function of CNR. We can conclude that the region growing method outperforms the deconvolution and the FCA method. To further clarify the increased sensitivity of the region growing method, we examine one particular case, with CNR at 0.6. The ROC curves of the three methods are depicted in Fig. 3a. The characteristics of the three models are well illustrated in this
460
Y. Lu et al. / NeuroImage 20 (2003) 455– 465 Table 1 Computation time in seconds for the three different methods in simulated dataset with different CNR Method
Deconvolution FCA Region growing
CNR 0.2
0.4
0.6
0.8
1.0
⬍30 4816 467
⬍30 3709 432
⬍30 2975 403
⬍30 3610 496
⬍30 3122 415
Note. FCA is calculated in the optimal cluster and the region growing method is calculated with T1 value that corresponds to the optimal classification point.
Fig. 2. Performance of the region growing method, the deconvolution method, and FCA method for the simulated data set.
figure. The ROC curves of the region growing method approached the top left corner primely, whereas the other two methods were far from the top left corner. In Fig. 3a, point x, y, and z are the points closest to the optimal classification (False Positive Rates ⫽ 0, True Positive Rates ⫽ 1) in the respective ROC curves. Figs. 3b, c, and d illustrate the activation maps obtained with the region grow-
ing method, the deconvolution method, and FCA, respectively. The corresponding thresholds were selected according to the value that corresponds to the optimal classification point x, y, and z. From this comparison, we can conclude that the region growing method is more sensitive than the FCA and the deconvolution method for activation detection. The computational complexities of the three methods are reported in Table 1. The computations were performed on a 1.8-GHz PC, and the region growing method, FCA, and the deconvolution method are all coded in Matlab (The MathWorks, Inc.). From this comparison, we can conclude that
Fig. 3. (a) ROC curves corresponding to the case of CNR ⫽ 0.6 for the simulated data. (b) Results of the region growing method. (c) Results of the deconvolution method. (d) Results of FCA.
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
461
Fig. 4. Averaged time series of the 18 clusters and the corresponding correlation coefficient between the reference and the averaged time series.
the region growing method has less computational complexity than the FCA. Auditory experiment To demonstrate how the region selection procedure affects the performance of the method, we set T1 to 0.95 TRG to 0.5 and T2 to 4. For simplicity, only the clusters in the 29 –32th slices of the total 64 slices are computed. There are a total of 133 clusters that survived from the region selection procedure. Thirty-two of the 133 clusters are taskrelated. The 32 clusters are merged into 2 larger clusters. The other 101 clusters are merged into 16 larger clusters. In Fig. 4 the averaged time series of the 18 clusters and the correlation coefficients with reference time series are shown. In Fig. 5, the two task-related clusters are shown in red; the other 16 clusters are shown in green. The region growing method provides, in addition to the activation clusters, a number of other clusters. Most of these are dominated by different sources of noise. In the context of fMRI data, how to interpret the huge amount of information provided by the data-driven methods is still a major concern. We refer to Fadili (2000) for an example on this streamline. We will not consider this issue here. The right auditory cortex in the 29th slice was set as the reference area. For the decovolution method, four thresh-
olds are selected corresponding to the P-value of 10⫺8, 10⫺7, 10⫺6, and 10⫺5. The corresponding threshold T1 of region growing method and membership value of FCA are listed in Table 2. Results of the three methods are shown in Fig. 6. All three methods can successfully detect the other activated areas such as Brodmann’s area 42 (primary auditory cortex) as well as the Brodmann’s area 22 (auditory association area). The detailed volume size of the activated region in each slice is listed in Table 2. Compared with the FCA method, the activated areas revealed by the region growing method are larger and more continuous. For example, the activation in the left auditory cortex is larger and continuous. In addition, small activated clusters scattered throughout the brain are fewer in the result obtained with the region growing method than that obtained with the FCA method. In Fig. 6d, for instance, there are some small false activated clusters with the FCA method. Compared with the deconvolution method, the cluster size detected with the region growing method is bigger and more continuous. In the left auditory cortex activated area, there are a total of 233 voxels that both activated in region growing and deconvolution method; the averaged time series of these 233 voxels is shown in Fig. 7a. There are 125 voxels that were activated in the region growing method but were not activated in the deconvlolution method. The averaged time series of these 125 voxels is shown in Fig. 7b.
462
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
Fig. 5. Clusters in the 29 –32th of the 64 axial slices. Red: task-related clusters; green: non-task-related clusters. The left side corresponds to the right hemisphere Fig. 6. Activation maps revealed by each of the three methods. Results of the three methods, i.e., the deconvolution method, the FCA, and the region growing method, are shown from the top row to the bottom of the figure. Four slices (29 –32th of the total 64 slices) involved in the auditory experiment are shown. The left side corresponds to the right hemisphere, a, b, c, and d are the results of the three methods that correspond to the deconvolution method thresholded at the voxelwise significance level of P ⬍ 10⫺8, P ⬍ 10⫺7, P ⬍ 10⫺6, and P ⬍ 10⫺5, respectively.
The correlation coefficient between this two averaged time series is 0.938. Therefore, we can conclude that the 125 voxels are all meaningful activated voxels. In the right auditory cortex activated region, there are totally 328 voxels that both activated in the deconvolution, and the region growing method; averaged time series of these 328 voxels is shown in Fig. 7c. There are 62 voxels that activated in region growing but nonactivated in the deconvolution mehod. Averaged time series of the 62 voxels is shown in Fig. 7d. The correlation coefficient between this two mean time series is 0.826. Therefore, we can conclude that the 62 voxels are all meaningful activated voxels. In this data set, these results indicate that the region growing method is more sensitive than the deconvolution method.
Discussion and conclusion Our ROC analysis is based on several assumptions, as other simulation studies in general. First, the property of true activated area is known a priori, as it is artificially simulated. This is the reason that we can get the optimal cluster number for FCA based on the area under the ROC curve. Although the optimal cluster number is used for FCA, the detection power is smaller than the region growing method. For the deconvolution method, though the linearity assumption is satisfied, the detection power is smaller than the region growing method. This indicates that the region growing method is not very sensitive to noise. Another issue is the fMRI noise. We used a computer-simulated Gaussian
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
463
Table 2 Volume size (voxel number) for different methods Method
Threshold
29(R)
30(R)
31(R)
32(R)
29(I)
30(I)
31(I)
32(I)
45 46 52
26 31 30
35 41 41
41 53 61
44 78 71
51 54 60
30 35 31
39 45 44
53 56 77
50 88 96
58 65 67
38 43 38
44 54 62
59 82 84
67 99 110
67 72 76
45 81 48
50 86 76
68 107 108
81 115 126
a Deconvolution FCA Region growing
10⫺8 0.19 0.895
97 98 96
66 77 77
52 53 54 b
Deconvolution FCA Region growing
⫺7
10 0.15 0.88
102 102 103
82 82 89
58 66 61 c
Deconvolution FCA Region growing
10⫺6 0.09 0.865
106 105 105
93 88 108
67 73 68 d
Deconvolution FCA Region growing
⫺5
10 0.02 0.85
117 115 116
102 94 116
78 89 82
Note. L: left; R: right. The threshold for the deconvolution, the FCA, and the region growing method is P value, membership, T1, respectively. Number with shading is the size of reference area. a, b, c and d correspond to the result of Figs. 6a, b, c, and d, respectively.
white noise. In practice, it is different from the experimental data sets. In addition, the shape of simulated activated clusters would be more complicated in reality. It is essential, therefore, to examine the effectiveness of our model using real experimental data. In this study, we employed the auditory data set. Compared with the other two methods, the region growing method has the three following advantages: (1) Volume size of the activated clusters is larger; (2) Activated clusters are more continuous; (3) Small scattered clusters are fewer. From the viewpoint of methodology, compared with the deconvolution method, the region growing method has three advantages: (1) It does not require any prior knowledge of the
HRF, which may not be generally known. (2) It does not need the linear assumptions. The evidence of linearity of the hemodynamic response should contrast with the work of Vazquez and Noll (1998), Friston et al. (1998b), and Glover (1999). (3) As a univariate technique, the deconvolution method ignores the relationship of neighboring voxels. Therefore, the deconvolution method is hindered from the feature of clustered activation. However, the region growing method can utilize the relationships of neighboring voxels. Compared with FCA, the region growing method has the following three advantages: (1) It does not need the predefined cluster number. (2) It is not sensitive to the ill-balanced data problem, so we minimize the risk of discarding possibly activated voxels. (3) FCA ignores
Fig. 7. (a) Averaged time series of the 233 voxels that activated in both region growing and the deconvolution in the left auditory cortex. (b) Averaged time series of the 125 voxels that activated in region growing but nonactivated in the deconvolution in the left auditory cortex. (c) Averaged time series of the 328 voxels that activated in both region growing and the deconvolution in the right auditory cortex. (d) Averaged time series of the 62 voxels that activated in region growing but nonactivated in deconvolution in the right auditory cortex.
464
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
the spatial relationships of neighboring voxels. However, the region growing method can utilize the relationships of neighboring voxels. In addition, region growing can serve as a data reduction preprocessing technique. As mentioned above, fMRI activation is more likely to occur in clusters of several spatial contiguous voxels than in a single voxel. The implementation of data reduction procedure by utilizing region growing is outlined in the following two steps: First, each voxel in the brain was set as the initial seed and regions are separately grown from the seed by adding neighboring voxels that satisfy the homogeneity criterion. Second, the voxels whose corresponding grown regions smaller than a predefined threshold, for instance, 4 – 6, are considered noise and will not participate in the subsequent analysis. This data reduction procedure can be applied to existing fMRI data analysis method, for example, GLM, the deconvolution, and FCA. A limitation of the region growing method is its computation complexity for larger activated regions. For instance, the computation time of FCA (recently, the application of FCA was expanded by introducing two additional stages (Jarmasz and Somorjai, 2002), whose effect is to speed up the computation and make the results more robust. In this article, we only utilized the typical FCA method for comparison) for the block data is 8676 s, but the computation time of the region growing method (with T1 ⫽ 0.85) is 15,442 s. The computation burden of the region growing method lies in the calculation of the region growing part. There are two feasible schemes to reduce the computation complexity. The first way is anatomical restriction. fMRI BOLD activation is thought to occur mainly in the gray matter. Thus, if only the voxels belonging to this issue are considered, the computation complexity will be reduced. The segmentation strategy proposed by Fadili (2000) can be applied for this purpose. The second way to reduce the computation complexity is to introduce parallel computation techniques. It is very easy to implement in a parallel way since each region can be grown separately and simultaneously on different CPUs or different computers. In summary, a new method has been presented for fMRI data analysis. An ROC analysis has demonstrated its effectiveness for fMRI activation detection. Its usefulness has been confirmed by applying it to a real experimental. The method can serve as a new reliable data-driven method suited to the nature of fMRI data analysis.
Acknowledgments The authors are grateful to the anonymous referees for significant and constructive comments and suggestions, which greatly improved the article. The corresponding author also thanks his colleague, Professor Ricardo Vilalta, Department of Computer Science, University of Houston,
for carefully checking the use of English in the manuscript. The authors give thanks to Professor Karl Friston for his kind permission to use the block trial data sets. This work was partially supported by Hundred Talents Programs of the Chinese Academy of Sciences, the Natural Science Foundation of China, Grant Nos. 60172056 and 697908001.
References Bandettini, P.A., Wong, E.C., Hinks, R.S., Tikofsky, R.S., Hyde, J.S., 1992. Time course EPI of human brain function during task activation. Magn. Reson. Med 25, 390 –397. Baumgartner, R., Ryner, L., Richter, W., Summers, R., Jarmasz, M., Somorjai, R., 2000. Comparison of two exploratory data analysis methods for fMRI: fuzzy clustering vs. principal component analysis. Magn. Res. Imaging 18, 89 –94. Baumgartner, R., Scarth, G., Teichtmeister, C., Somorjai, R., Moser, E., 1997. Fuzzy clustering of gradient-echo functional MRI in the human visual cortex: Part I: Reproducibility. J. Magn. Res. Imaging 7 (6), 1094 –1101. Calhoun, V.D., Adali, T., Pearlson, G.D., Pekar, J.J., 2001. Spatial and temporal independent component analysis of functional MRI data containing a pair of task-related waveforms. Hum. Brain Mapping 13, 43–53. Cox, R.W., 1996. AFNI: software for analysis and visualization of functional magnetic resonance neuoimages. Comp. Biomed. Res 19, 162– 173. Fadili, M.J., Ruan, S., Bloyet, D., Mazoyer, B., 2000. A multistep unsupervised fuzzy clustering analysis of fMRI time series. Hum. Brain mapping 10, 160 –178. Filzmoser, P., Baumgartner, R., Moser, E., 1999. A hierarchical clustering method for analyzing functional MR images. Mag. Res. Imaging 17, 817– 826. Friston, K.J., Fletcher, P., Josephs, O., Holmes, A., Rugg, M.D., Turner, R., 1998a. Event-related fMRI: characterizing differential responses. Neuroimage 7 (1), 30 – 40. Friston, K.J., Holmes, A.P., Worsley, K.J., Poline, J.P., Frith, C.D., Frackowiak, R.S.J., 1995. Statistical parametric maps in functional imaging: a general linear approach. Hum. Brain Mapping 2, 189 –210. Friston, K.J., Josephs, O., Ress, G., Turner, R., 1998b. Nonlinear eventrelated responses in fMRI. Magn. Res. Med 39, 41–52. Glover, G.H., 1999. Deconvolution of impulse response in event-related BOLD fMRI. Neuroimage 9, 416 – 429. Golay, X., Kollias, S., Stoll, G., Merier, D., Valavanis, A., Boesiger, P., 1998. A new correlation-based fuzzy logic clustering algorithm for fMRI. J Magn. Res. Med 40, 249 –260. Goutte, C., Toft, P., Rostrup, E., Nielsen, F.A., Hansen, L.K., 1999. On clustering fMRI time series. Neuroimage 9 (3), 298 –310. Hansen, L.K., Larsen, J., Nielsen, F.A., Strother, S.C., Rostrup, E., Savoy, R., Lange, N., Sidtis, J., Svarer, C., Paulson, O.B., 1999. Generalizable patterns in neuroimaging: how many principal components? Neuroimage 9, 534 –544. Jarmasz, M., Somorjai, R.L., 2002. Exploring regions of interest with cluster analysis (EROICA) using a spectral peak statistic for selecting and testing the significance of fMRI activation time-series. Artif. Intell. Med 25, 45– 67. Jung, T.P., Makeig, S., Mckeown, M.J., Bell, A.J., Lee, T.W., Sejnowski, T.J., 2001. Imaging brain-dynamics using independent component analysis. Proc. IEEE 89 (7), 1107–1122. Katanoda, K., Matsuda, Y., Sugishita, M., 2002. A spatial-temporal regression model for the analysis of functional MRI data. Neuroimage 17, 1415–1428. Kwong, K.K., Belliveau, J.W., Chesler, D.A., Goldberg, I.E., Weisskoff,
Y. Lu et al. / NeuroImage 20 (2003) 455– 465 R.M., Poncelet, B.P., Kennedy, D.N., Hoppel, B.E., Cohen, M.S., Turner, R., Cheng, H.M., Brady, T.J., Rosen, B.R., 1992. Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc. Natl. Acad. Sci. USA 89, 5675–5679. Lai, S.H., Fang, M., 1999. A novel local PCA-based method for detecting activation signals in fMRI. Magn. Res. Imaging 17 (6), 827– 836. Langer, N., 1996. Tutorial in biostatistics: statistical approaches to human brain mapping by functional magnetic resonance imaging. Stat. Med 15, 389 – 428. McKeown, M.J., Makeig, S., Brown, G.G., Jung, T.P., Kindermann, S.S., Bell, A.J., Sejnowski, T.J., 1998. Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapping 6 (3), 160 –188. Morel, J.M., Solimini, S., 1995. Variational Methods in Image Segmentation. Birkhauser Publishing, Inc., Basel. Moser, E., Diemling, M., Baumgartner, R., 1997. Fuzzy clustering of
465
gradient-echo functional MRI in the human visual cortex. Part II: Quantificaltion. J. Magn. Res. Imaging 7, 1102–1108. Ogawa, S., Lee, T.M., Ray, A.R., Tank, D.W., 1990. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc. Natl. Acad. Sci. USA 87, 9868 –9872. Scarth, G., Somorjai, R., Alexander, M., Wowk, B., Wennerberg, A., McIntyre, M., 1995. Cluster analysis of functional brain images, in: Proceedings of the First International Conference on Functional Mapping of the Human Brain, Vol. S1 of Hum. Brain Mapping 158. Tononi, G., McIntosh, A.R., Russell, D.P., Edelman, G.M., 1998. Functional clustering: identifying strongly interactive brain regions in neuroimaging data. Neuroimage 7 (2), 133–149. Vazquez, A.L., Noll, D.C., 1998. Nonlinear aspects of the BOLD response in functional MRI. Neuroimage 7, 108 –118. Ward, B.D., 2002. Deconvolution analysis of fMRI time series data, http://www.afni.nimh.gov/afni/edu.
www.elsevier.com/locate/ynimg
Region growing method for the analysis of functional MRI data Yingli Lu,a Tianzi Jiang,a,b,* and Yufeng Zanga a
National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, Beijing 100080, People’s Republic of China b Department of Computer Science, University of Houston, Houston, TX 77204-3010, USA Received 14 April 2003; revised 3 June 2003; accepted 5 June 2003
Abstract Existing analytical techniques for functional magnetic resonance imaging (fMRI) data always need some specific assumptions on the time series. In this article, we present a new approach for fMRI activation detection, which can be implemented without any assumptions on the time series. Our method is based on a region growing method, which is very popular for image segmentation. A comparison of performance on fMRI activation detection is made between the proposed method and the deconvolution method and the fuzzy clustering method with receiver operating characteristic (ROC) methodology. In addition, we examine the effectiveness and usefulness of our method on real experimental data. Experimental results show that our method outperforms over the deconvolution method and the fuzzy clustering method on a number of aspects. These results suggest that our region growing method can serve as a reliable analysis of fMRI data. © 2003 Elsevier Inc. All rights reserved.
Introduction Functional magnetic resonance imaging (fMRI) has become one powerful tool for research into human brain functions due to its noninvasiveness and high spatial-temporal resolution property. The signals detected with fMRI are mostly based on blood oxygenation level dependence (BOLD) (Ogawa et al., 1990). Therefore, fMRI can capture the spatial-temporal neural activity related to sensory stimulus. However, the functional signals are rather small and are contaminated by systematic and physiological noise. The complicated structure of the mixed signals presents a formidable challenge to analysis of fMRI data. A number of methods have been proposed for activation detection. Among them, we can classify two categories: model-driven methods and data-driven methods. The typical model-based method is the general linear model (GLM) (Friston et al., 1995). GLM is a framework that includes simple t-test, ANOVA, ANCOVA, and multiple regressions. In GLM, once the design matrix is specified, multiple regressions can be used to estimate the parameters. Statistics for the obtained parameters can be
* Corresponding author. Fax: ⫹86-10-62551993. E-mail address: [email protected] (T. Jiang). 1053-8119/03/$ – see front matter © 2003 Elsevier Inc. All rights reserved. doi:10.1016/S1053-8119(03)00352-5
considered to constitute a statistical parameter map. Then a hypothesis test is made under the null hypothesis that the voxel is nonactivated. The power of GLM is that the number of false positives can be limited through correctly assessing the P-value by utilizing the random field theory. On the other hand, limitations of the GLM are that it depends on a number of assumptions, as follows: (1) linearity is assumed (2) residuals are normally distributed, (3) the estimated parameters are normally distributed, and (4) the number of variables are chosen, i.e., how to specify the design matrix. Typical data-driven methods are clustering analysis (Baumgartner et al., 1997; Filzmoser et al., 1999; Goutte et al., 1999; Scarth et al., 1995), principal component analysis (PCA) (Hansen et al., 1999; Lai and Fang, 1999), and independent component analysis (ICA) (Calhoun et al., 2001; Jung et al., 2001; McKeown et al., 1998). These techniques can separate different types of hemodynamic responses without the hypothesis of paradigm or hemodynamic response function. The power of these methods is that they can identify nonanticipated or transient task-related components. These methods are beneficially complementary to model-based approaches. For clustering analysis, the underlying assumption is that the pattern of activation actually has a structure and can be classified into a few types of similar activations. Clustering methods have several
456
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
drawbacks; as follows: (1) a priori definition of the number of clusters, i.e., the “cluster validity” problem; (2) ill-balanced data problem, i.e., the activated regions represent a small proportion of the brain and can be embedded in the large amount of non-activated voxels (Fadili et al., 2000); (3) with larger cluster numbers, the computational complexity is high. PCA finds the orthogonal patterns capturing the greatest variance in the data. However, if task-related fMRI changes are only a small part of the total signal variance, orthogonal patterns with the greatest variance may lose the information of task-related activations or other interest signals. In addition, PCA is based only on second-order statistics. It cannot extract sources that are higher than secondorder dependencies. ICA is conceptually similar to PCA, and the only difference is that ICA assumes the produced patterns are mutually statistically independent components. In practice, for real experimental data, the orthogonality or statistical independent assumptions will not be fully satisfied. For fMRI data analysis, each model has its advantages and disadvantages. Therefore, it is important to realize that no single model can reflect all features of the data, and the assumptions underlying different models are different. Therefore, model selection for fMRI activation detection is an explorative procedure under the property of the actual data. In this article, we present a novel data-driven approach for fMRI data analysis based on the idea of region growing method. Region growing methods are applied on image segmentation literature (Morel and Solimini, 1995). These methods postulate that neighboring pixels within the same region have similar intensity values. For fMRI data, the inherent character is its tendency toward clustered activation. That is, fMRI activation detection is more likely to occur in clusters of several spatially contiguous voxels than in a single voxel (Katanota et al., 2002; Tononi et al., 1998). Therefore, it is natural to use the region growing method to capture the cluster feature of fMRI activation. The proposed method is suitable for the actual spatial and temporal feature of fMRI data and can realize the reliable and sensitive detection of clustered activation without any specific assumptions.
Materials and methods Data sets Null experiment Imaging was acquired on a 1.5 T scanner Siemens Sonata equipped with high-speed gradients. For functional imaging acquisitions, technical parameters were as follows: 2000/60 ms (TR/TE), 20 slices, 64 ⫻ 64 matrix, 90° flip angle, 22-cm FOV, 5-mm section thickness, and 2-mm gap. The null data set consisted of 180 volumes. During the acquisition, the subject was instructed to rest, be as motionless as
possible, and perform no specific cognitive task. We used actual baseline data from real fMRI scans instead of computer-simulated data, since the former more closely matches the noise structure of the real fMRI data. Simulated activations were added to the null datasets, as mentioned below. Auditory experiment This data set is from the Welcome Department of Cognitive Neurology and is used with the kind permission of Functional Imaging Laboratory (FIL) methods group. This experiment is a study about auditory bisyllabic stimulation. A modified 2 T Siemens scanner system was used to acquire the whole brain BOLD/EPI images with TR ⫽ 7 s. The volume size is 64 ⫻ 64 ⫻ 64, with the voxel size 3 ⫻ 3 ⫻ 3 mm. The paradigm design starts with rest and alternates between rest and stimulation. Each block is 42 s, and 84 scans are selected, i.e., there are 84/12 ⫽ 7 cycles. Structural images’ volume size is 256 ⫻ 256 ⫻ 54, with the voxel size 1 ⫻ 1 ⫻ 1 mm. Preprocessing Preprocessing was performed on both the null and the auditory experimental datasets. Images are first realigned with the least-squares approach and a six-parameter rigid transformation to reduce the effect of head motion. Then, a spatial smoothing is made to increase the signal-to-noise ratio with an isotropic Gaussian filter, in which the full width half-maximum (FWHM) was set to 6 mm. Last, time series were detrended using the first-order polynomial detrending method. The above-mentioned three procedures were all implemented within the SPM99 software. Simulated activation For the ROC analysis described below, simulated activations were added to the null data sets. The time series of the “activated” area was defined as a boxcar (starts with rest and alternates between rest and simulation, 10 sample points for both rest and stimulation, and totals 8 cycles) convolved with the hemodynamic response function (HRF). Here, the HRF was a combination of two ␥-functions (Friston et al., 1998a). Parameters of the two ␥-functions are the defaults of SPM99 software and the convolving of HRF with boxcar was implemented in the SPM99 software too. Then, the simulated activated time series were mixed with Gaussian white noise. Five simulated data sets were generated with the contrast-to-noise ratio (CNR) (Langer, 1996) varying among 0.2, 0.4, 0.6, 0.8, and 1.0, respectively. The corresponding amplitude of the activations was slightly smaller than a typical signal change observed in the BOLD contrast (Bandettini et al., 1992; Fadili et al. 2000; Filzmoser et al., 1999; Kwong et al., 1992). The volume size-of-the-presumed four activated regions corresponding to a 3 ⫻ 3 ⫻ 4, 4 ⫻ 4 ⫻ 4, 5 ⫻ 5 ⫻ 4, 6 ⫻ 6 ⫻ 4 cuboid, respectively (refer to Fig. 1).
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
457
the seeds by adding neighboring voxels that satisfy the homogeneity criterion. During the region selection process, some regions are selected from the obtained regions based on a predefined selection criterion. The rest of the regions are not considered.
Fig. 1. One slice of the simulated fMRI data. The white box represents the sites of artificial activated region. The top right-hand graph depicts the time series of an “activated” voxel, and the bottom right-hand graph displays the time series of a “nonactivated” voxel.
Region growing method The region growing method is a well-developed technique for image segmentation. It postulates that neighboring pixels within the same region have similar intensity values. The general idea of the region growing method is to group pixels with the same or similar intensities to one region according to a given homogeneity criterion. More precisely, the region growing method starts with a set of prespecified seed pixel(s) and grows from these seeds by merging neighboring pixels whose properties are most similar to the premerged region. Typically, the homogeneity criterion is defined as the difference between the intensity of the candidate pixel and the average intensity of the premerged region. If the homogeneity criterion is satisfied, the candidate pixel will be merged to the premerged region. The procedure is iterative: at each step, a pixel is merged according to the homogeneity criterion. This process is repeated until no more voxels are assigned to the region. In an in vivo fMRI experiment, the true neural activation typically tends to occur in a functional cluster (Katanota et al., 2002; Tononi et al., 1998). In other words, true fMRI activation is more likely to occur in clusters of some spatial connected voxels than a single voxel. Based on this feature, intuitively, we can utilize the idea of region growing method for fMRI activation detection. However, the following key questions remain: (1) How to set the homogeneity criterion? Which should be based on the spatial time series? (2) How to set the initial seed voxel(s)? Different from three-dimensional image, the fMRI data is typically a spatial-temporal signal, which is represented as a four-dimention array F ⫽ {fxyzt}MNOT, where (x, y, z) is the spatial 3D coordinates of brain voxels, and t is the time index. M ⫻ N ⫻ O denotes the total number of brain voxels in an image scan and T represents the total number of image scans. fxyzt indicates the image intensity of the brain voxel (x, y, z) at the tth instance of time. The proposed region growing method solves the abovementioned two questions successfully. It is divided into two procedures: the region growing and the region selection. In the region growing procedure, each voxel in the image was set as the initial seed and regions are separately grown from
Region growing During initialization, each voxel vi, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} is specified as the initial seed of region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. For each region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}, we can define a set E as the unassigned voxels which are boundary points of region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}, E ⫽ 兵q ⰻ A i 兩 N共q兲 艚 A i ⫽ 0其 where N(q) is the set of immediate neighbors of voxel q. For each voxel p in E, if the predefined homogeneity criterion is satisfied, voxel p is added to region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. This grown procedure is repeated until no more voxels are being assigned to region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. The final result is that for each voxel vi, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}, there is a corresponding grown region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. Region selection Set A ⫽ {A1, A2, . . . , AM⫻N⫻O}, S ⫽ null, i.e., S is defined as a null family of sets. Step 1: Find the maximal region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} in A, where maximal region indicates the region that has the largest voxel number. Reserve region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} in S, i.e., set S ⫽ S 艛 {Ai}, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. Then, regions derived from the seeds that belong to region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O} are removed from A. Mathematically, A is updated by set A ⫽ A ⫺ {Aj}J, J 傺 {1, 2, . . . , M ⫻ N ⫻ O}. Here, index J is defined as J ⫽ {j 兩 vj 僆 Ai}. Step 2: Repeat step 1 until the volume size of the remaining regions in A is smaller than a predefined threshold T2. This means that only clusters larger than T2 are significant. The rest are not considered. The final result of region selection procedure is saved in S. The number of regions in S is greatly smaller than M ⫻ N ⫻ O. During region selection procedure, why do we need to record the maximal region on each repeat step? It is known that in an in vivo fMRI experiment, the activation typically tends to constitute a cluster with a peak amplitude and a gradually descending scope around the peak. We can define the peak as the activation center, i.e., the point that has the strongest coherence with the simulation task. During the region growing procedure, if the initial seed is the activation center of activation cluster, the grown region will be the
458
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
maximal one. On the other hand, if the initial seed voxel deviates from the activation center, the grown region will be a smaller one. In each repeat step, it is obvious that only the maximal region can capture the feature of activated clusters.
Therefore, in each repeat step, only the maximal region was reserved. The implementation of region growing method for fMRI time series is summarized in the following pseudocode.
Region growing Specify each voxel as the initial seed of region A1, A2, . . . , AM⫻N⫻O For i ⫽ 1: M ⫻ N ⫻ O Repeat For each voxel p in the E set If the homogeneity criterion is satisfied Add p to region Ai End If End For Until no more voxels are being assigned to region Ai End For
Region selection Set A ⫽ {A1, A2, . . . , AM⫻N⫻O}, find the maximal region Ai in A, set S ⫽ S 艛 {Ai}, A ⫽ A ⫺ {Aj}J. Repeat. Find the maximal region Ai in A, set S ⫽ S 艛 {Ai}, A ⫽ A ⫺ {Aj}J. Until the volume size of remained regions in A is smaller than T2. Practical issues In the region growing procedure, it is critical to know how many and which neighboring voxels should be considered. For a 3D volume, 6-connected, 18-connected, or 26 connected neighborhood in three orthogonal directions can be considered. In this article, 26-connected neighborhood is selected. In the region growing procedure, the definition of homogeneity criterion is another critical problem. In this article it is defined as corr 共 f xyz:, M i兲 ⬎ T 1 where fxyz: is the time series of voxel (x, y, z), Mi is the averaged time series of voxels in the premerged region Ai, i 僆 {1, 2, . . . , M ⫻ N ⫻ O}. corr(,) is the Pearson correlation coefficient. T1 is the predefined threshold. The final result of the region growing method results in some clusters. Based on the clustered feature of fMRI activation, the task-related clusters are assuredly survived from the region selection procedure. That is to say, the result of region selection will not lose any task-related clusters. The problem is that besides task-related clusters, there are some other homogenous clusters reserved. To extract the task-related clusters, we have to distinguish among responses of the different anatomical structures. The procedure mentioned below will help us to this end. As a data-driven method, it is clear that region growing method cannot stand for a fMRI activation detection step. Some
additional treatments are necessary to extract the task-related clusters. If paradigm information can be utilized, one could test for the significant periodicity via Fourier transform methods (for periodic paradigm) or for the significant correlation to the predefined reference (Fadili et al., 2000). If the paradigm information cannot be utilized, some a priori neurophysiological and anatomical knowledge must be employed. In this article, a test of significant correlation was exploited. First, a reference is specified. Then, Pearson correlation coefficient was calculated between the reference and each cluster’s averaged time series. Clusters with a correlation coefficient greater than the predefined threshold TRG is considered task-related. In addition, the task-related clusters may be intersectant. Such intersectant clusters are merged into a single cluster. The final result of this procedure results in unconnected clusters (components), which are what we expect. ROC analysis on simulated data We conducted a receiver operating characteristic (ROC) analysis on the simulated data sets to evaluate the effectiveness of the region growing method and to compare its performance with that of other methods. The following three statistical models are compared: (1) the deconvolution method, (2) fuzzy c-mean clustering analysis (FCA), and (3) the region growing method. The reason the deconvolution and the FCA method were chosen for comparison is twofold. On one hand, the deconvolution method is a popular hypothesis-driven method for fMRI activation detection (Ward, 2002). On the other hand, FCA has been proven to be a successful data-driven analysis method in analyzing fMRI data (Baumgartner et al., 1997, 2000; Moser et al., 1997). The deconvolution method is an optimal program of AFNI software (Cox, 1996; Ward, 2002). It first estimates
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
the impulse response function(s). The estimated impulse response function(s) is then convolved with the stimulus paradigm to yield the estimated response. Then, various statistics are calculated to indicate the “goodness” of the fit. Finally, according to a predefined level of statistics, those activated voxels can be identified. FCA partitions the volumetric data into a number of fuzzy clusters. To obtain the task-related clusters, Pearson correlation coefficients were calculated between a predefined reference and the centroid of each cluster. Clusters with a correlation coefficient greater than the predefined threshold TFCA are considered task-related. For each cluster, FCA produces a cluster membership map for the volumetric data. By thresholding a membership map at a predefined level of threshold, those voxels that belong to the cluster can be identified. This is analogous to- but not equivalent to, selecting a significant level of a statistical test. To give an ROC analysis, four steps were processed: First, the three methods were applied to the simulated datasets. Second, according to the known number of simulated activation, the true positive ratio in the activated region and the false positive ratio in the nonactivated region were calculated. Third, an ROC curve was calculated by plotting the true positive ratio on the false positive ratio. The curve corresponding to a certain method closest to the upper left corner should be the best. Fourth, as the conventional ROC analysis, the area under the ROC curve was taken as the detectability measure of the different methods. In the implementation of the deconvolution method, the maximum time lag is set according to the known HRF, and the degree of the polynomial in the baseline model is set to 1. For FCA, the fuzziness exponent weight is set to 2 according to previous work of Fadili (2000). There is no theory to determine the number of clusters prior to FCA. In this article, for each of the five simulated datasets, the number ranging between 2 and 30 was set as the candidate number of clusters. The number that obtained the best results, i.e., the maximal area under the ROC curve, was recorded as the cluster number. In addition, the distance measure is the hyperbolic correlation measure (Golay et al., 1998),
459
Analysis of in vivo data sets The auditory dataset was analyzed with the three methods referred to above. For the decovolution method, the maximum time lag is set to 5 and the degree of the polynomial in the baseline model is set to 1. For the FCA method, the fuzziness exponent weight was set to 2, and the number of clusters was set to 30 (this is a typically used value for FCA) (Baumgartner et al., 2000). The hyperbolic correlation measure is utilized as the distance measure. To overcome the ill-balanced data problem, prior to FCA, noisy pixels were removed by preprocessing (for all voxel time series, we calculated their 1-lag-shifted autocorrelation function, and P ⫽ 0.01 was selected for the significance level) (Fadili et al., 2000). For the region growing method, T2 is set to 4. TRG and TFCA are all set to 0.5. To extract the task-related clusters of FCA and the region growing method, the reference time series is the paradigm convolved with the HRF (combinations of two ␥-functions, SPM99). For all of the three methods, only the brain voxels are considered, which are obtained from a low-level threshold estimated from histogram. In addition, for deconvolution and FCA method, only the clusters with a cluster size larger than four voxels were considered significant. The rest of the clusters were regarded as noise. The aim of this is to stay consistent with the T2 threshold of the region growing method. To make a comparison among the three methods, a quantitative method was implemented. The focus of quantitative comparison is to select the threshold of different activation maps, as the threshold determines whether a voxel is activated or not. Moreover, the threshold determines the fairness of comparison. The scheme of the proposed quantitative comparison method is as follows: (1) some or other activated area obtained with the deconvolution method is selected as the reference area; (2) thresholds of the deconvolution (P value), the FCA (membership value), and the region growing (T1) method are adjusted by keeping the same cluster size of the reference area; (3) the volume size and the averaged time series of the nonreference area are taken as measures of effectiveness and sensitiveness of different methods.
d共 x k,v i兲 ⫽ 共1 ⫺ c ik兲/共1 ⫹ c ik兲, Results where cik denotes the Pearson correlation coefficient between the time series xk and cluster centroid vi. The reference time series is the simulated activated time series without mixing Gaussian noise. TFCA is set to 0.5. For the region growing method, the reference time series is same with that utilized in FCA. TRG is set to 0.5 (this corresponding to TFCA ⫽ 0.5). In addition, T2 is set to 2. Accordingly, for the other two methods, only clusters with a cluster size larger than two voxels were considered significant. The rest clusters were regarded as noise and hence will not be considered.
Simulated data Fig. 2 shows the area under the ROC curve for each of the three statistical models plotted as a function of CNR. We can conclude that the region growing method outperforms the deconvolution and the FCA method. To further clarify the increased sensitivity of the region growing method, we examine one particular case, with CNR at 0.6. The ROC curves of the three methods are depicted in Fig. 3a. The characteristics of the three models are well illustrated in this
460
Y. Lu et al. / NeuroImage 20 (2003) 455– 465 Table 1 Computation time in seconds for the three different methods in simulated dataset with different CNR Method
Deconvolution FCA Region growing
CNR 0.2
0.4
0.6
0.8
1.0
⬍30 4816 467
⬍30 3709 432
⬍30 2975 403
⬍30 3610 496
⬍30 3122 415
Note. FCA is calculated in the optimal cluster and the region growing method is calculated with T1 value that corresponds to the optimal classification point.
Fig. 2. Performance of the region growing method, the deconvolution method, and FCA method for the simulated data set.
figure. The ROC curves of the region growing method approached the top left corner primely, whereas the other two methods were far from the top left corner. In Fig. 3a, point x, y, and z are the points closest to the optimal classification (False Positive Rates ⫽ 0, True Positive Rates ⫽ 1) in the respective ROC curves. Figs. 3b, c, and d illustrate the activation maps obtained with the region grow-
ing method, the deconvolution method, and FCA, respectively. The corresponding thresholds were selected according to the value that corresponds to the optimal classification point x, y, and z. From this comparison, we can conclude that the region growing method is more sensitive than the FCA and the deconvolution method for activation detection. The computational complexities of the three methods are reported in Table 1. The computations were performed on a 1.8-GHz PC, and the region growing method, FCA, and the deconvolution method are all coded in Matlab (The MathWorks, Inc.). From this comparison, we can conclude that
Fig. 3. (a) ROC curves corresponding to the case of CNR ⫽ 0.6 for the simulated data. (b) Results of the region growing method. (c) Results of the deconvolution method. (d) Results of FCA.
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
461
Fig. 4. Averaged time series of the 18 clusters and the corresponding correlation coefficient between the reference and the averaged time series.
the region growing method has less computational complexity than the FCA. Auditory experiment To demonstrate how the region selection procedure affects the performance of the method, we set T1 to 0.95 TRG to 0.5 and T2 to 4. For simplicity, only the clusters in the 29 –32th slices of the total 64 slices are computed. There are a total of 133 clusters that survived from the region selection procedure. Thirty-two of the 133 clusters are taskrelated. The 32 clusters are merged into 2 larger clusters. The other 101 clusters are merged into 16 larger clusters. In Fig. 4 the averaged time series of the 18 clusters and the correlation coefficients with reference time series are shown. In Fig. 5, the two task-related clusters are shown in red; the other 16 clusters are shown in green. The region growing method provides, in addition to the activation clusters, a number of other clusters. Most of these are dominated by different sources of noise. In the context of fMRI data, how to interpret the huge amount of information provided by the data-driven methods is still a major concern. We refer to Fadili (2000) for an example on this streamline. We will not consider this issue here. The right auditory cortex in the 29th slice was set as the reference area. For the decovolution method, four thresh-
olds are selected corresponding to the P-value of 10⫺8, 10⫺7, 10⫺6, and 10⫺5. The corresponding threshold T1 of region growing method and membership value of FCA are listed in Table 2. Results of the three methods are shown in Fig. 6. All three methods can successfully detect the other activated areas such as Brodmann’s area 42 (primary auditory cortex) as well as the Brodmann’s area 22 (auditory association area). The detailed volume size of the activated region in each slice is listed in Table 2. Compared with the FCA method, the activated areas revealed by the region growing method are larger and more continuous. For example, the activation in the left auditory cortex is larger and continuous. In addition, small activated clusters scattered throughout the brain are fewer in the result obtained with the region growing method than that obtained with the FCA method. In Fig. 6d, for instance, there are some small false activated clusters with the FCA method. Compared with the deconvolution method, the cluster size detected with the region growing method is bigger and more continuous. In the left auditory cortex activated area, there are a total of 233 voxels that both activated in region growing and deconvolution method; the averaged time series of these 233 voxels is shown in Fig. 7a. There are 125 voxels that were activated in the region growing method but were not activated in the deconvlolution method. The averaged time series of these 125 voxels is shown in Fig. 7b.
462
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
Fig. 5. Clusters in the 29 –32th of the 64 axial slices. Red: task-related clusters; green: non-task-related clusters. The left side corresponds to the right hemisphere Fig. 6. Activation maps revealed by each of the three methods. Results of the three methods, i.e., the deconvolution method, the FCA, and the region growing method, are shown from the top row to the bottom of the figure. Four slices (29 –32th of the total 64 slices) involved in the auditory experiment are shown. The left side corresponds to the right hemisphere, a, b, c, and d are the results of the three methods that correspond to the deconvolution method thresholded at the voxelwise significance level of P ⬍ 10⫺8, P ⬍ 10⫺7, P ⬍ 10⫺6, and P ⬍ 10⫺5, respectively.
The correlation coefficient between this two averaged time series is 0.938. Therefore, we can conclude that the 125 voxels are all meaningful activated voxels. In the right auditory cortex activated region, there are totally 328 voxels that both activated in the deconvolution, and the region growing method; averaged time series of these 328 voxels is shown in Fig. 7c. There are 62 voxels that activated in region growing but nonactivated in the deconvolution mehod. Averaged time series of the 62 voxels is shown in Fig. 7d. The correlation coefficient between this two mean time series is 0.826. Therefore, we can conclude that the 62 voxels are all meaningful activated voxels. In this data set, these results indicate that the region growing method is more sensitive than the deconvolution method.
Discussion and conclusion Our ROC analysis is based on several assumptions, as other simulation studies in general. First, the property of true activated area is known a priori, as it is artificially simulated. This is the reason that we can get the optimal cluster number for FCA based on the area under the ROC curve. Although the optimal cluster number is used for FCA, the detection power is smaller than the region growing method. For the deconvolution method, though the linearity assumption is satisfied, the detection power is smaller than the region growing method. This indicates that the region growing method is not very sensitive to noise. Another issue is the fMRI noise. We used a computer-simulated Gaussian
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
463
Table 2 Volume size (voxel number) for different methods Method
Threshold
29(R)
30(R)
31(R)
32(R)
29(I)
30(I)
31(I)
32(I)
45 46 52
26 31 30
35 41 41
41 53 61
44 78 71
51 54 60
30 35 31
39 45 44
53 56 77
50 88 96
58 65 67
38 43 38
44 54 62
59 82 84
67 99 110
67 72 76
45 81 48
50 86 76
68 107 108
81 115 126
a Deconvolution FCA Region growing
10⫺8 0.19 0.895
97 98 96
66 77 77
52 53 54 b
Deconvolution FCA Region growing
⫺7
10 0.15 0.88
102 102 103
82 82 89
58 66 61 c
Deconvolution FCA Region growing
10⫺6 0.09 0.865
106 105 105
93 88 108
67 73 68 d
Deconvolution FCA Region growing
⫺5
10 0.02 0.85
117 115 116
102 94 116
78 89 82
Note. L: left; R: right. The threshold for the deconvolution, the FCA, and the region growing method is P value, membership, T1, respectively. Number with shading is the size of reference area. a, b, c and d correspond to the result of Figs. 6a, b, c, and d, respectively.
white noise. In practice, it is different from the experimental data sets. In addition, the shape of simulated activated clusters would be more complicated in reality. It is essential, therefore, to examine the effectiveness of our model using real experimental data. In this study, we employed the auditory data set. Compared with the other two methods, the region growing method has the three following advantages: (1) Volume size of the activated clusters is larger; (2) Activated clusters are more continuous; (3) Small scattered clusters are fewer. From the viewpoint of methodology, compared with the deconvolution method, the region growing method has three advantages: (1) It does not require any prior knowledge of the
HRF, which may not be generally known. (2) It does not need the linear assumptions. The evidence of linearity of the hemodynamic response should contrast with the work of Vazquez and Noll (1998), Friston et al. (1998b), and Glover (1999). (3) As a univariate technique, the deconvolution method ignores the relationship of neighboring voxels. Therefore, the deconvolution method is hindered from the feature of clustered activation. However, the region growing method can utilize the relationships of neighboring voxels. Compared with FCA, the region growing method has the following three advantages: (1) It does not need the predefined cluster number. (2) It is not sensitive to the ill-balanced data problem, so we minimize the risk of discarding possibly activated voxels. (3) FCA ignores
Fig. 7. (a) Averaged time series of the 233 voxels that activated in both region growing and the deconvolution in the left auditory cortex. (b) Averaged time series of the 125 voxels that activated in region growing but nonactivated in the deconvolution in the left auditory cortex. (c) Averaged time series of the 328 voxels that activated in both region growing and the deconvolution in the right auditory cortex. (d) Averaged time series of the 62 voxels that activated in region growing but nonactivated in deconvolution in the right auditory cortex.
464
Y. Lu et al. / NeuroImage 20 (2003) 455– 465
the spatial relationships of neighboring voxels. However, the region growing method can utilize the relationships of neighboring voxels. In addition, region growing can serve as a data reduction preprocessing technique. As mentioned above, fMRI activation is more likely to occur in clusters of several spatial contiguous voxels than in a single voxel. The implementation of data reduction procedure by utilizing region growing is outlined in the following two steps: First, each voxel in the brain was set as the initial seed and regions are separately grown from the seed by adding neighboring voxels that satisfy the homogeneity criterion. Second, the voxels whose corresponding grown regions smaller than a predefined threshold, for instance, 4 – 6, are considered noise and will not participate in the subsequent analysis. This data reduction procedure can be applied to existing fMRI data analysis method, for example, GLM, the deconvolution, and FCA. A limitation of the region growing method is its computation complexity for larger activated regions. For instance, the computation time of FCA (recently, the application of FCA was expanded by introducing two additional stages (Jarmasz and Somorjai, 2002), whose effect is to speed up the computation and make the results more robust. In this article, we only utilized the typical FCA method for comparison) for the block data is 8676 s, but the computation time of the region growing method (with T1 ⫽ 0.85) is 15,442 s. The computation burden of the region growing method lies in the calculation of the region growing part. There are two feasible schemes to reduce the computation complexity. The first way is anatomical restriction. fMRI BOLD activation is thought to occur mainly in the gray matter. Thus, if only the voxels belonging to this issue are considered, the computation complexity will be reduced. The segmentation strategy proposed by Fadili (2000) can be applied for this purpose. The second way to reduce the computation complexity is to introduce parallel computation techniques. It is very easy to implement in a parallel way since each region can be grown separately and simultaneously on different CPUs or different computers. In summary, a new method has been presented for fMRI data analysis. An ROC analysis has demonstrated its effectiveness for fMRI activation detection. Its usefulness has been confirmed by applying it to a real experimental. The method can serve as a new reliable data-driven method suited to the nature of fMRI data analysis.
Acknowledgments The authors are grateful to the anonymous referees for significant and constructive comments and suggestions, which greatly improved the article. The corresponding author also thanks his colleague, Professor Ricardo Vilalta, Department of Computer Science, University of Houston,
for carefully checking the use of English in the manuscript. The authors give thanks to Professor Karl Friston for his kind permission to use the block trial data sets. This work was partially supported by Hundred Talents Programs of the Chinese Academy of Sciences, the Natural Science Foundation of China, Grant Nos. 60172056 and 697908001.
References Bandettini, P.A., Wong, E.C., Hinks, R.S., Tikofsky, R.S., Hyde, J.S., 1992. Time course EPI of human brain function during task activation. Magn. Reson. Med 25, 390 –397. Baumgartner, R., Ryner, L., Richter, W., Summers, R., Jarmasz, M., Somorjai, R., 2000. Comparison of two exploratory data analysis methods for fMRI: fuzzy clustering vs. principal component analysis. Magn. Res. Imaging 18, 89 –94. Baumgartner, R., Scarth, G., Teichtmeister, C., Somorjai, R., Moser, E., 1997. Fuzzy clustering of gradient-echo functional MRI in the human visual cortex: Part I: Reproducibility. J. Magn. Res. Imaging 7 (6), 1094 –1101. Calhoun, V.D., Adali, T., Pearlson, G.D., Pekar, J.J., 2001. Spatial and temporal independent component analysis of functional MRI data containing a pair of task-related waveforms. Hum. Brain Mapping 13, 43–53. Cox, R.W., 1996. AFNI: software for analysis and visualization of functional magnetic resonance neuoimages. Comp. Biomed. Res 19, 162– 173. Fadili, M.J., Ruan, S., Bloyet, D., Mazoyer, B., 2000. A multistep unsupervised fuzzy clustering analysis of fMRI time series. Hum. Brain mapping 10, 160 –178. Filzmoser, P., Baumgartner, R., Moser, E., 1999. A hierarchical clustering method for analyzing functional MR images. Mag. Res. Imaging 17, 817– 826. Friston, K.J., Fletcher, P., Josephs, O., Holmes, A., Rugg, M.D., Turner, R., 1998a. Event-related fMRI: characterizing differential responses. Neuroimage 7 (1), 30 – 40. Friston, K.J., Holmes, A.P., Worsley, K.J., Poline, J.P., Frith, C.D., Frackowiak, R.S.J., 1995. Statistical parametric maps in functional imaging: a general linear approach. Hum. Brain Mapping 2, 189 –210. Friston, K.J., Josephs, O., Ress, G., Turner, R., 1998b. Nonlinear eventrelated responses in fMRI. Magn. Res. Med 39, 41–52. Glover, G.H., 1999. Deconvolution of impulse response in event-related BOLD fMRI. Neuroimage 9, 416 – 429. Golay, X., Kollias, S., Stoll, G., Merier, D., Valavanis, A., Boesiger, P., 1998. A new correlation-based fuzzy logic clustering algorithm for fMRI. J Magn. Res. Med 40, 249 –260. Goutte, C., Toft, P., Rostrup, E., Nielsen, F.A., Hansen, L.K., 1999. On clustering fMRI time series. Neuroimage 9 (3), 298 –310. Hansen, L.K., Larsen, J., Nielsen, F.A., Strother, S.C., Rostrup, E., Savoy, R., Lange, N., Sidtis, J., Svarer, C., Paulson, O.B., 1999. Generalizable patterns in neuroimaging: how many principal components? Neuroimage 9, 534 –544. Jarmasz, M., Somorjai, R.L., 2002. Exploring regions of interest with cluster analysis (EROICA) using a spectral peak statistic for selecting and testing the significance of fMRI activation time-series. Artif. Intell. Med 25, 45– 67. Jung, T.P., Makeig, S., Mckeown, M.J., Bell, A.J., Lee, T.W., Sejnowski, T.J., 2001. Imaging brain-dynamics using independent component analysis. Proc. IEEE 89 (7), 1107–1122. Katanoda, K., Matsuda, Y., Sugishita, M., 2002. A spatial-temporal regression model for the analysis of functional MRI data. Neuroimage 17, 1415–1428. Kwong, K.K., Belliveau, J.W., Chesler, D.A., Goldberg, I.E., Weisskoff,
Y. Lu et al. / NeuroImage 20 (2003) 455– 465 R.M., Poncelet, B.P., Kennedy, D.N., Hoppel, B.E., Cohen, M.S., Turner, R., Cheng, H.M., Brady, T.J., Rosen, B.R., 1992. Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc. Natl. Acad. Sci. USA 89, 5675–5679. Lai, S.H., Fang, M., 1999. A novel local PCA-based method for detecting activation signals in fMRI. Magn. Res. Imaging 17 (6), 827– 836. Langer, N., 1996. Tutorial in biostatistics: statistical approaches to human brain mapping by functional magnetic resonance imaging. Stat. Med 15, 389 – 428. McKeown, M.J., Makeig, S., Brown, G.G., Jung, T.P., Kindermann, S.S., Bell, A.J., Sejnowski, T.J., 1998. Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapping 6 (3), 160 –188. Morel, J.M., Solimini, S., 1995. Variational Methods in Image Segmentation. Birkhauser Publishing, Inc., Basel. Moser, E., Diemling, M., Baumgartner, R., 1997. Fuzzy clustering of
465
gradient-echo functional MRI in the human visual cortex. Part II: Quantificaltion. J. Magn. Res. Imaging 7, 1102–1108. Ogawa, S., Lee, T.M., Ray, A.R., Tank, D.W., 1990. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc. Natl. Acad. Sci. USA 87, 9868 –9872. Scarth, G., Somorjai, R., Alexander, M., Wowk, B., Wennerberg, A., McIntyre, M., 1995. Cluster analysis of functional brain images, in: Proceedings of the First International Conference on Functional Mapping of the Human Brain, Vol. S1 of Hum. Brain Mapping 158. Tononi, G., McIntosh, A.R., Russell, D.P., Edelman, G.M., 1998. Functional clustering: identifying strongly interactive brain regions in neuroimaging data. Neuroimage 7 (2), 133–149. Vazquez, A.L., Noll, D.C., 1998. Nonlinear aspects of the BOLD response in functional MRI. Neuroimage 7, 108 –118. Ward, B.D., 2002. Deconvolution analysis of fMRI time series data, http://www.afni.nimh.gov/afni/edu.