The results of the clinical evaluations are summarized in Table 2. The UPDRS motor symptoms were significantly higher in the PDsevere group (n = 79) compared to the PDmild group (n = 48) for the total score (PDmild: 17.7 ± 9.9, PDsevere: 30.1 ± 12.4; W = 816, p < 0.01), rigidity (PDmild: 4.4 ± 3.2, PDsevere: 6.5 ± 3.5; W = 1231, p < 0.01), PIGD (PDmild: 3.5 ± 2.5, PDsevere: 6.8 ± 2.7; W = 687, p < 0.01), and bradykinesia (PDmild: 7.4 ± 5.6, PDsevere: 12.5 ± 6.0; W = 953.5, p < 0.01). There was no significant between-group difference in tremor score (PDmild: 1.2 ± 1.5, PDsevere: 2.2 ± 2.9; W = 1588, p = 0.11). The patient's age and the duration of disease were nearly the same in the PDmild and PDsevere groups (age: PDmild: 70.8 ± 8.4, PDsevere: 71.4 ± 8.1, W = 1784, p = 0.58; duration of disease: PDmild: 6.1 ± 3.8, PDsevere: 7.0 ± 5.1, W = 1776, p = 0.63).

Table 2 Characteristics of the patients with PD and the healthy controlsExploratory factor analysis

To extract distinct components of the postural instabilities in the PD patients, we performed an exploratory factor analysis (EFA) with 23 variables calculated from the patients' CoP time-series data. Four variables and a single variable were removed in the correlation and KMO criteria, respectively (MeanPos_AP, RMS, AAMV, MD3, and MV_CV). In addition, two variables were removed in the factor loading criteria (MeanPos_ML and MT3). A total of five factors were extracted by the EFA with 23 accepted variables (Fig. 1).

Figure 1's middle panel depicts the loading score of each variable on the detected factors. Each factor can be interpreted as follows: F1 was interpreted as a "magnitude of sway factor," with large contributions such as Area, MV, and Ds in the SDA. F2 was interpreted as a "ML frequency component factor," with large contributions from MFreq_ML and SLPL_ML. In contrast, F3 was interpreted as an "AP frequency component factor" with large contributions from MFreq_AP and SLPL_AP. F4 was interpreted as a "high-frequency component factor" with large contributions from SLPH_AP and ML and FD. F5 was interpreted as a "closed-loop control factor," with large contributions from Dl and CriT in SDA.

The relationship between factor scores and UPDRS scores and the duration of disease are shown as a color matrix in Fig. 2A. F1 showed a significant positive correlation with the disease duration (ρ = 0.29, p < 0.01) but a significant negative correlation with UPDRS rigidity (ρ =  − 0.25, p < 0.01). F4 was positively correlated with the UPDRS total score (ρ = 0.31, p < 0.01), rigidity (ρ = 0.30, p < 0.01), bradykinesia (ρ = 0.20, p = 0.02), and PIGD (ρ = 0.31, p < 0.01). F2, F3, and F5 showed no significant correlations.

Fig. 2

Clinical characteristics. A Color matrix: indicating the correlation coefficient between factor scores and UPDRS scores and the duration of disease. *p < 0.05, **p < 0.01 by Spearman's rank correlation coefficient. Plot graphs: the results of the comparison of clinical evaluations in each cluster. Data are mean ± SD. *p < 0.05, **p < 0.01 by post hoc test (Steel–Dwass test) and Kruskal–Wallis test. B Percentages of members. Left panel: The percentages of Elderly, PDmild and PDsevere in each cluster. Percentages of the clusters in the clinical classification. Right panel: The percentages of clinical classification in each cluster

Gaussian mixture model-based cluster analysis

Since the comparison of factor scores confirmed that the Young healthy group had postural control characteristics that were distinctly different from those of the Elderly healthy and PD groups, we performed GMM-based clustering on the data from the healthy Elderly and PD patients (n = 174). The GMM-based cluster analysis classified the six clusters.

