The study was approved by the Medical Ethics Committee of the Affiliated Brain Hospital of Nanjing Medical University. Written informed consent was obtained from all participants after a full explanation of the whole procedure. Patients were recruited from the Movement Disorder Clinic at the Department of Neurology of Brain Hospital Affiliated with Nanjing Medical University. PD was diagnosed by an experienced neurologist according to the UK Parkinson’s Disease Society Brain Bank diagnostic criteria (Hughes et al., 1992). Patients were excluded if they had (1) anti-Parkinsonism medications prior to enrollment; (2) history of cerebrovascular disorders, head injury, seizure, hydrocephalus, intracranial mass, previous neurological surgery, and other neurologic diseases; (3) meeting the diagnostic criteria for PD with dementia; (4) history of psychiatric diseases other than depression; (5) antidepressant treatment or other psychiatric therapy; (6) other major systemic comorbidities; (7) history of alcohol or drug dependency or abuse; and (8) excessive head motion and poor MRI image quality (see below). All the patients were followed up for at least one year after enrollment to confirm the diagnosis according to the disease evolution and response to the dopaminergic therapy. Healthy controls (HCs) were matched to patients with PD for age, gender, and education. This study finally included 55 de novo patients with PD and 43 HCs. All the participants completed the clinical assessments and MRI scans during the baseline visit.

Patients with dPD were diagnosed using the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-V) criteria by an experienced psychiatrist. The severity of depression was quantified with the 17-item Hamilton Depression Rating Scale (HAMD-17). A HAMD-17 score higher than 14 was required for a more accurate diagnosis of dPD (Leentjens et al., 2000). Twenty-one patients with PD were diagnosed with dPD at last. In addition, motor symptoms and disease severity were evaluated by the motor section of the Unified Parkinson’s Disease Rating Scale (UPDRS-III) and Hoehn and Yahr (H&Y) staging scale. Cognitive function was assessed using the Mini-Mental State Examination (MMSE). All the assessments were conducted immediately before the MRI scan.

Magnetic resonance (MR) images were required using a 3T MRI scanner (Siemens, Verio, Germany). All participants laid supine with their heads fixed by foam pads with a standard birdcage head coil to minimize head movement. The participants were instructed to remain as still as possible and to close their eyes while remaining awake without thinking of anything. Axial anatomical images were acquired using a T1 fluid-attenuated inversion recovery sequence with the following parameters: repetition time (TR) = 2,530 ms; echo time (TE) = 3.34 ms; field of view (FOV) = 256 mm × 256 mm; matrix = 256 × 192; slice thickness/gap = 1.33/0.5 mm; flip angle (FA) = 7 degrees; bandwidth = 180 HZ/PX; 128 slices covered the whole brain, for image registration and functional localization. Functional images were subsequently collected in the same slice orientation with a gradient-recalled echo-planar imaging pulse sequence. A total of 240 volumes were obtained (TR = 2,000 ms; TE = 30 ms; FOV = 220 mm × 220 mm; matrix = 64 × 64; thickness/gap = 3.5/0.63 mm; FA = 90 degrees; bandwidth = 2,232 HZ/PX; slice numbers = 31).

Preprocessing of the functional MRI (fMRI) data was carried out using data processing assistant for resting-state fMRI (DPARSF)¹ (Chao-Gan and Yu-Feng, 2010) based on SPM12². The first 10 volumes of each rest functional section were removed for signal equilibrium and participant adaptation to the scanning environment, resulting in a total of 230 volumes. The remaining volumes were then corrected for slice timing using the middle slice as a reference, realigned for head motion correction, segmented into white matter, gray matter, and cerebrospinal fluid (CSF) using tissue probability maps derived from the T1 images, normalized into the standard Montreal Neurological Institute space using diffeomorphic anatomical registration through exponentiated Lie algebra (DARTEL), resliced as a 3 mm × 3 mm × 3 mm voxel size, and spatially smoothed with a 6 mm full-width at half-maximum (FWHM) Gaussian kernel. Aiming to further minimize the potential effects of head motion, we excluded patients with a mean framewise displacement (FD) > 0.5 mm or whose head motion exceeded a maximum translation of 3 mm or rotation of 3° from the further analysis (1 patient with dPD, 3 patients with ndPD, and 2 HCs), referring to a previous study (Navalpotro-Gomez et al., 2020). There were no significant group differences in head motion between dPD (mean FD = 0.08 ± 0.04 mm), ndPD (mean FD = 0.10 ± 0.06 mm), and HCs (mean FD = 0.11 ± 0.06 mm) groups (ANOVA, P = 0.12). Additionally, considering the possible deleterious influence of head motion, we added a scrubbing step in the preprocessing procedures and repeated the following analyses for further verification (Power et al., 2012; Yan et al., 2016). we found nearly consistent results with regard to temporal properties (Supplementary Figure 2).

