Advances in functional magnetic resonance imaging data analysis methods using Empirical Mode Decomposition to investigate temporal changes in early Parkinson's disease.

INTRODUCTION: Previous neuroimaging studies of Parkinson's disease (PD) patients have shown changes in whole-brain functional connectivity networks. Whether connectivity changes can be detected in the early stages (first 3 years) of PD by resting-state functional magnetic resonance imaging (fMRI) remains elusive. Research infrastructure including MRI and analytic capabilities is required to investigate this issue. The National Institutes of Health/National Institute of General Medical Sciences Center for Biomedical Research Excellence awards support infrastructure to advance research goals.
METHODS: Static and dynamic functional connectivity analyses were conducted on early stage never-medicated PD subjects (N = 18) and matched healthy controls (N = 18) from the Parkinson's Progression Markers Initiative.
RESULTS: Altered static and altered dynamic functional connectivity patterns were found in early PD resting-state fMRI data. Most static networks (with the exception of the default mode network) had a reduction in frequency and energy in specific low-frequency bands. Changes in dynamic networks in PD were associated with a decreased switching rate of brain states. DISCUSSION: This study demonstrates that in early PD, resting-state fMRI networks show spatial and temporal differences of fMRI signal characteristics. However, the default mode network was not associated with any measurable changes. Furthermore, by incorporating an optimum window size in a dynamic functional connectivity analysis, we found altered whole-brain temporal features in early PD, showing that PD subjects spend significantly more time than healthy controls in a specific brain state. These findings may help in improving diagnosis of early never-medicated PD patients. These key observations emerged in a Center for Biomedical Research Excellence-supported research environment.

Introduction

Functionally related regions of the resting brain show a high degree of temporal correlation in blood-flow fluctuations, as measured by the blood-oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) signal [1]. Using either seed-based methods or data-driven approaches such as independent component analysis (ICA), brain regions that fluctuate in synchrony and constitute reliable and reproducible functional networks in the human resting brain can be identified [2], [3], [4], [5], [6], named resting-state networks. Resting-state networks were assumed to be static in nature in the past, and approximately a dozen of such static networks have been discovered and investigated in relation to how they are impacted by neurodegenerative disorders. However, more recently, it has been shown that resting-state networks are dynamic in character and change on a time scale of several seconds to a minute [7], [8], [9], [10]. Analyzing the temporal dynamics of resting-state connectivity provides a more accurate picture of the working brain and can help in the early detection of neurological disorders and in monitoring effects of potential treatments.

Both static and dynamic analysis methods have been applied to study resting-state functional networks in major neurodegenerative diseases, for example, Alzheimer's disease (AD). One of the major brain networks affected in AD is the so-called default mode network (DMN), which is heavily involved in memory formation and retrieval [11]. In normal subjects, the DMN shows functional connections between the posterior cingulate cortex, angular gyrus, hippocampus, and the medial prefrontal cortex. In AD patients, amyloid-beta (Aβ) protein has been found to accumulate in DMN and other regions, which may disrupt connections and lead to the symptoms of memory and cognitive impairment [12], [13], [14]. Early Aβ accumulation is associated with reduced static functional connectivity within the DMN and between the DMN and the frontoparietal network (FPN), a network that is involved in attention-demanding tasks [15], [16]. The dynamic aspect of the DMN shows significant changes in AD as well [17]. It has been reported that AD subjects spend less time in brain states with strong posterior DMN contributions and more time in states with dorsal medial prefrontal cortex contributions.

