Lipidomics Prediction of Parkinson’s Disease Severity: A Machine-Learning Analysis

Participants and clinical evaluation

The cohort analyzed was as described in an earlier study [4]. In brief, participants in the ‘Spot’ study [22, 23] were recruited between 2010–2016 from the Center for Parkinson’s Disease at Columbia University Irving Medical Center in New York, NY. Among the participants, we randomly selected PD patients (n = 150). All participants were non-carriers of SNCA, LRRK2, and GBA mutations (with the exception of one PD case who was a carrier of the GBA variant E326K, n = 1). Evaluation of all participants included the MoCA, and the UPDRS in the “on” state. Medical history, current medications, and demographics were also collected. The original study also included 100 controls, which were not analyzed here, as our primary aim was to test the link between lipids and disease severity. Data from these controls are presented only in Fig. 2 for illustrative purposes. Plasma collection for lipidomics was as previously described [4]. In brief, blood samples were collected in single EDTA tubes (10 cc in volume) and centrifuged, after which 1 cc plasma aliquots were extracted and stored in a –80°C freezer within one hour of collection. Table 1 includes all the lipids measured in the study and lipid class membership. Informed consent was signed by all participants and the Columbia University IRB approved all study procedures. Since the UPDRS III and UPDRS total score were missing for a single patient, in this analysis, we used the data of 149 participants.

Analysis of lipid classes and PD severity

The first phase of the study focused on identifying the most influential lipid classes that predicted PD severity as determined by the UPDRS and MoCA scores. We predicted both the UPDRS III and UPDRS total scores, although the latter composes of UPDRS I, II, and III and thus depends on the former. The 37 lipid classes, two clinical variables [age at onset (AAO) and disease duration], demographic variables [i.e., age, gender, height, body mass index (BMI), marital status], and education were the predicting variables.

The root mean square error (RMSE) [25], which is the most common measure to evaluate prediction performance, was calculated as an average over a Monte Carlo cross validation, where in each of 200 datasets, 80% of the (149) observations were sampled randomly (without replacement) for training the algorithm, and the remaining 20% used for testing it. This resampling method was chosen to reinforce the significance of the research results when experimenting with the cohort.

In this phase, we wanted first to validate our ML algorithm, the RF, in predicting PD severity. The RF is a popular accurate classifier, which was shown in many previous studies to be efficient and beneficial in various clinical [21, 26–28] and non-clinical [29, 30] domains. The RF classifier [31, 32] holds no assumptions about the data distribution, can cope with very complex problems with minimum overfitting [31, 33–35], and ranks variables by their contribution to accurate (or informative) prediction [36–38]. In this phase, we selected variables that were ranked in the top 30% as contributory by RF. For further explanation about the RF, see the Supplementary Material.

In our case, we expected the RF to help interpret the prediction results by indicating interrelations among lipid classes and species with different stages in the disease. Since the RF has never been used in predicting PD severity, we first validated it through conventional statistical methods. Each statistical method is limited by the assumptions it makes, hence we employed three of them. Each statistical method explored the data from a different angle and selected contributing variables, as specified below. The RF together with the statistical methods create an ensemble of classifiers (though not independent) that is used here as “a panel of experts”.

The statistical models that were used to validate and then to enforce the RF algorithm are: (1) Univariate linear regression (UR), where a variable is selected if considered statistically significant with p < 0.05; (2) Multivariate linear regression (MR), implemented by considering variables using sequential backward selection, sequential forward selection, or sequential floating forward selection (i.e., stepwise), where the model that achieves the lowest Akaike information criterion (AIC) [39] is chosen, and all variables that are statistically significant (p < 0.05) are selected; (3) Lasso regression [40], which is based on ordinary least squares with a penalty factor for the number of variables, forces some of the coefficients to equal to zero, and thereby selects the variables with coefficients greater than zero.

We let each of the four methods identify the most contributing variables among the predictors and kept only variables that at least three methods indicated as contributing highly to prediction of the PD severity.

Analysis of lipid species and PD severity