As shown in Fig. 2B, 31 patients (17.8%) were classified as Cluster 1, 32 patients (18.4%) as Cluster 2, 41 patients (23.6%) as Cluster 3, 39 patients (22.4%) as Cluster 4, 10 patients (5.7%) as Cluster 5, and 21 patients (12.1%) as Cluster 6. The percentage of Elderly and PD patients in the clusters are shown in the right panel of Fig. 2B. Clusters 4, 5, and 6 were considered to be PD-specific clusters due to the small number of Elderly belonging to them. There was no significant difference in the patient's age of each cluster (χ2 = 7.988, df = 5, p = 0.16). The duration of disease was significantly longer in Cluster 4 and Cluster 6 than Cluster 1 (χ2 = 13.353, df = 5, p = 0.02, Fig. 2B, lower right panel).

The UPDRS results for each cluster are shown in Fig. 2A right panel. The total score was significantly lower in Cluster 6 than Cluster 3 (χ2 = 11.37, df = 5, p = 0.04). Regarding the sub-scores, rigidity was significantly higher in Cluster 3 than Clusters 2 and 6, and the rigidity in Cluster 6 was significantly lower than that in Cluster 4 (χ2 = 22.384, df = 5, p < 0.01). There was no significant between-cluster difference in bradykinesia (χ2 = 8.587, df = 5, p = 0.13), PIGD (χ2 = 6.378, df = 5, p = 0.27), or tremor (χ2 = 4.249, df = 5, p = 0.51).

The results of a detailed comparison of the factor scores and the variables representing factors in each clinical classification (A) and cluster-based categorization (B) are summarized in Fig. 3. The radar chart at the figure's upper right shows the results of the comparison of the variables with the largest contribution of each factor in each group. In the case of clinical classification, the factor scores are largely overlapped among groups. F1 and F4 showed significantly higher values in Elderly, PDmild, and PDsevere compared to Young (F1: χ2 = 45.913, df = 3, p < 0.01, F4: χ2 = 43.685, df = 3, p < 0.01). The values were not significantly different in F2 (χ2 = 5.224, df = 3, p = 0.16), F3 (χ2 = 3.153, df = 3, p = 0.37), and F5 (χ2 = 6.908, df = 3, p = 0.07).

Fig. 3

Summary of the results of the exploratory factor analysis. A Comparison of factor scores by clinical classification. Data are mean ± SD. *p < 0.05, **p < 0.01 by post-hoc test (Steel–Dwass test) and Kruskal–Wallis test. Scatterplots: the FA in each cluster with 95% confidence ellipsoids and ellipses. Upper right radar chart: the z-scores of the main variables for each factor. The clinical classification shows that each variable has a large overlap in PD severity. B Comparison of factor scores by clusters (as in the right part of Fig. 1). Upper right radar chart: the z-scores of the main variables for each factor. For the categories by clusters, the different shapes of the radar chart show the independence and characteristics of each cluster

In contrast to the clinical classification, the plot of each cluster distributed independently. F1 was significantly lower for Clusters 1 and 3 than for Clusters 4 and 5 (χ2 = 85.622, df = 5, p < 0.01). Cluster 1 was also significantly lower than Cluster 2. In contrast, cluster 6 was significantly higher than all clusters excluding cluster 5. In F2, Cluster 5 was significantly higher than the other clusters, and Cluster 3 was significantly higher than Cluster 2 (χ2 = 40.304, df = 5, p < 0.01). In F3, Clusters 3 and 5 were significantly higher than Clusters 1, 2, 4, and 6 (χ2 = 73.777, df = 5, p < 0.01). In F4, Cluster 5 was significantly higher than the other clusters, and Clusters 3 and 4 were significantly higher than Clusters 1, 2, and 6 (χ2 = 47.148, df = 5, p < 0.01). In F5, Clusters 2 and 6 were significantly lower than the other clusters, and Cluster 4 was significantly higher than Clusters 1 and 3 (χ2 = 103.9, df = 5, p < 0.01). The detailed statistical results are provided in Table 2.

Comparison of CoP parameters

