Sequence of clinical and neurodegeneration events in Parkinson’s disease progression

See Le Heron et al. (doi:10.1093/brain/awab060) for a scientific commentary on this article. Using event-based modelling, Oxtoby et al. reveal the fine-grained probabilistic sequence of changes in clinical, cognitive and neuroimaging measures during Parkinson’s disease progression. Performance on visual tasks is impaired early on, typically before neurodegeneration and retinal thinning can be detected.


Cohort comparison
Supplementary Table 1 is a summary comparison of our discovery and external cohorts. We compare the cohorts statistically within groups: controls, PDD-LR and PDD-HR. Highlights include the following: • The PPMI cohort is younger on average (but we adjusted for age effects).
• PPMI patients had a shorter disease duration (and lower UPDRS scores), but the patient groups in both cohorts had comparable age at onset of symptoms. • Our discovery cohort controls had more females (but we adjusted for gender effects).
• PPMI patients had slightly lower MoCA scores.

Cross-cohort comparison, within groups.
Each p value shown is for a Mann-Whitney U-test (median) or χ 2 -test (proportion) of the null hypothesis that there is no statistical difference between the cohorts, within each clinical group: Controls, PDD-HR, PDD-LR. Comparisons of olfaction between the cohorts is made possible by converting UPSIT (PPMI) to a Sniffin' Sticks equivalent score using an equipercentile method (Lawton et al., 2016). Abbreviations: PD -Parkinson's disease; PDD-HR -PD dementia high-risk; PDD-LR -PD dementia lowrisk; UPDRS -Unified PD Rating Scale; F/M -female/male; RBDSQ -REM Sleep Behaviour Disorder Screening Questionnaire; MoCA -Montreal Cognitive Assessment; GNT -Graded Naming Test. MAPT +/refers

Mixture models for event probabilities
Patients have different marker profiles that reflect disease severity, e.g., different scores on clinical tests and different severity/extent of neurodegeneration. The event-based model exploits this fact to infer a longitudinal sequence of events from a cross-sectional sample by assessing marker severity across individuals. Severity is quantified using a mixture model that assigns "normal" and abnormal" labels to each individual's marker data (be they patient or control) separately. This is done in a data-driven fashion, usually starting with all data from controls assigned the "normal" ("pre-event") label and data from patients assigned the "abnormal" ("post-event") label, with labels swapped for outliers in each group until convergence. In the original event-based model demonstrated in sporadic Alzheimer's disease (Young et al., 2014), the mixture modelling was implemented by allowing only patient labels to be swapped, under the assumption that controls are extremely unlikely to be prodromal patients due to stringent inclusion criteria. In this work we take the same "fixedcontrols" mixture modelling approach but allow for outliers to swap labels -marker data from controls that is in the 90 th percentile or above (patients and controls together) is allowed to participate in the label swapping of the mixture modelling. We also inform the mixture modelling with prior information on disease direction for each biomarker. We encode these innovations within the nonparametric mixture modelling approach introduced in (Firth et al., 2020) and available at https://github.com/noxtoby/kde_ebm_open. Together, our innovations add clear interpretability to our models: they represent Parkinson's disease progression in patients at elevated risk of dementia as disease-specific deviations from normality.
Supplementary Figure 1 shows histograms and representative mixture models for selected features included in our model ( Figure 2, main manuscript). We include at least one of each  (Firth et al., 2020). Data shown has been detrended for gender, age and years of formal education as described in the main manuscript.

Combining event-based models from cross-validation
An event-based model consists of a sequence and uncertainty in that sequence. This is represented by a square probability matrix, with the th row being the posterior density for biomarker over sequence positions . We call this a positional density map, or a positional variance diagram (Fonteijn et al., 2012). Typically, the sequence is determined using maximum-likelihood estimation, and the density is determined by multiply initialised MCMC samples from the posterior.
In this paper we used 10-repeated 5-fold cross-validation which produces 50 event-based models, each trained on different partitions of 80% of the data. We average the posteriors of these 50 models to form our cross-validated posterior density map.