The second phase was done only for the UPDRS since the first phase results for the MoCA showed poor relation to the lipidome. The same methodology of the first phase was implemented here. However, in this phase, we faced another challenge since the number of variables exceeded the number of examples (patients); that is, the 517 lipid species and eight demographical and clinical variables were much greater than the 149 patients, which made variable selection even more challenging than in the first phase. To accommodate this, we selected only lipid classes that had at least two supporting methods (three methods were employed in the first phase). Then for the lipid species of these classes, we added the eight demographical and clinical variables and let the RF rank these variables as in the first phase. Finally, we sequentially removed lower-ranked variables as long as it did not worsen the RMSE, while keeping at least 10% of the variables to allow for retention of a sufficient number of them.

First in our research plan, we identified and ranked the most contributing lipid classes and species for UPDRS total and UPDRS III scores in a two-stage approach (first classes and then species of the most contributing classes) using the RF classifier that is validated using statistical methods. Second, we examined values of plasma concentrations of the most contributing lipid classes and species in different PD severity levels and considered their sensitivity to gender and age. Third, we identified lipid signatures (combinations of classes, species, and demographics) for different severity levels, and evaluated their interrelations with age, AAO, and disease duration and the impact of the combined signatures on the prediction model accuracy.

Data description

Table 2 presents descriptive statistics for the demographical and clinical data of the 149 PD patients’ cohort. The AAO is platykurtic (kurtosis = –0.52, skewness = 0.029), the disease duration is leptokurtic and positively skewed (kurtosis = 1.99, skewness = 1.35), with 82% of the patients with disease durations shorter than 10 years. Gender, by design, is evenly distributed (75 males and 74 females). The medication (Levodopa) dosage distribution for different UPDRS III scores is presented in Supplementary Figure 1, showing, as expected, that the dosage increases with the severity.

Lipid classes for prediction of disease severity

Table 3 shows the most contributing variables (as defined above and based on at least three out of four algorithms) from all lipid classes and the demographical and clinical variables for the UPDRS total and UPDRS III scores. As expected, age and disease duration were found among the most influential variables especially for UPDRS total score. Our data shows that dihydroglobotriaosylceramide (dhGB3), dihydrosphingomyelin (dhSM), glucosylceramide (GlcCer), and plasmalogen phosphatidylethanolamine (PEp) are the most contributory lipid classes in prediction of both UPDRS total and UPDRS III scores. While the UR and MR sometimes gave conflicting results regarding the significance of a variable, RF and Lasso were more consistent. Similar to age and disease duration, higher GlcCer, dhGB3, and globotriaosylceramide (GB3) values were associated with higher scores (the first two classes for both UPDRS total and III scores). In contrast, dhSM and PEp (for both UPDRS total and III scores) as well as monoacylglycerol (MG), phosphatidylethanolamine (PE), and dihydro GM3 (dhGM3) (only for UPDRS III score) affected UPDRS negatively. In addition, lipid classes sulfatide (Sulf), free cholesterol (FC), acyl phosphatidylglycerol (APG), lysophosphatidylcholine (LPC), phosphatidylglycerol (PG), and dhGM3, as well as lipid classes triacylglycerol (TG), dihydrolactosylceramide (dhLacCer), lysophosphatidylinositol (LPI), APG, and Sulf had only two supporting methods for UPDRS total and UPDRS III scores, respectively, and thus are not presented in Table 3. Additionally, BMI was not found to be a contributing variable.

Average plasma concentrations divided by tertiles of the UPDRS III and UPDRS total scores for the cohort are presented in Fig. 1 for dhSM, GlcCer, PEp, and dhGB3. Figure 1 shows that in most of the cases (e.g., dhSM and PEp), a clear trend can be observed (either increasing or decreasing monotonicity) between the scores in the lower and higher cohorts for both UPDRS total and UPDRS III scores. The figure also demonstrates that while female patients have higher concentrations than male patients in the lower score ranges (both UPDRS total and UPDRS III scores), the difference almost always becomes (very) small in the higher severity level.

OPEN IN VIEWER

Fig. 1

