Spatiotemporal Empirical Mode Decomposition of Resting-State fMRI Signals: Application to Global Signal Regression

Moradi, Narges; Dousty, Mehdy; Sotero, Roberto C.

doi:10.3389/fnins.2019.00736

ORIGINAL RESEARCH article

Front. Neurosci., 23 July 2019

Sec. Brain Imaging Methods

Volume 13 - 2019 | https://doi.org/10.3389/fnins.2019.00736

This article is part of the Research Topic Origins of the Resting-State fMRI Signal View all 14 articles

Spatiotemporal Empirical Mode Decomposition of Resting-State fMRI Signals: Application to Global Signal Regression

$\nNarges Moradi,,$ Narges Moradi^1,2,3^*

Mehdy Dousty^4,5

Roberto C. Sotero^1,2,3

¹Biomedical Engineering Graduate Program, University of Calgary, Calgary, AB, Canada
²Hotchkiss Brain Institute, University of Calgary, Calgary, AB, Canada
³Computational Neurophysics Lab, Department of Radiology, University of Calgary, Calgary, AB, Canada
⁴Institute of Biomaterials and Biomedical Engineering, University of Toronto, Toronto, ON, Canada
⁵KITE, Toronto Rehab, University Health Network, Toronto, ON, Canada

Resting-state functional connectivity MRI (rs-fcMRI) is a common method for mapping functional brain networks. However, estimation of these networks is affected by the presence of a common global systemic noise, or global signal (GS). Previous studies have shown that the common preprocessing steps of removing the GS may create spurious correlations between brain regions. In this paper, we decompose fMRI signals into 5 spatial and 3 temporal intrinsic mode functions (SIMF and TIMF, respectively) by means of the empirical mode decomposition (EMD), which is an adaptive data-driven method widely used to analyze non-linear and non-stationary phenomena. For each SIMF, functional connectivity matrices were computed by means of Pearson correlation between TIMFs of different brain areas. Thus, instead of a single connectivity matrix, we obtained 5 × 3 = 15 functional connectivity matrices. Given the high correlation and global efficiency values of the connectivity matrices related to the low spatial maps (SIMF3, SIMF4, and SIMF5), our results suggest that these maps can be considered as spatial global signal masks. Thus, by summing up the first two SIMFs extracted from the fMRI signals, we have automatically excluded the GS which is now voxel-specific. We compared the performance of our method with the conventional GS regression and to the results when the GS was not removed. While the correlation pattern identified by the other methods suffers from a low level of precision in identifying the correct brain network connectivity, our approach demonstrated expected connectivity patterns for the default mode network and task-positive network.

1. Introduction

Resting-state functional connectivity MRI (rs-fcMRI) has considerable potential for mapping functional brain networks (Biswal et al., 1995; Kandel et al., 2000; De Luca et al., 2006; Fox et al., 2006; Shmuel and Leopold, 2008; Friston, 2011). This mapping, which reveals the brain's functional architecture and operational principles (Kandel et al., 2000; Friston, 2011), can be used for early detection of brain connectivity pathologies in neuropsychiatric patients (Erdoğan et al., 2016). However, the presence of broadly shared synchronous fluctuations, termed as the global signal (GS) in Blood Oxygen Level Dependent (BOLD) responses, is a significant problem for fcMRI analysis. Its presence is problematic as it is of unknown origin (Damoiseaux et al., 2006; Fox et al., 2009; Erdoğan et al., 2016). Therefore, effective removal of GS has become an important step in data preprocessing and must be done prior to fcMRI analysis. GS is generally defined as the average of the BOLD signals over the whole brain (Zarahn et al., 1997; Fox et al., 2009; Liu et al., 2017) and can be computed from the raw images or after some preprocessing steps (Liu et al., 2017). The average-based GS is typically called conventional GS (or static GS (SGS) Erdoğan et al., 2016).

Application of SGS regression (SGSR) was at first just limited to task-related fMRI imaging (Zarahn et al., 1997; Macey et al., 2004). More recently, SGSR usage has received more attention in the analysis of resting-state fMRI than in task-related fMRI studies (Liu et al., 2017). Some studies suggest that application of SGSR improves the functional specificity of resting-state correlation maps and helps to remove non-neuronal sources of global variance like respiration and movement (Fox and Raichle, 2007; Fox et al., 2009; Liu et al., 2017). However, other studies found that these improvements are limited to systems that would exhibit only positive correlations with the specific selected seeds (Fox et al., 2009; Weissenbacher et al., 2009). On the other hand, many studies have shown that the common preprocessing steps of removing GS via a general linear model can create correlations between regions that may never have existed (Murphy et al., 2009; Anderson et al., 2010; Saad et al., 2012; Murphy and Fox, 2017), which results in spurious fcMRI values. Moreover, it has been shown that SGSR do not consider the brain's spatial heterogeneities and biases correlations in different regions of the brain (Saad et al., 2012). Accordingly, the extracted correlation maps are known to present artifacts and do not reflect underlying neurological properties (Murphy et al., 2009; Anderson et al., 2010; Saad et al., 2012; Murphy and Fox, 2017). Therefore, regressing out GS is under debate as its removal by applying current approaches may introduce artifacts into the fMRI data or cause the loss of important neuronal components (Murphy et al., 2009; Anderson et al., 2010; Saad et al., 2012; Murphy and Fox, 2017). These concerns about the GSR methods and the need for accurate brain functional connectivity maps motivate the need to develop new methods for dealing with GS. Moreover, it has been shown that GS has a variety of sources with different spatial distributions which are voxel-specific. Accordingly, it is desirable to use a new method that works selectively and is able to identify and remove the spatially specific GS for each voxel or region (Saad et al., 2012; Tong and Frederick, 2014; Chang et al., 2016; Power et al., 2017; Turchi et al., 2018), and also produce known connectivity patterns in networks such as the default mode network and task-positive network (Fox et al., 2009; Erdoğan et al., 2016), thus avoiding the creation of spurious correlations.

