Energy-Period Profiles of Brain Networks in Group fMRI Resting-State Data: A Comparison of Empirical Mode Decomposition With the Short-Time Fourier Transform and the Discrete Wavelet Transform

Traditionally, functional networks in resting-state data were investigated with linear Fourier and wavelet-related methods to characterize their frequency content by relying on pre-specified frequency bands. In this study, Empirical Mode Decomposition (EMD), an adaptive time-frequency method, is used to investigate the naturally occurring frequency bands of resting-state data obtained by Group Independent Component Analysis. Specifically, energy-period profiles of Intrinsic Mode Functions (IMFs) obtained by EMD are created and compared for different resting-state networks. These profiles have a characteristic distribution for many resting-state networks and are related to the frequency content of each network. A comparison with the linear Short-Time Fourier Transform (STFT) and the Maximal Overlap Discrete Wavelet Transform (MODWT) shows that EMD provides a more frequency-adaptive representation of different types of resting-state networks. Clustering of resting-state networks based on the energy-period profiles leads to clusters of resting-state networks that have a monotone relationship with frequency and energy. This relationship is strongest with EMD, intermediate with MODWT, and weakest with STFT. The identification of these relationships suggests that EMD has significant advantages in characterizing brain networks compared to STFT and MODWT. In a clinical application to early Parkinson’s disease (PD) vs. normal controls (NC), energy and period content were studied for several common resting-state networks. Compared to STFT and MODWT, EMD showed the largest differences in energy and period between PD and NC subjects. Using a support vector machine, EMD achieved the highest prediction accuracy in classifying NC and PD subjects among STFT, MODWT, and EMD.


SUPPLEMENTARY MATERIAL
Suppl. S1. Peak frequency, full-width-at-half-maximum (FWHM), frequency range, mean energy density, and period of the first nine intrinsic mode functions (IMFs) for the default mode network (DMN). Thalamus, Cerebellum_6, Cerebellum_Crus1 Note: Blank entries are unknown networks that have not been reported or investigated in the literature.

S4.1 Detrending
IMFs are also useful for local de-trending of time series. For example, the IMF s for k > 9 have frequencies less than 0.01 Hz (using our acquisition protocol). The local trend of data can be computed by ∑ IMF and is certainly more accurate than a trend estimated by Fourier or other stationary filtering methods, as was recently shown (Kaleem and Cordes, 2016). However, using EMD to estimate IMFs with frequencies less than 0.01Hz is more involved because different voxel time series may have frequency components less than 0.01Hz decomposed into IMFs with different indices. This makes detrending with EMD more complicated than a simple high-pass filtering using regression with cosine basis functions (as we have done in this study). Alternatively, the continuous wavelet transform could also be used for high-pass filtering with cut-off frequency of 0.01Hz. We refrained from studying these very low drift frequencies with EMD or the continuous wavelet transform in this project and defer to such a study in a future application.

S4.2 Global Signal Regression
Conventionally, global signal regression uses the average of whole-brain time series as a nuisance regressor in a linear model to remove the effect of the global signal. However, this regression method does not allow for spatial variations of the global signal and treats every voxel the same, i.e. every voxel has the same regressor. Moradi et al. (2019) recently applied a spatiallyadaptive EMD approach called Fast and Adaptive Tridimensional (3D) EMD (FATEMD) (Riffi et al., 2015) for removing the global signal by decomposing the fMRI data first into 5 spatial IMFs (SIMFs) and associating the spatial components SIMFs with indices 3-5 as spatially adaptive global signal masks. Then, by summing up the maps SIMF 1 and SIMF 2 , the global signal can be excluded from the data.

S4.3 Single-Scale Time-Dependent (SSTD) Window-Sizes in Sliding-Window Dynamic Functional Connectivity (dFC) Analysis
In sliding-window dFC analysis, a fixed window-size is usually used and heuristically selected since no consensus exists yet on the choice of the optimal window-size. Recently, Zhuang et al. (2020) proposed the SSTD window-sizes computed from all EMD IMFs as optimal window-sizes in the sliding-window dFC analysis and applied this method to a large group of professional fighters with and without cognitive impairment, and healthy controls. Specifically, at every time point, an SSTD window-size is computed as the average of the instantaneous period over all IMFs weighted by the corresponding energy. In this case, SSTD window-sizes are based on the frequency content at every time point and are able to capture more temporal dynamic information. Both a higher classification accuracy in predicting cognitive impairment status in fighters and a larger explained behavioral variance in normal controls were found when using dynamic FC matrices computed with SSTD window-sizes as features, as compared to using dynamic FC matrices computed with conventional fixed window-sizes.

