REM Sleep Behaviour Disorder Dominates Heterogeneity in Longitudinal Analysis of Parkinson's Disease

Introduction¶

Translational research mandates from the National Institutes of Health emphasize the urgent need to characterize the natural history of Parkinson’s disease (PD) and develop objective stratification tools Sieber et al., 2014. Large-scale cohort studies, such as the Parkinson’s Precision Medicine Initiative (PPMI), now provide the infrastructure to discover data-driven progression subtypes and accelerate targeted therapeutic trials Marek et al., 2018. Historically, stratification relied on baseline motor features such as tremor-dominant and PIGD phenotypes Jankovic et al., 1990. However, both clinical and data-driven subtypes are predominantly early-stage phenomena, often coalescing into a more uniform clinical presentation as the disease advances Sauerbier et al., 2016.

Comprehensive multidomain clustering studies have significantly advanced our understanding of PD heterogeneity. By defining phenotypic profiles at a static baseline, standard distance-based clustering approaches have successfully identified subgroups with distinctly divergent clinical outcomes, such as the fast-progressing “diffuse malignant” phenotype Fereshtehnejad et al., 2015Fereshtehnejad et al., 2017Velucci et al., 2025. Notably, this same temporal approach is frequently mirrored even in studies utilizing advanced machine learning frameworks Markello et al., 2021: phenotypic classes are established cross-sectionally, and longitudinal follow-up is only used post hoc to observe the progression of these fixed groups. Conversely, approaches that explicitly subtype patients based on progression rates have demonstrated the immense prognostic value of temporal data Faghri et al., 2018. Yet, these machine learning methods often compress longitudinal follow-up into static summary vectors, obscuring the dynamic shape of the disease course.

To capture actual symptom evolution, univariate latent class models have mapped individual domains over time, such as cognition Pourzinal et al., 2024, autonomic function Chen et al., 2021, and motor severity He et al., 2023. Evaluating these axes in isolation limits our understanding of PD as a multi-system disorder. To bridge this gap, our research undertakes an exploratory investigation using a multivariate longitudinal Latent Class Mixed Model (LCMM). Rather than assuming equal contribution across all symptom domains, this multi-dimensional approach allows the natural variance of the cohort to dictate the clustering, revealing which clinical scales predominantly drive longitudinal heterogeneity.

Finally, we relate these emergent clinical phenotypes to targeted biological metrics. Specifically, further relate these phenotypes to structural MRI, as atrophy patterns track both clinical severity and trans-neuronal spread of PD pathology Zeighami et al., 2015. As highlighted in a recent review Filidei et al., 2025, bridging data-driven subtypes and biological correlates is a critical priority to ensure clinical classifications reflect true pathophysiological differences

Results¶

Latent Class Identification and Trajectories¶

Multivariate LCMM identified a three-class solution as optimal based on the lowest BIC, mean posterior probabilities, class size, and relative entropy criteria (Table 1), see Supp.Methodology for details : a Stable High-Burden class (Class 1, n = 173, 20.2%), a Low-Burden class (Class 2, n = 568, 66.4%), and an Increasing-Burden class (Class 3, n = 114, 13.3%). Observed mean trajectories are shown in Figure 1 (A–D). Class 1 started with the highest REM Sleep Behavior Disorder Screening Questionnaire (RBDSQ) scores, which showed a slight decline by Year 5. Class 2 remained consistently below the cutoff throughout follow-up. Class 3 crossed the cutoff around Year 2 and reached levels comparable to Class 1 by Year 5. For orthostatic systolic blood pressure drop (ΔSBP), Classes 1 and 3 showed increases in later years, while Class 2 remained low. Movement Disorder Society–Unified Parkinson’s Disease Rating Scale Part III (UPDRS Part III) scores increased in all classes, with the steepest rise observed in Class 1. Montreal Cognitive Assessment (MoCA) scores declined over time in Classes 1 and 3, whereas Class 2 remained stable. RBDSQ accounted for the largest variance proportion (39.1%), substantially exceeding ΔSBP (0.74%), MoCA (0.50%), and UPDRS-III (0.20%); the multivariate class structure showed strong agreement with the RBDSQ-only solution (ARI = 0.96; Cramér’s V = 0.95) (Table S1; Table S2). Baseline characteristics are described in Supp.Baseline.