In addition to GS, in fMRI studies, BOLD signal is low-pass filtered (<0.1 Hz) during the preprocessing procedure to be sure that the effects of the high frequency physiological noises are removed from the data (Boubela et al., 2013; Brooks et al., 2013; Liu et al., 2017). This is because, physiological noises which are mainly cardiac and respiratory, are spatially widespread and have cycles located prominently at the frequency range of 0.1–2.5 Hz. It is indicated that, among different noise-removal methods (such as band-pass filtering and Independent component analysis), EMD based methods facilitate the noise removal from fMRI data. In EMD-based methods, IMFs with specific frequency bands are identified and removed from fMRI data to enhance the functional sensitivity of the data (Typically the first and second IMFs which have the highest frequency bands among all IMFs are considered as a noise) (Lin et al., 2016). However, removing the whole high-frequency data from fMRI time series is controversial, as smoothing the signals via low-pass filtering decreases the signal to noise ratio by smoothing the peaks and amplifying the noise (Brooks et al., 2013). In fact, it has been shown that filtering high frequency modes may also remove the signal of interest that contains similar frequencies. The main reason is that the TR time for sampling fMRI data is too low to distinguish the high frequency components and causes signal's frequencies being aliased that can not be separated by temporal filtering (Brooks et al., 2013). Furthermore, even using very high sampling rate (TR < 0.4 s) to detect the high frequency modes may cause losing information of neuronal activation in high frequencies by filtering high frequency modes (Tagliazucchi et al., 2011, 2012; Boubela et al., 2013). Accordingly, in resting-state studies, we cannot do the band-pass filtering through previous methods as the brain dynamics in all frequency bands needs to be investigated. Therefore, we need a method that can remove the physiological noises more specifically from BOLD signal.

There are several signal processing methods, such as Fourier transform (Gallagher et al., 2008), Wavelet transform (Yves, 1993), spatial and temporal Blind source separation (Comon and Jutten, 2010), and the EMD (Huang et al., 1998). However, all of the mentioned method except EMD require predefined basis function or some prior knowledge to decompose the signal. Considering the fact that real-world signals including fMRI signals are non-linear and non-stationary data and do not perfectly obey our assumption, EMD method would be the best method to apply, as it does not need any basis functions and parameters that need to be adjusted such as wavelet type in wavelet transform or informed separation ideas in Blind source separation method (Liutkus et al., 2013; Riffi et al., 2014; He et al., 2017). EMD is a computationally efficient method that can adaptively decompose any non-linear and non-stationary signals into Intrinsic mode functions (IMF) and obtain meaningful frequencies estimation (Huang et al., 1998; Mandic et al., 2008; Riffi et al., 2014; He et al., 2017).

In this paper, we define an adaptive global signal regression (AGSR) by performing a spatiotemporal decomposition of the fMRI time series through EMD-based methods. The GS which is computed using this method is voxel-specific and depends on brain regions' heterogeneity.

Additionally, we show that by applying AGSR, we do not need the traditional low-pass filtering methods as the proposed method exhibits the potential to adaptively remove the physiological noises from high temporal frequency modes of fMRI time series, that are shared in whole brain regions. Therefore, AGSR method, besides removing the GS, helps to eliminate the high frequency physiological noises without needing to perform the low-pass filtering step separately.

In AGSR method, We do not use the Multidimensional EMD approach as it requires great runtime and cannot decompose a multidimensional signal into multidimensional components (Wu et al., 2009; Riffi et al., 2014; He et al., 2017). Consequently, in this paper, two EMD-based methods are used sequentially to decompose the fMRI signals adaptatively and voxel-specifically. We acquired the Spatial and Temporal Intrinsic Mode Functions (SIMF and TIMF, respectively) of fMRI data by applying FATEMD (Riffi et al., 2014) and ICEEMDAN (Colominas et al., 2014) methods, respectively (Huang et al., 1998; Mandic et al., 2008). It has been shown that applying EMD-based methods on fMRI data separate inherent brain oscillations and fundamental modes embedded in BOLD signal. Each of these oscillations occupies a unique frequency band and can be used to investigate the frequency characteristics in resting-state brain networks (McGonigle et al., 2010; Zheng et al., 2010; Niazy et al., 2011; Song et al., 2014, 2015; Qian et al., 2015; Lin et al., 2016; Cordes et al., 2018).

To explore the frequency characteristics of the brain networks, first, we obtain the average functional connectivity matrices for different TIMFs of each SIMFs over all subjects. Functional connectivity was computed using pearsons' coefficient between the peak voxels of each brain regions included in the AAL 116 atlas (Tzourio-Mazoyer et al., 2002).

