Data used in the preparation of this article were obtained from the Parkinson’s Progression Markers Initiative (PPMI) database,¹ which is a multicenter international study. For currently updated information on this study visit.² All data used in our study were downloaded before January 2021. For detailed inclusion criteria of the PPMI dataset, please visit (see text footnote 2).

A total of 101 PD patients scanned fMRI were recruited in PPMI dataset. We selected patients with below inclusion criteria to our study: (1) visited for at least 4 years; (2) scanned fMRI at baseline (PD-BL) and followed-up (PD-FU); (3) older than 50 years old, as pathogenesis of late-onset PD and early onset PD is different (Schrag and Schott, 2006; Fereshtehnejad and Posteuma, 2017) and this study focused on late-onset PD. We focused on the late-onset PD, because late-onset PD population was larger than early onset PD patients. Late- and early onset PD were separated in many studies (Rango et al., 2021; Sigirli et al., 2021). The patients with low image quality or diagnosed with other diseases, such as Alzheimer’s disease, multiple system atrophy, dementia with Lewy bodies, during visiting period were excluded from this study. Imaging data and information of Healthy control (HC) subjects were also downloaded from PPMI dataset, who were scanned fMRI. In total, data of 20 PD patients and 19 healthy control subjects were utilized in our study.

The disease severity of PD patients was assessed using the Movement Disorder Society Unified Parkinson Disease Rating Scale (MDS-UPDRS) (Goetz et al., 2008), Part I of the UPDRS was used to evaluate non-motor symptoms of PD, and parts II and III evaluated motor symptoms. Total UPDRS was sum of part I-IV. Hoehn and Yahr (H and Y) staging (Hoehn and Yahr, 2001) were used to assess severity of PD. The UPDRS-III and H and Y staging scores used in this study were assessed at “ON” state, which was defined as the last dose of levodopa or dopamine agonist were taken <6 h before assessment. The “OFF” state meant the last dose of levodopa or dopamine agonist were taken ≥6 h before assessment. The drug treatment dosage was calculated as the levodopa-equivalent daily dose (LEDD) (Tomlinson et al., 2010). Patients’ cognitive states were assessed using the Montreal Cognitive Assessment (MoCA) scale, Benton Judgment of Line Orientation Score (BJLOT), Letter Number Sequencing Score (LNS), Semantic Fluency Total Score (SFT), Symbol Digit Modalities Score (SDM), and Hopkins Verbal Learning Test–Revised (HVLT) scales. The total SFT score and subtest such as animal, fruit and vegetable scores were calculated in this study. HVLT contains several aspects, such as discrimination scores, immediate/recall scores, retention score, false alarms, delayed recall and delayed recognition. The ability of daily living (ADL) scale was used to evaluate life quality. Depression and anxiety were assessed by the Geriatric Depression Scale (GDS) score and State-Trait Anxiety Inventory (STAI). The University of Pennsylvania Smell Identification Test (UPSIT) was used to evaluate olfactory nerve dysfunction. Sleepiness quality and rapid eye movement (REM) sleep were assessed by the Epworth Sleepiness Scale (ESS) score and REM Sleep Behavior Disorder Questionnaire score. Autonomic function was assessed by the Scales for Outcomes in Parkinson’s Disease-Autonomic dysfunction (SCOPA-AUT).

T1 MRI and rs-fMRI data were collected on Prisma fit 3-tesla Siemens MR scanner (Siemens, Erlangen, Germany). All participants were requested to remain calm and not think during scanning. T1 MRI data were acquired with the following parameters: repetition time (TR) = 2,300 ms, echo time = 3.0 ms, flip angle = 9°, matrix size = 176 × 240 × 170, and voxel size = 1 mm × 1 mm × 1 mm. The rs-fMRI data, consisting of 210 volumes per patient, were acquired with TR = 2.4 s, echo time = 25 ms, flip angle = 80°, matrix size = 68 × 66 × 40, and voxel size = 3.3 mm × 3.3 mm × 3 mm.