S4.4 Analytic Extension of the EMD Method
Variational Mode Decomposition (VMD) is an analytic extension of EMD (Dragomiretskiy and Zosso, 2014). VMD can be thought of as a more formal EMD-type decomposition technique, but as with all such techniques that utilize such variations of mathematical formulation of EMD-type behavior, certain parameters always must be adjusted which in our view defeats the purpose of a completely model-free and data-adaptive decomposition method. Resting-state fMRI applications have been investigated with VMD by Yuen et al. (2019) for TR=2s data and for TR=323ms data after low-pass filtering with cutoff 0.25Hz so that IMFs can be compared for different TRs. Reproducibility of IMFs and corresponding frequency clusters were found across fMRI trials. Similar to Song et al. (2014) and our studies, Yuen et al. (2019) found cortical restingstate networks to be best represented by low frequency oscillations. VMD has also been used to avoid potential mode mixings (Yuen et al., 2019). A comparison of VMD with EMD is beyond the scope of the current study.

Suppl.S5. Mode Mixing in EMD and Applying EMD to ICA time series
Mode mixing in EMD refers to an intermittent signal, for example high frequency oscillations that occur in a limited time segment of otherwise low-frequency oscillations throughout the IMF. An IMF that shows mode mixing gives rise to significant non-orthogonality to the other IMFs. Mode mixings have been reported in task-activated data (Lin et al., 2016), but in resting-state fMRI data, mode mixing may not be easily observable, since no a priori knowledge of the correct signal is known. In our analysis, we have not observed mode mixing in our ICA time series. Without ICA, however, EMD can be more difficult to apply directly because of possible mode mixings that may occur for those voxel time series in which more than one network is active (Fosso (2019), Lin et al. (2016) ). For example, if a voxel time series contains signals belonging to two different networks that operate at different frequencies in a dynamic fashion, a decomposition by EMD at a specific time point could result in an IMF order that is different for different time points (for example IMF ( becomes IMF ( ).
Without ICA, applying EMD (or related methods) directly on resting-state fMRI data can be difficult for the purpose of extracting meaningful intrinsic components and it has been speculated that the energy-period distributions would be close to white noise distributions (Lin et al., 2016). In this scenario, Ensemble EMD (EEMD) and related methods (complementary EEMD) can lead to improvements in computing IMFs by avoiding mode mixings (Qian et al., 2015(Qian et al., , 2018. Using EEMD, noise is injected (added) into the data and multiple EMD runs are carried out for different injected noise realizations. The IMFs obtained from each run are ensemble averaged to eliminate the effect of noise. It has been shown that for a large number of noise realizations (≅60 or more), robust IMFs can be obtained that are more orthogonal to each other. Non-orthogonal contributions are evidence of under-sampled processes, thus, if non-orthogonal IMFs are found, the data segment may be too small and EEMD may provide a better estimation of IMFs. EEMD is computationally more expensive than EMD, however and furthermore, the amount of noise that should be injected for fMRI data is unknown and may have to be estimated by applying different models to the data. As a consequence, EEMD cannot be categorized as a model-free method (for detailed information on EEMD, see Huang and Shen (2014)). For completeness, we would like to mention that EEMD has been further improved by Complementary EEMD (CEEMD) (Yeh et al., 2010), then by CEEMD with Adaptive Noise (CEEMDAN) (Torres et al., 2011) and later on by Improved Complementary EEMD with Adaptive Noise (ICEEMDAN) (Colominas et al., 2014;Moradi et al., 2019). A comparison of these advanced EMD methods with our method is beyond the scope of the current study.
Different temporal characteristics of Fourier (STFT) components, wavelet (MODWT) components, and EMD components (IMFs) as measured by the log vs. log relationship for six resting-state networks. The given numbers indicate the mean difference (over subjects) in log between PD and NC. Only values in magnitude of at least 0.1 are shown (which was only possible in EMD analysis (20 instances) and MODWT (2 instances) but not in STFT (0 instances)). A "*" indicates a significant difference between PD and NC with p<0.05 (uncorrected). The horizontal and vertical bars indicate the standard deviation in log and log , respectively. Note that for most IMFs of the PAR, CCN, PFC, lFPN and rFPN, the mean period is larger for PD (circles) versus NC (squares), indicating that the decomposition of these networks by EMD leads to IMFs containing lower frequencies in PD than NC. Furthermore, most low frequency IMFs (i.e. IMFs with index >1) have lower energy for PD than NC.