We then compute the efficiency (Fair et al., 2007; Rubinov and Sporns, 2009; Cohen and D'Esposito, 2016) of the brain network of different spatiotemporal IMFs, which is used as a measure of integration. Integration values are used to identify the GS, since GS is defined as a synchronous fluctuation which is shared among all brain regions that makes it being highly integrated in the whole brain. Given the high values of efficiency in the low spatial maps (SIMF3, SIMF4, and SIMF5), our results suggest that these maps can be considered as spatial global signal masks. The performance of the proposed method is compared with the SGSR method, and also with the results when GS is not removed. This is done by investigating the functional connections within an extracted peak voxel of the known network's regions and the selected seed region. While the correlation pattern identified by the other methods suffers from a low level of precision, our method demonstrates a high level of accuracy due to its data-driven adaptive nature.

2. Methods

2.1. fMRI Data Acquisition

The resting-state fMRI preprocessed data-set of 21 subjects from the NIH Human Connectome Project (HCP) (https://db.humanconnectome.org) (Essen et al., 2013) is used in this research. Each subject was involved in 4 runs of 15 min each using a 3 T Siemens, while their eyes were open and had a relaxed fixation on a projected bright cross-hair on a dark background. The data were acquired with 2.0 mm isotropic voxels for 72 slices, TR = 0.72 s, TE = 33.1 ms, 1,200 frames per run, 0.58 ms Echo spacing, and 2,290 Hz/Px Bandwidth (Moeller et al., 2010). Therefore, the fMRI data were acquired with a spatial resolution of 2 × 2 × 2 mm and a temporal resolution of 0.72 s, using multibands accelerated echo-planar imaging to generate a high quality and the most robust fMRI data (Moeller et al., 2010). The fMRI data were preprocessed to remove spatial artifacts produced by head motion, B0 distortions, and gradient non-linearities (Jovicich et al., 2006). As comparison of fMRI images across subjects and studies is possible when the images have been transformed from the subject's native volume space to the MNI space (Evans et al., 1993; Ashburner and Friston, 1999), fMRI images were wrapped and aligned into the MNI space with FSL's FLIRT 12 DOF affine and then a FNIRT non-linear registration (Jenkinson and Smith, 2001; Jenkinson et al., 2002; Jahanshad et al., 2013). In this study, the MNI-152-2 mm atlas (Mazziotta et al., 1995, 2001a,b) was utilized for fMRI data registration.

2.2. Estimation of the Temporal IMFs (TIMFs)

EMD is an adaptive data-driven signal processing method, which does not need any prior functional basis such as the wavelet transform (Mandic et al., 2008). The basic functions are derived adaptively from the data by the EMD sifting procedure. The EMD method developed and established by Huang et al. (1998) decomposes non-linear and non-stationary time series into their fundamental oscillatory components, called Intrinsic Mode Functions (IMFs). There are two criteria defining an IMF during the sifting process: 1) the number of extrema and zero crossings must be either equal or differ at most by one, and, 2) at any instant in time, the mean value of the envelope defined by the local maximum and the envelope of the local minimum is zero. The EMD algorithm for estimating the IMFs of the time series x(t) is:

1. r₀(t) = x(t), j = 1.

2. For extracting the j-th IMF:

(a) h₀(t) = r_j(t), k = 1,

(b) Locate local maximum and minimum of h_k−1(t),

(c) Identify the average envelope using cubic spline interpolation to define upper and lower envelope of h_k−1(t),

(d) Calculate the mean value m_k−1(t),

(e) Put h_k(t) = h_k−1(t) − m_k−1(t),

(f) Check the stopping criteria. The stopping criteria determines the number of sifting steps to decompose an IMF Huang et al. (1998). If stopping criteria is satisfied then h_j(t) = h_k(t) otherwise, go to (a) to extract next IMF with k = k + 1.

3. r_j(t) = r_j−1(t) − h_j(t).

4. If at least two extrema were in the r_j(t), the next IMF is extracted otherwise the EMD algorithm is finished and r_j(t) is the residue of x(t). Accordingly, x(t) is defined as:

\begin{array}{l} x (t) = \sum_{j = 1}^{n} h_{j} (t) + r_{n} (t), & (1) \end{array}

where h_j(t) is the j-th IMF, n is the number of IMFs, and r_n(t) is the residue of the signal. Thus, the EMD method adaptively decomposes a time series into a set of IMFs and a residue where the first IMF (IMF1) corresponds to the fastest oscillatory mode and the last IMF (IMFn) to the slowest one, the sum of these components yields the original signal (Huang et al., 1998; Hassan and John, 2005). However, frequent occurrences of the mode-mixing phenomenon in analyzing real signals using EMD algorithm is problematic. To address this problem and improve the spectral separation of modes, the ensemble empirical mode decomposition (EEMD) method was proposed (Wu and Huang, 2009). This method extracts modes by performing the decomposition over an ensemble of noisy copies of the original signal combined with white Gaussian noises, and taking the average of all IMFs in the ensemble (Colominas et al., 2014).

The EEMD method solves the mode mixing problem, but certain issues remain. First, the number of IMFs extracted from each of the noisy signal copies is different, and this creates a problem when averaging the IMFs. The second problem is a reconstruction error in the EEMD method (Wu and Huang, 2009; Colominas et al., 2014). To fix this error the complementary EEMD (CEEMD) was proposed (Yeh et al., 2010). In the CEEMD algorithm, pairs of positive and negative white noise processes are added to the original signal to make two sets of ensemble IMFs. Accordingly, the CEEMD effectively eliminates residual noise in the IMFs which alleviate the reconstruction problem. Nonetheless, the problem of the different number of modes when averaging still persists. To overcome this problem, the CEEMD with adaptive noise (CEEMDAN) was developed (Torres et al., 2011; Colominas et al., 2014). In this approach, the first mode is computed exactly as in EEMD. Then, for the next modes, IMFs are computed by estimating the local means of the residual signal plus different modes extracted from the white noise realizations. CEEMDAN decomposition can create some spurious modes with high-frequency and low-amplitude due to overlapping in the scales. Additionally, some residual noise is still present in the modes. As a consequence, the new optimization algorithm, Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (ICEEMDAN), was proposed (Colominas et al., 2014).