The preprocessing of the fMRI data was conducted using the functional connectivity (CONN) toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012),³ based on MATLAB (R2018a, MathWorks, Inc., Natick, MA, USA). The first ten scans were removed to maintain the stability of fMRI data. Thus, each participant provided 200 volumes for further analysis.

A slice-timing correction step was performed to correct the image acquisition time between slices. The fMRI data were realigned to the first volume to correct for head movement and then co-registered to the T1 MRI data using an affine transformation. The T1 MRI data were normalized to standard Montreal Neurological Institute (MNI) space and segmented into gray matter, white matter, and cerebrospinal fluid using the tissue probability maps. A non-linear transformation was applied to the functional data with the same parameters as the T1 MRI data. The functional data were smoothed with a Gaussian kernel of 6 mm full width half maximum, to reduce the noise signal.

Independent component analysis (ICA) grouping was performed to obtain specific data components using Group ICA of the functional MRI Toolbox (GIFT v4.0 c).⁴ Voxel-level variance normalization was performed on all data. The number of independent components was estimated using a minimum description length approach. Functional data were decomposed using two data reduction steps (subject-specific and group-level PCA). The Infomax method was used to calculate independent components, and the ICASSO method implemented in GIFT was used to maintain reliability.

ICASSO was run 100 times, and 80% similarity was selected. Significant components from all independent components were chosen according to the location of the peak coordinate and time course characteristics (Allen et al., 2014; Kim et al., 2017). The selected components were identified as seven resting state networks (RSN): the default modal network (DMN), attention network (ATT), basal ganglia network (BG), visual network (VIS), sensorimotor network (SMN), auditory network (AUD), and frontoparietal network (FP) (Yeo et al., 2011).

The time courses of the selected components were extracted. Correlations between every pair of selected components formed the subject-specific static FC matrix. Static FC group comparisons were calculated among HC, PD-BL, and PD-FU using the multivariate analysis of covariance (MANCOVA) toolbox contained in GIFT with a corrected false discovery rate (FDR). Post hoc analyses were performed between HC and PD-BL and between PD-BL and PD-FU.

Dynamic FC analysis was conducted using the dynamic functional network connectivity toolbox (dFNC) contained in the GIFT toolbox. The time courses of selected components were linear detrended and 3D-despiked. The time courses were also filtered using a low-pass filter with a cutoff of 0.15 Hz, as recommended in the GIFT manual.

A sliding window approach was used to analyze dynamic functional connectivity changes. The sliding window method is the most commonly used method for dynamic FC analysis, and a 30-60 s window size is appropriate, according to previous studies (Allen et al., 2014; Kim et al., 2017). The sliding window size was set to 22 TR (52.8 s) and the step was set to 1 TR (2.4 s), based on previous research (Allen et al., 2014; Kim et al., 2017). As a result, each participant had 178 windows across the entire scan. Correlations between the time course of each selected component pair in each window formed dynamic FC matrices. The K-means clustering method was performed to sort all dynamic FC matrices into different clusters, also known as states. A penalty on the L1 norm was imposed with 100 repetitions, and the elbow method was used to decide the clustering number.

Temporal properties of the dynamic FC analysis, such as fractional windows, mean dwell times, and numbers of transitions, were calculated for further statistical analysis. The fractional window is the proportion of time spent in each state, and the mean dwell time is the average time spent by each participant in each state. The number of transitions is calculated as the number of times switching between states.

Group differences in general information between the HC and PD groups, such as age, gender, years of education, and clinical assessments, were calculated using the two-sample t-test and the chi-square test. Differences in temporal properties among HC and PD subgroups were calculated using two sample t-test or rank sum test. Statistical differences of clinical assessments and temporal properties between PD-BL and PD-FU were tested using paired t-test or Wilcoxon signed rank test. Selecting t-test or non-parametric test was depended on the result of normality test. Correlations between temporal properties and clinical assessment scores were analyzed using partial correlation analysis, which controlled effect of age, sex, and education years. Above statistical analysis was performed in SPSS Statistics 22.0 (IBM Corporation, Armonk, NY, USA), and the threshold for statistical significance was p < 0.05. Multiple comparisons of p-values were corrected with FDR correction.