Parkinson's disease (PD) is, after AD, the second most common neurodegenerative disorder in the elderly and is characterized by degeneration of dopaminergic neurons in the substantia nigra pars compacta with resulting striatal dopaminergic deficiency [18]. Previous neuroimaging studies of PD patients have shown that whole-brain functional networks such as the DMN and networks involving the motor pathway are affected, leading to different functional connectivity patterns when compared to those found in normal controls (NC) [19]. Studies of the temporal characteristics of fMRI resting-state brain networks have also shown abnormal spontaneous low-frequency content in PD [20]. The dynamic aspects of brain networks have been widely studied using electrophysiological recordings. Intraoperative electrophysiological data have shown that the occurrence of motor symptoms in PD is associated with changes in synchronizations within and between brain regions and changes in phase-amplitude coupling between brain regions [21], [22]. However, it is not clear if static changes in resting-state networks are present in the very early stages (first 3 years) in drug-naïve never-medicated patients with PD. Furthermore, whether changes in temporal dynamics occur in resting-state functional networks in de novo PD subjects is unknown.

In the present study, as part of the National Institutes of Health/National Institute of General Medical Sciences Centers of Biomedical Research Excellence grant to the Center for Neurodegeneration and Translational Neuroscience, we investigated low-frequency BOLD fluctuations of major resting-state networks in early PD using data from the Parkinson's Progression Marker Initiative (www.ppmi-info.org). Previously, frequency-specific analysis of resting-state networks has been carried out using bandpass filtering in which the frequency intervals were specified using information from electrophysiological data [23] or simply by dividing the possible frequency range into equal intervals that were specified by the user [24], [25]. An alternative approach toward finding frequency intervals in resting-state data is by Empirical Mode Decomposition (EMD) [26], [27]. EMD is a data-adaptive analysis method for studying the naturally occurring frequency bands in time series [28]. EMD can be used, in particular, for nonstationary signals and allows the decomposition of time series into nearly orthogonal modes spanning narrow frequency bands. The oscillatory modes are called intrinsic mode functions (IMFs) and are obtained by a sifting algorithm. The novelty of our EMD approach lies in the adaptive decomposition of fMRI data using EMD and identification of resting-state networks based on energy and period (inverse of frequency) characteristics of IMFs. These novel energy-period relationships of resting-state networks in PD may allow use of imaging biomarkers in characterizing or detecting PD in the early stages of the disease. Early stage identification of PD may improve diagnostic accuracy, enrollment in clinical trials of disease-modifying agents, and allow for more effective treatments.

In a second investigation, we explored the dynamic aspects of functional connectivity using the same data set. In previous research studies, dynamic functional connectivity analysis was carried out mainly by using a sliding-window method, in which pairwise linear correlations among network components are captured in subsequent temporal windows with a fixed window size and further clustered into multiple dynamic functional brain states [29], [30]. To find an appropriate window size is challenging because the windows size should be small enough to capture existing temporal transients and large enough to produce stable results [10]. The EMD method, however, provides us an alternative way to compute a time-dependent optimum window size in the sliding-window analysis. IMFs obtained from EMD track local periodic changes of nonstationary time series and an optimum window size can be determined at each time point. We incorporated the optimum window size in the sliding-window method to explore dynamic functional connectivity within and between the major resting-state networks in early PD.

Materials and methods

Empirical mode decomposition and intrinsic mode functions

EMD is a method defined by an algorithm to decompose any time series, whether nonstationary or nonlinear, into a set of IMFs. This decomposition is based on local characteristics of the time series, which transforms instantaneous amplitude and instantaneous frequency information into meaningful quantities to be computed by the Hilbert transform of the IMFs. Although Fourier and wavelet transforms use preassigned basis functions, the EMD basis functions are data-derived IMFs. The EMD method operates at the scale of one oscillation and is fully data-driven. An IMF represents a simple oscillatory mode but is more general than a harmonic function of one frequency component. An IMF can have variable amplitude and frequency along the time axis.

In general, an IMF is a function that must satisfy two conditions: (1) For the entire time series, the number of extrema and the number of zero crossings must be equal to or can differ at most by one and (2) the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero at every time point. A signal x(t) can be decomposed in terms of its K IMFs f(t) bywhere K is the number of IMFs, f(t) is the k-th IMF, and r(t) is a small monotone residual (trend) function. The basic algorithm for obtaining the decomposition is an iterative sifting algorithm described in detail in the study by Huang et al. [28]. This iterative process sequentially explores the natural constitutive scales of a time series. The IMF with index 1 (IMF1) contains the highest frequencies, and IMF with index K contains the lowest frequency components. It has been shown that the frequency arrangement in IMFs mimics that of a dyadic filter bank [31]. Instantaneous frequency and amplitude of the IMFs can be computed by extending the signal into the complex plane with the Hilbert transform. In the literature, EMD combined with the Hilbert transform is referred to as the Hilbert-Huang Transform.