During the sifting process using ICEEMDAN method the local mean of realizations is estimated, instead of using the average of modes from the first step. This change in the algorithm reduces the amount of noise present in the final computed modes. To deal with the issue of creation of spurious modes in the final results, ICEEMDAN method proceeds differently than the EEMD and CEEMDAN methods. In ICEEMDAN, white noise is not added directly; instead EMD modes of white noise are added to the original signal and to the IMFs during the sifting process (Wu and Huang, 2009; Colominas et al., 2014). Furthermore, in this method as in CEEMDAN, a constant coefficient is added to the noise that makes the desired signal to noise ratio between the added noise and the residue to which the noise is added. This coefficient is computed based on the standard deviation of the residue at each step of the sifting process. Therefore, the IMFs computed with ICEEMDAN have less noise and more physical content than IMFs obtained with other methods (Colominas et al., 2014) (More detailed description of ICEEMDAN method can be found at Colominas et al., 2014). The high accuracy rate, reduction in the amount of noise contained in the modes, and the alleviation of mode mixing phenomenon qualify this method to effectively decompose biological signals. In this paper the ICEEMDAN method with 300 ensembles and a level of noise of 0.2 (Wu and Huang, 2009) is used to extract the Temporal Intrinsic Mode Functions (TIMFs) from the fMRI data.

2.3. Estimation of the Spatial IMFs (SIMFs)

A fast, time efficient, and effective method is essential for processing real images that have a large size. Previous EMD-based methods were limited to small size images as the extrema detection, interpolation at each iteration, and the large number of iterations make their processing time consuming and complicated (Bhuiyan et al., 2008; Riffi et al., 2013, 2014; He et al., 2017). Therefore, those methods were just applicable to reduced size images, which resulted in losing some information during their process. Fast and Adaptive Tridimensional (3D) EMD, abbreviated as FATEMD, is a recent extension of the EMD method to three dimensions (Riffi et al., 2014). The FATEMD method is able to estimate volume components called tridimensional Intrinsic Mode Functions (3D-IMFs) quickly and accurately by limiting the number of iterations per 3D-IMF to one, and changing the process of computing upper and lower envelopes, which reduce the computation time for each iteration (Bhuiyan et al., 2008; Riffi et al., 2014; He et al., 2017). In the FATEMD method, the steps of extracting 3D-IMFs are almost the same as the previous EMD based methods, except for the number of iterations and the estimations of the maximum and minimum envelopes. The steps for decomposing a volume V(m, n, p) with dimensions m, n, and p using the FATEMD approach are as follows (Bhuiyan et al., 2008; Riffi et al., 2014):

1. Set i = 1 and R_i(m, n, p) = V(m, n, p).

2. Determine the local maximum and minimum values by browsing R_i(m, n, p) using a 3D window (cube) with a size of 3 × 3 × 3 which results in an optimum extrema maps (Map_max(m, n, p) and Map_min(m, n, p)). These local maximum (or minimum) values are strictly higher (or lower) than all of their neighborhoods contained in the cube.

3. Calculate the size of the Max and the Min filters which will be used in making extrema envelopes and their smoothness. The maximum and minimum filters are made by computing the nearest Euclidean distances between the maximum (d_adj.max) (minimum (d_adj.min)) points. The cubic window width (w_en) then is determined by using one of the following four formulae for both maximum and minimum filters. Here, we used the 4-th formula as outlined below, although using the other formulas will result in approximately the same decomposition result:

\begin{array}{l} w_{en} = min {min {d_{adj . max}}, min {d_{adj . min}}}, \\ w_{en} = min {max {d_{adj . max}}, max {d_{adj . min}}}, \\ w_{en} = max {min {d_{adj . max}}, min {d_{adj . min}}}, \\ w_{en} = max {max {d_{adj . max}}, max {d_{adj . min}}} . & (2) \end{array}

4. Create the envelopes of maxima and minima (Env_max(m, n, p) and Env_min(m, n, p)) of size (w_en).

5. Use the mean filter to compute the smoothed envelopes:Env_max−s(m, n, p) and Env_min−s(m, n, p).

6. Calculate the mean filter by averaging the smoothed upper and lower envelopes (Env_A(m, n, p)).

7. Calculate the i-th 3D-IMF: IMF_i(m, n, p) = R_i(m, n, p)−Env_A(m, n, p).

8. Calculate R_i+1(m, n, p) = R_i(m, n, p)−IMF_i(m, n, p).

9. If R_i+1(m, n, p) contains more than two extrema then

Go to the step 2 and set i = i + 1,

Else

The FATEMD decomposition is completed.

Therefore, FATEMD is an adaptive approach as all of the processes for computing filters and making the maximum, minimum, and the mean envelops are data driven. FATEMD decomposes a volume into a set of 3D-IMFs (Riffi et al., 2014). In general, a volume V can be reconstructed from the summation of the K 3D-IMFs and the residue as follows:

\begin{array}{l} V (m, n, p) = \sum_{i = 1}^{K} {IMF}_{i} (m, n, p) + R_{K + 1} (m, n, p) . & (3) \end{array}

K is the number of IMFs, and R(m, n, p) is the residue of the signal.

In this paper, we apply the FATEMD method at each time instant to decompose the resting-state fMRI data into tridimensional IMFs called Spatial Intrinsic Mode Functions (SIMF). Figure 1 shows the spatial decomposition results of a sample resting-state fMRI image. The ICEEMDAN method is then utilized to decompose each SIMF into its corresponding TIMFs.

FIGURE 1

Figure 1. Spatial decomposition of a sample fMRI image using FATEMD method. The original fMRI image at one TR time is decomposed into 5 SIMFs (SIMF1 to SIMF5) and a residue.

2.4. Spatiotemporal Pattern Analysis of the fMRI Data

