Depression recognition using resting-state and event-related fMRI signals

Available online at www.sciencedirect.com

Magnetic Resonance Imaging 30 (2012) 347 – 355

Depression recognition using resting-state and event-related fMRI signals☆ Qing Lua,⁎, Gang Liua , Jing Zhaoa , Guoping Luoa , Zhijian Yaob,⁎ a Research Centre for Learning Science, Southeast University, Nanjing, 210096, China Academic Department of Psychiatry, Nanjing Brain Hospital, Nanjing Medical University, Nanjing, 210029, China Received 6 July 2011; revised 10 October 2011; accepted 4 December 2011

b

Abstract Purpose: This paper aimed to develop a method for depression detection using blood-oxygen-level-dependent (BOLD) response estimated from event-related signals and resting-state functional magnetic resonance imaging (fMRI) signals together. Materials and Methods: Thirteen patients with unipolar depression and matched healthy subjects were recruited. Resting state data of each subject were collected. Thereafter, event-related paradigm was undertaken using sad facial stimuli. The resting-state fMRI signal was deemed as the baseline of each subject's activity. Coefficient marks were designed to sort and select temporal independent components of eventrelated signals. Thereafter, stimulus-evoked BOLD response components inside event-related signal were extracted and taken as features to discriminate depressive patients from healthy controls. Results: Accuracy rate for depression recognition was 77.27% with P value of .017 for whole-brain analysis and 81.82% with P value of .009 for region-of-interest analysis. The effectiveness and the superiority of the proposed method for disease recognition were demonstrated via the performance comparison with three other typical methods. Conclusions: The proposed model was effective in depression recognition. © 2012 Elsevier Inc. All rights reserved. Keywords: Depression; Classification; Resting state; Event related; Functional magnetic resonance imaging (fMRI)

1. Introduction Functional magnetic resonance imaging (fMRI) can detect in real time the local hemodynamic changes which are associated with cerebral cortex neural activation [1]. Comparisons of the activation in brain regions have revealed significant differences between depressive patients and healthy controls. Based on group-level statistic discriminative analysis, the pattern recognition problem of ‘automatically separating the depressive from the healthy’ was proposed as an essential step toward a clinical application. Typical work can be found in Ref. [2], where the bloodoxygen-level-dependent (BOLD) signals evoked via sad facial stimuli were utilized to classify depressive patients ☆ Competing interests statement: The authors declare that they have no competing financial interests. ⁎ Corresponding authors. E-mail addresses: [email protected] (Q. Lu), [email protected] (Z. Yao).

0730-725X/$ – see front matter © 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2011.12.016

from healthy volunteers. Other explorations trying to recognize various types of psychiatric disorders have been undertaken as well, such as Refs. [3–5]. In order to estimate BOLD signal for disease recognition, the most commonly used approach is the general linear model, which requires the specification of an expected hemodynamic response function (HRF). This function is defined as the expected time course of the BOLD signal after the onset of a stimulus. Related work can be found in Ref. [2]. The authors utilized a convolution model and predesigned HRFs to estimate BOLD response for each event. Thereafter, the estimated time series were employed to provide a classification. Other related works using the statistically calculated contrast images were reported as well. They involved the following features: the number of activated voxels, the value of primary peak, the size of cluster where the statistical significant voxel belongs to and the number of significant clusters. However, in case the task is cognitively complex or in cases of diseases, the misspecification of the HRF may reduce the reliability of results.

348

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

