Introduction

Parkinson’s disease (PD) is a common progressive neurodegenerative disorder, characterized by significant clinical heterogeneity and prognosis1,2, which makes investigating the clinical and biological subtypes a crucial breakthrough for exploring its underlying mechanisms. Previous studies have demonstrated that patients with dominant postural instability and gait difficulty or with significant non-motor symptoms would have a poor quality of life and suffer a fast disease progression3,4,5,6. However, it is not possible to accurately identify the subtype based on patient symptoms alone until the disease progresses to an advanced stage and presents with significant clinical symptoms. The brain change underlying clinical symptoms is expected to provide a more characteristic and objective characterization of the disease7. Unfortunately, the heterogeneity of brain degeneration patterns in the early stages of PD is also largely unknown.

The hallmark pathophysiological change in PD is the widespread aggregation of α-synuclein pathology throughout the brain1, which leads to the progressive and heterogeneous brain alterations with different aspects of clinical relevance8. In the recent years, magnetic resonance imaging (MRI) has proven to be a valuable non-invasive technique for in vivo quantifying brain tissues, thereby providing insights into disclosing the clinical presentations and disease prognoses of PD9. Unsupervised clustering enables the discovery of latent patterns underlying the systematical brain information characterized by multimodal MRI, which would enhance the identification of PD subtypes that share common degeneration mechanisms in a data-driven way3,10. Based on these techniques and algorithms, previous studies have demonstrated the spatial and temporal brain degeneration of PD, successfully reflecting their pathophysiological information and revealing significant heterogeneity. Uribe et al. discovered distinct patterns of cortical atrophy associated with cognitive performance11, they also focused on de novo PD patients12, revealing longitudinal atrophy differences across these distinct cortical degeneration patterns over a 4-year follow-up13; Inguanzo et al. examined both gray and white matter, providing a more comprehensive view of brain degeneration patterns in PD14; Guo et al. integrated clinical symptoms with fiber connectivity, revealing distinct patterns of brain damage associated with the depression-dominant and motor-dominant subtypes15. All these studies strongly supported the existence of differentiated brain degeneration in PD. However, the systemic characterization of brain degeneration patterns in early-stage PD remains insufficient, with few studies incorporating degeneration information from the core nigrostriatal regions. The disruption of iron homeostasis in these key regions could exacerbate oxidative reactions between dopamine and iron, thereby playing a crucial role in the pathophysiological mechanisms of PD16, which could be non-invasively measured by Quantitative Susceptibility Mapping (QSM)17. While previous studies have revealed the iron heterogeneity in the early disease18, data-driven clustering research in this field remains limited. Incorporating measurements of iron could enrich the quantification of brain tissue, offering a broader perspective on brain degeneration in PD. Therefore, we hypothesized that, the brain degeneration has already differentiated when PD patients are at early-stage and might exhibit different clinical profiles and the longitudinal progression.

This study aimed to identify the existence of different brain degeneration patterns among early-stage PD patients by employing multimodal MRI and investigate its clinical relevance according to the following steps: (1) utilizing an unsupervised clustering method to classify early-stage PD patients into subtypes that were sharing different degeneration of brain morphometry, white matter (WM) microstructure and iron metabolism; (2) demonstrating the differential clinical profiles and their longitudinal progression between the early-stage PD patients with different brain degeneration patterns.

Results

Demographic information of PD and HC

The demographic information of the discovery dataset was shown in Table 1. There was no significant difference in age (p = 0.930) and sex (p = 0.362) between PD patients and healthy controls (HC). However, HC had received a longer duration of education compared to PD patients (p < 0.001).

Table 1 Demographics of healthy controls and Parkinson’s disease patients in discovery dataset
Full size table

Two PD subtypes identified by the multimodal MRI parameters

We ranked the independent principal component (PC) generated during the reduction of principal component analysis (PCA) dimensionality based on decreasing variance explainability. Subsequently, we selected the initial 27 PCs with variance explainability >1%, retaining 70.57% of the effective information. Details of the loadings of each imaging feature for the first PC and the cumulative contributions for the first 27 PCs can be found in Supplementary Fig. 1. The hierarchical clustering analysis applied to the 27 PC scores identified two PD subtypes as the optimal solution. The overall changes in their imaging variables of 126 brain regions were illustrated in Fig. 1.

Fig. 1: The imaging differences between two subtypes and HC.
figure 1

All significant p values were corrected using FDR correction. Using warm and cool color tones to indicate the relative magnitude of imaging features between groups. a Comparison of cortical thickness; b Comparison of Fractional Anisotropy; c Comparison of magnetic susceptibility. HC Healthy Control. CN caudate nucleus, GP globus pallidus, SN substantia nigra, RN red nucleus. All comparisons of imaging features above were adjusted for age, sex, and education, and the comparisons of cortical thickness were adjusted for eTIV and mean cortical thickness additionally. The results were corrected for multiple comparisons using FDR correction.

