Frequency Clustering Analysis for Resting State Functional Magnetic Resonance Imaging Based on Hilbert-Huang Transform

Objective: Exploring resting-state functional networks using functional magnetic resonance imaging (fMRI) is a hot topic in the field of brain functions. Previous studies suggested that the frequency dependence between blood oxygen level dependent (BOLD) signals may convey meaningful information regarding interactions between brain regions. Methods: In this article, we introduced a novel frequency clustering analysis method based on Hilbert-Huang Transform (HHT) and a label-replacement procedure. First, the time series from multiple predefined regions of interest (ROIs) were extracted. Second, each time series was decomposed into several intrinsic mode functions (IMFs) by using HHT. Third, the improved k-means clustering method using a label-replacement method was applied to the data of each subject to classify the ROIs into different classes. Results: Two independent resting-state fMRI dataset of healthy subjects were analyzed to test the efficacy of method. The results show almost identical clusters when applied to different runs of a dataset or to different datasets, indicating a stable performance of our framework. Conclusions and Significance: Our framework provided a novel measure for functional segregation of the brain according to time-frequency characteristics of resting state BOLD activities.

The time-based methods are usually convenient and effective in examining the point-to-point relationship between regional blood oxygen level dependent (BOLD) signals. Recent studies suggested that the frequency dependence between BOLD signals may also convey meaningful information regarding interactions between brain regions (Wen et al., 2012;Yu et al., 2013;Wei et al., 2014). In a recent fMRI study (Song et al., 2014), a ReHo based frequency clustering analysis framework was introduced for resting-state fMRI analysis. The BOLD time series of each voxel was decomposed into several frequency components using empirical mode decomposition (EMD), and the ReHo values of the components were used as features for clustering the voxels based on similar frequency-specific ReHo signature. The forging studies indicated that analyzing time-frequency characteristics is equally important for comprehensively exploring how different brain systems/sub-systems coordinate.
One challenge of further extracting the time-frequency characteristic in frequency clustering analysis is time-frequency representation of fMRI signals. In many previous studies, the time-frequency characteristics of fMRI time series were usually measured using short-time Fourier transform (Mezer et al., 2009) or wavelet transform (Bullmore et al., 2001;Shimizu et al., 2004) which always assume the linearity or stationarity of input signals (Huang and Shen, 2005). However, BOLD time series may not conform to these assumptions (Lange and Zeger, 1997). Furthermore, constrained by the Uncertainty Principle (Robertson, 1929), most of the traditional timefrequency methods are limited in providing both high temporal resolution and high frequency resolution at the same time.
HHT is a novel time-frequency method suitable for both non-linear and non-stationary signals. Its application to electrophysiological studies has demonstrated its efficay in providing fine expressions of instantaneous frequency (Huang and Shen, 2005;Peng et al., 2005;Donnelly, 2006;. For example, HHT has been successfully applied in EEG-based seizure classification (Oweis and Abdulhay, 2011), detection of spindles in sleep EEGs (Yang et al., 2006), and ECG de-noising (Tang et al., 2007). However, its application in fMRI studies is rare.
Other challenges in time-frequency analysis based frequency clustering analysis voxel-wised analysis at different frequency bands may demand great amount of calculation. Furthermore, in the stage of clustering analysis, the labels of the clusters change randomly across analyses (Mezer et al., 2009), causing difficulty in cross-condition/datasets comparisons.
In this article, we introduced a novel frequency clustering analysis method based on HHT and an improved k-mean clustering method using label-replacement procedure. In our framework, first, the time series from multiple predefined regions of interest (ROIs) [i.e., 90 ROIs defined by the Automated Anatomical Labeling (AAL) template Tzourio-Mazoyer et al., 2002] were extracted. Second, each time series was decomposed into several intrinsic mode functions (IMFs) of which the instantaneous frequency characteristics were subsequently calculated using HHT. Third, the improved k-means clustering method using a label-replacement method was applied to the data of each subject to classify the ROIs into different classes. To test the efficacy of our frequency clustering analysis method, two independent resting-state fMRI data sets of healthy subjects (198 subjects in Dataset I; 88 subjects in Dataset II) were analyzed. The results demonstrated that for different dataset, our method generated stable clusters of the brain regions according to time-frequency characteristics of their resting state BOLD activities.

MATERIALS AND METHODS fMRI Data Acquisition
In this study, we used a resting-state fMRI dataset (Dataset I) provided by the open source website of "1,000 Functional Connectomes' Project" (http://www.nitrc.org/projects/fcon_ 1000/). The dataset included functional and structural MR images recorded from 198 healthy subjects (18-26 years old, 122 females) by Yu-Feng Zang (Song et al., 2014). No subject had a history of neurological, psychiatric or medical conditions. Written informed consents were given to all subjects in accordance with Institutional Review Board guidelines and in compliance with the Declaration of Helsinki.
The scanning was performed using a 3.0-Tesla scanner (Siemens TRIO TIM, Munich, Germany). The subjects were instructed to rest with their eyes closed, keep their heads still, and not to fall asleep. A gradient echo T2 * -weighted EPI sequence was used for acquiring resting state functional images with the following parameters: TR = 2,000 ms; TE = 30 ms; 33 slices; matrix size = 64 × 64; FOV = 240 × 240 mm 2 acquisition voxel size = 3.75 × 3.75 × 3.50 mm 3 ; 225 volumes.
In order to verify the stability of the results, an independent resting-state data set (Dataset II) was also analyzed. Dataset II was collected from 88 healthy young right-handed college students (19-26 years old, 44 females) performing eyes-closed resting state task. The scanning was performed using a 3.0-Tesla Siemens whole-body MRI system in Brain Imaging Center of Beijing Normal University. All subjects were given the written informed consents before scanning. No subject had a history of neurological, psychiatric, or medical conditions. The scan was performed during a resting-state condition. The detailed parameters used were as follows: TR = 2,000 ms; TE = 30 ms; 33 slices; matrix size = 64 × 64; FOV = 240 × 240 mm 2 acquisition voxel size = 3.13 × 3.13 × 3.60 mm 3 ; 145 volumes. The experiment was approved by the Institutional Review Board of the Beijing Normal University.

Image Preprocessing
In current research, the images were analyzed using SPM8 (http://www.fil.ion.ucl.ac.uk/spm). For each subject, the original first five functional volumes were removed to avoid the possible disturbance caused by non-equilibrium effects of magnetization. The remaining functional images (220 in Dataset I; 140 in Dataset II) were corrected for slice timing, motion corrected, and spatially normalized into a Montreal Neurological Institute (MNI) space using the standard EPI template (Evans et al., 1993). The normalized image had 61 slices, a matrix size of 61 × 73, and a voxel size of 3 mm × 3 mm × 3 mm. No translation or rotation movement in any data set exceeded ±2 mm or ±2 degree. The data had originally been "cleaned" through the use of confound regressors derived from CSF and white matter masks, as well as head motion parameters. The linear trend was regressed out for each voxel's time course to remove signal drifts caused by scanner instability or other factors. The time course of each voxel was normalized by subtracting the temporal mean and dividing by the temporal standard deviation. After preprocessing, for each subject and for each of the 90 regions of interest (ROIs) defined using the AAL template (Tzourio-Mazoyer et al., 2002), an ROI time course was extracted by averaging the time courses of all voxels in the ROI.
Using Hilbert-Huang Transform (HHT) to Acquire Instantaneous Frequency and Hilbert Weighted Frequency (HWF) The forgoing ROI time courses were fed into HHT to acquire instantaneous frequency and HWF feature. The HHT consists of three main processes. First, major IMFs are extracted from the input signal based on empirical mode decomposition (EMD). Second, Hilbert transform is applied to each IMF to obtain the analytic transform of the original signal. Last, the instantaneous frequency is calculated according to the analytical transform of each IMF (Huang and Shen, 2005;Huang and Pan, 2006;Ding et al., 2007) and the Hilbert weighted frequency (Xie and Wang, 2006) of each IMF is calculated according to the instantaneous frequency of the IMF. The detailed descriptions of each step are as follows: (1) Empirical mode decomposition The EMD method (Huang and Shen, 2005) decomposes an input signal into a finite set of intrinsic oscillatory components, namely, the IMFs. Mathematically, for fMRI time series, EMD generates a set of IMFs and a monotonic residue signal r(t): where N is the number of the IMFs. Each IMF must satisfy two conditions: (i) Along the time course of the IMF, the number of the local extrema and the number of zero crossings are either equal or differ by one; (ii) The sum of the envelope defined by the local maxima and the envelope defined by the local minima is constantly zero.
To extract IMFs using EMD, an iterative method known as the sifting algorithm is used as follows: Step 1: Find the local extrema of the input signal; Step 2: Use interpolation to generate the lower envelope elow(t) and the upper envelope eup(t) of the current signal according to the local minima and local maxima respectively; Step 3: Calculate the local mean time course emean(t): Step 4: Obtain the "oscillatory-mode" r(t) = x(t) − emean(t); Step 5: If r(t) satisfies the stopping criteria (the two conditions of IMF), IMF i = r(t) becomes an IMF, otherwise set x(t) = r(t) and repeat the process from Step 1.
To obtain remaining IMFs, the same procedure is applied iteratively to the residual r(t) = x(t) − IMF i (t) until r(t) is monotonic.
Hilbert transform was used to extract the instantaneous frequency of each IMF. For signal x(t), its Hilbert transform y(t) is defined as: where P is the Cauchy principal value (Surhone et al., 2013). Hilbert transform is capable of describe the local properties of x(t) (Peng et al., 2005). The analytic transform of z(t) x(t) is defined as: where a(t) is the instantaneous amplitude, and ϕ (t) is the instantaneous phase. The instantaneous frequency ω(t) is defined as the time derivative of ϕ (t): (3) Hilbert weighted frequency (HWF) based on instantaneous frequency The Hilbert weighted frequency (Xie and Wang, 2006) of each IMF is also calculated based on the instantaneous amplitude and phase to reflect the mean oscillation frequency of the IMF. The HWF(j) of the jth IMF is defined as: where ω j (i) is the instantaneous frequency, a j (t) is the instantaneous amplitude, and m is the number of time point.

Identify the Brain Networks Using k-means Clustering Analysis Based on HWF Characteristics
In order to identify the brain regions sharing common instantaneous frequency characteristics, we employed k-means clustering analysis to the two resting-state datasets respectively.
In each analysis, the 90 ROIs were classified into different clusters based on the HWF feature vectors each of which comprised the first five HWFs. For each feature vector, each HWF form all subjects was concatenated to yield a group feature vector for the following analysis. The clustering analysis was performed for different k from 1 to 90. The squared Euclidean distance index (Mezer et al., 2009) for different k values was plotted as a function to determine the appropriate k.
Label-Replacement Method to Improve k-means Clustering Analysis In conventional k-means clustering analysis, the labels of each cluster are unordered due to random initialization of the algorithm which hinders the compare between conditions and datasets. The first impede caused by the randomization is that the label of a cluster, of which the spatial structure changes little, may change dramatically from run to run even for the same condition in the same dataset. For example, some brain regions are classified into a cluster labeled as "1" in one run, and into the same cluster but labeled as "2" in another run. In this study, we proposed a method for sorting the label of clusters. The method composed of two steps: (1) label-sorting and (2) label-matching. In label-sorting, a hash table (Maurer and Lewis, 1975) was used to improve the computing efficiency. The detailed process was as follows: Step 1: Obtain an unordered label table (ULT(i)) according to the raw results of k-means clustering.
Step 2: A hash table (Hash(i)) is constructed to record the labels in turn.
Step 3: If the label is not found in Hash(i); add the label to Hash(i) and record the order by Lab If the label already exists in Hash(i); replace the label of ULT(i) by Lab; update Lab and return to Step2.
A sorted label table (SLT(n)) is constructed by the follow process (shown in Figure 1): In this part, the hash table is used to record a new label of each ULT(n). In each repetition, the algorithm searches the hash table to decide whether the unordered label needs to be replaced. After label-sorting, the ULT(i) was sorted in the order of brain regions (AAL).
The second impede caused by the randomization is that label may change across different conditions/datasets. That is, for the same cluster, its label within one condition/dataset could be different from that in another condition/dataset, making it hard to compare between different groups or between different conditions. This could be even worse when the cluster changes slightly across conditions/dataset. Therefore, we developed a label-matching method besides the forgoing label sorting method. The basic idea is that use the label setting of one of the condition/dataset as a reference, then go through the target clusters of another condition/dataset, when a target cluster share a similar spatial pattern with a reference cluster, the label of this reference cluster is assigned to the target cluster.
The detail of matching the target and the reference cluster is as below (shown in Figure 2).

RESULTS
For each dataset, BOLD time series from the 90 ROIs defined by AAL template were extracted. First, EMD was applied to decompose the BOLD time series into different frequency components. The EMD outcome of both datasets showed that the BOLD signals could be decomposed into five major IMFs (please see examples in Figure 3).
In order to determine the number of IMFs which were selfadaptive decomposed from fMRI signal, the EMD of each vowel in each subjects were calculated.
According to the results (Figure 3), most of the voxels have at least 5 IMFs. Therefore, in our study, the first 5 IMFs were selected for the following analysis (an example of EMD was shown in Figure 4, left panel).
Second, the instantaneous frequency of each major IMF was calculated using Hilbert transform and its corresponding power spectrum was calculated. The results demonstrate that the instantaneous frequency of the IMFs varies across frequency bands centered by different dominant frequency from low (around 0.01 Hz) to high (around 0.12 Hz) (please see the example in Figure 4).
Third, the HWF of each major IMF was calculated. The histograms of HWF distributions in the whole brain across all subjects showed that the major IMFs occupy certain frequency bands: IMF1, 0-0.01 Hz; IMF2, 0.005-0.015 Hz; IMF3 0.01-0.03 Hz; IMF4, 0.03-0.07 Hz; and IMF5, 0.08-0.18 Hz (Figure 5).  The first five HWF of each ROI were selected as the features for k-means clustering analysis. For selecting the appropriate parameter (k), k-means clustering analysis was repeated for different k from 1 to 90 (shown in Figure 6). Then the clusters were evaluated using a squared Euclidean distance index (Mezer et al., 2009) and the appropriate parameter, k = 20, was selected for each subject.
A comparison of the stability of before and after label-sorting method was shown in Figure 7. After label-sorting, the results of cluster were sorted in the order of brain regions.
The results of the two datasets show almost identical clusters of the ROIs and a comparison of the results of before and after label-replacement (shown in Figure 8).

DISCUSSION
In this study, we presented a new frame work for brain region clustering based on HHT method and an improved k-mean clustering method adopting cluster label sorting, and applied the frame work to analyzing two independent resting-state fMRI data sets of healthy subjects (198 subjects in Dataset I; 88 subjects in Dataset II). The results showed that our method is efficacy in functional segregation based on time-frequency characteristics of resting state BOLD activities.

Applying Hilbert-Huang Transform to BOLD Signal Analysis
In our study, HHT was introduced to frequency clustering analysis of resting-state fMRI. The advantage of using HHT is mainly based on two features of the fMRI data. First, previous studies have shown that the fMRI data may not strictly conform to the assumptions of linearity and stationarity (Lange and Zeger, 1997). Compared with traditional time-frequency analysis method such as wavelet transform, short-time Fourier transform, HHT can be directly applied to the non-linear or non-stationary signals. Second, the fMRI signals mainly record the slow BOLD change in low sampling intervals around 1-3 s. Traditional timefrequency methods constrained by the Uncertainty Principle (Robertson, 1929) are limited in providing both high temporal resolution and high frequency resolution at the same time. While many previous studies have shown that HHT does not suffer from the trade-off between frequency resolution and temporal resolution (Peng et al., 2005;Donnelly, 2006; and thus may be a suitable candidate for time-frequency representation of the fMRI signals. Our results show that HHT can represent the fMRI signals in both high temporal resolution and high frequency resolution.

Number of IMFs in Empirical Mode Decomposition
EMD method, as an important part of HHT, is a local, fully data driven and self-adaptive analysis approach. The results of EMD show that different brain voxel/regions may contain multiple IMFs, the number of IMF voxels may affect the efficacy of clustering in the features space. On one hand, too small number of IMFs may quicken the processing but provides inadequate information to describe the functional significance of the clusters. On the other hand, too larger number of IMFs may impede the clustering progress. In our method, the number of IMFs was determined by picking up the number that present in most of the brain voxels (Figure 3). The HWF of each IMF was determined using instantaneous amplitude and phase information to reflect the mean oscillation frequency of the IMF. Previous research has shown that HWF has clear physical meaning and has low variability in terms of robustness (Xie and Wang, 2006). Therefore, the HWF of the first five IMFs were calculated to construct the feature vector. And then, the histogram of HWF was also calculated (shown in Figure 5) to show the mean frequency distribution of the five IMFs. Each of the five histograms is a statistic of the whole-brain voxels. Since the frequency content of different voxels at different sites Frontiers in Human Neuroscience | www.frontiersin.org  of the brain (and subjects) are roughly similar, the same IMF (IMF(j), j = 1, 2, 3, 4, 5) from any voxel will roughly fall into the same frequency band. These results suggest that EMD works well in adaptively decomposing the fMRI signals into different IMFs that fall into distinctive frequency bands and is a promising method for non-stationary and non-linear neurological signal processing.

Label-Sorting for k-means Clustering
In previous frequency clustering analysis, k-means clustering method has been applied to resting-state fMRI network analysis (Song et al., 2014) and generates meaningful results. However, in the previous work, the labels of the clusters were randomly assigned and changed from analysis to analysis, making it hard to compare between conditions/datasets. Our study presents a label-sorting method which uses Hash table to obtain an ordered and stable clustering result across different runs of analysis within a condition/dataset, and further a cluster-label matching method to deal with cluster matching and label assignment across conditions/datasets. The verification results showed almost identical clusters no matter when the method was applied to different runs of a dataset or to different datasets, indicating a stable performance of our framework (Figures 7, 8). It is worth noting that when condition or dataset changes, the spatial representation of an underlying brain cluster may also change in some extends according to the real scenario. A careful visual inspection for potential unmatched cluster and FIGURE 7 | A comparison of the stability of before and after label-sorting method in each run.
FIGURE 8 | The k-means frequency clustering analysis of two independent fMRI data sets and a comparison of clusters of before and after label-replacement. label caused by dramatic brain change between conditions or subject groups is recommended besides our frame work. The change itself, if significant, could deliver meaningful clinical, neurological and psychophysiological information.

Selecting the Regions of Interest
It is worth noting that using AAL template or for ROI selection is not part of the major line but an alternative module of our frame work. In the current study we used the classical AAL template for ROI selection to demonstrate the performance of our method. However, AAL template is defined anatomically. The current work only clustered the known anatomical structures function similarly on frequency domain. Indeed, the AAL ROI template can be replaced by other ROIs or voxels according to the interest of the researchers. The ROIs can be a set of task-activated sites with their intrinsic relationship to be clarified, or a set of anatomically defined structures. The analysis can also be performed in a whole-brain or partial-brain voxel-wise fashion.

Limitation of the Current Work
In the current work, the data driven process introduce five IMFs referring to different frequency bands. The IMF 5 was corresponding to a frequency band of 0.08∼0.18 Hz which had been excluded in most of the previous resting state functional connectivity analysis. It is a nontrivial question that what information the higher frequency bands of BOLD change provide. Although, faster neural electrophysiological activities have been found in higher order regions such as the frontal lobe and were proposed to carry important cognitive meanings (Lang et al., 1986), the neurocognitive meaning of the higher frequency components of BOLD which were usually considered as noise in many previous functional connectivity studies is under debate (Michels et al., 2010;Boubela et al., 2013). Understanding the neurocognitive meanings of the clustering results requires further careful works in the future.

CONCLUSION
In this study, a novel frequency clustering analysis method based on HHT and a label-replacement procedure was introduced. First, the ROI time series were extracted. Second, each time series was decomposed into several intrinsic mode functions (IMFs) by using HHT. Third, the improved k-means clustering method using a label-replacement method was applied to the data of each subject to classify the ROIs into different classes. Two independent resting-state fMRI dataset of healthy subjects were analyzed to test the efficacy of method. The results showed that for different dataset, our method can stably cluster the brain regions according to the time-frequency characteristics of their resting state BOLD activities.