The present study was approved by the Medical Research Ethics Committee of the Third Military Medical University (Chongqing, China), and written informed consent was obtained from all subjects. Fifty right-handed PD patients [25 AR patients (M:9/F:16, Age(y): 65.12 ± 10.11)] and 25 TD patients [M:9/F:17, Age(y): 60.28 ± 11.14] were recruited from the movement disorders outpatient clinic of Southwest Hospital. Each PD subject was diagnosed by an experienced movement disorder specialist according to the United Kingdom Parkinson’s Disease Society Brain Bank criteria (Hughes et al., 1992). Importantly, all PD patients were screened to exclude patients with visual hallucination symptoms.

In the current study, only those patients who could tolerate long periods of immobility in the MRI scanner were included. Subsequently, motor performance was quantitatively and accurately assessed in detail by using the motor section (part III) of the Unified Parkinson’s Disease Rating Scale (UPDRS) (Fahn et al., 1987) and used to classify PD subjects by subtype. Specifically, the PD subjects were classified as either TD or AR based on the ratio of the mean tremor score (TD score) to the mean akinetic-rigidity score (AR score) (Eggers et al., 2011; Lewis et al., 2011). Just like the previously published study (Zhang et al., 2015b), for each patient, we computed a mean tremor score for 9 tremor items (right and left arm tremor from patient history, and lips or chin tremor, tremor in all 4 limbs, and activity in both arms or postural tremor, on examination) and a mean score for 5 AR items (falling, freezing, and walking difficulty, from patient history, and gait and postural instability, on exam-ination). Briefly, patients were classified as having TD-PD when the ratio of the mean tremor score to the mean AR score was ≥1.5 and as having AR-PD when this ratio was ≤1 (Eggers et al., 2011; Lewis et al., 2011). In addition, a Mini Mental State Examination (MMSE) score > 28 and a Montreal Cognitive Assessment (MoCA) score > 26 were required for all subjects to ensure that no patients met the criteria for dementia (Dubois et al., 2007). The patient exclusion criteria included secondary parkinsonism, atypical parkinsonian disease, and pre-existing neurological or psychiatric disorders (including seizures, aphasia, neglect, substantial sensory disturbances, severe depression, or claustrophobia). Meanwhile, patients with severe tremor which might disturb the MRI procedure were excluded.

All patients were substituted with dopaminergic treatment but none were treated with antidepressants. Three of 25 TD patients were also receiving low-dose anticholinergics (benzhexol, 0.5 mg/3 times daily).

In addition, 26 healthy controls (HCs) [M:9/F:17, Age(y): 59.84 ± 7.01] were recruited from among local individuals who volunteered to participate in scientific studies. All of the control subjects had normal neurological examination results and underwent neuropsychological testing prior to the MRI scan. The exclusion criteria were the same as those applied to the patients.

Functional images were acquired using a 3.0 T Siemens Tim Trio whole-body MRI system (Siemens Medical Solutions, Erlangen, Germany). Subjects were instructed to stay awake and close their eyes, and to try not to think of anything. All subjects confirmed they did not fall asleep during the MRI procedures. Foam padding and earplugs were used to reduce head motion and scanner noise. Imaging data were collected transversely by using an echo-planar imaging (EPI) sequence with the following settings: TR = 2000 ms, TE = 30 ms, flip angle = 90°, FOV = 192 mm × 192 mm, in-plane matrix = 64 × 64, thickness = 3 mm, voxel size = 3.0 mm × 3.0 mm × 3.0 mm. For each subject, a total of 240 volumes were acquired with scan time of 480 s. To minimize the impact of dopaminergic medication, the MRI scans were performed during a relatively hypodopaminergic state (12 h after a last dose of dopaminergic treatment). Other medications (anticholinergics) were withdrawn for 12 h prior to scanning (ALI–MELKKILÄ et al., 1993; Khor and Hsu, 2007).

Data preprocessing was carried out using Neuroscience Information Toolbox (NIT¹). A series of conventional preprocessing steps were performed, which include (1) discarding first five volumes; (2) slice-timing correction; (3) head-motion correction; (4) normalization with an EPI template in the Montreal Neurological Institute (MNI) space; (5) linear detrending; (6) temporal filtering (band-pass 0.01–0.08 Hz); and (7) Regressing out six motion parameters, white matter and cerebrospinal fluid. Participants were excluded from further analyses if head motion exceeded 2.0 mm or 2.0°. We further evaluated the mean absolute displacement of each volume as compared with the previous volume (Zhang et al., 2015b). The average/standard deviation of frame-wise displacement (FD) for each group is present at Supplementary Table 1. All the three groups showed very slight head motion (mean FD < 0.1). In addition, there is no differences among the three groups (F = 1.302, p = 0.278, ANOVA). Considering that the censoring may lead to fewer degrees of freedom and introduce noisy (Power et al., 2015), we did not apply the volume censoring in the preprocessing (Supplementary Data).