The preprocessed resting-state data of all subjects were decomposed into functional networks by applying spatial group independent component analysis (ICA) implemented in the Group ICA of the fMRI Toolbox (GIFT v4.0b)³. Subject-specific data were reduced to 120 independent components (ICs) with the principal component analysis (PCA), and further decomposed into 100 ICs using the expectation-maximization algorithm at the group level as previously done (Kim et al., 2017). The Infomax ICA algorithm in ICASSO (Bell and Sejnowski, 1995) was repeated 20 times to ensure stability and validity (Himberg et al., 2004). Subject-specific spatial maps and time courses were obtained using the GICA back-reconstruction approach (Calhoun et al., 2001).

Out of the 100 independent components, 37 ICs were identified as meaningful according to the following criteria by Allen et al. (2014): (1) peak coordinates of spatial maps located primarily in the gray matter; (2) low spatial overlap with known vascular, motion, and susceptibility artifacts; (3) time courses dominated by low-frequency signals (ratio of the integral of spectral power < 0.10 Hz to.15–0.25 Hz); and (4) time courses characterized by a high dynamic range (a range difference between the minimum and maximum power frequencies). These 37 ICs were then sorted into six functional networks based on the spatial correlation values between ICs and the template and visual inspection (Shirer et al., 2012). The six functional networks were: basal ganglia (BG), auditory (AUD), visual (VIS), sensorimotor (SMN), cognitive executive (CEN), and default mode (DMN) (Figure 1A and Supplementary Table 1).

The following postprocessing steps were performed on the time courses of 37 ICs to remove remaining noise sources, including detrending linear, quadratic, and cubic trends, regressing out six realignment parameters and their temporal derivatives, despiking detected outliers by 3DDESPIKE, and low-pass filtering with a high cutoff frequency of.15 Hz (Allen et al., 2014). To obtain the static functional connectivity matrix, pair-wise Pearson’s correlations between ICs were computed and converted to z-scores via Fisher’s z-transformation using the post-processed time courses (Kim et al., 2017; Figure 1B).

Dynamic functional connectivity analysis was conducted using the GIFT toolbox through two steps: a sliding window approach and k-means clustering. First, as in previous studies (Kim et al., 2017; Fiorenzato et al., 2019; Chen et al., 2021), we applied a sliding time window of 22 TR (44 s) to segment the resting-state data, a setting that was proved to provide a good trade-off between the ability to resolve the dynamics of FC and the quality of the correlation matrix estimation (Allen et al., 2014), with a Gaussian window alpha value of 3 and a step of 1 TR, resulting in the analysis of 208 windows. To assess the robustness of the results regarding different window sizes, we repeated all the analyses for window sizes ranging from 20 to 28 TR (Supplementary Figure 3). For each window, a 37 × 37 pairwise covariance matrix was calculated by the regularized inverse covariance matrix (Varoquaux et al., 2010). Further, a penalty on the L1 norm of the precision matrix was imposed in the graphical LASSO framework with 100 repetitions to promote sparsity (Friedman et al., 2008). The resulting functional connectivity matrices were stabilized via Fisher’s z transformation and then residualized with age, gender, education, and MMSE scores. We finally obtained 208 functional connectivity matrices for each subject representing the dynamic changes of functional connectivity during the entire scan time for further analysis. Second, we applied the k-means clustering method to the windowed covariance matrices of all subjects to assess the reoccurring functional connectivity patterns (states), which were repeated 100 times (Allen et al., 2014). L1 (Manhattan) distance function was implemented considering its adaptation to high-dimensional data (Aggarwal et al., 2001). Based on the elbow criterion of the cluster validity index suggested by previous studies (Allen et al., 2014; Díez-Cirarda et al., 2018), the optimal number of clusters (k) was determined to be four (k = 4) (Figure 2B), and all matrices for each subject were then categorized into one of the four specific states. Thus, each subject can enter one to four of the defined states throughout the entire scanning (Figure 2A). Additionally, to assess the robustness of the results regarding the optimal value for k, we repeated the above analyses for k ranging from 2 to 5 (Supplementary Figure 4).