Full size image

Specifically, compared to HC, Subtype1 patients (n = 72) exhibited a neurodegeneration pattern characterized by reduced cortical thickness in frontal and temporal regions, including the right caudal anterior cingulate (p = 0.005) and left lingual (p = 0.022) gyrus, and increased cortical thickness in parietal regions, including the left inferior parietal (p < 0.001), bilateral superior parietal (p < 0.001), and right supramarginal (p = 0.001) gyrus. Additionally, these patients manifested as extensive decreases in fractional anisotropy (FA) values compared to HC, encompassing a wide range of 36 regions including 17 association fibers (p-range: <0.001 to 0.019), 3 commissural fibers (p-range: <0.001 to 0.014), 14 projection fibers (p-range: <0.001 to 0.002), and 2 brainstem fibers (p-range: <0.001 to 0.028), while increased FA value in the middle cerebellar peduncle (p < 0.001) was observed. In the consideration of regional iron deposition, Subtype1 patients exhibited significantly increased iron deposition in the left caudate nucleus (CN) (p = 0.030), bilateral putamen (p < 0.001), bilateral substantia nigra (SN) (left: p = 0.024; right: p < 0.001), and bilateral red nucleus (RN) (p < 0.001) compared to HC.

Subtype2 patients (n = 107) exhibited a different pattern characterized by reduced cortical thickness in the left rostral middle frontal gyrus (p = 0.002), and increased cortical thickness in parieto-occipital regions, including the left inferior parietal (p = 0.002), left lateral occipital (p = 0.020), and left supramarginal (p = 0.010) gyrus compared to HC. Additionally, there were 18 WM skeletons with increased FA values compared to HC, predominantly involving 9 projection fibers (p-range: <0.001 to 0.049), as well as some association fibers (p-range: <0.001 to 0.013) and brainstem fibers (p-range: <0.001 to 0.037). No significant differences were observed in iron deposition within the nucleus compared to HC (p-range: 0.109 to 1).

Compared to the Subtype2 patients, patients in Subtype1 showed cortical thinning in the right temporal pole (p = 0.003) and cortical thickening in the right superior parietal cortex (p = 0.036). Additionally, Subtype1 patients exhibited significant reductions in FA values in 44 fibers (p-range: <0.001 to 0.010). Regarding regional iron deposition, Subtype1 patients exhibited significantly increased iron deposition in the left CN (p = 0.005), bilateral putamen (left: p = 0.006; right: p = 0.008), and RN (left: p = 0.016; right: p = 0.004) than Subtype2. Detailed results were provided in Supplementary Table 1.

Demographic and clinical characteristics of PD subtypes

Significant differences were observed in age (p < 0.001), with patients of Subtype1 (63.86 ± 7.44 years) being older than those of Subtype2 (51.99 ± 9.59 years); sex (p < 0.001), with Subtype2 (68 males, 63.55%) contained more males than Subtype1 (32 males, 44.44%); and education (p < 0.001), with Subtype1 (6.00 (6.80) years) having a lower level of education compared to Subtype2 (9.00 (5.00) years). No significant difference in disease duration (p = 0.365) and Hoehn Yahr (HY) stage (p = 0.535) was observed between the two subtypes.

Patients of Subtype1 exhibited more severe clinical manifestations, encompassing both motor and non-motor aspects. Specifically, Subtype1 patients had higher scores in the following assessments and tests: Unified Parkinson’s Disease Rating Scale (UPDRS)-II scores (p = 0.047), UPDRS-III scores (p = 0.047), and bradykinesia scores in UPDRS-III (p = 0.047) for motor symptom and Geriatric Depression Scale (GDS) scores (p = 0.047) for emotional function; and lower scores including Mini-Mental State Examination (MMSE) scores (p = 0.047) and number of correct responses on the Symbol Digit Modalities Test (SDMT) (p = 0.023) for cognitive ability. No significant differences were observed in tremor (p = 0.412) or rigidity (p = 0.075) scores within UPDRS-III, performance in the Auditory Verbal Learning Test (AVLT) (p = 0.309), Semantic Fluency Test (SFT) (p = 0.247), Digit Span Test (DST) (p = 0.095), rapid eye movement (REM) sleep behavior disorder questionnaire-Hong Kong (RBDQ-HK) (p = 0.847), Parkinson’s Disease Questionnaire-39 items (PDQ-39) (p = 0.110), and Scale for Outcomes in Parkinson’s Disease for Autonomic Symptoms (SCOPA-AUT) scores (p = 0. 075). All details were summarized in Table 2 and Fig. 2.