All evaluated parameters were shown in both the clinical classification (Fig. 4, left panel) and the cluster-based categorization (Fig. 4, right panel). Although the data among the four subject groups (Young, Elderly, PDmild, and PDsevere) widely overlapped, the data in the cluster categorization showed clear separation, which we speculate was due to different contributions of elements comprising the postural instability. Figure 4A shows the data plots of CoP area and mean velocity (MV), which are general spatial and temporal evaluation parameters, respectively. Both the CoP area and MV were significantly higher in the Elderly, PDmild, and PDsevere groups was compared to the Young group (CoP area: χ2 = 46.937, df = 3, p < 0.01, MV: χ2 = 77.058, df = 3, p < 0.01). PDsevere was also significantly higher than Elderly (p = 0.04). There was no significant difference in either the CoP area or the MV between PDmild and PDsevere groups. In the cluster categorization, the CoP area was significantly lower in Clusters 1 and 3 than in other clusters and significantly higher in Clusters 5 and 6 versus the other clusters (χ2 = 80.672, df = 5, p < 0.01). The MV was significantly lower in Cluster 1 compared to the other clusters and significantly higher in Cluster 5 than in the other clusters. Cluster 2 was significantly lower than Cluster 4. Cluster 6 was significantly higher than Clusters 1, 2, and 3 (χ2 = 85.183, df = 5, p < 0.01).

Fig. 4

Comparison of CoP parameters by clinical classification and cluster categorization. Left panel: Comparison by clinical classification. Right panel: Comparison of factor scores by clusters. All data are mean ± standard deviation (SD). *p < 0.05, **p < 0.01 by post-hoc test (Steel–Dwass test) in the Kruskal–Wallis test. A Comparison of spatiotemporal parameters (Area and MV). B Comparison of power-spectrum parameters by anterior–posterior direction (SLPL_AP and SLPH_AP). Each left part: PSD log plot from 0.15 to 10 Hz. The blue line is the regression line for 0.15–1 Hz, and the red line is the regression line for 1–5 Hz. C Comparison of results of the stabilogram diffusion analysis (SDA). Each left part: Linear stabilogram diffusion plots in CoP displacements in the planar direction

The logarithmic plots of power spectral density (PSD) by the anterior–posterior direction for each group are demonstrated in Fig. 4B. The clinical classification showed no significant differences in power slopes in the low-frequency band down to 1 Hz (SLPL_AP: χ2 = 7.806, df = 3, p = 0.05). In contrast, the slope of power in the high-frequency band above 1 Hz was significantly higher in the PD group (SLPH_AP: χ2 = 34.484, df = 3, p < 0.01). In the cluster classification, SLPL_AP was significantly higher in Clusters 3 and 5 than in the other clusters (χ2 = 66.017, df = 5, p < 0.01). SLPH_AP was significantly higher in Cluster 4 than in clusters 2, 3, and 6 (χ2 = 20.194, df = 5, p < 0.01).

The SDA results for each group and each cluster are summarized in Fig. 4C. As shown by the stabilogram diffusion plot in the figure's left panel, each group and cluster showed clearly different characteristics in both short- and long-term regions. In the clinical classification, Ds, which indicates the expansion of the short-time region, was significantly higher in the Elderly and PD groups than in the Young group (χ2 = 46.308, df = 3, p < 0.01). CriD was also significantly higher in the Elderly and PD groups compared to the Young group (χ2 = 43.289, df = 3, p < 0.01). In addition, Dl, which indicates the long-time region, was significantly higher in the PDsevere group versus the Young group (χ2 = 13.41, df = 3, p < 0.01). There was no significant difference in Crit (χ2 = 1.852, df = 3, p = 0.60).