Baseline Neuroimaging and Biomarker Associations¶

Baseline MRI ORs are reported per standard deviation increase to facilitate interpretation, as eTIV-normalized volumes are measured on very small scales. Clinical and biofluid markers remained on their raw scales.

Baseline UPSIT was significantly lower in both the high-burden ($p$ = 0.005) and increasing-burden ($p$ = 0.002) classes (Table 2), validating the model’s capacity to capture external clinical heterogeneity Sun et al., 2025. Whole striatum SBR was also significantly lower in both classes ($p$ = 0.020 and $p$ = 0.012, respectively), indicating reduced overall DAT binding which has often been found for early PD patients with pRBD Xu et al., 2025. However, when caudate and putamen SBR were evaluated jointly (Multivariate DAT), neither showed an independent association with class membership. In multivariate MRI regression, smaller baseline thalamic ($p$ = 0.026) and putaminal ($p$ = 0.023) volumes were independently associated with high-burden class membership, whereas larger baseline pallidal volume ($p$ = 0.033) was independently associated with increasing-burden membership.

Longitudinal Structural Change¶

Longitudinal Linear Mixed Models revealed divergent temporal dynamics between classes (Table 3). The high-burden class showed unadjusted trends of accelerated atrophy across five regions: the amygdala, inferior temporal gyrus, parahippocampal gyrus, superior parietal lobule, and caudal middle frontal gyrus. Conversely, the increasing-burden class exhibited four unadjusted longitudinal trends: slower caudate atrophy alongside accelerated ventricular and choroid plexus expansion.

No regional structural slopes survived False Discovery Rate (FDR) correction across the 20 evaluated regions. Consequently, these longitudinal variations represent exploratory trends rather than definitive trajectory markers.

Early-Stage Subtype Prediction¶

The XGBoost model achieved an AUC of 0.88 and cross-validated balanced accuracy of 0.73. Baseline RBDSQ was the single most discriminative predictor of class membership by both split frequency and information gain, tripling the next-ranked feature. Class-specific SHAP analysis (Figure 1 (F)) revealed that Class 1 was overwhelmingly driven by RBD severity; Class 2 showed a more distributed profile led by RBD, olfactory function, autonomic dysfunction, and anxiety; and Class 3 had the broadest SHAP profile with postural instability and dopaminergic imaging as notable contributors, consistent with its lower classification accuracy (F1 = 0.36). CSF α-synuclein SAA ranked prominently in native tree metrics but not in SHAP values, reflecting population-level rather than patient-level discriminative utility.

Discussion¶

While RBD is a well-established prognostic marker, it is typically deployed as a static or binary baseline feature Velucci et al., 2025Liu et al., 2021. Our findings highlight its dynamic evolution: longitudinal RBDSQ trajectories contributed the majority of variance, driving the latent class structure and demonstrating that RBD-related heterogeneity persists well beyond the prodromal stage.

This longitudinal approach yields markedly different prognostic insights compared to baseline clustering. For instance, the highly cited diffuse/malignant phenotype Fereshtehnejad et al., 2015 couples a wide breadth of severe baseline symptoms with rapid global progression. While our Class 1 shares this severe baseline multi-domain profile (Supp.Baseline) and resembles the RBD+ cluster Velucci et al., 2025, its symptom progression was not exclusively the most rapid. Instead, Class 3 demonstrated the steepest multi-domain deterioration despite much milder baseline symptoms. Because our multivariate LCMM models domain-specific trajectories rather than collapsing metrics into a single composite score, these results suggests that initial symptom burden and progression rate represent distinct, decoupled dimensions of PD heterogeneity.