Fig. 2: The differences of clinical features between two PD subtypes at baseline.
figure 2

The differential clinical features between the two subtypes are represented using a raincloud plot. In this plot, box plots with center line indicating median, bounds of boxes showing upper and lower quartile, whiskers illustrating 1.5*interquartile range, and the position of the dots and the width of the violin representing the distribution of raw data. All comparisons above were adjusted for age, sex, and education. The p values were FDR corrected. UPDRS Unified Parkinson’s disease rating scale, MMSE mini-mental state examination, SDMT symbol digit modalities test, GDS geriatric depression scale. *p < 0.05.

Full size image
Table 2 Demographics and clinical characteristics of Parkinson’s disease patients in two subtypes
Full size table

Comparison of clinical progression between PD subtypes

Forty-seven patients were followed up at least once (median: 1.00; interquartile range: 1.00 to 2.00), with a median duration of 1.58 (interquartile range: 1.17 to 2.33) years, enabling us to assess the discrepancy in clinical progression between the two subtypes. The comparisons of clinical characteristics and imaging features for these patients at baseline were presented in Supplementary Tables 2 and 3. Through linear mixed model (LMM), we found distinct disease progression trajectories for the two subtypes. Specifically, significant progression was exclusively observed in Subtype2 for the following clinical and neuropsychological assessments: UPDRS-II scores (Time: β = 1.061, p = 0.007), UPDRS-III scores (Time: β = 3.365, p < 0.001), rigidity scores in UPDRS-III (Time: β = 1.003, p = 0.003), bradykinesia scores in UPDRS-III (Time: β = 1.538, p = 0.003), and SCOPA-AUT scores (Time: β = 0.959, p = 0.020). Additionally, a comparison between the two subtypes revealed that Subtype2 exhibited a trend of a higher rate of change in bradykinesia scores in UPDRS-III (Subtype × Time: β = 1.572, uncorrected p = 0.042) compared to Subtype1. Conversely, Subtype1 exhibited a trend of a faster decline in AVLT long-delayed recall compared to Subtype2 (Subtype × Time: β = 1.363, uncorrected p = 0.030). The difference in progression between the two subtypes regarding clinical and neuropsychological assessments were demonstrated in Table 3 and Fig. 3. After additionally considering the influence of baseline conditions, significant progression was still observed in Subtype2 for UPDRS-III scores (Time: β = 3.358, p = 0.023), rigidity scores in UPDRS-III (Time: β = 0.961, p = 0.023), bradykinesia scores in UPDRS-III (Time: β = 1.424, p = 0.039), and SCOPA-AUT scores (Time: β = 1.137, p = 0.023) was still observed in Subtype2. Detailed information was provided in Supplementary Table 4.

Fig. 3: The differences of clinical progression between two PD subtypes.
figure 3

The red/blue dots and connected lines indicate the changes of variables during follow-up at the subject level for Subtype1/Subtype2. The corresponding bold lines represent the changes at the group level for each subtype. All comparisons were adjusted for age at baseline, sex, and education. The p values were FDR corrected. Significant p values are shown in bold. β1, p1: The time effect of variable in Subtype1; β2, p2: The time effect of variable in Subtype2; β3, p3: The interaction effect of subtype and time. UPDRS unified Parkinson’s disease rating scale, AVLT auditory verbal learning test, SCOPA-AUT scale for outcomes in PD for autonomic symptoms.

Full size image
Table 3 Linear mixed model analysis results in two subtypes
Full size table

Reproducibility of PD subtypes

We assessed the reproducibility of the subtypes in an independent validation dataset including 82 patients (age: 59.68 ± 8.59 years; sex: 42 males, 51.22%). Employing the same methodology, each patient was categorized into a subtype, enabling the calculation of the average z-score of each imaging feature across individuals within the same subtype, serving as a representation of subtype-specific degeneration signature. As a result, we identified two PD subtypes in the validation dataset and observed a high consistency of z-scores between the corresponding subtypes in both the discovery and validation datasets (r = 0.830, p < 0.001 in Subtype1; r = 0.739, p < 0.001 in Subtype2) (Supplementary Fig. 2).

Correlation analysis

Twenty-six imaging features were chosen to examine their correlation with clinical features across the entire cohort, Subtype1, and Subtype2. Additional details on the imaging feature selection process were provided in Supplementary Table 5.