In order to eliminate the influence from model, volume calculation, such as averaging volumes within the event minus a control volume, is another direct and simple substitute for BOLD response modeling [2]. Other methods, such as independent component analysis (ICA) [6,7], cluster analysis [8–10] and wavelet correlation analysis [11], provide more flexible alterations. However, limitations still exist among them. For example, component selection problem makes the ICA method more suitable for the paradigm of block design, where BOLD response has distinctive morphological features. Correlation analysis between separated components and stimuli series may be a possible way to find functional components in these cases. However, indeterminate time delay inherent in BOLD response might deteriorate the results. Wavelet correlation analysis aims to identify regions with common activation signals across subjects. Therefore, critical subject variability for disease detection is impaired. Similar problems are inevitable in the volume averaging approach. This paper aimed to develop a method for depression detection using BOLD response estimated from eventrelated signals and resting-state fMRI signals together. The signals under event-related design and resting-state fMRI signals were collected from the same subject. The restingstate fMRI signals were considered for the purpose of two aspects. Firstly, resting-state ICA components helped schedule event-related ICA components. Secondly but potentially more importantly, the integration of restingstate information can help eliminate the effect of individual variability, such as different brain activation patterns of different individuals. The work was based on an assumption that the event-related fMRI signals under stimuli were a linear combination of the evoked BOLD response and the underlying brain activities in the resting state. It was well founded via Ref. [12]. Thereafter, the selected BOLD response signals were utilized as input features for depression recognition purpose. Its classifying performance was compared with that of three other typical event-related temporal feature extraction methods. 2. Methods and materials 2.1. Subjects All depressed patients were recruited from inpatient facilities. Thirteen patients with unipolar depression were recruited. Eligibility screening procedures included the Structured Clinical Interview for the DSM-IV (SCID), the Brief Psychiatric Rating Scale, the 24-item Hamilton Depression Rating Scale and common clinical laboratory tests. The initial diagnoses of depression were made by the participants' treating psychiatrists and confirmed by an expert psychiatrist according to SCID. Patients were required to complete a minimum 2-week medication washout prior to MRI scans. Patients with other psychiatric illness and a history of electroconvulsive therapy were excluded.

Nine healthy subjects were recruited from a roll living in the same or nearby places with depressed patients by print advertisements. Each healthy subject was required to match with patients in gender, age and education level (Table 1). They all did not have any psychiatric illness presently or a history of psychiatric illness. All subjects were right-handed and satisfied the criteria to undergo an MRI scan. Exclusion criteria for all subjects included the following: subjects with serious medical and neurological disorders, a family history of serious psychiatric or neurological illness in first-degree relatives, abnormal clinical laboratory tests, acutely suicidal or homicidal, present or previous substance abuse and substance dependence, currently taking any prescription, currently pregnant or breastfeeding. After a complete description of the study to all subjects, written informed consents were obtained. The study was approved by the Research Ethics Review Board. 2.2. Implicit sad facial affect recognition task Sad, happy, food chewing and neural emotional faces were performed by several young students from the Academy of Art. Emotional stimuli and baseline trials (cross-hair fixation) were presented randomly in video form. Food chewing was embedded here in order to exclude the influence of muscular movement on the emotional recognition. Each type of stimuli was presented 20 times (10 stimuli were performed by female students, and the rest were performed by male students) for a total of 100 trials. Each trial was performed for 3 s, and the intertrial interval was randomly varied according to a Poisson distribution, with mean intertrial interval of 4 s. Total duration was 6 min and 39 s. When scanned, the subjects were asked to identify whether the stimuli was sad or not by button press. 2.3. Data acquisition Imaging data were acquired using a GE Signa 1.5-T MRI scanner. All subjects were placed in a birdcage head coil and fitted to a foam padding to reduce head motion. The fMRI sequences were as follows: 1. Anatomic scan: three-dimensional T1-weighted axial images: repetition time/echo time (TR/TE)=15 ms/4.4 ms, thickness=1.0 mm, flip angle=15°, matrix=256×256, field of view (FOV)=240×240 mm 2. 2. Resting-state and functional scans: TR/TE=3000 ms/40 ms, thickness=4.0 mm, flip angle=15°, matrix=128×128, Table 1 Demographic features

Age (years) HAMD Gender Education

Depressed patients

Healthy control subjects

40.1 (12.1) 40.9 (10.5) Female 12.8 (3.4)

42.5 (11.9) Female 13.6 (4.3)

Mean scores for each variable are presented with standard deviations in parentheses. HAMD: Hamilton Depression Rating Scale.

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

FOV=240×240 mm 2, 27 sequence slices for both type of scans, time: 6 min and 39 s.