Our trajectory-derived subtypes align more closely with longitudinal RBD progression studies than with traditional baseline clustering. Our data-driven classes mirror the a priori groups defined by Ye et al. (2022). Their largest group, the non-RBD-stable phenotype, matches our low-burden stable class (Class 2). Our Class 3 strongly aligns with their “late-RBD” group (12.1% of their cohort), showing a late-emerging probable RBD trajectory that crosses the clinical threshold around Year 2. This transitional, high-risk phenotype exhibits baseline olfactory impairment and orthostatic hypotension Y. Saitoh et al., 2018. Our high-burden Class 1 likely represents a combination of their pRBD-stable and pRBD-reversion phenotypes—supported by a slight reversion in Class 1’s RBDSQ during Years 4 and 5.

The independent associations identified in the multivariable MRI regression are not directly comparable with most previous neuroimaging studies, which primarily relied on univariate statistical tests. Kruskal-Wallis tests identified overall differences for the thalamus ($p$ = 0.005) and putamen ($p$ = 0.036), but not the pallidum ($p$ = 0.226) (Supp.Baseline), while post-hoc Dunn’s tests identified no significant pairwise differences. Together, these findings suggest that adjusting for demographic covariates and for the shared variance between regions improved the identification of pairwise class-specific associations. The thalamic and putaminal findings are consistent with their identification as early structural markers in pRBD Boucetta et al., 2016Ellmore et al., 2010Rahayel et al., 2019Salsone et al., 2014. Evidence for baseline pallidal volumetric differences is limited. However, pallidal shape contraction has been reported predominantly in PD patients without RBD Rahayel et al., 2019, while pallidal hypertrophy has also been observed in PD patients with excessive daytime sleepiness compared to patients without Gong et al., 2020, suggesting that larger pallidal volumes relative to the low-burden reference class are not necessarily implausible. Nevertheless, given the relatively small MRI sample for the increasing-burden class (n = 63), this finding requires validation in another cohort.

Although longitudinal structural alterations did not survive strict FDR correction, their unadjusted trends offer exploratory mechanistic insights. The high-burden class exhibited accelerated atrophy in the amygdala, parahippocampal gyrus, and superior parietal lobule. Amygdalar atrophy directly aligns with longitudinal pRBD findings Yoon & Monchi, 2021 and contextualizes this cohort’s significantly elevated baseline depression and anxiety (GDS/STAI). Furthermore, parahippocampal and superior parietal thinning mirror longitudinal structural progression patterns distinguishing RBD from non-RBD PD phenotypes Ye et al., 2022. Conversely, the increasing-burden class trended toward accelerated ventricular and choroid plexus expansion—radiological markers of widespread central atrophy and altered CSF dynamics linked to impaired glymphatic clearance He et al., 2023.

Finally, while we did not directly measure α-synuclein pathology, these distinct trajectories tentatively align with proposed models of Lewy body spatial progression. Class 1’s early concurrent triad of RBD, autonomic, and olfactory dysfunction resembles a “body-first” or brainstem-early trajectory Borghammer & Van Den Berge, 2019Mastenbroek et al., 2024. Conversely, Class 3’s post-motor RBD onset and rapid cognitive decline suggests a “brain-first” origin with delayed, steep brainstem involvement. Class 2, persistently lacking RBD, may represent a phenotype where pathology remains temporarily confined to olfactory regions. Ultimately, the longitudinal timing of RBD expression appears to carry critical, albeit interpretive, pathophysiological weight.

Supplementary material¶

List of Abbreviations¶

Methodology¶

Participants¶