In Subtype1, lower FA values in the right cingulum (hippocampus) (r = −0.389, p < 0.001) and the fornix (cres)/stria terminalis (r = −0.415, p < 0.001) were significantly correlated with higher UPDRS-III scores at baseline (Fig. 4a). However, within each subtype, we did not observe any correlation between imaging features and clinical progression rates. But in the entire cohort, higher FA values in the genu of corpus callosum (r = 0.514, p < 0.001) and the left posterior thalamic radiation (r = 0.470, p < 0.001) were significantly correlated with a greater increase in bradykinesia scores of UPDRS-III. Additionally, lower FA values in the left corticospinal tract (r = 0.554, p < 0.001) was significantly correlated with a faster decline in AVLT long-delayed recall scores (Fig. 4b). Non-significant results of the correlation analysis were presented in Supplementary Fig. 3.

Fig. 4: Correlation between imaging and clinical features.
figure 4

a Associations between imaging features and baseline clinical characteristics in Subtype1; b Associations between imaging features and the rates of clinical symptom progression in all PD patients. The heatmap provides an overview of the correlation between the 26 selected imaging features and clinical indicators that show differences at baseline or follow-up. All features included were adjusted for age, sex, and education, the cortical thickness was further adjusted for eTIV and mean cortical thickness. The scatter plots display results with p < 0.001. FA fractional anisotropy, UPDRS unified Parkinson’s disease rating scale, MMSE mini-mental state examination, SDMT symbol digit modalities test, GDS geriatric depression scale, AVLT auditory verbal learning test, SCOPA-AUT scale for outcomes in PD for autonomic symptoms.

Full size image

Lateralization analysis

An analysis of brain changes associated with the laterality of motor symptoms was presented below. Among the 179 PD patients, 166 showed lateralization of motor symptoms, with 71 exhibiting higher UPDRS-III scores on the left side and 95 on the right side. Among the 126 brain imaging features, there were 60 paired comparisons. None of the brain imaging features showed a significant difference associated with the laterality of motor symptoms (p-range: 0.080 to 0.976). More details were provided in Supplementary Table 6.

Alternative analysis using mean diffusivity (MD)

To further validate the reliability of the clustering results using FA as the primary WM metric, we conducted clustering analysis using MD as an alternative. The clustering results still identified two subtypes (Subtype1: n = 59; Subtype2: n = 120), exhibiting similar brain imaging features, demographic characteristics, and clinical assessments as those observed with FA. Detailed results were provided in Supplementary Tables 7 and 8.

Analysis of the suboptimal clustering solution

To explore potential insights from other clustering solutions, we conducted an analysis of the suboptimal clustering solution. The suboptimal solution identified by the ‘NbClust’ package classified early-stage PD patients into three subtypes (Supplementary Fig. 4). Compared to the two-cluster solution, the three-cluster solution further subdivided Subtype1 (n = 72) into Subtype1 (n = 37) and Subtype3 (n = 35), while Subtype2 remained unchanged (n = 107). To simplify the presentation of the new clustering results, we compared brain imaging features between Subtype1 and Subtype3 after incorporating the HC. Compared to HC, Subtype3 patients exhibited more widespread (38 WM fibers) reductions in FA values (p-range: <0.001 to 0.026) and normal subcortical iron deposition, while Subtype1 patients showed higher iron deposition in the bilateral putamen (left: p < 0.001; right: p = 0.002), left GP (p = 0.012), bilateral SN (p < 0.001), and bilateral RN (p < 0.001), with relatively mild (8 WM fibers) or non-significant WM degeneration. These findings revealed results that were additional to those observed in the two-cluster solution, which may suggest the existence of distinct degeneration sequences, especially in the original Subtype1. More details were provided in Supplementary Table 9.

Discussion

This unsupervised data-driven study revealed different neurodegeneration patterns in early-stage PD patients utilizing multimodal MRI information, and identified their clinical symptoms and progression trajectories. Our study’s findings and inferences were summarized in Fig. 5. Based on the clustering scheme, the primary brain degeneration patterns were labeled as “Subtype1: early-deterioration” and “Subtype2: early-compensatory”. The early-deterioration subtype exhibited atrophy of frontotemporal cortex, thickening of parietal lobes, widespread reduction in FA values, and increased subcortical iron content. They manifested severe motor and non-motor symptoms at baseline and demonstrated a trend of faster memory decline longitudinally. In this subtype, lower FA values of the right hippocampus part of cingulum and fornix were associated with more severe motor symptoms at baseline. The early-compensatory subtype showed atrophy of the left rostral middle frontal gyrus, thickening of the parieto-occipital cortex, extensively increased FA values, and normal subcortical iron content. This subtype presented with mild symptoms at baseline but exhibited accelerated progression in motor and autonomic nervous system symptoms. Across the entire cohort, the higher FA values of the genu of corpus callosum and the left posterior thalamic radiation were correlated with the faster bradykinesia progression, while the lower FA values of the corticospinal tract were correlated with the faster memory decline

Fig. 5: Synthesis of information for the two subtypes.
figure 5