Average plasma concentrations (pmol/μl) for different ranges of the UPDRS total score (left) and UPDRS III score (right) for four contributory lipid classes common to UPDRS total and III scores (Table 3). Numbers in parentheses indicate numbers of observations (patients) in different ranges of severities. The black solid, red dotted, and blue broken lines represent all patients, female patients, and male patients, respectively. The distributions of females and males over UPDRS III and UPDRS total scores ranges are given at the top of the figure. The confidence interval (CI) results for multi comparisons between the severities are based on a post-hoc Tukey HSD test and those between females and males on a two-sample t-test [both CIs are presented below each plot only for statistically significant (p < 0.05) differences].

Figure 2A demonstrates differences between PD patients in three tertiles of their UPDRS III scores with respect to values of their lipid classes and demographics, showing distinct demographic-lipidomic signatures of patients with different disease severities in comparison with controls.

OPEN IN VIEWER

Fig. 2

A) Radar plot representing UPDRS III scores and average values of the most contributing lipid classes (Table 3), age, AAO, and gender for three PD stages: early (low severity; yellow), intermediate (red), and late (high severity; blue) in comparison to control (green). The controls are 100 subjects genetically unrelated frequency-matched by gender and age (Avg. = 66.11 years, std. = 9.4) with average and std. BMI values of 25.7 and 4.4, respectively. AAO for the control was set to the value of the low severity. Numbers in parentheses indicate numbers of observations (patients). As we move along an axis from the center/origin of the outer polygon to any of its vertices, the corresponding variable increases its value, e.g., towards higher age, AAO, and lipid class concentrations (and for higher percentages of men than women for the gender variable). The figure demonstrates differences between PD patient groups in three tertiles of UPDRS III scores of 0–11 (lower tertile), 12–20 (middle tertile), and 21–48 (upper tertile) in comparison to control. Among the patients, early stage patients (lowest tertile), mostly women, are recognized by the lowest values of GlcCer (and, of course, of age and AAO values), with extreme values of dhSM, PEp, and PE; the latter value is extremely larger than those of the intermediate and late stage patients (and similar to that of control). Patients with intermediate UPDRS III scores (middle tertile) have dhSM and PEp values similar to those of early stage patients, but also the highest levels of GlcCer and lowest levels of dhGB3. The patients with the highest UPDRS III scores (upper tertile), as expected, are older and have higher AAO values, highest levels of dhGB3, intermediate levels of GlcCer, and PE, and the lowest values of PEp and dhSM. When considering also the control group, we notice two results: a) a perfect order of values from control to low, intermediate, and up to high severities for dhSM and almost perfect order for GlcCer and PE (the two highest severities swapped their order although are very similar in values), and b) for most lipid classes (except PEp and dhGB3), lipid levels of controls are at the extreme ends, i.e., having either the lowest or highest level. B) Similar to (A) for a combination of average (over the patients) values of dhGM3 and its most contributing species. The figure shows that patients in the lower UPDRS III score tertile have the highest values of dhGM3, dhGM3(22:0), and dhGM3(22:1) but the lowest values of dhGM3(16:0). The case is nearly opposite for patients in the highest UPDRS III score tertile, i.e., the lowest values of dhGM3, dhGM3(22:0), and dhGM3(22:1), and relatively high values of dhGM3(16:0). Interestingly, while the intermediate disease severity is characterized by intermediate values of dhGM3, dhGM3(22:0), and dhGM3(22:1), it is also characterized by the highest values of dhGM3(16:0). C) Similar to (A) for a combination of average (over the patients) values of GlcCer and its most contributing species together with age and AAO. The figure shows patients in the higher tertile of UPDRS III scores were older with later AAOs, but also had the highest values of GlcCer(16:0) and relatively high values of GlcCer and all of its other species. On the other hand, patients with the lowest UPDRS III scores are the youngest with the lowest AAO on average and have the lowest values of GlcCer and its species, except for GlcCer(24:1), where these patients have the highest plasma levels. D) The same as (C) including the control group for a reference (the small values of the control subjects changed the scale drastically so previous differences among PD severities were diminished), showing that none of the species of the GlcCer class is expressed in the control subjects, and thereby emphasizing the role of this lipid class in PD.

Finally, the results for the MoCA indicate four variables that had three or four supporting methods: age, education, disease duration, and PE. Age and disease duration were negatively correlated with the MoCA score, the older the patient or the longer the disease duration, the lower the score. Education was positively correlated with MoCA score; the more educated the patient, the higher the score. Note that the MoCA analysis is the sole place where education came up as contributable. PE has a strong monotonic negative impact on the MoCA score (data not presented).