Energy versus period relationship of intrinsic mode functions

The time series in fMRI data are known to contain structured as well as white noise sources. Because the IMFs are basis functions that are derived from the data rather than functions that satisfy given analytic expressions, it is important from a statistical perspective to understand the IMFs of data that contain only noise sources so that IMFs of noisy signals can be compared with IMFs of pure noise data. Comparison with artificial noise data provides a reference standard of results obtained by EMD and allows a statistical significance to be associated to IMFs. Of particular importance is the relationship of the mean energy as a function of the mean inverse of the frequency (mean period) for each IMF. The mean energy per unit time, E, of the k-th IMF, f(t), is defined by the mean instantaneous squared amplitude of the IMF. This definition leads towhere N is the number of data points and Δt is the sampling time which is equal to the TR in fMRI. The mean period, T, is defined by the mean value for the inverse of the instantaneous frequency obtained from the Hilbert transform, that is,

However, owing to outliers in the estimation of the instantaneous frequency spectrum, especially for frequencies close to zero, Eq. (3) does not provide robust values of T. Instead, we determine the density of v(t) using kernel density estimation (with a Gaussian kernel). For white Gaussian noise, it has been shown that [32]

Thus, y = log(E) as a function of x = log(T) is distributed approximately along the diagonal line y = −x for all IMFs of white noise data.

Time-dependent window size in dynamic functional connectivity analysis

The time-dependent window size at each time point of two fMRI time series x1 (t) and x2 (t) can be computed from the instantaneous period p(t) and average energy density E of the k-th IMF. The instantaneous periods p(t) capture the local nonstationarity of the original signal [33], and the average energy densities E summarize the energy contributions of each IMF to the original signal. A time-dependent period, T(t), for each time course is then determined as an average of p (t) weighted by E, that is,

The final time-dependent window size, T(t), of x1 and x2 to obtain an optimal sliding-window correlation is chosen to be the maximum of T(t) and T(t) (where T(t) and T(t) are the time-dependent window sizes of x1 and x2, respectively) to ensure that the data in one instantaneous period are included in calculating the correlation coefficient [34]. Thus, T(t) captures the local nonstationarity of original time courses, summarizes different contributions of each IMF to original signals, ensures that the data in one period are included in calculating the correlation coefficient, and therefore is optimum in capturing temporal dynamics between fMRI time series, as compared to a fixed window size. Dynamic functional connectivity analysis is then carried out using a sliding-window approach with this time-dependent window size T(t).

Simulation

The simulation aimed at demonstrating that the sliding-window correlation computed with a time-dependent window size can capture local transients and avoid unstable fluctuations, as compared to the correlation values computed using a fixed window size. To provide a specific example, two nonstationary time seriesandwere simulated with a sample rate of 1 second (TR = 1 s) with a static correlation coefficient of −0.02. The dynamic correlations between y1 and y2 were calculated and compared using the sliding-window method with two fixed window sizes, 10 TR and 50 TR, as well as the time-dependent window size T(t).

Participants

The data used in this study were obtained from the publicly available anonymized Parkinson's Progression Marker Initiative database [Marek et al., 2011]. We included 18 NCs (14 male (M); age: 64.25 ± 9.78 years (mean ± SD); years of education: 16.72 ± 2.67 years) and 18 newly diagnosed, early stage, and never-medicated PD subjects (10 M; age: 57.11 ± 11.63 years; years of education: 17.00 ± 2.77 years; disease duration: 0.83 ± 0.84 years) in our analysis. A chi-square test was performed to check statistical significance for gender difference between the two groups, and Wilcoxon rank-sum test was performed to check for differences of age and year of education. Differences in age (P = .06), gender (P = .15), and years of education (P = .76) were not significant between the two groups.