Subtype1 (early-deterioration) was characterized by an older age of onset and significant brain degeneration at baseline (green box), along with more severe clinical symptoms and a trend of faster memory decline during follow-up (red lines and dots). In contrast, Subtype2 (early-compensatory) was characterized by a younger age of onset and widespread compensatory brain changes at baseline (yellow box), with milder clinical symptoms initially, but accelerated progression of both motor and non-motor symptoms during follow-up (blue lines and dots).

Full size image

The degeneration pattern of early-deterioration subtype was primarily characterized by frontotemporal atrophy, widespread reduction in FA values across WM tracts and elevated subcortical iron content in the left CN, bilateral putamen, SN, and RN. The morphometry atrophy can result from the neuron loss and dendrite degeneration19. The FA reduction indicates the microstructural damage of WM, like axonal atrophy and myelin degradation20, while the excessive iron deposition can increase regional oxidative stress and injury, and facilitate a-syn aggregation17. Here, we reported a number of impaired hub regions that may closely associate with PD. Specifically, we identified widespread subcortical iron accumulation in bilateral SN, RN, and putamen, which was the phenomenon previously observed in the PD patients at late-stage21,22, indicating a high iron-related oxidative stress in these patients. Similarly, we detected widespread structural damage in gray and white matter regions responsible for motor control and spatial navigation, including the pontine crossing tract, superior longitudinal fasciculus, posterior corona radiata, and lingual gyrus, as well as in areas associated with cognitive function and emotional regulation, such as the cingulum, fornix, uncinate fasciculus, and the caudal anterior cingulate gyrus23,24,25. In contrast, parietal thickening in these patients was observed. From a functional perspective, the parietal cortex is known to play a critical role in integrating sensory and motor functions26. We speculated this alteration may be related to the dendrite regeneration, enhanced synaptic plasticity, and morphometry reorganization, possibly disclosing a complementarily enhanced brain function from the higher cortical areas27. Therefore, these findings offered a new insight into the existence of pathological heterogeneity in early-stage PD, thus a malignant brain deterioration pattern was suggested in the early disease3. Consistently, underlying the poor intrinsic brain architecture, these patients were observed with worse motor, cognitive, and emotional performance at baseline, along with a trend of faster memory decline during follow-up. It is worth noting that, these patients were older, with a similar disease duration compared to those of another subtype, indicating that they had a significantly later age of disease onset. As well established, general brain organization would get fragile when the brain chronological age gets old, with a reduced capacity for adaptation and increased susceptibility when facing pathological conditions28,29; this may explain one of the reasons of the occurrence of malignant brain degeneration pattern in these early PD patients. A recent study identified a subtype characterized by predominant degeneration of neocortical regions, which was associated with older age of onset, worse cognitive performance, and emotional disturbances30. We identified an early-deterioration subtype characterized by widespread WM damage connecting these neocortical regions, accompanied with similar clinical features. This further supports the empirically proposed subtypes based on age of onset31. Additionally, males exhibited a higher incidence of PD in our cohort, but a higher proportion of female patients was observed in Subtype1. This finding is consistent with demographic studies which reported a delayed age of onset of PD in female patients32. In summary, these findings supported the existence of an early-deterioration subtype characterized by an older age of onset and an early malignant brain deterioration pattern, accompanied by overt clinical symptoms.

The cortical alterations of early-compensatory subtype patients were generally similar to those of another subtype, characterized by frontal lobe atrophy and thickening of the parieto-occipital regions. We speculated that cortical morphological change is mild and display analogous trends in the early-stage. Here, we reported a differential alteration pattern in WM and subcortical nuclei. Specifically, increased FA values across multiple tracts including corticospinal tract, superior cerebellar peduncle, and internal capsule were observed. Wen et al.33 reported increased FA values in multiple WM tracts in the early disease, with these changes diminishing as the disease progresses, which further supported our observations. This may indicate the myelin remodeling and axonal regeneration could serve as a compensatory reorganization of neural circuits to counteract dopaminergic deficits, particularly in the early stages of the disease34. Abnormal regional iron deposition was not detected, suggesting high iron-related oxidative stress was not the case for these patients. Considering the earlier age of onset, this mild pattern of brain alterations may be attributable to an initial neural adaptation driven by the better neuroplasticity35,36. These findings were highly aligned with the mild clinical symptoms at baseline; however, after a median of 1.75 years’ follow-up, a rapid progression of motor symptoms occurred, along with worse autonomic nervous system symptoms; thus a nonlinear progression with later acceleration was suggested. Consistently, some researchers had proposed a nonlinear progression pattern characterized by an initial period of symptom stability, followed by an accelerated progression in the later-stage37,38,39. Besides, higher FA values in the genu of the corpus callosum and the left posterior thalamic radiation were associated with faster progression of bradykinesia in the entire cohort. This may suggest that the enhanced integrity of WM connecting the bilateral prefrontal motor cortices and those responsible for integrating sensory information could be a potential compensatory mechanism in the early motor dysfunction40,41. Even though our findings had explained the initial symptom stability with mild brain degeneration in these patients, the source of the later acceleration of symptoms are largely unknown. Future researches are warranted to clarify this important clinical phenomenon. In summary, here we also suggested an early-compensatory subtype that were with an earlier age of onset and mild brain degeneration at baseline, but with significantly accelerated clinical progression during follow-up.