Lipid species for prediction of disease severity

We started this analysis with 517 species and 8 demographical and clinical variables and then minimized this list using the feature selection ML methodology (described in the Material and Methods section), identifying the species that contributed the most for PD severity. This was implemented for UPDRS total and UPDRS III scores, as presented in Tables 4 and 5, respectively. As in the lipid class analysis, age and disease duration are highly ranked in predicting PD severity in the lipid species analysis of UPDRS total score (Table 4). Predictors of the UPDRS total score included 18 species from 11 classes, but the majority of the influential species (n = 15; 83%, Table 4) are from five classes: PEp (5), GB3 (3), dhSM (3), GlcCer (2), and dhGM3 (2). In addition to age and AAO that are influential on both UPDRS total and UPDRS III scores, half of the species (9/18) are influential on both scores (with the same directionality). Most of the classes (9/12) identified for UPDRS III in the first phase have representative species in Table 5. There are three species whose directions are not corresponding with those of their class directions. Two of them, PEp (42:7) and dhGM3(16:0), also appear influential on the UPDRS total score with no corresponding directions with their classes (Table 4).

Average plasma concentrations of the cohort divided by tertiles based on the UPDRS III score are presented in Fig. 3 for most of the contributing lipid species in Table 5. In most of the cases [e.g., GlcCer(16:0), dhSM(20:0), dhGM3(22:0), and dhGM3(22:1)], the cohort concentration increases/decreases linearly with the disease severity, whereas in other cases [e.g., GlcCer(24:1) and dhGM3(16:0)], a point in an intermediate severity shows a change in the trend of the concentration. We believe this change point is responsible for the mismatch in directionality of impact between a species and its class, as was observed before [e.g., dhGM3(16:0) in Tables 4 and 5]. We also see a similar pattern to that of lipid classes (Fig. 1), where female patients usually have higher concentrations in the low and medium severities, but the difference to male patients usually vanishes for the highest UPDRS III scores. Similar to the lipid class analysis, we show in Figs. 2B-2D radar plots presenting distinct demographic-lipidomic signatures of patients with different disease severities in comparison with controls.

OPEN IN VIEWER

Fig. 3

Average plasma concentrations (pmol/μl) for different ranges of the UPDRS III score for species that were identified among the eighteenth important ones (Table 4). Numbers in parentheses indicate numbers of observations (patients) in different ranges of severities. The black solid, red dotted, and blue broken lines represent all patients, female patients, and male patients, respectively. The distributions of females and males over UPDRS score ranges are given at the top of the figure. The CI results for the comparisons between the severities are based on a post hoc Tukey HSD test and those between females and males on a two-sample [both Cis t-test, are presented below each plot only for statistically significant (p < 0.05) differences.]. Note that dhGB3(16:0) is the only species in the dhGB3 class, what explains the identity between the corresponding figures here and in Fig. 1.

Prediction model accuracy for UPDRS III and UPDRS total scores as measured by the RMSE is presented in Fig. 4, with an average of 6.81 and 9.42, respectively. We see that our models are most accurate in predicting UPDRS III and UPDRS total scores in those 50% of the patients whose scores range between 9 and 22, with RMSE values of 2.85–5.72 and 4.48–7.56, respectively, and less accurate in those extreme scores (with UPDRS III score under 5 or over 27). As Supplementary Figure 2 shows, this difference in accuracy is because of few extreme prediction errors due to outliers and high variability in the patient data in the extreme severity ranges.

OPEN IN VIEWER

Fig. 4