According to the approach introduced by Tomasi and Volkow (2010) we used custom-written software in Neuroscience Information Toolbox (NIT¹) to evaluate the individual voxel-wise FCD maps. The global FCD of a given voxel was defined as the number of functional connections between that voxel and all other voxels. The local FCD was obtained using a “growing” algorithm. In this algorithm, given a voxel χ0, an additional voxel χj was added to the list of neighbors of χ0 if it was adjacent to a voxel that was linked to χ0 by a continuous path of functionally connected voxels and the correlation coefficient between χ0 and χj was >0.6. This calculation was repeated for all voxels that were adjacent to the neighbors of χ0 in an iterative manner until no new neighbors could be added to the list. The local FCD (intraregional) of χ0 was defined as the number of elements in the list of neighbors. Therefore, the global FCD reflected the whole brain FC, whereas the local FCD reflected local FC. The threshold of the correlation coefficient is a key parameter in FCD analysis. Based on prior knowledge, the threshold was set at 0.6 (Supplementary Data and Figure 1; Tomasi and Volkow, 2010, 2011). For each subject, a global FCD map and a local FCD map were computed. For standardization purposes, individual FCD maps were divided by the mean FCD value. Prior to the group-level analyses, we applied 6-mm FWHM spatial smoothing to the normalized FCD maps.

The brain sites identified by FCD group comparisons between the two subtypes were used as seeds in subsequent FC analyses. Specifically, a sphere ROI was generated for each of the brain sites, the center of which corresponded to the peak voxel (radius = 3 mm); the mean BOLD time series was then extracted and correlated to the time series of all of the voxels in the brain (Li et al., 2015). In this way, FC maps of each seed were produced; these maps were transformed by Fisher’s r-to-z transformation to improve normality.

One-way analyses of variance and Chi-squared tests were used to analyze demographic characteristics among the three groups, and two-sample t-tests were used to perform post hoc analysis. All tests were two-tailed, and p < 0.05 was considered statistically significant. These analyses were performed with SPSS software (version 20.0; SPSS, Chicago, IL, United States).

To determine the effect of group, a voxel-based one-way analysis of covariance (ANCOVA) was used with age, gender, and education level as covariates in both the FCD and FC analyses using REST software. The mean FCD and FC values were extracted for post hoc analyses using the least significant difference (LSD) t-test. Multiple comparisons correction was performed using a height threshold (p < 0.005) of individual voxel and a cluster size based on Monte Carlo simulations (Ledberg et al., 1998), which corresponds to cluster-level p < 0.05 by AlphaSim correction. The AFNI AlphaSim program² was used to correct for multiple comparisons.

We evaluated the clinical data for normality and then performed Spearman rank correlation analyses to investigate the relationships between each of the clinical variables (duration, age of onset, UPDRS-III, H-Y scale score, TD score, and AR score) and the altered functional brain properties separately (the normalized FCD values of the ROIs and the average z-score of the altered FC). ROC curves were computed using the mean global FCD differences between the TD patients and the AR patients. The discriminatory performance of each index was evaluated separately. Additionally, the indices were combined using a binary logistic regression model, and the performance of this model was also assessed. In each case, the area under the ROC curve (AUC) and the 95% confidence interval (CI) were calculated using non-parametric methods. Spearman rank correlation analyses and ROC analysis were performed with SPSS software.

There were no significant differences among the three groups in terms of gender, age, or educational level and no significant differences between the TD and AR patients in terms of the duration of disease, Hoehn and Yahr scale score or UPDRS-III score (Table 1).

One-way analysis of variance did not reveal any significant differences in head motion when measured by either the mean or maximum relative displacement [F(2,69), translational mean: F = 0.039, p = 0.961; translational max: F = 0.493, p = 0.613; rotational mean: F = 2.93, p = 0.061; rotational max: F = 1.012, p = 0.369]. In addition, there is no differences among three groups in term of FD [F = 1.302, p = 0.278, ANOVA, F(2,69)], indicating that the PD patients and controls exhibited similar head motion characteristics.