Cross-validation sequence
Aggregating the 50 model sequences is not straightforward because rank aggregation is not a solved problem and depends on the context, e.g., voting (Arrow, 1950), and event-based modelling (Huang and Alexander, 2012;Venkatraghavan et al., 2019). Rather than averaging the 50 model sequences, we introduce the concept of "event onset" to define the CV (cross validation) sequence. This is motivated by the original hypothetical models of dynamic biomarkers in Alzheimer's disease Frisoni et al., 2010; which inspired the original event-based model. The hypothesis is that biomarker change/dynamics occurs in a sequential manner. We quantify biomarker change using the derivative of the cross-validated cumulative positional density map (row-by-row for each event separately), then use this to define event onset as the point of maximum rate-ofchange, subject to reaching a minimum cumulative abnormality of > 0.25. We also use the value of cumulative abnormality as a tiebreaker, with higher cumulative abnormality ranking an event earlier. If the tie remains, then we rank higher the earliest marker to reach 50% cumulative abnormality using min� ( ) − 0.5�, where is the sequence position.
Mathematically, event onset for marker is the sequence position that maximises the derivative of posterior cumulative positional density We illustrate this in Supplementary Figure 2 using our full model from Figure 2 in the main manuscript. Supplementary Figure 2A is a representation of the right panel of Figure 2: the cumulative positional density map consisting of curves for events (curves smoothed using kernel density estimation for visualisation purposes). Supplementary Figure 2B shows the corresponding derivatives / with red dots showing event onset as defined above. Supplementary Figure 2C is a 2D scatter plot of cross-validated sequence position against event onset, which reemphasises the uncertainty found later in the model sequence (see discussion around Figure

3.1.2.
Similarity reference: randomised models As described in the main manuscript, we use the Bhattacharyya coefficient as a similarity measure between event-based models across the CV folds. This summarises the statistical overlap of the posterior positional density maps as a single number between zero and one. We calculate pairwise between all 50 CV models, then calculate the average similarity of each model (fold) with the other 49. Our overall CV similarity is reported as the mean and standard deviation of these 50 similarity scores (see Supplementary Section 3.2 below).
To provide some context, we calculate a reference distribution for by repeating the above procedure on 100 randomisations of the posterior positional density maps. As shown below in Supplementary Figure 3 (right, inset), this gives a reference similarity score of 0 = 0.37 ± 0.02 for randomised posteriors.

Cross-validation: supplementary results
Supplementary Figure 3 shows CV similarity of our posteriors by event (left panel) and across folds (right panel). Events having high similarity across CV folds showed consistent posterior positional density. Events having lower similarity across folds may be due to heterogeneity in Parkinson's disease progression, e.g., smell-early vs smell-late patients. See Discussion in the main manuscript for planned future work to unravel this heterogeneity.

Robustness to PDD-HR threshold (age at onset)
The average age at diagnosis of Parkinson's dementia is not known, but the literature suggests that Parkinson's dementia is diagnosed on average 10 years after Parkinson's symptom onset and is most common after age 70 (Lewis et al., 2016;Williams-Gray et al., 2013). Thus, Parkinson's patients are at elevated risk of progressing to Parkinson's dementia if their age at onset (of Parkinson's diagnosis, not dementia) exceeds 60 years old.
To fit our models of progression to Parkinson's dementia, we focussed on high-risk patients ("PDD-HR") whose age at onset was ≥ 65. We verified the robustness of our results by repeating the main experiment using PDD-HR age at onset thresholds of 60 and 70 years, then comparing: 1) the probability of abnormality for each event/feature included in the model -the more this score changes with the threshold, the more likely it is that the event-based model sequence will differ; and 2) the event-based model sequences using Kendall's tau rank correlation.
Supplementary Table 2 shows the effect of the PDD-HR age at onset threshold on the probability of abnormality. We summarise this in a single number by reporting the marker value for which ( |Event) = 0.5, as well as the relative ratio / 65 to show the magnitude of the change with respect to the threshold used in our main experiments. We report only those markers affected by the PDD-HR threshold by at least 5%. We report covariateadjusted values (age, gender, years of education: see main manuscript for details).
We found that increasing the PDD-HR age at onset threshold tended to accentuate the value of for which ( |Event) = 0.5, for most markers. That is, the event "occurs" earlier -at a less severe score or biomarker value -because these older patients have more severe pathology and symptoms. For example, QSM markers decrease, and symptoms get worse. This supports the existing notion that older age at onset bestows a higher risk of progressing to dementia. However, this tendency was not universal, which we suspect is due to disease heterogeneity that we hope to unravel in future work.  Figure 4 is a scatter plot of event-based model sequences using alternative PDD-HR age at onset thresholds. We found high (and statistically significant) rank correlation of the estimated sequences for age at onset thresholds of 60 (PDD-HR n=50, blue dots) and 70 (PDD-HR n=18, orange stars) with our main model that used a threshold of 65 (PDD-HR n=36). Kendall's tau correlations were = 0.65 ( < 2 × 10 −9 ) and = 0.58 ( < 7 × 10 −8 ), respectively.