Average prediction RMSE values for UPDRS III and UDPRS total scores and relatively equally populated UPDRS III score ranges. RMSE = 1 n ∑ j = 1 n ( y j - y ˆ j ) 2 , where n is the number of observations (patients), y _j is the actual UPDRS total and UPDRS III scores for the j ^th patient, and y ˆ j is the corresponding predicted value. Numbers in parentheses indicate numbers of patients for different severity ranges. The figure presents the RMSE of the RF trained using the most contributing lipid species (Tables 4 and 5) in eight ranges of UPDRS III score severities. While the average RMSE for UPDRS III and UPDRS total scores for all 149 patients is 6.81 and 9.42, respectively (not shown in the figure), the figure shows that for 50% (75) of the patients, those with UPDRS III scores between 9 to 22 (the four middle severity ranges), the average RMSE in predicting UPDRS III score is between 2.85 and 5.72, and that in predicting UPDRS total score is between 4.48 and 7.56. These RMSE values are much smaller than the 9.3 and 13.8 standard deviations of UPDRS III and UPDRS total scores in our data, respectively (Table 2). For another ~25% (38) of the patients, in severity ranges of 6–8 and 23–26, the average error is between 7.07 and 7.57 for UPDRS III score and between 9.79 and 10.11 for UPDRS total score, which are also smaller compared with the corresponding 9.3/13.8 standard deviations. However, for extreme (low and high) UPDRS III scores, 24% (36) of the patients have average error higher than the standard deviation, which is attributed to the high variance of, and outliers in, RMSE in these severity ranges (Supplementary Figure 2). Also note from comparing the graphs for UPDRS total and UPDRS III scores, that errors due to UPDRS I and UPDRS II scores (the difference between UPDRS total and III scores) are greater for extreme severity values.

Supplementary Figures 3A and 3B show the result of excluding age, AAO, and disease duration (the variables and their interactions with the lipids) from the models, relying only on the prediction capability of the contributing lipid species (Tables 4 and 5). While excluding age and AAO (duration was not found contributory to prediction in Table 5 probably because it is a derivation of the former two) did not affect UPDRS III score (Supplementary Figure 3A; RMSE increase of only 3% –4% for all severities), our data show that the absence of these two variables together with duration increased the prediction error for UPDRS total score for most of the severities (Supplementary Figure 3B) between 26% and 198%. We believe that age, AAO, and disease duration mostly affect UPDRS II score (Activities of Daily Living) and especially UPDRS I score (Mentation, Behavior, and Mood) and to a lesser degree the motor examination in the test (UPDRS III score). Supplementary Table 1 shows ranked by RF values of variable importance (ordered by the Gini index) for UPDRS total and UPDRS III scores, reinforcing that age and duration (much more than AAO) contribute to more accurate prediction of UPDRS total score, whereas age and AAO are less influential compared to lipid species (and duration is not even ranked) for accurate prediction of UPDRS III score (compare Gini index values of these variables between UPDRS total and UPDRS III scores). However, the lowest (under 10) and highest (above 36) severities for UPDRS total score (Supplementary Figure 3B) show the opposite picture. The improvement in the average error while excluding age, AAO, and duration is explained by a skewed disease duration in the lowest range (the 25^th and 48^th percentiles are almost identical) and high variability of this variable in the highest range that together undermine the accuracy of the model that uses all variables.

Strengths and limitations

One strength of our study is the novel machine-learning technologies applied to the data. These offer a solid prediction framework around the RF by: 1) designing a careful, statistically-supported training-validation-test ML methodology; 2) initially applying the classifier to select lipid classes (and demographics) that contributed to the smallest prediction error; 3) repeating this application with only highly ranked classes (and demographics) to select the most contributing lipid species; 4) validating and reinforcing this selection using statistical methods; and 5) profiling gender-age lipid signatures in relation to PD severity, establishing lipid-based markers of disease progression.

A second strength of the study is the extensive unbiased lipidomic analysis we applied using our ML methodology to validate our findings. The cohort is carefully phenotyped, and SNCA, GBA (with the exception of one), and LRRK2 mutation carriers were not included. The main limitation of the study is that the participants did not fast at time of the blood extraction. Fasting status could affect the levels of blood sphingolipids [62]. Therefore, ideally, our findings should be confirmed in larger studies, ideally including fasting participants with longitudinal data points, like in the Parkinson’s Progression Markers Initiative (PPMI) study.

In conclusion, we demonstrated lipid signatures for motor disease severity in PD. The links between PD status, disease severity, and lipidomics require further investigation, but we believe ML adds an extra tool that improves the level of interpretation for large and interconnected lipidomics data.