To define an adaptive and voxel-specific GS, the spectral information of fMRI data is investigated by constructing the functional connectivity matrices using extracted TIMFs and SIMFs data. To fulfill this aim, first, the SIMFs of the fMRI data at each TR time are computed by applying the FATEMD method, then, all spatial components are merged together in time to construct the time series of each SIMF. Second, the peak voxel at each region, that is, the voxel of maximal activation, is selected by computing the Root Mean Square (RMS) for each voxel's signal over all time. It has been shown that peak voxel provides the best effect of any voxel in the ROI (Sharot et al., 2005). Additionally, the peak voxel activity correlates better with evoked scalp electrical potentials than approaches that average activity across the ROI. This means that the peak voxel represents the ROI's activity better than other choices (Arthurs and J Boniface, 2003). The peak voxel in each region is determined using previously published Talairach coordinates (after conversion to MNI coordinates and using AAL 116 atlas) (Fox et al., 2005). After determining the peak voxels of each region, the ICEEMDAN method is applied to its time series to compute the TIMFs. Thus, the TIMFs of all regions for each SIMF are computed.

We then compare the predefined distinct frequency bands presented in fMRI studies (slow5 [0.01–0.027 Hz], slow4 [0.027–0.073 Hz], slow3 [0.073–0.198 Hz], slow2 [0.198–0.25 Hz], and slow1 [0.5–0.75 Hz]) (Penttonen and Buzsáki, 2003; Zhan et al., 2014), to the frequency content of the extracted TIMFs. In all subjects, TIMFs consistently corresponded to the same frequency bands. As seen in the Figure 2, the frequency range comprised in TIMF1 to TIMF3 is approximately the same as the frequency range of the sum of slow1 to slow3. The frequency range of TIMF4 is the same as slow4, and the frequency range of the sum of TIMF5 to TIMF9 has the same frequency range as the slow5 frequency band. Accordingly, we label the summation of TIMF1 to TIMF3 as TIMF1, TIMF4 as TIMF2, and the summation of TIMF5 to TIMF9 as TIMF3. Figure 3 represents the pipeline used in computing SIMFs and TIMFs for each resting-state fMRI data. Accordingly, the functional connectivity matrices are constructed by computing the average of correlation coefficients between all possible pairs of TIMFs correspond to different Spatial domains for all brain regions comprised in the AAL 116 atlas over all 21 subjects. Consequently, instead of the classical functional connectivity matrix, the decomposition presented here produces 5 × 3 = 15 connectivity matrices (each with size 116 × 116), 3 temporal domains and 5 spatial domains, encompassing the rich spatiotemporal dynamics of brain activity.

FIGURE 2

Figure 2. Temporal IMFs and their corresponding frequency spectrum of the sample SIMF time series. (A) Sample of SIMF time series before temporal decomposition (B) 9 decomposed TIMFs of a sample SIMF by applying the ICEEMDAN method with 300 ensembles and a level of noise of 0.2. (C) The 9 decomposed TIMFs are divided into three different frequency bands. According to slow1 to slow3 and slow5 frequency bands defined in the literature, TIMFs1 to 3 and 5 to 9 are combined, respectively. (D) Represents the frequency spectrum of the 9 TIMFs. (E) The frequency spectrum of TIMFs in part (C) that correspond to frequency bands used in the literature for slow1 to slow5 Penttonen and Buzsáki (2003), Zhan et al. (2014). TIMF, Temporal Intrinsic Mode Function; ICEEMDAN, Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise.

FIGURE 3

Figure 3. Pipeline for computing spatial and temporal IMFs (SIMF and TIMF) of the fMRI data. (A) A sample of fMRI data. (B) Splitting each fMRI data in time at each TR time. (C) SIMFs at each TR time which are computed by applying FATEMD approach. (D) Shows the AAL 116 atlas used after merging SIMFs in time to select the peak voxel of each region. (E) Time series of all brain ROIs for each SIMF. (F) The TIMFs' time series of a sample SIMF for one ROI computed by using ICEEMDAN approach. (G) Summation of time series of the TIMFs in (F) based on frequency bands of slow1 to slow5 defined in the literature. The summation of the TIMF1 to TIMF3, TIMF4, and the combination of TIMF5 to TIMF9 are labeled as TIMF1, TIMF2, and TIMF3 in the rest of the paper, respectively. rfMRI, resting-state fMRI; TIMF, Temporal Intrinsic Mode Function; SIMF, Spatial Intrinsic Mode Function; ICEEMDAN, Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise; FATEMD, Fast and Adaptive Empirical Mode Decomposition; ROI, Region of Interest.

2.5. Topological Properties of the Brain Network