Twenty PD patients (13 male and 7 female) and 19 healthy subjects (15 male and 4 female) met the inclusion and exclusion criteria described in (Section “2.1. Participants”). There were no significant differences in age and gender between the PD and HC groups. The years of education of the PD group were statistically shorter than those of the HC group. After years of follow-up, the H and Y stage (ON state) of the PD group increased. The LEDD of PD-FU was significantly higher than PD-BL. These results indicated that the severity of PD increased after years of follow-up. The UPDRS-I scores of PD-FU increased compared to PD-BL. The HVLT-False alarms and SDM scores of PD-FU were smaller than HC group. The SCOPA-AUT scores of PD-FU were bigger than HC group. The UPSIT scores of PD-BL were smaller than HC group. Above results indicated that the non-motor symptoms of PD patients were worsen after years of follow-up. Other clinical assessments, such as UPDRS-III (ON state), total UPDRS, MoCA, BJLOT, ESS, GDS, SFT, and STAI scores, did not significantly change between PD-BL and PD-FU. Details of these results are shown in Table 1.

The minimum description length approach estimated an average of 47 components. Sixty independent components were selected for greater accuracy. Twenty-nine components were identified as meaningful according to the selection criteria described in (Section “2.5. ICA analysis”), and the selected components were grouped into seven resting-state intrinsic networks. The selected components and networks are shown in Figure 1.

We performed two steps analysis to reveal static FC changes of PD in different stages. Static FC differences between the HC and PD-BL group were calculated with age, sex and education years as covariates, using two sample t-test. The difference between HC and PD-BL did not survive the FDR correction, which was used for multiple p-value correction. Differences between PD-BL and PD-FU were calculated using paired t-test, with FDR correction method to correct p-values. Results showed that network averaged connections between FP and SMN of PD-FU were lower than PD-BL (p < 0.05, t = −3.64) (Figure 2).

All dynamic FC matrices are clustered into four states according to the elbow method. 19% of the matrices clustered into state 1, characterizing positive coupling within the VIS network. State 2 contained 22% of the matrices, which characterized positive coupling within SMN, within VIS, between SMN and VIS, between DMN and SMN, and between DMN and VIS. State 3 had 56% of the matrices, showing hypo-coupling within and between all networks. State 4 was the smallest, with only 3% of matrices (Figure 3).

Temporal characteristics such as fractional windows, mean dwell time and transition number of HC group and PD subgroups were calculated. Age, sex, and education years may affect cognitive function. The effect of age, sex and education years were controlled when calculating difference between HC and PD-BL group. Results showed that there was no statistical difference between HC and PD-BL group of above temporal characteristics. Changes between PD-BL and PD-FU were performed sign rank test, because data wasn’t normality distributed. The temporal properties analysis showed that the fractional windows and mean dwell time of PD-FU at state 2 were statistically lower than PD-BL. Fractional windows and mean dwell time of PD-FU at state 3 were statistically higher than PD-BL (Figure 4).

Correlations between clinical scale scores and measurements of rs-fMRI were computed. Clinical scales included UPDRS-I, UPDRS-II(ON), UPDRS-III(ON), SCOPA-AUT, and ADL. Rs-fMRI measures included static FC of brain areas, fractional window, mean dwell time and transition number. Age, sex and education years may affect correlation result of abovementioned score. The effect of age, sex, and education years were controlled as variances, using partial correlation method. Results showed that SCOPA-AUT scores positively correlated with the mean dwell time of state 3 of PD-FU (Figure 5). There was no correlation between clinical scores and temporal characteristics of PD-BL.