This study had several limitations. First, this study merely included single-center data and a relatively small sample size, with only a subset of patients having a median follow-up of 1.58 (interquartile range: 1.17 to 2.33) years. Since PD is a progressive neurodegenerative disease with high heterogeneity as this study preliminarily disclosed, future studies with larger longitudinal sample size and longer follow-up duration are warranted to further validate and clarify the disease trajectory. Second, even though PD patients discontinued anti-parkinsonian medications for more than 12 h before data acquisition, the underlying effects of drugs cannot be completely eliminated especially for patients at re-visit.

In conclusion, two distinct early-stage brain degeneration patterns of PD were identified using unsupervised clustering methods based on multimodal MRI data. These subtypes exhibited differential clinical symptoms and progression characteristics, supporting the existence of high disease heterogeneity in early PD population.

Methods

Participants

This research was approved by the Medical Ethics Committee of the Second Affiliated Hospital of Zhejiang University School of Medicine and informed consent forms were obtained from all participants. We prospectively recruited 376 subjects: the discovery dataset included 179 early-stage PD patients (disease duration ≤5 years)42 in the study from our institution from February 2019 to November 2023, and 115 HC from the social community. In addition, a validation dataset consisting of 82 PD patients was enrolled from August 2014 to January 2019 at the same institute. The diagnosis of PD was made by experienced neurologists according to UK Parkinson’s Disease Society Brain Bank criteria43 before 2015 and Movement Disorder Society diagnostic criteria44 thereafter. Exclusion criteria for all subjects were as follows: (i) a history of stroke; (ii) a history of head injury; (iii) severe WM hyperintensity or massive cerebral infarction; (iv) obvious cerebral atrophy; (v) Left-handed or double-handed; (vi) severe metal dentures; (vii) motion artifact of imaging data; (viii) incomplete T1, diffusion tensor imaging (DTI) or enhanced susceptibility weighted angiography (ESWAN) data. Among these included patients, 47 of them had attended clinical follow-up visits, with a median follow-up duration of 1.58 years. The study flowchart was shown in Fig. 6.

Fig. 6: Study flowchart of data processing and analysis.
figure 6

(Step1) A total of 376 individuals (294 in discovery dataset and 82 in validation dataset) were enrolled in this study. Imaging data from three modalities and both motor and non-motor clinical features at baseline and follow-up were collected; (Step2) All imaging features were first corrected for covariates using linear regression, then subjected to dimensionality reduction through principal component analysis for hierarchical clustering analysis, thereby obtaining imaging subtypes of PD; (Step3) The differences of clinical/imaging variables, and longitudinal progression were assessed between subtypes; the correlations between imaging features and the clinical severity were analyzed. MRI magnetic resonance imaging, QSM quantitative susceptibility mapping, REM rapid eye movement, FA fractional anisotropy, PCA principal component analysis, UPDRS unified Parkinson’s disease rating scale, MMSE mini-mental state examination.

Full size image

Clinical and neuropsychological assessments

A battery of clinical questionnaire assessments was obtained from 179 PD patients during the “OFF state” (a period at least 12 h after withholding PD medications). Each patient was assessed using the UPDRS. The subscores of tremor, rigidity, and bradykinesia symptoms were calculated: subscore for tremor was obtained by adding the UPDRS-III items 20–21; subscore for rigidity was equal to the UPDRS-III item 22; subscore for bradykinesia was obtained by adding the UPDRS-III items 23–26 and 31. Cognitive assessments included MMSE for global cognition, Chinese Version of the AVLT for verbal learning and memory, SFT for executive function, DST, and SDMT for attention and working memory45,46. Depressive mood was measured using GDS, REM sleep behavior was measured using RBDQ-HK, quality of daily living was measured using PDQ-39 and the autonomic nervous function was calculated using SCOPA-AUT. The same evaluation instruments were used at all follow-up visits for the PD patients.

Image acquisition and analysis