Additional information about the distribution of participants with hypercholesterolemia, hypertension, and diabetes was also obtained from the Parkinson's Progression Marker Initiative database and listed in Table 1. A chi-square test was performed to check statistical significant difference for each disease distribution between the PD and NC groups. Differences in the distribution of participants with hypercholesterolemia (P = .72), hypertension (P = .07), and diabetes (one in the NC group only) were not significant between the two groups.

Table 1

The distribution of participants with hypercholesterolemia, hypertension, and diabetes between the two study groups

Number of subjects has the condition	PD	NC	Differences (P-value)
Hypertension	3	8	.07
Hypercholesterolemia	5	6	.72
Diabetes	0	1	N/A

Abbreviations: PD, Parkinson's disease; NC, normal controls.

MRI data acquisition and preprocessing steps

All subjects underwent resting-state fMRI scans on 3T Siemens scanners. The resting-state fMRI involved an 8 minutes and 24 seconds echo-planar acquisition with 210 time points (TR = 2,400 ms, TE = 25 ms, field of view = 22.4 cm, flip angle = 80°, resolution = 3.3 × 3.3 × 3.3 mm3, and 40 axial slice). In addition, a T1-weighted structural image was also acquired for each subject (TR = 2,300 ms, TE = 2.98 ms, flip angle = 9°, and voxel size = 1 × 1 × 1 mm3).