Using the averaged local and global FCD maps, we found that those regions with high local and global FCD (Hubs) were mainly distributed bilaterally, with maximal magnitude in the posterior cingulate, occipital and prefrontal cortices (Figure 1A and Table 2).

For comparison among three groups of global FCD, the brain regions that showed a significant main effect (one-way ANCOVA, p < 0.005, AlphaSim corrected) predominantly included prefrontal regions (superior/middle frontal gyrus), the temporal cortex (middle temporal gyrus), and the cerebellum (Inferior and posterior lobe) (Figure 1B and Table 2).

Post hoc analyses were applied to identify significant pairwise differences between three pairs of groups (post hoc two-sample t-test, p < 0.005, AlphaSim corrected). (1) TD vs. AR. The TD patients had significantly decreased global FCD in the left inferior frontal gyrus, right middle frontal gyrus and right superior frontal gyrus, and increased global FCD in the cerebellum anterior lobe relative to the AR patients (Figure 2 and Table 2). (2) TD vs. HC. Compared with HCs, the TD patients had significantly increased global FCD in the right cerebellum posterior lobe, left superior frontal gyrus and right fusiform gyrus (Figure 2 and Table 2). (3) AR vs. HC. Compared with the HCs, the AR patients exhibited significantly increased global FCD in the left posterior lobe of the cerebellum, the right middle temporal gyrus, the bilateral middle frontal gyrus, and the right inferior frontal gyrus (Figure 2 and Table 2).

No significant differences in local FCD were found among groups. Thus, the four regions that exhibited significant differences in global FCD between TD and AR (the left cerebellum anterior lobe, right superior frontal gyrus, right middle frontal gyrus, and left inferior frontal gyrus) were defined as seeds in the following FC analysis.

We correlated the clinical variables with the normalized FCD values of the four regions just mentioned (the left cerebellum anterior lobe, right superior frontal gyrus, right middle frontal gyrus, and left inferior frontal gyrus) for the TD and AR patients separately. In the PD group (TD and AR patients), a significant negative correlation was identified between H-Y scale score and the FCD value of the cerebellum anterior lobe (Figure 3A).

The results of the ROC analysis indicated that the discriminatory performance of each FCD index was far lower than the inclusion of the four FCD indices (Table 3) and the inclusion of the four FCD indices in a binary logistic regression model had the highest power to discriminate the AR patients from the TD patients, with an AUC of 0.96 (95% CI = 0.91–1.00; Figure 3B and Table 3). At a cutoff of 0.23, the sensitivity and specificity were 95.2 and 80.9%, respectively (Table 3).

For the seed of the left cerebellum anterior lobe, TD patients showed significantly increased FCs in the left ventral and dorsal post-central gyrus, left superior frontal gyrus, right insula and orbital inferior frontal gyrus compared with AR patients (Figure 4 and Table 4).

For the remaining seeds, no significant difference in FC was found between TD and AR patients (Table 4).

Some methodological limitations of this study should be mentioned. First, the small sample size and the relatively weak correction strategy employed (AlphaSim program) limited the conclusions that could be drawn; thus, a large number of participants should be included in future studies. Second, most of the PD patients were administered dopaminergic medications and a few with anticholinergics which might affect brain condition (Ehrt et al., 2010). Studies of drug-naïve individuals are warranted to exclude the effects of dopaminergic and anticholinergics medications on the findings. Longitudinal data are also required to further study the neuroimaging changes observed in the two subtypes. Third, the threshold (Tc) of the correlation coefficient used to define the functional connection was set at 0.6 in the FCD analysis based on prior knowledge, and this fixed value may lead to some false positive or negative findings. Flexible thresholds, such as a set of continuous thresholds, might be an alternative option for overcoming this limitation. We also evaluated the difference in FCD among groups using five thresholds (from 0.3 to 0.7 at intervals of 0.1), and similar findings were observed across several calculations using different thresholds (Supplementary Data and Figure 1). Lastly, considering other literature on PD pathophysiology, there are certainly other ways to subtype PD patients, i.e., according to different genetic backgrounds and/or clinical presentations. Our study, however, focused on TD versus AR subtype, one of the first clinical sub-types described in Jankovic et al. (1990). Future studies, therefore, are warranted to understand the relevance of resting state FC differences between PD subtypes other than TD versus AR.