Data acquired from Parkinson’s Progression Markers Initiative (PPMI) dataset (https://www.ppmi-info.org/), on 21 March 2026. PPMI is a multi-center, longitudinal, and observational study that was launched in 2010. Each PPMI site was approved by the appropriate institutional review board before study initiation, and they all fully adhere to the principles set forth in the Declaration of Helsinki. All subjects provided written informed consent prior participation.

Participants were included if they met the following criteria at baseline: (1) drug-naïve with a levodopa equivalent daily dose (LEDD) of 0; (2) disease duration within 2 years; (3) early-stage disease defined by Hoehn-Yahr stage < 3; (4) no dementia; and (5) age onset ≥ 50 years to exclude early-onset Parkinson’s disease. Participants were followed for up to 5 years, and only those with two or more follow-up visits were included, resulting in a total of 855 participants with Parkinson’s disease.

Trajectory analysis¶

Analyses were performed in R (v4.5.3) and Python (v3.12.13). multlcmm function in the R package lcmm Proust-Lima et al., 2017 was applied for trajectory analysis. This approach follows the rationale of group-based trajectory modeling Nagin & Odgers, 2010, allowing several longitudinal markers measured on different clinical scales, to inform a common underlying latent disease process while accounting for marker-specific measurement relationships. Latent classes and individual membership probabilities were estimated within a maximum-likelihood framework, providing asymptotically unbiased parameter estimates under a missing-at-random (MAR) assumption. Follow-up time since baseline, measured in years, was used as the time indicator. The following steps were performed to optimize the analysis:

1. Literature-informed scale selection: Candidate scales were chosen to represent major PD progression domains based on prior subtyping studies. The initial set included RBDSQ, SCOPA-AUT, STAI, age- and education-adjusted SDMT T-scores, and the MDS-UPDRS Part III OFF-medication score (UPDRS Part III) Velucci et al., 2025He et al., 2023. However, modeling above five scales together failed to meet minimum criteria: mean posterior probabilities > 0.7 or minimum class proportions > 5%.

2. Indicator refinement: To ensure model constrction from scales with longitudinal signals and optimal class seperation, we evaluated candidate scales in univariate LCMM, and prioritized scales that had been studied in univariate model. The final set is RBDSQ, MoCAWang et al., 2025, UPDRS Part III, ΔSBP Chen et al., 2021.

3. Link function & distributioanl consideration: Although nonlinear link functions better accommodate ceiling/floor effects and curvilinearity of psychometric scales Proust-Lima et al., 2011, their application in our multivariate framework resulted in reduced classification quality (relative entropy <0.7 or OCC <5). Therefore, a linear link function was adopted for all indicators. To address MoCA’s known curvilinearity and ceiling effect, square root transformation was applied prior to modeling Wang et al., 2025.

4. Random effects specification: Models with random intercept-slope and random intercept only were both evaluated. Both identified a 3-class solution as optimal under linear link; however, the random intercept-slope model did not converge. The final model therefore specified random intercept only, with fixed and mixture components including both intercept and slope terms.

5. Model selection: Models with 1 to 4 classes were fitted. The final model was selected based on lowest BIC, mean posterior probabilities >70%, minimum class size >5%, and relative entropy >0.7 Lennon et al., 2018.

6. Scale contribution assessment: Residual standard error and variance explained proportion were examined to evaluate each indicator’s contribution to the multivariate model.

7. Class assignment validation: To assess the consistency of class solutions, comparisons were performed using confusion matrices, Adjusted Rand Index (ARI), and Cramér’s V.

8. Sensitivity analyses: Models were initially estimated using raw/pre-transformed scores. As z-standardized scores yielded identical class solutions (ARI = 1, Cramér’s V = 1) while facilitating convergence in downstream analyses, z-standardized scores were adopted as the primary model specification.

Missingness and attrition¶

The missing rates for RBDSQ, MoCA, ΔSBP, and UPDRS Part III were 0.88%, 1.09%, 2.99%, and 15.99%, respectively. LCMM accommodates incomplete longitudinal data, so no additional missingness handling was performed. Little’s MCAR test was significant (χ² = 208, df = 28, p < .001), indicating that data were not missing completely at random. Differential attrition was observed across classes, with Year 5 completion rates of 20.2%, 28.5%, and 41.2% for Classes 1, 2, and 3, respectively. As Class 1 also exhibited the overall highest baseline disease burden, attrition was likely associated with observed disease severity, supporting MAR as a reasonable assumption. Although LCMM is expected to limit the impact of differential attrition under MAR, later trajectory estimates for Class 1 are based on a smaller and potentially less severely affected subsample, which may limit their representativeness.

Clinical assessments¶

The above selected four input clinical scales, each representing a core clinical domain (sleep, cognitive, autonomic, and motor):

1. REM sleep behavior disorder (RBD): Defined by a score ≥5 on the REM Sleep Behavior Disorder Screening Questionnaire (RBDSQ), a 10-item self-report instrument (maximum total score 13 points) designed to screen for RBD Stiasny-Kolster et al., 2007.

2. Global cognitive function: Assessed through the Montreal Cognitive Assessment (MoCA), adjusted for education. A score below 26 was used as the cutoff for cognitive impairment Nasreddine et al., 2005.

3. Orthostatic hypotension: Quantified as the orthostatic change in systolic blood pressure (ΔSBP, supine SBP minus standing SBP upon standing). A ΔSBP ≥20 mmHg within 3 min of standing was considered indicative of clinically significant orthostatic hypotension Freeman et al., 2011.

4. Motor severity: Evaluated using the Movement Disorder Society – Unified Parkinson’s Disease Rating Scale (MDS-UPDRS) Part III. A score between 33 to 58 was considered moderate motor impairment Martínez-Martín et al., 2015.

MRI processing¶

Baseline morphological and quality control data were obtained directly from the PPMI repository, derived from the FreeSurfer (v7.3.2) and MRIQC (v23.1.0) pipelines within the nipoppy framework Bhagwat et al., 2023. To extend this to a longitudinal framework while maintaining computational efficiency, we independently processed all participants with available structural MRI at two or more visits (n = 382 of the total N = 855 inclusion cohort) using FastSurfer (v2.4.2) Henschel et al., 2020. A total of 1036 scans were segmented, with volume statistics collated for regions defined by the Desikan-Killiany Atlas Desikan et al., 2006. One participant was subsequently excluded due to technical issues involving missing entries and extreme hemispheric asymmetry. The specific processing parameters and code for this pipeline are documented in the Colab notebooks within our GitHub repository.

Quality control for the Freesurfer data involved a rigorous outlier detection process using the Gap Statistic algorithm Tibshirani et al., 2001 via the GapStatistics Python package Loehr & Musab, 2025. We focused on three key MRIQC metrics: the coefficient of joint variation (CJV), contrast-to-noise ratio (CNR), and entropy focus criterion.

Secondary Analysis¶

Relating latent class models to external variables requires careful handling of estimation bias. The traditional one-step method where covariates and the latent class model are estimated simultaneously often suffers from model instability, as the inclusion of predictors can shift the latent structure itself. To avoid this, many studies fall into the trap of the “naive” three-step method (assigning participants to classes before regression), which produces biased parameter estimates by ignoring classification uncertainty. We instead employed the improved three-step and two-step frameworks developed to account for this uncertainty Bolck et al., 2004Vermunt, 2010Bakk & Kuha, 2018Nylund-Gibson & Masyn, 2016. Specifically, we utilized the externVar function in the lcmm package Proust-Lima et al., n.d., opting for the two-step method over the three-step bootstrap to maintain computational efficiency while achieving comparable bias reduction.

Because the two-step regression function explicitly accounts for latent class assignment uncertainty, it requires substantially more parameters than a conventional regression model. Attempting to fit all predictors into a single unified model led to non-convergence (the algorithm destabilized with >15 predictors) and would have caused severe sample reduction due to varying missingness profiles across modalities. Consequently, a tiered modelling strategy was implemented. Non-MRI variables (UPSIT, fluid biomarkers, and DAT imaging) were evaluated in a complete-case overlap cohort (N = 240, denoted Multivariate 1 in the Results section), while MRI-derived regional volumes and cortical thicknesses were evaluated in a dedicated parallel model (N = 474, Multivariate 2 in the Results section).

To avoid unnecessary loss of information, UPSIT and the SBRs from DAT imaging were additionally evaluated in available-case regressions using their full respective cohorts (N = 834 and N = 846), as restricting these variables to the overlap cohort would have reduced the available sample by more than threefold. Whole striatum SBR was analysed in a univariate regression to provide an overall picture of DAT binding differences between classes, whereas caudate and putamen SBRs were evaluated jointly for any associations with class membership independent of each other (Multivariate DAT in the Results section). The available-case results for UPSIT and DAT imaging are reported alongside the overlap results for these variables, to determine whether restricting the sample primarily affected statistical power or materially altered the estimated associations.

Due to extreme data sparsity of the categorical CSF SAA variable, it was analysed only in its independent available-case model (N = 796) to prevent sparse-data bias from corrupting concurrent covariates, as this effect has been noted to potentially be more detrimental than missing-variable bias Greenland et al., 2016. Furthermore, SAA was binarized into “LBD-like” and “other” (combining MSA-like, inconclusive, and negative results), as these alternate categories were too sparse to be modelled independently even within the larger available-case cohort. For the fluid biomarkers, red blood cells represent a significant source of interference in α-synuclein assays Barbour et al., 2008. To account for this, we leveraged the PPMI hemoglobin (Hb) threshold indicators. Comparative analysis confirmed that α-synuclein levels did not differ significantly in median (Mann-Whitney U, $p$ = 0.627) or distribution (Kolmogorov-Smirnov, $p$ = 0.495) between samples with detectable Hb (N = 60) and those without (N = 265); consequently, the full sample was retained to maximize statistical power. Associations between monogenic PD variants and class membership were not evaluated due to the high prevalence of sporadic cases (N = 808, 94.5%). Similarly, APOE $\epsilon$4 status was binarized (carrier vs. non-carrier) because homozygous cases were too infrequent for independent analysis.

For the cross-sectional multinomial logistic regression, we utilized a mix of volume and thickness metrics derived from FreeSurfer, depending on the specific region. Subcortical volumes were normalized by estimated Total Intracranial Volume (eTIV) to account for head size Voevodskaya et al., 2014, while cortical thicknesses were normalized by mean cortical thickness. For the longitudinal Linear Mixed Models (LMM) investigating atrophy and ventricular expansion rates, we utilized FastSurfer outputs. Because FastSurfer does not provide an eTIV estimate, these volumes were instead normalized using MaskVol; this shift was a necessary adaptation to the respective software pipelines rather than a change in statistical strategy. Notably, FastSurfer segmentation directly provides volumes for cortical regions, allowing us to include them in the longitudinal analysis without the need to run time-intensive surface reconstructions across multiple timepoints.

Given the strict parameter limits of the uncertainty-adjusted regression, structural predictor selection required targeted refinement. Because our derived latent classes were predominantly characterized by their evolving REM sleep behaviour disorder (RBD) phenotypes, we restricted a priori region of interest (ROI) selection to 13 specific cortical and subcortical structures demonstrated in the literature to differ morphologically between PD patients with and without RBD Boucetta et al., 2016Ellmore et al., 2010Rahayel et al., 2019Salsone et al., 2014Lim et al., 2016Ye et al., 2022Yoon & Monchi, 2021. Comparisons to healthy controls were omitted from this rationale as our models exclusively evaluated intra-disease phenotypic progression. To prevent multicollinearity in the regression, left and right hemisphere outputs were averaged across all regions. Additionally, we consolidated the caudal and rostral sub-regions of the anterior cingulate and middle frontal cortices by summing their volumes, as the broader neuroimaging literature rarely distinguishes between these specific Desikan-Killiany parcellation subdivisions. Multicollinearity among the final selected predictors was assessed using Variance Inflation Factors (VIF) via the car library John Fox & Sanford Weisberg, 2019, with all predictors yielding acceptable values (VIF < 4).

The relationship between class membership and longitudinal atrophy was implemented using the hlme function via LMM. We acknowledge that this longitudinal analysis utilized a modal (naive) class assignment, as a corrected bias-adjustment method for LMMs was not available in the current lcmm implementation. Because the LMMs were treated as an exploratory longitudinal extension, the caudal and rostral parts of the cortical regions were analysed independently rather than summed. We also expanded this longitudinal analysis to include several additional whole-brain metrics: the ventricles (as indicators of general central atrophy), the choroid plexus (a marker of altered glymphatic clearance) He et al., 2023, total cerebral white matter (a marker of general structural connectivity), and white matter hypointensities, which served as a marker of white matter lesion burden Wei et al., 2019.

XGBoost¶

We employed XGBoost, a gradient boosting framework optimised for tabular data, to predict LCMM-derived latent classes from baseline features with the aim of developing a lightweight predictive model. The dataset included PATNO, LCMM class assignment, age, sex, race, baseline clinical scales, DaTScan features, genetics, biofluid markers. The dataset was split into training, validation, and test sets at a 70:15:15 ratio. Hyperparameters were optimised using random search over 100 configurations with 3-fold cross-validation. Sample weights and balanced accuracy were used to address class imbalance. Model performance was evaluated using AUC, and SHAP was used to interpret the XGBoost results.

Limitations, Strengths, and Future Directions¶

Cohort and Clinical Measurement Constraints¶

Several limitations should be considered when interpreting these findings. First, while the PPMI dataset provides an unprecedented de novo PD cohort, it represents a highly educated, predominantly white demographic that excludes atypical parkinsonism or early dementia, limiting immediate generalizability to broader clinical populations. Clinically, REM sleep behavior disorder was assessed via the RBDSQ rather than polysomnography, reflecting probable symptom trajectories (pRBD) rather than confirmed diagnoses. Additionally, simplifying the multisystem complexity of PD by representing each clinical domain with a single scale, alongside collecting data post-diagnosis, prevents us from inferring the exact temporal order of early pathological events or $\alpha$-synuclein propagation.

Statistical Modeling, Missing Data, and Distal Outcomes¶

Second, the use of Latent Class Mixed Models (LCMM) inherently assumes discrete categorical subpopulations within a continuous neurodegenerative spectrum. While LCMM natively handles missing longitudinal data, our secondary covariate analysis evaluates baseline variables as class predictors using a bias-adjusted 2-step method (externVar). This framework highlights associations rather than strict causality and is currently restricted to a complete-case subset with overlapping biomarker data. To eliminate missing variable bias and leverage the full $N = 855$ cohort for secondary inference, future work will implement multiple imputation for these baseline biomarkers. Furthermore, while externVar was deployed here exclusively for baseline predictors, it can also be used in future iterations of this work to relate trajectory classes directly to distal clinical outcomes, such as reaching Hoehn & Yahr Stage 3, mild cognitive impairment (MCI), and dementia.

Neuroimaging Feature Selection and Pipeline Consistency¶

Third, our structural MRI analysis was restricted to a targeted subset of mostly subcortical regional volumes. This restricted feature selection was a strict statistical necessity: the parameter-heavy, bias-adjusted 2-step method introduces severe convergence issues if overloaded with variables, constraining our initial model and leaving broader cortical alterations unexplored. To resolve this and expand our feature space into cortical thickness and broader regional volumes, future work will integrate a wider array of cortical metrics. To ensure strict within-subject consistency across these expanded longitudinal metrics, future iterations will transition from standard automated segmentation to a dedicated longitudinal pipeline utilizing FastSurfer with Surf-Recon, contingent on sufficient computational resources.

Study Strengths¶

Despite these limitations, our longitudinal multidomain approach successfully identifies progression phenotypes based on within-person change rather than static baseline severity alone. The robust associations demonstrating clear alignment between these latent trajectory classes and underlying clinical, biomarker, and MRI features strongly support their physiological relevance to dissecting PD heterogeneity.

The 4-class model yielded an empty class (0.00%) and undefined posterior probability (NaN), indicating a degenerate solution.

ARI = 0.956; Cramér’s V = 0.950.

Baseline Characteristics¶

Baseline characteristics across the three classes were compared using the Kruskal–Wallis test for continuous variables. For categorical variables, the chi-square test was used by default; when expected cell frequencies were small (any expected count < 1, or more than 20% of cells with an expected count < 5), an exact test was applied instead. Specifically, Fisher’s exact test was used for 2×2 tables and a Monte Carlo–based Fisher’s exact test for larger R×C tables. Pairwise comparisons between classes were performed using Dunn’s test for continuous variables and the same chi-square/exact-test procedure for categorical variables, with Benjamini–Hochberg false discovery rate (FDR) correction applied across the pairwise comparisons within each variable. Continuous variables are expressed as mean ± standard deviation (SD), and categorical variables as number (percentage).

Baseline clinical characteristics showed that the three classes did not differ significantly in age at diagnosis or symptom onset, disease duration, education, Hoehn-Yahr stage, dominant side of symptom onset, motor severity (UPDRS Part III), global cognition (MoCA), or the battery of neuropsychological tests. Class 1 carried the heaviest overall disease burden, with the highest RBDSQ (8.65±1.83), SCOPA-AUT (14.33±7.49), STAI (66.69±18.46), GDS (2.78±2.58), ESS (6.42±3.83), QUIP (0.37±0.80), and MDS-UPDRS Part I (8.28±5.12), all significantly exceeding Class 2 (p < 0.05). Olfaction was impaired (UPSIT 20.33±7.76, p < 0.001 vs Class 2), and orthostatic blood pressure drop was elevated (ΔSBP 6.61±13.64, p = 0.006 vs Class 2), however these two features were most severe in Class 3. This group also had the greatest motor and functional impact, with the highest MDS-UPDRS Part II, PIGD score, lowest ADL score, and a marked male predominance (80.9% vs 61.4% in Class 2, p < 0.0001). Class 2 repesented the low-burden group, with overall mildest non-motor, motor and functional symptoms, serving as the reference against which the other two groups were defined.

Class 3 occupied an intermediate position overall, with a baseline profile that showed early signs of divergence from Class 2. Compared with Class 2, Class 3 had significantly worse olfactory function (UPSIT 19.85 ± 6.24, p < 0.001), higher autonomic burden (SCOPA-AUT 11.43 ± 6.43, p = 0.002), higher PIGD scores, lower ADL scores, and a greater male predominance (all p < 0.05). Although baseline RBDSQ was significantly higher than Class 2 (3.82 vs 2.77), it remained below the diagnostic cutoff. ΔSBP was numerically highest in Class 3 (6.82 ± 15.07) but did not differ significantly from Class 2. Compared with Class 1, Class 3 remained less globally impaired, with significantly lower RBDSQ, SCOPA-AUT, PIGD, and MDS-UPDRS Part I/II scores. Overall, Class 3 was characterized by selective early abnormalities, most notably hyposmia and autonomic dysfunction.

Baseline biomarker and neuroimaging profiles were largely comparable across classes. CSF Aβ, total α-synuclein, phosphorylated tau, total tau, serum NfL, and APOE ε4 burden showed no significant differences. Although serum urate differed modestly among classes (p = 0.034), with Class 2 showing the lowest levels, this class also had the lowest proportion of male participants, so the difference may partly reflect sex imbalance. Pairwise comparisons between classes were not statistically significant. CSF α-synuclein seed amplification assay (SAA) status differed among classes (overall p = 0.014), with Class 2 showing a higher proportion of SAA-negative cases (11.7%) than Classes 1 (5.0%) and 3 (4.6%). Age-/sex-expected lowest putamen ratio was similar across groups, indicating comparable baseline nigrostriatal dopaminergic deficits. MRI revealed several nominal overall differences, including thalamus (p = 0.005), putamen (p = 0.036), nucleus accumbens (p = 0.009), third ventricle (p = 0.024), total gray matter (p = 0.008), subcortical gray matter (p = 0.022), and cortex volume (p = 0.009), though pairwise comparisons did not survive FDR correction. Class 2 generally showed numerically larger brain volumes. Cortical thickness measures did not differ significantly across classes.

Detailed values for each variable are provided in the table below.