In order to explore the temporal properties of dynamic functional connectivity states, we calculated three different indices for each subject, including: (1) fraction time, defined as the percentage of total time one subject spent in a given state; (2) mean dwell time, defined as the time a subject remained in a certain state; and (3) number of transitions, defined as the number of translations between states. The differences among groups (HCs, ndPD, and dPD) were investigated using a one-way ANOVA test, followed by post hoc two-sample t-tests. P-values were corrected by false discovery rate (FDR) for multiple comparisons. The threshold for statistical significance was set at P < 0.05.

Additionally, the subject-specific connectivity pattern of each state was represented by the median value of all functional connectivity matrices assigned to that state. Group differences among patients with HCs, ndPD, and dPD in dynamic functional connectivity pairs (666 pairings) within each state were calculated using one-way ANOVA, followed by post hoc two-sample t-tests (P < 0.05, FDR-corrected).

We applied a graph theory approach to analyze the dynamics of graph characteristics of the network using GRETNA software⁴. Thirty-seven ICs were defined as nodes and the functional connectivity between each two ICs was constructed as edges for each time window. The resulting 208 functional connectivity matrices per subject were binarized with respect to a fixed sparsity threshold (the existing number of edges divided by the maximum possible number of edges), set as 0.1 to 0.34 in 0.01 increments based on prior studies (Kim et al., 2017; Navalpotro-Gomez et al., 2020). Only positive correlations were considered in our analysis. We calculated the following parameters by integrating over all the defined threshold ranges within each time window, including small-worldness, global efficiency, betweenness centrality, clustering coefficient, and local efficiency of the network. Then the variance of these parameters over the time windows was computed to evaluate the dynamic graph properties for each subject. Two-sample t-tests (or Mann–Whitney U test) and ANOVA test (or Kruskal–Wallis test) were respectively conducted to assess whether there were differences in topological variability between PD and HCs groups or among the HCs, ndPD, and dPD groups (P < 0.05).

Statistical analyses were calculated using SPSS Statistic 24.0 (Chicago, IL, United States). A two-sample t-test or Mann–Whitney U test was applied to compare two groups. One-way ANOVA or Kruskal–Wallis test was performed to compare the HCs, ndPD, and dPD groups. A chi-squared test was used to compare categorical variables such as gender. Spearman’s correlation analyses were conducted to assess the correlations between the detected temporal properties and HAMD-17 scores. P-values were corrected for multiple comparisons and P < 0.05 was set as a threshold for statistical significance.

A total of 34 patients with ndPD, 21 patients with dPD, and 43 HCs were eventually included in the further analysis. There were no significant differences among the three groups in terms of age, gender, education, and MMSE scores (P = 0.963, P = 0.780, P = 0.150, P = 0.056). No significant differences between patients with dPD and ndPD were found in disease duration, H&Y stage, UPDRS-III scores (P = 0.823, P = 0.199, P = 0.097). In addition, the HAMD-17 score of the patients with dPD was significantly higher than that of patients with ndPD or HC (P < 0.001) (Table 1).

Spatial maps of all 37 independent components defined by using group ICA are illustrated in Figure 1A. These 37 ICs made up the following six networks: BG (IC 19), AUD (IC 55, 91), VIS (IC 10, 23, 34, 44, 45, 58, 64, 72, 82), SMN (IC 1, 3, 5, 6, 8, 21, 36), CEN (IC 20, 22, 33, 57, 61, 63, 69, 86, 89, 93), and DMN (IC 7,12, 24, 31, 32, 41, 51, 67). The specific distribution of the 37 ICs is showed in Supplementary Table 1.