In the cluster categorization, Ds was significantly higher in Clusters 5 and 6 than in other clusters, and significantly lower in Cluster 1 than in the other clusters (χ2 = 71.736, df = 5, p < 0.01). Dl was significantly higher in Cluster 5 than in Clusters 1, 2, and 3 (χ2 = 28.359, df = 5, p < 0.01). CriT was significantly lower in Clusters 3, 4, and 5 versus the other clusters, and significantly higher in Cluster 6 than in Cluster 1 (χ2 = 74.211, df = 5, p < 0.01). CriD was significantly higher in Cluster 6 than in the other clusters, and significantly higher in Clusters 2 and 4 than in Clusters 1, 3, and 5 (χ2 = 101.03, df = 5, p < 0.01). Cluster 1 was also significantly lower than Cluster 3.

Comparison of four representative cases

To validate the interpretation of each cluster, we carried out the characteristics of four representative patients (Fig. 5). Each patient had undergone electromyography (EMG) of the ankle joint (tibialis anterior and gastrocnemius medialis) measured simultaneously with the posturography. Although all four patients were at HY stage 3, they belonged to different clusters and showed different postural control characteristics.

Fig. 5

Representative cases. Left panel: The CoP Lissajous figure of CoP, the EMG activity in the gastrocnemius (blue) and tibialis anterior muscles (red), and the power spectral density (PSD) log plot in four representative cases are shown. The blue line in the CoP Lissajous figures shows the trajectory of the CoP, and the green line shows the 95% confidence ellipse. The PSD log plot showed a range from 0.15 to10 Hz. The blue line is the regression line for 0.15–1 Hz, and the red line is the regression line for 1–5 Hz. Right panel: Radar chart of key variables for each factor in four representative cases. The cases are shown as individual Z-scores and Young as mean Z-scores

Patient #1 (Cluster 1)

The patient was a 61-year-old man who had been diagnosed for 3 years. His UPDRS motor score was 35 points, with high scores mainly for tremor, rigidity, and bradykinesia. His postural alignment showed a dropped head, freezing of gait, and easy falling backward in daily life. The postural sway characteristics showed a narrow sway range (F1) and low-frequency characteristics (F2, F3, F4). The tibialis anterior muscle was dominant over the medial gastrocnemius muscle in muscle activity. In this patient's case, the postural sway was smaller than that of healthy young adults, but this may have been due to compensatory postural control for the easy falling backward.

Patient #2 (Cluster 4)

The patient was a 75-year-old man who had been diagnosed for 3 years. His UPDRS motor score was 49 points, which was high mainly in the scores of rigidity, bradykinesia, and PIGD. He showed decreased flexibility mainly in his trunk. The postural sway characteristics showed a rather wide sway range (F1) but low-frequency characteristics (F2, F3, F4). F5 was high compared to other patients' cases. This patient's muscle activity showed phasic activity of medial gastrocnemius and tibialis anterior muscles with an anterior–posterior CoP shift. Gradual swaying in the low-frequency band was increased due to decreased postural localization by the closed-loop control (F5).

Patient #3 (Cluster 5)

The patient was a 63-year-old man who had been diagnosed for 9 years. The UPDRS motor score was 29 points, which was high mainly in rigidity, bradykinesia, and PIGD. In daily life, the patient showed easy falling due to freezing of gait. The postural sway characteristics showed a narrow range of sway range (F1) as in Patient #1 but high-frequency characteristics (especially F4). The muscle activity showed co-contraction of the medial gastrocnemius and tibialis anterior muscles. In this patient, the postural sway was narrowed due to the high ankle joint stiffness caused by the co-contraction (shown in F4), suggesting that the postural control was continuous rather than intermittent [40].

Patient #4 (Cluster 6)

The patient was a 78-year-old man who had been diagnosed for 16 years. His UPDRS motor score was 41 points, and it was high mainly in bradykinesia and PIGD scores. In addition, dyskinesia was observed in the on-medication state. In daily life, the patient was highly active by walking outdoors and engaging in hobby activities, but he experienced frequent falls. The postural sway characteristics showed a wide sway range (F1), high-frequency characteristics (F2, F3, F4), and low F5. The muscle activity showed phasic activity of medial gastrocnemius and tibialis anterior muscles with an anterior–posterior CoP shift. This patient's dyskinesia caused fast and large anterior–posterior swaying, which may have indicated a time delay in the functioning of the closed-loop control.