All participants were scanned on a GE Discovery MR750 3.0 T MRI scanner. Earplugs and foam pads were used to reduce noise and head motion, respectively. High-resolution 3D T1-weighted imaging, DTI, and ESWAN were performed. T1WI were acquired using a fast-spoiled gradient recalled sequence: repetition time (TR) = 7.336 ms; echo time (TE) = 3.036 ms; inversion time = 450 ms; flip angle = 11°; field of view (FOV) = 260 × 260 mm2; matrix = 256 × 256; slice thickness = 1.2 mm; 196 continuous sagittal slices. DTI images were acquired using a spin echo-echo planar imaging sequence: TR = 8000 ms; TE = 80 ms; flip angle = 90°; FOV = 256 × 256 mm2; matrix = 128 × 128; slice thickness = 2 mm; slice gap = 0 mm; number of slices = 67 (axial). Diffusion images were acquired from 30 gradient directions (b = 1000 s/mm2), and included five acquisitions without diffusion weighting (b = 0). ESWAN images were acquired using gradient recalled echo sequence: TR = 33.7 ms; first echo time/spacing/eighth echo time = 4.556 ms/3.648 ms/30.092 ms; flip angle = 20°; FOV = 240 × 240 mm2; matrix = 416 × 384; slice thickness = 2 mm; slice gap = 0 mm; 64 continuous axial slices.

Image processing