Through the k-means clustering analysis, a total of four highly structured functional connectivity states were identified by the elbow criterion. We displayed the 5% of the functional connectivity network with the strongest positive or negative connections to observe the different patterns of each state (Figure 3B). Combined with the visualized cluster centroids (Figure 3A), state 1, which was the most frequent one (47%), was characterized by sparse connections both within and between networks, with relatively weak connectivity mainly located within DMN or VIS. State 2 (25%) was characterized by DMN-dominated and CEN-disconnected patterns, in which strong correlations were mainly exhibited within the DMN and between DMN and other networks (CEN/SMN/AUD), while the CEN showed sparse intra-network connections and inter-network connections with SMN/VIS/AUD. State 3 (17%) was characterized by relatively stronger positive connections within and between each state. State 4 (11%) was characterized by a regionally densely connected pattern, which showed positive correlations located mainly within VIS/SMN and between these two networks.

One-way ANOVA indicated that there were significant differences between groups in terms of fraction time and mean dwell time in state 2 (ANOVA, fraction time: P = 0.004; mean dwell time: P = 0.048, FDR-corrected). Post hoc t-tests demonstrated that patients with dPD had a significantly higher fraction time in state 2 than that of patients with ndPD (P = 0.030 FDR-corrected) and healthy controls (P = 0.009, FDR-corrected). Additionally, patients with dPD spent significantly more dwell time in state 2. This is consistent with fraction time when compared to patients with ndPD (P = 0.048, FDR-corrected) and HC (P = 0.015, FDR-corrected). We found that fraction time (P = 0.006, r = 0.277) and mean dwell time (P = 0.035, r = 0.215) in state 2 were slightly correlated with HAMD scores in all participants when regressing out a group as a covariate. No differences were found in the number of state transitions or temporal properties in other states between the three groups (ANOVA, P > 0.05) (Figure 4 and Supplementary Table 2).

Validation analyses showed that the main results of temporal properties in state 2 were consistent for different window sizes (Supplementary Figure 3). Besides, with regard to different numbers of clusters, the main results also persisted, but there was a redundant state in k = 5, and the clustered states were incomplete for k less than 4 (Supplementary Figure 4).

We subsequently examined between-group differences in the strength of connections for each state. ANOVA results (P < 0.001, FDR-corrected) showed a significant overall difference between the HC, patients with ndPD, and patients with dPD. Post hoc t-tests indicated that patients with dPD and patients with ndPD both had significantly decreased connection within the SMN (IC1/IC3, left postcentral cortex, and right precentral cortex) in state 4 compared to healthy controls (P < 0.001, FDR-corrected). In contrast, no differences were found between the ndPD and dPD groups (P = 0.58) (Figure 2C).

When comparing the variance of global efficiency, neither HCs group and PD group (Mann–Whitney U test, P = 0.755) nor HCs group and PD sub-groups (Kruskal-Wallis test, P = 0.552) differed significantly in this aspect. No significant differences were found between the HCs and PD groups (two-sample t-test, P = 0.501) or between the dPD, ndPD, and HCs groups (ANOVA, P = 0.286) with respect to betweenness centrality. In terms of the local efficiency of the network, there were significant differences between the dPD, ndPD, and HCs groups (ANOVA, P = 0.049). Post hoc t-tests revealed that the dPD group exhibited smaller variance than the HCs group did (P = 0.010), but no significant differences were found between the dPD and ndPD groups (P = 0.143) or between the ndPD and HCs groups (P = 0.319). There was a trend for less variability in the PD group than in the HCs group (two-sample t-test, P = 0.075). Further, the analysis of small-worldness variability demonstrated that the PD group was less variable than the HC group (two-sample t-test, P = 0.034), and there may be a trend of differences between the dPD, ndPD, and HC groups (ANOVA, P = 0.097), consistent with the results of clustering coefficient (two-sample t-test, P = 0.041; ANOVA, P = 0.076) (Figure 5 and Supplementary Table 2).