Cortical thinning model
Supplementary Figure 5 is a visualisation of the model-based pattern of cumulative abnormality in cortical thickness from Figure 4 (main manuscript), using Brain Painter (Marinescu et al., 2019). Cortical thinning is a late event in our data-driven model of Parkinson's disease progression, with occipital and parietal lobes being the earliest regions affected.
Visual representation of cortical thinning model

Supplementary Figure 5. Estimated sequence of cortical thinning (discovery cohort).
Visualisation of cortical thinning generated using Brain Painter (Marinescu et al., 2019). Lower: right hemisphere. Upper: left hemisphere (reflected orientation for comparison). Regional cumulative abnormality (0 to 1) is proportional to colour intensity.

Cognitive decline in patients enriched for dementia risk
We analysed cognitive decline on the MoCA in both cohorts. Supplementary Table 3 shows the results of comparing the PDD-HR and PDD-LR groups in each cohort at baseline and followup. Supplementary

Genetic variation
We investigated dementia-specific genetic variation in both cohorts. We excluded a small number of GBA-positive participants (discovery: n=7 patients; external: n=46 patients, n=10 controls), which is a risk factor for more rapid progression to dementia, especially in younger PD patients. We excluded these cases as they are likely to have a more rapid progression to dementia and may show a divergent sequence of events (Blauwendraat et al., 2020). The low number in our cohort precludes a statistical analysis, but we note that our analyses are essentially unchanged when these patients are included: the estimated sequences with/without GBA cases (not shown) have high rank correlation with Kendall's tau of = 0.75 ( = 2 × 10 −12 ).
For included patients with available genetic information we examined genetic variation due to MAPT (H1/H1 haplotype) and APOE4 (at least one e4 allele), both of which increase risk for dementia. Specifically, we performed Mann-Whitney U tests of group differences in two key Parkinson's disease dementia measures: UPDRS-3 (motor) scores and MoCA scores.  (14) The data can support comparisons between genetic subgroups (carriers vs non-carriers) of patients (PDD-HR and PDD-LR combined) for APOE4 and MAPT. The genetic breakdown for both cohorts is as follows:  (222) Mann-Whitney U tests of carriers versus non-carriers in UPDRS-3 and MoCA all produce ≥ 0.19 (discovery cohort) and ≥ 0.13 (PPMI) -see below. From this we conclude that these dementia-specific risk factors do not influence PD progression in our study and that larger numbers may be required to detect differences in sequences between these genetic variations.

Carriers
Non

Cross-cohort Model Comparison
Supplementary Figure 6 shows comparable event-based model posterior density maps for the discovery and external cohorts. For comparison, individual events have been averaged into modality-based features. For example, all diffusion weighted imaging features of white matter neurodegeneration are averaged into a single row labelled DWI -likewise for cortical neurodegeneration from MRI, and verbal fluency. The ordering (vertical axis) was determined by the median position in the full sequence ( Figure 2, main manuscript). The models show high concordance with Kendall's tau rank correlation of = 0.87 ( = 0.017) and Bhattacharyya coefficient of = 0.96.
Supplementary Figure 6. Cross-cohort comparison. Event-based model posterior density maps for discovery cohort (left) and external cohort (right), with imaging events (DWI neurodegeneration in the substantia nigra, MRI cortical thinning) and verbal fluency events averaged together into individual rows.

Dysautonomia in Parkinson's progression: SCOPA-Aut
Baseline measures of autonomic function are not available in our discovery cohort but are available in the external data set (PPMI) via the SCOPA-Aut score (Scales for Outcomes in Parkinson's Disease -Autonomic Dysfunction). Supplementary Figure 7 shows the corresponding event-based model with SCOPA-Aut included. The positional density map shows that, as expected, dysautonomia is an early event along with olfactory and sleep abnormalities.

Supplementary Figure 7. Event-based model of Parkinson's progression including autonomic function.
The data and model estimate that autonomic dysfunction (SCOPA-Aut score) is an early event in Parkinson's disease progression.