NeuroImages software [13] (http://afni.nimh.nih.gov), images were bandpass filtered (0.01–0.08 Hz) to reduce low-frequency drift and physiological high-frequency noise. Linear trend was removed to reduce the influence of the rising temperature of the MRI equipment.

2.4. Data preprocessing 2.4.1. Event-related data preprocessing Data were realigned to remove residual motion effects, spatially normalized into standard space and smoothed with a Gaussian kernel of 6 mm full width at a half maximum using SPM2 (Wellcome Trust Centre for Neuroimaging, London, United Kingdom).

2.5. BOLD response estimation for disease recognition The main purpose of this part is to model the evoked BOLD response inside event-related signals, depending on the baseline provided via resting-state signals. Temporal ICA was applied separately to decompose each type of signals first. After that, the ICA components of these two types of signals were put together to sort out the signal components, which were associated with the functional activations. A diagram description was shown in Fig. 1.

2.4.2. Resting-state data preprocessing Slice timing, head motion correction and spatial normalization were performed using SPM2. Through above processes, images were normalized to the standard SPM2 echo planar imaging template. Using Analysis of Functional

Subject

Event-related signal

4

4

2

2

0

0

-2

-2

-4 0

Resting-state signal

-4 40

80

120

0

40

80

120

PCA

PCA

ICA Sei(t) (i=1,2,…,25)

ICA Sri(t) (i=1,2,…,25)

Calculate coefficient CCi mark for Sei T

∑ ( Se (t ) − Se )( Sr (t ) − Sr ) i

d ij =

i

j

j

t =1

T

∑

349

T

( Sei (t ) − Sei ) 2

t =1

∑ ( Sr (t ) − Sr ) j

2

j

t =1

CCi = arg max(d ij ) Sr j

Select component Sei with small CCi, determine the BOLD response component Fig. 1. Diagram for one subject's BOLD response component extraction.

350

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

The preprocessed fMRI volumes were divided into 116 regions according to the anatomically labeled template via [14]. The data representation process was carried out in the regional level. That is, only the ICA components of resting-state signals and event-related signals from the same regions can be taken together for further comparison. Before ICA analysis, principal component analysis (PCA) was performed on the signals from each anatomical region of each subject separately. PCA helps find the information that explains the greatest amount of variance of the data. Data reduction was performed ensuring that more than 90% signal information was retained. Data reduction helps ensure the efficiency and stability of ICA process. Thereafter, temporal ICA was applied separately to discover temporal independent components of two types of data from each subject. Let a matrix X denote the observed event-related or resting-state data, where the row dimension was the number of principal components selected via PCA from the related region and the column dimension was the number of sampling points. The observed signal X was written as: X ¼ AS

ð1Þ

A was the mixing matrix, and S was the source matrix. Taking the minimization of mutual information as the goal, a matrix W was found, so that 1

S ¼ WX

ð2Þ

1

S was the estimate of the source signal distribution in the brain. Thereafter, the temporal independent components over the brain could be estimated. Each anatomical region had 25 components kept for further analysis. The ICA process was implemented via fastICA [15] (http://www.cs. helsinki.fi/u/ahyvarin/papers/fastica.shtml). Correlation analysis between every event-related independent component and every resting-state independent component was carried out then. Let dij (i=1,2, ..., M; j=1,2, ..., N) record the absolute value of the correlation coefficient between event-related component Sei(t) and resting-state component Srj(t) from the same region. Here, M was the number of event-related independent components, and N was the number of resting-state independent components. M and N were equal and set as 25 in this study. Absolute correlation coefficient dij was calculated via following equation: T  P Z  Z Sei ðt Þ− Sei Sr j ðt Þ− Sr j t¼1 ffi d ij ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T  T  P Z 2 P Z 2 Sei ðt Þ− Sei Sr j ðt Þ− Sr j t¼1

t¼1

ð3Þ

― ― where Sei and Srj were the average signals over the time course of Sei(t) and Srj(t). Set the coefficient mark for each event-related component Sei(t) as:   ð4Þ C C i ¼ arg max d ij Sr j

That is, for a certain event-related component, the coefficient mark was set via the maximum value of the correlation coefficients between this event-related component and all resting-state components. We sorted all eventrelated independent components Sei according to their coefficient mark CCi in an incremental order. It was more likely that the top event-related components, which had relatively low coefficient marks, were the components containing dominating evoked BOLD response information. The related reasons were analyzed as follows. Early fMRI studies have already described signal fluctuations that were closely related to physiological motions like respiration and cardiac action, as well as lowfrequency oscillations. Image-to-image fluctuation due to physiological motion is a major limitation to an accurate detection of neuronal activity with fMRI. The respiration and heart beat during the acquisition of imaging data are the primary source of such physiological noise inherent in the fMRI signals. Although regions significance and magnitude of cardiac- and respiratory-related signal variance exist, similar patters can be found between the event-related scanning and the resting-state scanning for an identical subject. It was also well founded via Ref. [12] that measured event-related fMRI signal is a combination of task-related signal in addition to low-frequency (‘resting state’) intrinsic fluctuations. As a result, if event-related signals and restingstate signals were decomposed by ICA, respectively, the independent components for physiological motions and the underlying brain activates in the event-related signals were able to find their corresponding components in the restingstate signals. Under such a situation, the correlation coefficient between these similar components were relatively large and assigned as coefficient marks CC for the corresponding event-related components, according to Eq. (4). Due to the similar pattern of these components in both types of signals, it can be inferred that these coefficient marks would be larger than those of the independent components for BOLD response and random noise. As a result, if the event-related components were sorted increasingly according to their coefficient marks, the top components might contain the dominating evoked BOLD response information or relate to uncorrelated random noise. In order to ensure that the selected components described BOLD response mainly but not random noise, further spectrum analysis was carried out on these event-related components. Generally speaking, the distributed energy is concentrated within certain frequency range for BOLD response signals, while even distribution is common for random noise. Hence, the components of BOLD response and random noise can be easily distinguished via their

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

spectral characteristics. Fig. 2 gave an illustration of spectral characteristics for BOLD response signal and random noise, respectively. As shown in Fig. 2(A), the BOLD response signal of a typical activated voxel was reconstructed after the individual statistical analysis inside SPM software using the event-related scanning of an arbitrarily selected subject under our paradigm. Different from the even energy distribution in Fig. 2(B) for random signal, the typical BOLD response signal under emotional stimuli had its energy concentrated within the range less than 0.1 Hz. Such phenomenon was exactly consistent with the characteristics of BOLD signals. Hence, the procedure for the BOLD response estimation can be summarized as follows: Step 1: Apply PCA method on event-related data and resting-state data, respectively, in each region to reduce the data dimension from brain voxel level to eigenvalue level. Step2: Calculate the temporal independent components Sei(t) and Srj(t) for event-related data and resting-state data, respectively. Step 3: Calculate coefficient mark CCi for each eventrelated components using Eq. (4) and sort the components in an increasing order of CCi. Step 4: Study the spectrum analysis results for the top event-related components and select the first component, which shows uneven energy distribution over its frequency domain, as the major BOLD response component for the related region. We have also investigated the signal characteristics of the frequency spectrogram after second-order polynomial fitting. The quadratic coefficient of the polynomial function represents the energy distribution pattern. The larger the coefficient, the more concentrated the signal energy is. If the

351

coefficient for the quadratic term is close to zero, it means an even distribution of signal energy. A random test was carried out on 10 000 time series of BOLD response signals (calculated via standard SPM procedure on various subjects under various stimuli) and 10 000 time series of noisy signals with equal length. There was significant statistical difference between them with Pb.0001. Hence, we can use the quadratic term's coefficient as an assistant numeric indicator in Step 4. A region-of-interest (ROI) analysis was also conducted in this study. The selected regions have shown differential activation in depressed patients and control subjects [16– 17]. Anatomical ROI masks for bilateral inferior frontal cortices and anterior cingulate gyrus were obtained from the Anatomical Automated Labeling library in the Marsbar toolbox (http://marsbar.sourceforge.net/). Data for these areas were represented following the same procedure for whole-brain analysis stated as above. 2.6. Classifier design The naive Bayes classifier technique was applied to recognize each subject, indicating whether he/she is depressed or healthy. The naive Bayes classifier is based on the so-called Bayesian theorem. It is particularly suited when the dimensionality of the inputs is high while the sample number is comparatively small. The pattern recognition toolbox PRtools for MATLAB (http://www.prtools. org/) was used to perform the classifications. In our study, the whole time series of the selected BOLD response ICA components in event-related signals were taken as the classification features. Since the ICA components comparison was employed in the regional level, there were 116 BOLD response ICA components from different regions according to the anatomically labeled template via Ref. [14]. We used PCA again in order to reduce the size of the classification features for Bayesian analysis. Dimension reduction was performed ensuring that n-2 eigenvectors with the largest eigenvalues were reserved. Here, n was the total number of training samples. Given a set of variables X=[x1, x2, …, xd] produced by PCA, we want to construct the posterior probability for one subject among a set of possible classes C. Since naive Bayes assumes that the conditional probabilities of the independent variables are statistically independent, we got the posterior probability of each class membership as: d     Y pðC j jX Þ~p C j p x k C j

ð5Þ

k¼1

Fig. 2. Spectral characteristics for BOLD response signal and random noise. (A) Amplitude frequency diagram for a typical BOLD response reconstructed via SPM. (B) Amplitude frequency diagram for random noise.

where p(Cj|X) was the posterior probability of class membership, i.e., the probability that this subject belonged to Cj. Here, equal prior probabilities were set for both classes. Naive Bayes was modeled via normal functions to estimate p(xk|Cj). Lastly, we labeled a new subject into the class that achieved the highest posterior probability level.

352

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

In order to give more insight into the discriminative characters of the classifier's features, three more methods for model feature design were adopted for comparison. In the first method, BOLD response was created by averaging the volumes within the sad event minus a control volume. In the second method [2], BOLD response was modeled by an events-related wave convolved with the canonical HRF and its temporal derivative. The whole procedure can be undertaken via SPM2. In the third method, activation of the voxel with highest t-value in each region was calculated as the model feature. The t-value was calculated via individual-level statistical analysis. It refers to the statistical significance of the subject's activation under a certain condition. The whole procedure can be employed via software SPM2 with its ‘Model specification and parameter estimation’ function. Both whole-brain analysis and ROI analysis were carried out using these approaches. As depression is associated with a negative bias, the neural activations engaged by processing of sad facial expressions were utilized here. The criterion for sad stimuli selection rested on the argument that increased brain activity has been observed in acutely depressed patients with aversive and negative visual stimuli.

3. Results In this study, since behavioral performance was above 90% correct, we included both the correct and incorrect trials for analysis. After data representation on these trials' time series as depicted above, the event-related ICA components were sorted. Fig. 3 provided sequentially the spectrum diagrams of the top three components in an arbitrarily selected brain region for an arbitrarily selected subject. We can find that component C1 had similar characteristics with typical BOLD response as shown in Fig. 2(A). Therefore, component C1 can be deemed as a component that contained dominating evoked BOLD response information, while

component C2 was a noisy component. As for component C3, it might be composed by some information of BOLD response and some other possible information. Exploring the quadratic coefficient after second-order polynomial fitting, we found that the value relationship was C1NC3NC2. As a result, component C1 was taken as the model feature for further classification analysis. Moreover, an overview of the spectrum diagrams for all ICA components of these subjects showed that BOLD response information always appeared via the top two components. This result was consistent with the analysis in Section 2.5. Leave-one-out cross-validation test was applied to evaluate classifier's performance. In each trial, data from all but one subject were utilized to train the classifier. Subsequently, the class assignment of the remaining subject was calculated during the test phase. Generalization rate, sensitivity and specificity were utilized to quantify the performance. Sensitivity indicates the proportion of patients correctly predicted, and specificity indicates the proportion of normal controls correctly predicted, while generalization rate is the overall proportion of samples correctly predicted. Using the event-related BOLD response extracted via the above procedures, classification of subjects was carried out. From the whole-brain analysis, the accuracy was 77.27%. Among them, 83.33% (sensitivity) of the patients were correctly classified into the patient group, and 70.00% (specificity) of the healthy subjects were correctly assigned to the control group. The accuracy for ROI analysis was 81.82%, with the sensitivity being 84.62% and the specificity being 77.78%. In order to study whether the observed classification accuracy was statistically significant, permutation test was carried out to compute the probability of a classification result at least as accurate as that observed level if the subjects were allocated at random to either group. The calculated probabilities were P=.017 and P=.009 for whole-brain analysis and ROI analysis, respectively. Hence, the performance of the classification method was statistically significant. The top 10 discriminative regions were provided in Fig. 4 and Table 2 for each group. We found that the regions with the highest contribution for group discrimination were those emotionally related regions, such as the cingulate gyrus, the inferior frontal gyrus, the thalamus, the medial frontal gyrus and the hippocampal gyrus. In order to give more insight into the classifier performance, three more methods for model feature design (see Section 2.6 for method details) were adopted for comparison. However, none of them yielded classification performance as good as our proposed method did. All the results were summarized in Table 3.

4. Discussions

Fig. 3. Spectrum diagrams of the top three components in an arbitrarily selected brain region for an arbitrarily selected subject.

Depressive disorder is a common mental disorder, which is characterized by persistent depressed mood, anxiety and dysphoria, psychomotor changes, alterations of motivation

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

353

Fig. 4. Weight map (absolute values) overlaid on structural images illustrating the highly discriminate regions, with Z value indicating slice position. (A) Top 10 discriminate regions for patient group. (B) Top 10 discriminative regions for healthy control group.

and social behavior, and sleep abnormalities. Current clinical practice relies on psychiatrists' verbal communication with patients or accompaniers, trying to seek disease-related symptoms, and finally reaching a diagnosis and treatment decision. Such approaches depend on psychiatrists' experience and patients' cooperation greatly. Therefore, neuroimage-based objective diagnosis or measurement is required to support clinical practice. In such exploration, the features extracted from scanning images determine the diagnosis performance greatly. This paper has compared the proposed classification feature design method with three other typical ways of constructing the event-related features. The validity of these extracted features was tested via depression detection based on our samples. Among them, one general way of estimating BOLD response signal was using the convolution model in SPM2 to estimate the temporal dynamics (method 2 in Table 3); thereafter, the related features were extracted for disease detection. The classification accuracy can only reach 63.64% based on our samples. It can be due to the fact that the intra- and intersubject variability was averaged since stimuli responses were modeled by standard hemodynamic response functions.

In general, the BOLD response is associated with a given impulse of neural activity enduring for many seconds. In such a prolonged temporal response, variability is aggravated due to the complex mechanism inherent in the brain, including intersubject variability and intrasubject variability. In order to deal with the limits arising from model dependency, several efforts were also reported and taken as comparison in this study as well. Averaging volumes within the event minus a control volume was a general and simple way (method 1 in Table 3) to estimate BOLD response. The fair accuracy of 59.09% in this study demonstrated such limits. As a substitute way, activation of the most active voxel in each region was calculated as model features in method 3. It helped reduce the data size depending on prior knowledge. However, similar fair disease detection performance was achieved again using this method. One major reason for such phenomenon was that subtle regional differences may not be revealed from representative voxels. Different from all above investigations, we combined event-related and resting-state fMRI signal together for BOLD response detection (method 4 in Table 3). Here, the resting-state fMRI signal was deemed as the baseline of each

354

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

Table 2 Top 10 highly discriminative regions Regions Depression Right superior temporal gyrus Left medial frontal gyrus Right inferior temporal gyrus Right medial frontal gyrus Left superior frontal gyrus Left cingulate gyrus Right parahippocampal gyrus Right thalamus Right inferior frontal gyrus Left inferior frontal gyrus Controls Right superior temporal gyrus Left medial frontal gyrus Right inferior temporal gyrus Right cingulate gyrus Right medial frontal gyrus Left superior frontal gyrus Right parahippocampal gyrus Left superior temporal gyrus Right middle temporal gyrus Right thalamus

Discriminative weights 0.4024 0.2799 0.2601 0.2242 0.1750 0.1205 0.1091 0.0799 0.0777 0.0749 0.3646 0.3256 0.2785 0.2454 0.1728 0.1080 0.0982 0.0976 0.0911 0.0840

subject's activities. Gusnard and Raichle [18] proposed that there exists a physiological baseline of the human brain, which can be observed when subjects are awake and resting with their eyes closed. Based on positron emission tomography studies, Raichle and colleagues [18–20] have defined a default network of brain activity consisting of areas that are metabolically active in the resting state. The uniformity of the oxygen extraction fraction (OEF) at rest indicated that equilibrium has been reached between the local metabolic requirements, which were necessary to sustain a long-term continuing level of neural activity, and the level of blood flow in that region. Consequently, those areas with a reduced/ increased OEF relative to the mean OEF of the brain were defined as ‘activated’/‘deactivated.’ Based on these studies, the scanned signal from event-related experiments can be deemed as the final time series of BOLD response inspired by stimuli based on an equilibrium level at resting state. Comparison of the ICA components from resting-state signals and eventrelated signals helped find the time series which represented the major evoked BOLD information without eliminating the variability between subjects. A temporal correlation will require a phase coherence of two periodic signals. It could be zero or any other values depending on the relative phase. As a result, an accidental consequence of phase offsets between two periodic signals would lead to a small correlation value as well, no matter how interrelated actually they were. Our method was demonstrated to be free from such phenomenon when it tried to calculate the coefficient mark. If the correlation of one component in the resting state to the event-related component is small due to phase offset, it was excluded from further consideration as coefficient marks for this eventrelated component since we tried to find the maximum

coefficient as the mark according to Eq. (4). Thereafter, the sort mechanism did not consider the small correlation value arising from phase offset, although it aimed to find the component with the smallest coefficient mark. Another important point is that event-related data and resting-state data were not acquired simultaneously. It may raise a question whether we can utilize the resting-state data as the baseline for the event-related signals. It was supposed that correlation analysis based on the scanning length of 6 min 39 s was long enough to eliminate the effect from the variation between sessions. Another support for our work can be found in the related work of De Luca et al. [21], who stated that distinct resting state network patterns are reproducible. Accordingly, we simplified our study based on the assumption that critical signal dynamics under resting state would not vary too much between trials. Classification of groups with ROI analysis method, which relied on experienced results observing significant differences between patients and healthy control groups, was deemed insufficient to reveal subtle regional differences in previous discussion based on some typical methods [2]. However, it is worth noting that the accuracy for ROI analysis using the proposed method was 81.82%, with the sensitivity being 84.62% and the specificity being 77.78%. The related reason may be behind the possibility that the sorting methods inside the proposed approach helped reduce the irrelevant information and sharpen the subtle group differences. The weight map shown via Fig. 4 demonstrated the importance of brain regions for illness estimation in the proposed discriminative mechanism. The regions with the highest weight included the right temporal lobe, the right inferior frontal gyrus, the right cingulate gyrus, the hippocampus and so on. Such results were consistent with the previous study on depression. For example, Davidson Table 3 Group classification of depressed and healthy individuals

Method 1 ROI Whole brain regions Method 2 ROI Whole brain regions Method 3 ROI Whole brain regions Method 4 ROI Whole brain regions

Accuracy

Sensitivity

Specificity

59.09% 45.45%

64.29% 55.56%

50.00% 38.46%

59.09% 63.64%

64.29% 69.23%

50.00% 55.56%

45.45% 59.09%

53.85% 70.00%

33.33% 50.00%

81.82% 77.27%

84.62% 83.33%

77.78% 70.00%

Method 1: BOLD response was created by averaging the volumes within the sad event minus a control volume. Method 2: BOLD response was convoluted by a standard hemodynamic response function. Method 3: Activity of the most active voxel in each region was calculated as model features. Method 4: BOLD response was extracted via the proposed method in this paper.

Q. Lu et al. / Magnetic Resonance Imaging 30 (2012) 347–355

et al. (2003) [16] underscored the importance of the neural circuitry activated by negative affect in depression. The cingulate cortex was the core region inside such a neural circuitry, which included the lateral prefrontal cortex, the amygdala, the insular cortex hypothalamus and the hippocampus. It was supported by Drevets [22] and Mayberg et al. [23] as well, who stated that the affective component of depressive psychopathology may be primarily related to functional alterations in the brain regions such as the anterior cingulate cortex (ACC). It was also demonstrated that brain areas involved in executive function include the dorsal lateral frontal cortex (DLPFC), the orbitofrontal cortex (OFC), the ACC and the basal ganglia. Those areas form a frontalstriatum loop, with the prefrontal cortex as the center of the loop. In the process of attention and inhibition, increased activation showed in the ACC; in the process of response forbidden and task management, the main regions of increased activation showed in the DLPFC, especially in the right DLPFC. Considering the stimuli in the study, experience of sad emotion was related with activation areas involved in the temporal lobe and right DLPFC. Emotion inhibition experience was related with activation areas involved in the right DLPFC and right OFC. Sad emotion was related with cingulate cortex activation. The above analysis demonstrated that signal components involved in emotion and recognition can be extracted effectively by the proposed method in this paper. Despite other well-developed complex classifiers, Bayes classification method was applied in this paper. Although the assumption that the predictor variables are independent is not always accurate, it does simplify the classification task dramatically. In effect, naive Bayes reduces a highdimensional density estimation task to a one-dimensional kernel density estimation. Furthermore, the assumption does not seem to greatly affect the posterior probabilities, especially in regions near decision boundaries, thus leaving the classification task unaffected. Acknowledgment The work was supported by grants of The National Natural Science Foundation of China (30900356, 81071135) and Nanjing Medical Technology Development Project (ZKX09037). References [1] Ogawa S, Lee TM, Kay AR, Tank DW. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc Natl Acad Sci USA 1990;87:9868–72. [2] Fu CHY, Mourao-Miranda J, Costafreda SG, Khanna A, Marguand AF, Williams SCR, et al. Pattern classification of sad facial processing: toward the development of neurobiological markers in depression. Biol Psychiatry 2008;63:656–62.

355

[3] Demirci O, Calhoun VD. Detection of schizophrenia using FMRI data via projection pursuit. IEEE Workshop Mach Learn Signal Proc 2007: 199–204. [4] Kawasaki Y, Suzuki M, Kherif F, Takahashi T, Zhou SY, Nakamura K, et al. Multivariate voxel-based morphometry successfully differentiates schizophrenia patients from health controls. NeuroImage 2007;34:235–42. [5] Zhu CZ, Zang YF, Cao QJ, Yan CG, He Y, Jiang TZ, et al. Fisher discriminative analysis of resting-state brain function for attentiondeficit hyperactivity disorder. NeuroImage 2008;40:110–20. [6] Calhoun VD, Adali T, Pearlson GD. Independent components analysis applied to fMRI data: a generative model for validating results. J VLSI Signal Process Syst Signal Image VideoTechnol 2004;37:281–91. [7] McKeown MJ, Hansen LK, Sejnowski TJ. Independent component analysis of functional MRI: what is signal and what is noise? Curr Opin Neurobiol 2003;13:620–9. [8] Dimitriadou E, Barth M, Windischberger C, Hornik K, Moser E. A quantitative comparison of functional MRI cluster analysis. Artif Intell Med 2004;31:57–71. [9] Windischberger C, Barth M, Lamm C, Schroeder L, Bauer H, Gur RC, et al. Fuzzy cluster analysis of high-field functional MRI data. Artif Intell Med 2003;29:203–23. [10] Sato JR, Fujita A, Amaro E, Miranda JM, Morettin PA, Brammer MJ. DWT-CEM: an algorithm for scale-temporal clustering in fMRI. Biol Cybern 2007;97:33–45. [11] Lessa PS, Sato JR, Cardoso EF, Neto CG, Valadares AP, Amaro Jr E. Wavelet correlation between subjects: a time-scale data driven analysis for brain mapping using fMRI. J Neurosci Methods 2011;194:350–7. [12] Fox MD, Snyder AZ, Vincent JL, Raichle ME. Intrinsic fluctuations within cortical systems account for intertrial variability in human behavior. Neuron 2007;56:171–84. [13] Cow RW. AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput Biomed Res 1996;29: 162–73. [14] Tzourio-Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, Delcroix N, et al. Automated anatomical labelling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single subject brain. Neuroimage 2001;15:273–89. [15] Hyvärinen A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans Neural Netw 1999;10:626–34. [16] Davidson RJ, Irwin W, Anderle MJ, Kalin NH. The neural substrates of affective processing in depressed patients treated with venlafaxine. Am J Psychiatry 2003;160:64–75. [17] Yao ZJ, Wang L, Lu Q, Liu HY, Teng GJ. Regional homogeneity in depression and its relationship with separate depressive symptom clusters: a resting-state fMRI study. J Affect Disord 2009;115:430–8. [18] Gusnard DA, Raichle ME. Searching for a baseline: functional imaging and the resting human brain. Nat Rev Neurosci 2001;2:685–94. [19] Raichle ME, MacLeod AM, Snyder AZ, Powers WJ, Gusnard DA, Shulman GL. A default mode of brain function. Proc Natl Acad Sci USA 2001;98:676–82. [20] Raichle ME, Mintun MA. Brain work and brain imaging. Annu Rev Neurosci 2006;29:449–76. [21] De Luca M, Beckmann CF, De Stefano N, Matthews PM, Smith SM. fMRI resting state networks define distinct modes of long-distance interactions in the human brain. Neuroimage 2006;29:1359–67. [22] Drevets WC. Neuroimaging and neuropathological studies of depression: implications for the cognitive-emotional features of mood disorders. Curr Opin Neurobiol 2001;11:240–9. [23] Mayberg HS, Brannan SK, Tekell JL, Silva JA, Mahurin RK, McGinnis S. Regional metabolic effects of fluoxetine in major depression: serial changes and relationship to clinical response. Biol Psychiatry 2000;48:830–43.