T1 structural MRI processing: all parcellations were performed with fully automated segmentation software FreeSurfer (version 6.0.0, http://surfer.nmr.mgh.harvard.edu/) using “recon-all” pipeline which included an automated procedure of motion correction, skull stripping, spatial normalization, registration, cortical parcellation, and volumetric segmentation. To ensure quality control, we visually inspected the segmentations of 68 cortical regions based on the Desikan-Killiany atlas47. The cortical thickness features were defined and extracted according to the Desikan-Killiany atlas, which contained 34 cortical ROIs in each hemisphere. Estimated total intracranial volume (eTIV) was calculated simultaneously.

Diffusion MRI processing: DTI data were processed using the FMRIB Software Library (FSL, http://www.fmrib.ox.ac.uk/fsl). The preprocessing procedures included the following steps: (1) correction of head movement and eddy-current-induced handling motion with the “eddy” tool; (2) skulls were stripped from the DTI data for each participant with the “BET” tool; and (3) fitting of FA with the “DTIFIT” tool. Tract-based spatial statistics (TBSS) were applied for skeletonized WM analysis to mitigate the impact of local misregistration48. The specific steps of TBSS were as follows: (1) DTI metric maps from all subjects were nonlinearly registered to an FMRIB58 FA template and subsequently projected to the mean FA images. (2) A threshold of 0.2 was applied to the mean FA images to generate a mean FA skeleton. (3) The registered FA data from all subjects were then projected to the mean FA skeleton. (4) Nonlinear warping and projection were also applied to the non-FA data, such as MD. The JHU atlas, with its 48 WM ROIs, was used for further WM metrics extraction and analysis. FA serves as an indicator of cellular structures within WM fiber tracts, including axon diameter and fiber directionality, MD is considered as an indicator of the average magnitude of water diffusion49.

QSM data processing: Susceptibility Tensor Imaging Suite V3.0 software package (https://people.eecs.berkeley.edu/~chunlei.liu/software.html) was used to calculate the susceptibility maps from the phase images50. Specifically, the phase images were unwrapped using the Laplacian approach, which relies solely on the sine and cosine functions of the phase angle51. The background phase was eliminated using the V-SHARP method with a spherical kernel radius rising from 0.6 mm at the periphery of the brain to 25 mm toward the center of the brain51. QSM images were calculated using the streaking artifact reduction for QSM (STAR-QSM) method52. The mean signal from the individual brain was used as a susceptibility reference.

To automatically extract subcortical nuclei, we used a deep learning-based end-to-end tool (DeepQSMSeg)22,53, which was developed by setting the ground truth as the previous manually and semiautomatically segmented data21,54, to segment bilateral CN, putamen, globus pallidus (GP), SN, and RN. All the automatically segmented data were checked and manually revised by an experienced neuroradiologist using the ITK-SNAP (www.itksnap.org). Then, the regional magnetic susceptibility of the bilateral CN, putamen, GP, SN, and RN were extracted.

Clustering analysis

As described above, we acquired a dataset consisting of 126 imaging features, encompassing individual cortical thickness from 68 cortical regions defined in the Desikan-Killiany atlas, mean FA values from the skeletons of 48 WM tracts defined in the JHU atlas, and the magnetic susceptibility of 10 segmented subcortical nuclei. To further validate the reliability of the clustering results using FA as the primary WM metric, we also performed an extended analysis using MD as an alternative metric to FA. All the following clustering analyses were performed in R 4.3.1.

To eliminate the influence of covariates such as age, sex, and education, we used a linear regression model to examine the relationship between 126 imaging features and the covariates. Firstly, we identified 115 HC as the robust normative control group. For each imaging feature, we regressed values on age, sex, and education to get the regression weights (β values). Then, we utilized these β values derived from the normative control group to calculate the predicted value for each imaging feature for all participants. While assessing regional cortical thickness, eTIV and mean cortical thickness were additionally considered as covariates.

Subsequently, we conducted dimensionless processing through standardization. To further exclude undesired background noise and improve the clustering performance, PCA was used to extract a lower dimensional component space of the selected features. Specifically, we took the variance explainability >1% as the selection criterion for calculated PC55.

In hierarchical clustering analysis, the Euclidean distance between every pair of patients was calculated. Then, we used Ward’s clustering linkage method to combine pairs of clusters at each step while minimizing the sum of square errors from the cluster mean12,56. Next, we used the ‘NbClust’ package to determine the best number of clusters in our analysis. To explore potential insights from other clustering solutions, we conducted an analysis of the suboptimal clustering solution.

Reproducibility of PD subtypes

In order to validate the reproducibility, we used our validation dataset consisting of 82 PD patients. Following the same methodology for image processing and clustering analysis described before, we extracted 126 imaging features for each patient and classify PD patients into different subtypes. Subsequently, we separately calculated the average z-score of each imaging feature for the discovery and validation datasets after combining each with HC. Spearman correlation coefficient was used as a quantitative coefficient to evaluate the consistency of z-scores within corresponding subtypes across the discovery and validation datasets57.

Other statistical analyses

Statistical analyses comparing demographic, clinical, neuropsychological, and imaging variables at baseline between PD subtypes and HC participants were conducted using IBM SPSS 26.0 statistical software for Windows. The Shapiro-Wilk test was used to assess the normality of the distribution of continuous data. In comparing demographic information between PD and HC, or between the two PD subtypes, the Student’s t-test for independent samples was used for variables that followed a normal distribution; Mann–Whitney U test was employed for variables that did not conform to a normal distribution; and the Chi-square test was used for the assessment of sex distribution difference. Subsequently, general linear model (GLM) was employed to compare the other clinical, neuropsychological, and imaging variables. The comparisons of clinical and neuropsychological variables, as well as imaging variables, were adjusted for age, sex, and education. For cortical thickness comparisons, additional corrections were made for eTIV and mean cortical thickness. The results except for demographic variables were corrected for multiple comparisons using FDR correction, and the significance level was set at p < 0.05. The same statistical methods were applied to the sub-sample of patients with follow-up data, aiming to test whether these information of this sub-sample could adequately represent the entire cohort.

To evaluate the differential progressions of PD subtypes, a LMM was used through the ‘lme4’ package in R 4.3.1. The following formulas were applied to the models:

Model 1: lmer (variable ~ Subtype*Time + age_bl + sex + education +(1 + Time | Subject_ID))

Model 2: lmer (variable ~ Time + age_bl + sex + education + (1 + Time | Subject_ID))

Subtype*Time represented the interaction between subtype and time, while age_bl meant age at baseline. Model 1 was designed to calculate the interaction effect between subtype and time for each clinical and neuropsychological variable. Meanwhile, Model 2 was devised to assess the time effects within each subtype separately, and to extract the progression rates for all individuals within each subtype’s model for further correlation analysis. To further examine the influence of baseline values of the variables, we included them as covariates and applied the same method for calculating longitudinal progression. The results were corrected for multiple comparisons using FDR correction, statistical significance was determined at p < 0.05.

To explore the clinical significance of imaging features, Spearman correlation analysis was conducted between the adjusted brain metrics and clinical assessments using a linear regression model, as described previously in the “Clustering analysis” section. To mitigate any redundant features, we employed PD subtypes as outcomes and subjected all imaging features showing disparities after comparison using GLM to Lasso regression with 1000 random seeds. Then, we retained the top n features (n = the average number of retained features per regression iteration) as imaging features of interest. Lasso regression was performed using the ‘glmnet’ package in R 4.3.1. Subsequently, we selected clinical features that differed at baseline and those with different interaction or time effect between subtypes as features of interest for correlation analysis with the aforementioned imaging features. Statistical significance was determined at p < 0.001, to reduce the false-positive rate58.

To investigate whether the brain degeneration in our study was associated with the lateralization of motor symptoms, we calculated UPDRS-III scores for the left and right sides separately. Based on the laterality of motor symptoms, we classified the brain imaging features of patients into the “Contralateral side of the less affected limb” and “Contralateral side of the more affected limb”. Paired t-tests were conducted using IBM SPSS 26.0 to compare the brain imaging differences between the two sides, and the results of multiple comparisons were corrected using FDR correction.