The first 5 time points (12 seconds) were removed to allow the MR signal to achieve T1 equilibrium. Echo-planar data were slice-timing corrected and realigned to the mean echo-planar image in Statistical Parametric Mapping 12 (http://www.fil.ion.ucl.ac.uk/spm/), further coregistered to the subject T1 space, and then normalized to the standard Montreal Neurological Institute-152 2-mm template using Advanced Normalization Tools software (http://stnava.github.io/ANTs/). Six head motion parameters, signals extracted from subjects' white matter and cerebrospinal fluid (3-mm cubes centered at Montreal Neurological Institute [26, −12, 35] and [19, −33, 18]), were regressed out from each data set. fMRI data were further spatially smoothed using an 8-mm 3D Gaussian filter.

Static analysis of resting-state networks

To obtain the spatial resting-state networks, group ICA [35] (based on the FastICA algorithm [36]) was performed by stacking all data in the temporal domain to obtain 30 resting-state networks. Then, spatial regression was used on the networks of the group time series data to obtain the time series signatures for NC and PD. To get more detail on the time signatures of the resting-state networks, we decomposed the corresponding time signatures using EMD into the first 5 IMFs for each spatial resting-state network. These 5 IMFs covered a frequency range from 0.01 Hz to the Nyquist frequency (0.5/TR) of the data. For each IMF, the average instantaneous energy, period, and their standard deviations were computed for NC and PD.

Dynamic functional connectivity analysis

Dynamic functional connectivity analysis was carried out using a sliding-window approach with the optimum time-dependent window size as determined by EMD. Specifically, whole-brain dynamic functional connectivity was captured by computing the windowed correlation between time courses from every pair of nodes, in which each node is a network component. To obtain network components, another group ICA with 100 components was performed by stacking data from both PD and NC subjects. Seventy-two ICA components were visually identified as network-related components, and the corresponding subject-specific ICA maps and time courses were calculated using dual regression [37]. A voxel-wise comparison was conducted for each of the 72 ICA components. Two sample t-tests were performed with age and gender as covariates, and each ICA component was spatially masked with the thresholded group ICA component map.

The connectivity matrices for each subject were then calculated using sliding-window correlations between each pair of nodes with the time-dependent window size T(t). The connectivity matrices in each sliding window (size, 72 × 72) were concatenated in time for each subject and further stacked for both PD and NC subjects. Standard k-means clustering was performed in MATLAB (www.mathworks.com) on the concatenated connectivity matrix from all subjects to estimate dynamic functional states for both the groups. The optimum cluster number K was determined by a leave-one-out cross-validation. Finally, the time spent in each state and the frequency-of-state alternations were calculated for every subject separately and used to compare the temporal dynamics between the PD and NC groups. Two sample t-tests were carried out with age, gender, and the distribution of hypertension as covariates for these comparisons. To compare results with traditional methods which are previously published, the same analyses were also repeated with fixed window sizes of 13 TRs (∼30 seconds), 25 TRs (∼60 seconds), and 50 TRs (∼120 seconds), as suggested in other studies [10], [29].

Results

Static functional connectivity analysis in early PD

Fig.1 shows the spatial maps and corresponding temporal IMFs computed for the DMN. A t-test showed no significant spatial differences in this network between NC and PD. IMF1 of PD patients shows some variation in the amplitude of the time series signal. All other IMFs have similar characteristics for the same index in NC and PD participants. We calculated the energy and period for each IMF and plotted this information using group-specific markers in a log (energy) versus log (period) diagram. Standard deviations of the markers are indicated by horizontal and vertical lines for log (period) and log (energy), respectively. We found no significant difference in energy or period for any of the IMFs of the DMN for NC versus PD.

Fig. 1

Top: The spatial DMN determined by group ICA for NC and subjects with Parkinson's disease in the early stages (PD). There are no significant differences in the spatial distribution of this network between NC and PD. Middle: Time signatures of the DMN for NC and for patients with PD were decomposed into five IMFs. Bottom: For each IMF, the average energy and period were calculated and displayed in a log (energy) versus log (period) diagram. The diagonal line indicates the expected mean of Gaussian white noise. Different markers specify the IMF properties for NC (squares, connected by blue dotted line) and for PD (solid circles, connected by red dashed line). The horizontal and vertical bars through the markers indicate the standard deviation in log (period) and log (energy), respectively. Abbreviations: ICA, independent component analysis; DMN, default mode network; NC, normal controls; PD, Parkinson's disease; IMF, intrinsic mode function.

In Figs. 2 and 3, we show all six resting-state networks (out of 30) in which the period of the IMF (with the same index) differed by a large effect size (Cohen's d > 0.8 [38]) between NC and PD. The obtained networks are the executive control network (ECN), the parietal network (PAR), the cognitive control network (CCN), the prefrontal cortex network (PFC), and the left/right frontoparietal network (lFPN and rFPN). All these networks show spatial and temporal differences. The ECN, CCN, and PFC have reduced activations in PD, whereas the PAR, lFPN, and rFPN have spatially extended activations in PD.

Fig. 2

Six ICA resting-state networks show different spatial patterns for NC and early PD. Note that the ECN and PFC networks show reduced spatial activity, whereas the PAR and FPN show increased activity in PD. The CCN has decreased activity in the hippocampus but increased activity in the inferior temporal lobes in PD. Abbreviations: ICA, independent component analysis; FPN, frontoparietal network; NC, normal controls; ECN, executive control network; PAR, parietal network; CCN, cognitive control network; PFC, prefrontal cortex network; PD, Parkinson's disease.

Fig. 3

Different temporal characteristics of IMFs as measured by the log (energy) versus log (period) relationship for the six resting-state networks in Fig. 2. All relationships shown with a “*” indicate a large effect size (Cohen's d > 0.8) (either for energy or for period) for specific IMFs. The ECN was the only network found where the period for PD was reduced (in IMF5). For all other networks that show a large effect size between log (period) of PD and NC (i.e., PAR, CCN, PFC, lFPN, rFPN), the mean period is always larger for PD versus NC, indicating that these networks operate at lower frequencies in PD. The letter T in the relational statements indicates the period. Abbreviations: PD, Parkinson's disease; lFPN, left frontoparietal network; rFPN, right frontoparietal network; NC, normal controls; IMFs, intrinsic mode functions; ECN, executive control network; PAR, parietal network; CCN, cognitive control network; PFC, prefrontal cortex network.

The temporal characteristic of these networks differ; the ECN has increased frequency content (less period) for the very low–frequency band in IMF5, which is in the drift range (f < 0.01 Hz), whereas all other networks show a decrease in low frequencies for some of the higher bands (IMFs with index ≤ 4). We found that the period for the same indexed IMF is always larger for PD, for all networks except the ECN, irrespective of the effect size being small (Cohen's d = 0.2) or large (Cohen's d = 0.8). The amplitude of oscillations of the IMFs as measured by log (energy) is generally smaller for PD for most of the IMFs.

Dynamic functional connectivity analysis in early PD

To illustrate the advantage of using a time-dependent sliding-window size, Fig. 4A shows simulated nonstationary time series y1 (blue) and y2 (red). Instantaneous periods of y1 (dashed blue) and y2 (dashed red) and the time-dependent window size T(t) between y1 and y2 (solid green) are plotted in Fig. 4B at every time point. As shown in Fig. 4C, dynamic correlations between y1 and y2 calculated with the time-dependent window size (solid green) capture existing local transients without creating unstable fluctuations, as compared to the correlation computed with the fixed window sizes of 10 TRs (10 seconds, dashed blue) or 50 TRs (50 seconds, dashed red line).

Fig. 4

Simulation. (A) Two nonstationary simulated time series: y1 (blue) and y2 (red). (B). Instantaneous period of y1 (dashed blue) and y2 (dashed red) and the time-dependent window size (green) at each time point. (C). The dynamic correlation between two synthetic time series was calculated using sliding-window method with different window sizes: 10 TR (dashed blue), 50 TR (dashed red), and time-dependent window size (solid green). The static correlation between these two time series is indicated by the dashed purple line.

Real fMRI data

A flow chart of the dynamic connectivity analysis is shown in Fig. 5. Eight major resting-state networks, formed by 72 network-related ICA components, were investigated, including the subcortical network, auditory network, sensorimotor network, visual network, CCN, DMN, medial temporal network, and cerebellum network (Fig. 6A). Both the ECN and FPN obtained from the other ICA run with 30 components are combined to be the CCN in the dynamic functional connectivity analysis. A group comparison was conducted for each of the 72 ICA components, and no significant spatial difference was found at a family-wise corrected error rate of P < .05. A static correlation matrix between each pair of nodes was converted to Fisher's z statistics and shown in Fig. 6B. Three dynamic functional states are determined from the leave-one-out cross-validation in k-means clustering for both PD and NC subjects (Fig. 7A) with the time-dependent window size. The average window size is listed in Table 2. Fig. 7B and C show that NC subjects spend significantly more time in state II, which has stronger connections, both between and within networks, whereas PD patients tend to stay longer in the more weakly connected functional states I and III. Furthermore, a significant reduced frequency-of-state alternation (P = .006) is found in the PD group (Fig. 7D). The same analysis was repeated with a fixed window size of 30 seconds, 60 seconds, and 120 seconds. The number of dynamic states determined by cross-validation and the between-group comparison results of frequency-of-state alternation are listed in Table 2.

Fig. 5

Whole-brain dynamic functional connectivity analysis flow chart. Abbreviation: ICA, independent component analysis.

Fig. 6

(A) Spatial components from 72 network-related ICA components. (B) Whole-brain static functional connectivity matrix. Abbreviations: ICA, independent component analysis; AUD, auditory network; SMN, sensorimotor network; Vis, visual network; CCN, cognitive control network; DMN, default mode network; MTN, medial temporal network; CBN, cerebellum network.

Fig. 7

(A) Three dynamic states were obtained from whole-brain dynamic functional connectivity analysis with the time-dependent window size. (B) Average window fractions of time spent in each state (window in each state/total number of windows) for PD (orange) and NC (blue) groups. (C). Statistical comparisons of window fractions spent in each state between PD and NC groups; *P < .05. (D). Statistical comparison of frequency-of-state alternations between PD and NC groups. Abbreviations: PD, Parkinson's disease; NC, normal controls.

Table 2

Number of dynamic states computed from cross-validation and the difference of frequency-of-state alternation between the two groups obtained with the time-dependent window size and fixed window size used in previous studies

Window size	Number of dynamic states	Difference of frequency-of-state alternation between NC and PD (P value)
30 s	11	.05
60 s	7	.10
120 s	3	.24
Time-dependent window size (97.75 s ± 41.36 s)	3	.006

Abbreviations: PD, Parkinson's disease; NC, normal controls.

Discussion

Static resting-state analysis

In this study, we developed a novel method for the identification of abnormal temporal signatures associated with brain states in early PD using resting-state fMRI data. We used EMD as a data-adaptive method to determine energy and period characteristics of temporal signatures of major cortical networks that were obtained with ICA. The advantage of using EMD is that the temporal signatures can be decomposed into basic modes (namely the IMFs) that are subject-specific, and each basic mode can be characterized by energy density and period content. This analysis is different from a Fourier or wavelet analysis because there are no parameters that need to be adjusted (such as predefined frequency intervals or wavelet types), and the analysis is completely data-driven. We have obtained consistent features of energy and period for all subjects as shown by the small standard deviations about the mean values for average energy and period content. EMD shares a frequency decomposition feature with the discrete wavelet transform in which both methods exhibit a dyadic filter bank decomposition. However, the discrete wavelet transform has a frequency decomposition in the sense of components at different but fixed nonadaptive frequency scales. Using a discrete wavelet transform, energy can be calculated based on wavelet coefficients, and period information from the fixed frequency ranges can be obtained. However, the period information is nonadaptive to individual subject data, and one would obtain the same period information for both PD and NC. An alternate approach would be the use of a continuous wavelet transform, but a relationship between its decomposition level and frequency is not directly defined. For these reasons, EMD is a superior time series analysis method because it determines subject-specific energy densities and periods of fundamental modes.

Comparison with other studies

FMRI resting-state data have several advantages over other modalities such as fluorodeoxyglucose (FDG) positron emission tomography (PET) imaging to detect characteristic features in early PD. In previous studies using FDG-PET imaging in early PD, three disease-specific spatial covariance patterns were found [39], namely the PD motor-related pattern, the PD cognition-related pattern, and the PD tremor-related pattern. Of particular importance is the cognition-related pattern that showed metabolic reductions in preSMA, medial prefrontal cortex, precuneus, and metabolic increases in the cerebellum. Similar results with fMRI have been shown in resting-state data [40]. However, PD-specific resting-state networks were found using a spatial analysis method and not by characteristics of temporal signatures of resting-state as we have proposed here. Furthermore, the spatial covariance pattern was not obtained for never-medicated first-3-year PD patients. For this reason, previous spatial fMRI results are difficult to compare with those of our temporal analysis.

Compared with [15O] H2O PET, fMRI has greater temporal resolution, slightly greater spatial resolution, and greater sensitivity. However, the BOLD signal is more difficult to relate to dopamine deficiency, and how blood flow, blood volume, and oxygen consumption (which lead to the BOLD signal) are related to neurotransmitters such as dopamine is currently unknown. In a recent study, static spatial brain networks obtained by FDG-PET were compared with resting-state fMRI in a larger cohort of healthy nondemented subjects from the Alzheimer's Disease Neuroimaging Initiative database [41]. Most networks obtained showed similar spatial covariance patterns in FDG-PET and fMRI. However, discrepancies were observed for some important networks. For example, some of the anterior-posterior networks (e.g., the DMN and lFPN) could only be partially obtained by FDG-PET. Furthermore, reduced correlations were observed in anterior-posterior correlations in FDG-PET when compared to BOLD networks, which may indicate different signal mechanisms of metabolism for FDG-PET and BOLD in fMRI. This discrepancy could also arrive due to the different temporal scales of the imaging modalities (seconds for fMRI and minutes for FDG-PET), which lead to different inter-regional couplings [41]. The detection of temporal changes of brain states' dynamic functional connectivity using FDG-PET may not be feasible because of the lower temporal resolution of FDG-PET.

We demonstrated that in early never-medicated PD, there are no spatial or temporal changes detectable in the DMN. However, we found 6 other cortical brain networks that showed spatial and temporal differences of resting-state signal characteristics. These networks are the ECN, PAR, CCN, PFC, lFPN, and rFPN which show overlap with the cognitive resting-state pattern in previous PD studies [39], [40]. We have analyzed the first 5 IMFs of the associated temporal profiles in terms of content in average energy density and average period. We found that in early PD, the PAR, CCN, PFC, lFPN, and rFPN networks are driven by reduced frequencies (increased period) for all IMFs in the low-frequency range of less than 0.1 Hz (which is covered by IMFs with k = 2, 3, 4, 5). Several of the IMFs also showed a large effect size for reduced frequency content. In addition, most corresponding energy densities were lower in PD. There was only one network, namely the ECN, with different but not significant characteristics for frequencies above the drift range (f > 0.01 Hz). Overall, most networks in early PD were characterized by a reduction in frequency and energy in specific low-frequency bands of less than 0.1 Hz as determined by EMD. For future studies, it may be interesting to study early PD patients with mild cognitive impairment to see if DMN abnormalities can be detected because mild cognitive impairment occurs in approximately 1 of 4 PD patients.

To date, the sliding-window method with a fixed window size is most commonly used for examining dynamics in resting-state functional connectivity [9], [29], [30]. Using simulation, we have demonstrated that compared with the fixed window size, the time-dependent window size computed from instantaneous period of IMFs can more precisely capture the local periodic changes without creating unstable fluctuations. This advantage of using time-dependent window size in dynamic functional connectivity analysis is further demonstrated using real fMRI data from PD and NC groups (Table 2). Significant (P < .01) reduced frequency-of-state alternations in the PD group is found when a time-dependent window size is used, which is not observed with the fixed window size.

Altered functional dynamics in PD

In our analysis, three whole-brain functional dynamic states are found for both PD and NC subjects. As shown in Fig.7A, robust within-network functional connectivity is observed in all three states, whereas stronger between-network functional connectivity is observed in state II, as compared to state I and III. Our results indicate that NC subjects stay significantly longer in state II (Fig. 7B), which is consistent with previous findings of increased between-network functional connectivity in healthy aging [42]. The altered dynamics of basal ganglia-cortical circuits in PD subjects have been widely reported using electrophysiological data. Specifically, unmedicated PD subjects exhibit aberrant coherent activity patterns and excessive synchronization of neuronal activities in the basal ganglia-cortical loop, which will in turn affect the neuronal circuits' dynamics [21], [43], [44], [45], [46]. Using public neuroimaging data, we observe a significant reduced frequency-of-state alternation among the three dynamic functional states, which demonstrates the limited dynamic range of whole-brain functional connectivity in treatment-naive PD subjects.

Conclusions

We used EMD to study the energy and period content of IMFs for static resting-state networks in early PD and found a reduction in both IMF frequency content and energy for PAR, CCN, PFC, lFPN, and rFPN. In contrast, ECN showed increased low-frequency content. For the DMN, no spatial or temporal changes were observed. We further studied the dynamic functional connectivity in the same cohort using a sliding-window method with a time-dependent window size computed from IMFs and obtained using EMD. Altered temporal behaviors and reduced whole-brain temporal dynamics were found in early PD subjects.

Systematic review: The authors reviewed the literature using PubMed. The relevant citations are appropriately cited. Static and dynamic functional connectivity in Parkinson's disease subjects has not been systematically study with Empirical Mode Decomposition.

Interpretation: Our findings lead to a characterization of resting-state data in early never-medicated Parkinson's disease patients by providing static and dynamic imaging markers of functional connectivity using Empirical Mode Decomposition.

Future directions: The manuscript proposes new techniques to assess static and dynamic functional connectivity using Empirical Mode Decomposition which needs to be further validated in a larger cohort of PD subjects.