The GS is a synchronous fluctuation shared among all brain regions. Consequently, the GS component in the brain connectivity matrix should present a high integration value, where integration is the topological property of a network that describes how information from distributed brain regions is combined (Fair et al., 2007; Rubinov and Sporns, 2009; Cohen and D'Esposito, 2016). To compute the integration of the brain network at different spatiotemporal scales we use the global efficiency measure (Fair et al., 2007; Rubinov and Sporns, 2009). The global efficiency is computed as the average inverse shortest path length between all the node pairs of the network that is normalized by the maximal number of network's links. Therefore, the weighted global efficiency is computed via the following equation:

\begin{array}{l} E^{w} = \frac{1}{N (N - 1)} \sum_{j = 1}^{N} \sum_{i = 1, i \neq j}^{N} {({d_{i j}}^{w})}^{- 1}, & (4) \end{array}

where N is the number of nodes in the network and d_ij is the minimum path length between nodes i and j (Fair et al., 2007; Rubinov and Sporns, 2009). The shortest path length is computed by counting the smallest number of edges needed to get from node i to node j which is inversely related to node weight. The information needed to estimate the weight of all pairs of brain regions are provided by functional connectivity matrices (Rubinov and Sporns, 2009; Cohen and D'Esposito, 2016), strong association between regions has a large weight which leads to a shorter length. When two nodes are disconnected the length of that path would be infinite and correspondingly, the efficiency would be zero (Fair et al., 2007; Rubinov and Sporns, 2009).

3. Results

3.1. Defining Adaptive Global Signal (AGS)

We computed the functional connectivity matrices between all pairs of brain regions for different spatiotemporal domains extracted from fMRI data for each subject. Figure 4 shows the average connectivity matrices computed by Pearson's coefficient over the 21 subjects. As seen in the figures, SIMF1 and SIMF2 in all TIMFs showed low connectivity whereas SIMF3 to SIMF5 in all TIMFs showed high connectivity. Besides, they indicate that the magnitude of the correlation does not significantly depend on the TIMFs. Thus, based on the connectivity strength for different spatiotemporal domains, the summation of the SIMF1 to SIMF2 and the SIMF3 to SIMF5 including all TIMFs, were considered as two separate signals. We also averaged the six connectivity matrices resulting from the summation of TIMF1 to TIMF3 with SIMF1 and SIMF2 (Figure 4) and labeled it as AGSR (Figure 5A), and the nine connectivity matrices resulting when combining TIMF1 to TIMF3 with SIMF3 to SIMF5, which we labeled as AGS (Figure 5B).

FIGURE 4

Figure 4. Functional connectivity matrices of the whole brain regions using AAL 116 atlas for different spatial and temporal IMFs. Pearson's correlation coefficient (r) with P ≤ 0.01 is computed between all the brain regions' spatiotemporal domains extracted from fMRI data. Spatial domains are extracted by applying FATEMD method on fMRI signal. The three temporal domains including TIMF1, TIMF2, and TIMF3 are computed by applying ICEEMDAN on each SIMF. SIMF, Spatial Intrinsic Mode Function; TIMF, Temporal Intrinsic Mode Function; ICEEMDAN, Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise; FATEMD, Fast and Adaptive Empirical Mode Decomposition.

FIGURE 5

Figure 5. Average functional connectivity matrices of the whole brain regions using AAL 116 atlas over all subjects. (A) Average connectivity matrix of fMRI data applying AGSR which means the connectivity matrix of combination of SIMF1 and SIMF2 including all TIMFs of the fMRI data, (B) connectivity matrix of the AGS which is the combination of SIMF3 to SIMF5 including all TIMFs. AGSR, Adaptive Global Signal regression; AGS, Adaptive Global Signal; SIMF, Spatial Intrinsic Mode Function; TIMF, Temporal Intrinsic Mode Function.

We also computed the global efficiency (Figure 6A) for different spatial and temporal IMFs using Equation (4) and also based on functional connectivity results. Figure 6A shows that there are high values of efficiency in the low frequencies of spatial domains, SIMF3, SIMF4, and SIMF5, which indicate active shared connections between all the nodes in the brain, suggesting the existence of GS in the low-frequency spatial domains, called Adaptive Global Signal(AGS). Furthermore, SIMF3 to SIMF5 with high temporal frequency mode (TIMF1) which is included in the AGS can be considered as an adaptive filter to reduce the effects of the highly integrated physiological noises in high frequency modes instead of applying low-pass filtering (Shmueli et al., 2007; Boubela et al., 2013; Liu et al., 2017).

FIGURE 6

Figure 6. Integration of the brain network at different spatiotemporal scales. (A) Average efficiency of the whole brain network for different spatial and temporal IMFs defined in functional connectivity. (B) Comparing the magnitude of average efficiency of the brain network over all subjects when the AGS is removed from the fMRI time series (the sum of SIMF3 to SIMF5 in all TIMFS are removed and the sum of SIMF1 to SIMF2 including all TIMFs of the fMRI signal are considered to compute the connectivity), and the average efficiency of the AGS (summing up SIMF3 to SIMF5 in all TIMFs). High efficiency values in the SIMF3 to SIMF5 which represent the AGS in the fMRI data are seen in the figures. GS, Global Signal; AGS, Adaptive GS.

As seen in Figure 6B and Table 1, the high values of integration of AGS (summation of SIMF3 to SIMF5 including all TIMFs) confirm that they can be considered as a GS which has to be removed from the fMRI data to have more accurate brain connectivity results. In the last results' section (represented in Figures 8, 9) we show that, including low frequency spatial domains may cause spurious connectivity results between brain regions.

TABLE 1

Table 1. Integration of AGSR and AGS. The average efficiency of the brain network over all subjects, when the AGSR are performed, and the average efficiency of the AGS. AGS, Adaptive Global Signal; AGSR, AGS Regression.

3.2. Regressing Out the AGS and SGS From fMRI Data

According to the definition of AGS, for each brain voxel signal, there is a corresponding AGS while the SGS is common for the whole brain voxels. The AGS for each voxel is computed by summing up the SIMF3, SIMF4, and SIMF5 with all TIMFs while the SGS is computed by taking the average of all brain voxels' time series. It should be noted that in computing AGS, the residues of spatiotemporal decomposition of the fMRI data are added to the last TIMF and SIMF. The three time courses in Figures 7A–C correspond to the AGS, the fMRI sample time course [the peak voxel's time course in Medial Prefrontal cortex (MPF) ROI], and the conventional or Static GS (SGS), respectively. Figures 7D,E show resting-state fluctuations of the sample fMRI time series from MPF ROI after regressing out (subtracting) the AGS and the SGS. It also has to be mentioned that to be consistent with the previous fMRI studies, data are conventionally low-pass filtered except when the AGSR method is applied.

FIGURE 7

Figure 7. AGSR and SGSR of a sample fMRI data. (A) The voxel-specific AGS of Medial Prefrontal (MPF) cortex region and (B) the original fMRI time series of the peak voxel in MPF cortex region. (C) The SGS which is common for all region's voxels. (D,E) Show the time series with the SGSR and AGSR, respectively. These time series are computed by subtracting the AGS and SGS from the original time series. MPF, Medial Prefrontal cortex; AGS, Adaptive Global Signal; SGS, Static(conventional) Global Signal; AGSR, AGS Regression; SGSR, SGS Regression.

3.3. Connectivity Map of Task-Positive and Task-Negative Networks

The default mode network or Task Negative Network (TNN) is a state of brain activation whereby the individual is not attending to any external cues in the environment but certain brain regions are still activated and they are less active during task performance rather than during the resting-state. It has been shown that (Fox et al., 2005) the default mode network responses are significantly activated in three of the seeded regions: the Posterior Cingulate Cortex (PCC), Medial Prefrontal cortex (MPF), and Lateral Parietal cortex (LP). The efficacy of our approach is examined by computing the connectivity map. In computing functional connectivity maps, we computed Pearson's correlation which is popular in fMRI studies and also allows our findings to be comparable with other papers to test the validity of the proposed method. We computed the average connectivity between the time course of the PCC region as a seed region and the main regions of the Task Positive Network (TPN) which are the Middle Temporal (MT), right Frontal Eye Field (FEF), left Intraparietal Sulcus (IPS), Supplementary Motor Area (SMA), Inferior Parietal Lobule (IPL), Visual regions, and the left Auditory region and the TNN ROIs which are MPF, PCC, and left LP which includes the Angular Gyrus, Hippocampus, and Cerebellar tonsils ROIs (Fox et al., 2009; Erdoğan et al., 2016).

Considering the AGS definition, the combination of the SIMF1 and SIMF2 was used to compute the functional connectivity between PCC and TNN and TPN including visual ROIs by using Pearson's correlation coefficient (r), P ≤ 0.01. Figure 8 is functional connectivity brain map for different brain layers along the Z axis which show the mean connectivity over all subjects between brain regions and the PCC ROI as a seed region when the AGSR, NR, and the SGSR are performed.

FIGURE 8

Figure 8. Comparing the average functional connectivity between the PCC ROI as a seed region and the brain ROIs using the AAL 116 atlas for fMRI data of all subjects. The average functional connectivity applying (A) AGSR, (B) NR, and (C) SGSR. Slices shown in the maps are at Z = 09, 15, 25, 35, 45, 55, 65, 75, respectively. AGSR, Adaptive Global Signal Regression; NR, No Regression; SGSR, Static (conventional) Global signal regression.

Figure 9 shows expected average connectivity between the PCC ROI and different regions of the TPN and the TNN (positive correlation between the PCC and the TNN and negative correlation between the PCC and TPN) applying the new approach of GSR in resting-state fMRI data.

FIGURE 9

Figure 9. The average connectivity map between the PCC as a seed region and (A) the TPN and (B) TNN ROIs for fMRI data of all subjects. Connectivity results of applying AGSR and SGSR are shown in green and red, respectively, and the blue ones are the results of computing connectivity without applying any GSR (NR). Connectivity map is made by computing Pearson's correlation coefficient (r) with P ≤ 0.01 between the PCC region as a seed region and the main regions of the TPN and TNN. PCC, Posterior Cingulate Cortex; MPF, Medial Prefrontal cortex; LP, Lateral Parietal cortex; MT, Middle Temporal; FEF, Frontal Eye Field; IPS, Intraparietal Sulcus; SMA, Supplementary Motor Area; IPL, Inferior Parietal Lobule; ROI, Region Of Interest; NR, No Regression; AGS, Adaptive Global Signal Regression; SGS, Static (conventional) Global Signal Regression.

While the NR and SGSR (conventional GSR which is based on averaging) are unable to identify the expected connectivity in some regions for TPN and TNN ROIs, the AGSR approach obtains expected functional connectivity for all regions in TNN and TPN which confirms the effectiveness of the proposed method for GSR (Figure 9). As AGSR is an adaptive and voxel-specific method, we have a unique local signal for each voxel which by being removed from fMRI data augments the precision of the rsfc-MRI results.

4. Discussion

In contrast to previous works (Zarahn et al., 1997; Fox et al., 2009; Liu et al., 2017), the present study provides a new method for GSR, called AGSR, that works voxel-specifically and adaptively. It is believed that fMRI data are a superposition of the GS and network-specific fluctuations. However, the main reason for the controversy over the use of GSR in fMRI studies is that the average-based GS is a mixture of signals from multiple brain regions without considering the possibility of spatial heterogeneity in the GS (Fox et al., 2009; Murphy et al., 2009; Weissenbacher et al., 2009; Saad et al., 2012; Murphy and Fox, 2017). It has been shown that regressing out average-based GS results in negative correlations that do not have a biological basis and are artifacts in the voxels' time series which lead to distortion in the connectivity results or activation measures (Fox et al., 2009; Murphy et al., 2009; Murphy and Fox, 2017). In this paper, we showed that the AGSR method works voxel-specifically and can compute the neuronal correlations of the brain's networks more accurately. This is because using the FATEMD method in computing AGS maximizes the spatial contributions to the GS. In other words, decomposing fMRI data in space using the FATEMD approach, which is done by considering features of each voxel's neighbors, makes the computed AGS sensitive to brain regions' heterogeneity.

When assessing the efficiency for different spatiotemporal domains of the fMRI data, no large differences in different temporal IMFs at the same spatial IMF were obtained. Thus, we concluded that the variability of efficiency is just related to the spatial frequency domains. The high values of the efficiency in the low spatial frequencies demonstrated the existence of the GS. On the other hand, high spatial frequencies, SIMF1 and SIMF2, represented the most network-specific data. Accordingly, the low spatial frequencies, SIMF3 to SIMF5 including all TIMFs, were considered as the AGS.

Additionally, it has been shown that motion, cardiac, and respiratory noise components which have high frequency cycles and are spatially coherent, cause spatially widespread fluctuations in the BOLD signals that contribute to the global signal (Shmueli et al., 2007; Liu et al., 2017). Conventionally, filtering the high frequency components of the fMRI data to remove above mentioned physiological noises and the GSR are done separately as two preprocessing steps in fMRI studies (He and Liu, 2012; Caballero-Gaudes and Reynolds, 2017; Liu et al., 2017), however, common low-pass filtering methods through removing high frequency components cause missing a considerable amount of information on resting-state brain functional network (Tagliazucchi et al., 2011, 2012; Boubela et al., 2013; Turchi et al., 2018). In our proposed method, in addition to GSR, physiological noise components that are common across voxels and are mainly included in the high frequency modes are also removed from the data by removing the SIMF3 to SIMF5 of TIMF1 through AGSR. Thus, our proposed method, through AGSR, filters the highly connected part of high frequency modes adaptively without applying low-pass filter separately. It can help to provide more informative data by involving high frequency modes in the data.

We examined the efficacy of our method by computing the seed-based functional connectivity for the TPN and TNN regions. Our results in agreement with previous studies (Chang and Glover, 2009; Fox et al., 2009; Chai et al., 2012), show that the negative correlations are intrinsic to the brain and do not appear just as a result of the GSR. We found that the AGSR method identifies the connectivity between the TPN and TNN regions according with the expected results of prior studies (Fox et al., 2005, 2009). We compared the connectivity results of the AGSR with the SGSR and when there is NR in the fMRI data. Despite the connectivity results of the SGSR method and when there is NR, applying our proposed method resulted in an enhancement to the detection of network-specific fluctuations of the brain. Furthermore, although the strength of the correlations is related to cognitive function, in auditory regions, lower activity seen in the result of applying AGSR appears to be related to the better removal of the acoustic noise heard by subjects during fMRI. This shows that the acoustic noise of the fMRI device which is almost constant in all TR times and interferes with auditory system activity can be removed better through AGSR (Ravicz et al., 2000; Moelker and Pattynama, 2003). Thus, it is inferred from the results that AGSR method is able to remove physiological and remained systemic noises after preprocessing more correctly and without introducing artifactual correlations as confirmed by correlations between PCC and the reference regions.

In conclusion, AGS is a unique local signal for each voxel's BOLD signal. In the AGSR method, the first and second spatial IMFs of each fMRI data, decomposed by FATEMD method, are simply summed up to have a band-pass filtered fMRI data without GS. AGSR is a reliable method that works voxel-specifically for all subjects which leads to provide information about brain function with more accuracy. There are some limitations to the methods used in this study that should be noted. Although the FATEMD and ICEEMDAN are optimized approaches for finding the best IMF sets, they still need more improvement in the sifting procedure to yield better decomposition performance. For instance, finding the optimum values of added white noise and the ensemble number to overcome the mode mixing problem and speed up the calculation in ICEEMDAN approach are two drawbacks of this approach.

We computed the GS for each region of the AAL 116 atlas specifically, however, as this method has a “voxel-specific” nature, it can be applied to all voxels of the brain. Computing voxel-specific GS just needs more memory and computer power, such as a larger computer cluster but no additional changes to the underlying algorithm are needed. It is more feasible to compute the AGSR for all the voxels when we are interested in some specific regions of the brain and not the whole brain.

Therefore, the proposed method in this paper provides the opportunity to characterize the whole brain function and reflect the intrinsic property of the spatiotemporal nature of the fMRI data through removing the voxel-specific GS and not removing the whole high frequency modes. Future studies can be devoted to the application of our proposed method to the other image processing areas.

Author Contributions

NM and RS designed the research. NM analyzed the data, interpreted the results and wrote the main manuscript. MD assisted with analysis and interpretation of data. RS supervised the development of the work, helped in data interpretation and manuscript evaluation.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to Doug Phillips for generous assistance in computational work on Compute Canada cluster and Jordan Chad for proofreading the manuscript. This work was partially supported by grant RGPIN-2015-05966 from Natural Sciences and Engineering Research Council of Canada. Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. The computational work was enabled by support provided of WestGrid (www.westgrid.ca) through Compute Canada (www.computecanada.ca).

References

Anderson, J. S., Druzgal, T. J., Lopez Larson, M., Jeong, E. K., Desai, K., and Yurgelun-Todd, D. (2010). Network anticorrelations, global regression, and phase-shifted soft tissue correction. Hum. Brain Mapp. 32, 919–934. doi: 10.1002/hbm.21079

PubMed Abstract | CrossRef Full Text | Google Scholar

Arthurs, O., and J Boniface, S. (2003). What aspect of the fmri bold signal best reflects the underlying electrophysiology in human somatosensory cortex? Clin. Neurophysiol. 114, 1203–1209. doi: 10.1016/S1388-2457(03)00080-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, J., and Friston, K. J. (1999). Nonlinear spatial normalization using basis functions. Hum. Brain Mapp. 7, 254–266. doi: 10.1002/(SICI)1097-0193(1999)7:4<254::AID-HBM4>3.0.CO;2-G

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhuiyan, S. M. A., Adhami, R. R., and Khan, J. F. (2008). Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation. EURASIP J. Adv. Signal Process. 2008:728356. doi: 10.1155/2008/728356