Data-driven analysis of heterogeneous gait subgroups and ground reaction forces based on integrated center of pressure–center of mass dynamics in poststroke hemiparesis
Plos.org·July 20, 2026
AI Summary
Researchers analyzed gait patterns in stroke survivors with hemiparesis by examining the relationship between ground reaction forces and the dynamics of center of mass and center of pressure. The study identified distinct subgroups of patients based on data-driven analysis, revealing how abnormal force control and balance mechanisms vary across the post-stroke population.
Hemiparetic gait is characterized by abnormal ground reaction forces (GRFs) with impaired control of the center of mass (CoM) relative to the center of pressure (CoP). Although both anteroposterior and mediolateral gait control have been examined separately, how their integrated stance-phase–derived CoP–CoM interactions and local CoP-based features relate to GRF characteristics remains unclear.
This exploratory study aimed to explore the relationships between CoP–CoM and CoP-based dynamics and GRFs and descriptively identify candidate gait subgroups of individuals with poststroke hemiparesis using a data-driven clustering approach.
Seventy-eight community-dwelling individuals with poststroke hemiparesis participated in a three-dimensional gait analysis. Stance-phase–derived CoP–CoM parameters of the transverse plane and local CoP-based loading metrics were extracted during the paretic stance phase. Relationships between these gait metrics and GRFs, particularly early braking force, propulsion, and late braking force, were examined using nonparametric correlation analyses. K-means clustering was performed to explore candidate gait subgroups, and inter-cluster differences were exploratorily examined.
Across all participants, several CoP–CoM interaction metrics showed significant correlations with GRFs. In particular, insufficient forward progression of the CoM relative to the CoP during late stance showed a strong correlation with late braking force (rs = 0.79). Clustering suggested the presence of four candidate gait subgroups characterized by differing combinations of anteroposterior and mediolateral CoP–CoM dynamics and local CoP loading features. However, the results of clusters with small sample sizes warrant cautious interpretation.
Full story reconstructed from Plos.org. Formatting and media may differ from the original.
Integrated stance-phase–derived CoP–CoM dynamics and CoP-based parameters may highlight heterogeneity in hemiparetic gait and may have meaningful associations with GRF characteristics. The candidate gait subgroups should be interpreted as hypothesis-generating and may offer a descriptive framework for understanding diverse gait disorders.
Citation: Mori K, Teramae T, Wakida M, Mano N, Chujo Y, Kuwabara T, et al. (2026) Data-driven analysis of heterogeneous gait subgroups and ground reaction forces based on integrated center of pressure–center of mass dynamics in poststroke hemiparesis. PLoS One 21(7): e0354290. https://doi.org/10.1371/journal.pone.0354290
Editor: Anne E. Martin, Pennsylvania State University Main Campus: The Pennsylvania State University - University Park Campus, UNITED STATES OF AMERICA
Received: February 24, 2026; Accepted: July 6, 2026; Published: July 20, 2026
Data Availability: All relevant data are within the paper and its Supporting information files. The anonymized dataset used in this study is provided as Supporting information (S1 Dataset.csv).
Funding: This work was supported by JSPS KAKENHI (Grant Number JP23K16560) awarded to KM (https://www.jsps.go.jp/), by JSPS KAKENHI (Grant Number JP24K21325) awarded to TN (https://www.jsps.go.jp/), and by the JST Moonshot Research and Development Program (Grant Number JPMJMS2034) awarded to TN (https://www.jst.go.jp/moonshot/en/). Additional support was provided by the Strategic Project for Proofreading and Submission Support of International Academic Papers, Kansai Medical University (https://www.kmu.ac.jp/).
Competing interests: The authors have declared that no competing interests exist.
Hemiparetic gait after stroke is characterized by lower limb motor dysfunction and accompanied by heterogeneous impairments across individuals [1,2], including abnormal gait patterns [1], reduced walking speed [3], and impaired symmetry of ground reaction force (GRF) [4] and spatiotemporal parameters [5,6], resulting in limited walking efficiency and activities of daily living [7]. The GRF imparts acceleration to the center of mass (CoM) in coupled anteroposterior–mediolateral (AP–ML) directions; therefore, paretic limb control of the CoM should be understood as an integrated AP–ML interaction mediated by the GRF. Regarding AP control, the leading limb angle during early stance and the trailing limb angle, as defined by the relative position of the CoM to the center of pressure (CoP) during late stance, are related to the AP component of the GRF; however, the former has been related to early braking force (EBF) [8] and the latter has been related to propulsive force (PF) [9]. Regarding ML control, the foot progression angle (FPA) [3] and lateral foot placement relative to the CoM are associated with walking speed and stability [10,11], and cane use is considered compensation for ML instability [12]. Furthermore, CoP–CoM control in the transverse plane plays a crucial role in stabilizing gait and adapting to neurological deficits [13,14]. Therefore, both AP and ML controls should be understood as a coupled factor rather than as isolated factors.
While the global CoP–CoM relationship captures overall balance and propulsion strategies, local features of the CoP trajectory on the paretic side provide additional insight into the ankle control mechanisms that underlie these dynamics. In particular, difficulty with hindfoot loading during stance and a concomitant reduction in anterior displacement of the CoP are associated with slower walking speed [15] and the occurrence of late braking force (LBF) [16], with lateral CoP deviation further necessitating cane use [12]. Data of healthy individuals and individuals with subacute stroke have shown that AP and ML shifts in the CoP are linked to ankle inversion control and eversion control as well as trunk control [17,18], suggesting the need to integrate local CoP trajectories arising from ankle control with global CoP–CoM dynamics to capture the dynamic characteristics of hemiparetic gait. Moreover, characteristics of hemiparetic gait should be interpreted across the stance phase rather than at a single point, as events at one phase can influence those at subsequent phases; excessive paretic braking forces during early stance often hinder forward CoM progression and reduce backward foot placement during late stance [19], highlighting the interdependence of early and late stance dynamics. Nevertheless, how stance-phase–derived features of AP–ML interactions, which integrate the global trajectory of the CoP–CoM with the local CoP trajectory and loading in the transverse plane, are associated with AP components of the GRF, particularly PF and LBF, is unclear. This gap is clinically important and particularly relevant given that these integrated AP–ML CoP–CoM interactions are likely to vary across individuals depending on their walking ability [19] and compensatory strategies related to the GRFs [16,20].
Therefore, combined stance-phase–derived CoP–CoM and CoP-based parameters may help capture the heterogeneous biomechanical characteristics of hemiparetic gait across individuals. However, because each parameter carries distinct clinical relevance, it remains unclear which combination of parameters most meaningfully relates to GRF characteristics across individuals. Data-driven clustering offers a systematic approach to address this question, and has been increasingly applied to characterize the heterogeneity of post-stroke gait. Prior studies using spatiotemporal, kinematic [21,22], and limb kinematic variables [23] have identified biomechanically distinct subgroups associated with differences in motor impairment severity, compensatory strategies, and walking ability. However, as highlighted in a recent narrative review [24], cluster differentiation has often been dominated by walking speed, and approaches integrating AP–ML CoP–CoM dynamics with GRF-related kinetic features remain unexplored. Exploring candidate gait subgroups based on these integrated dynamics could advance understanding of the diverse GRF characteristics in poststroke hemiparesis and inform targeted rehabilitation strategies.
This study aimed to characterize the heterogeneous biomechanical characteristics of hemiparetic gait by extracting stance-phase–derived CoP–CoM and CoP-based parameters, examining their associations with GRF components from an integrated AP–ML perspective, and applying data-driven clustering to explore candidate gait subgroups. Based on the biomechanical relationships, we hypothesized that clustering would reveal candidate subgroups characterized by differing AP–ML CoP–CoM control patterns, in which insufficient anterior CoP displacement combined with lateral foot placement compensation would contribute to reduced PF and increased EBF, while lateral CoP deviation and decreased local anterior CoP displacement would be associated with increased LBF [16].
We performed a retrospective analysis of data from a database at our institution. This study included 78 community-dwelling individuals with poststroke hemiparesis who were capable of walking independently on a level surface with or without the use of a T-handle cane. The exclusion criteria were the inability to walk independently without assistance from a therapist and the presence of bilateral lesions, lower extremity joint pain, or other neurological or musculoskeletal disorders that could affect walking. We measured Fugl–Meyer Assessment of Lower Extremity synergy (item E, Stages II–IV; maximum score = 22) [25] and sensory scores (item H; lower extremity items only; maximum score = 12) and assessed the daily use of canes among individuals after stroke. Data were collected at three university hospitals as part of routine medical care or other studies between March 1, 2018 and March 1, 2024. This study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the University Human Research Ethics Committee (Kansai Medical University, no. 2022150). Data were first accessed for research purposes on April 1, 2024. Authors did not have access to information that could directly identify individual participants during or after data collection, and all data were anonymized prior to analysis. Informed consent was obtained through an opt-out process.
Gait performance was evaluated using an optical three-dimensional motion capture system (Locus 3D MA-3000; Anima Corp., Tokyo, Japan) and two force plates (MG-1190; Anima Corp.). The system comprised 12 infrared cameras and 1.2-m force plates with sampling rates of 100 Hz and 1000 Hz, respectively, that were synchronized. The 22 reflective markers had a diameter of 12 mm and were placed on the body surface of the acromion, anterosuperior iliac spine, posterosuperior iliac spine, greater trochanter, medial and lateral femoral epicondyles, medial and lateral malleolus, first metatarsal head, fifth metatarsal head, and heel, bilaterally.
GRF data were collected as participants walked at a self-selected speed along a 6-m walkway. To avoid unnatural gait adaptations, participants were instructed to walk naturally without attempting to adjust their steps to target the force plates. The use of a specialized cane (MK-1000A; Anima Corp.) that enables separation of GRF components from the cane and foot when both contact the same force plate was permitted to the minimum extent necessary to prevent falls caused by lateral instability during walking, particularly for individuals with lower walking ability. However, all data used in the present study were measured under the condition without an ankle–foot orthosis.
Marker trajectories and GRFs were filtered using second-order Butterworth low-pass filters at 10 Hz and 20 Hz, respectively. Gait events were identified based on the vertical GRF thresholds as described previously [26]. Specifically, initial contact was defined as the time when the vertical GRF exceeded 5% of body weight (BW), and foot-off was defined as the point when it decreased below 5% BW. The stance phase from the paretic initial contact to foot-off was normalized from 0% to 100%. The whole-body CoM and paretic CoP were calculated using marker locations and force plate signals, respectively. The CoM was estimated using a simplified seven-segment rigid body model including the bilateral feet, shanks, thighs, and trunk, as described by Yamaguchi et al. [27]. The results of at least three paretic stance trials of all individuals except one were averaged; for that one individual, the results of two paretic stance trials were averaged. Walking speed and stance duration relative to the gait cycle were calculated using at least four strides. The AP GRFs, such as the peak and mean values of the EBF during early stance, PF during late stance, and subsequent LBF, were calculated and normalized according to BW (%BW) (Fig 1).
Early braking force (EBF) was defined as the negative AP GRF phase following the heel strike transient until the onset of the propulsive phase. Propulsive force (PF) was defined as the subsequent positive AP GRF phase until the transition to negative values or foot-off. Late braking force (LBF) was defined as the negative AP GRF phase following the propulsive phase until foot-off. Peak values are indicated by circles. Mean values of each phase are shown as dashed lines. Positive values represent the anterior component and negative values indicate the posterior component.
https://doi.org/10.1371/journal.pone.0354290.g001
To ensure consistency during the analysis, the ML coordinates of the markers and CoP were standardized to account for the side of hemiparesis. Thus, the X-axis and Y-axis in the global coordinate system were defined as the anterior and medial directions, respectively, with positive values corresponding to the paretic limb; however, in the local coordinate system, the X-axis and Y-axis represented the anterior and medial directions, respectively, which were defined based on the foot (Fig 2A and 2B). In the local coordinate system, the X-axis was defined as the line extending from the midpoint between the first and fifth metatarsal head markers to the heel marker. The Y-axis (medial direction) was defined as the perpendicular direction to the X-axis in the transverse plane. The origin of this coordinate system was set at the heel marker. This framework was applied to all relevant markers and CoP measurements to allow direct comparisons across participants.
(A) Schematic illustration of the relative position of the CoM and CoP during gait and the ground reaction force (GRF). (B) Trajectories of the CoM and CoP in the global and local coordinate systems in the transverse plane. In the local coordinate system, the foot is divided into three anatomical regions (hindfoot, midfoot, and forefoot), which correspond to the regions used to define hindfoot_CoP_duration, midfoot_CoP_duration, and forefoot_CoP_duration (see Table 1 for definitions). The anteroposterior CoP position at foot-off (AP_CoP_FO), the mediolateral CoP position at foot-off (ML_CoP_FO), and the foot progression angle (FPA) are also indicated. (C) Parameters derived from the trajectory of the CoP relative to the CoM in the transverse plane, including anteroposterior parameters (AP_CoP–CoM_max, AP_CoP–CoM_min, AP_CoP–CoM_FO, and AP_CoP–CoM_min-FO) and mediolateral parameters (ML_CoP–CoM_max, ML_CoP–CoM_max-FO, ML_CoP–CoM_FO, and ML_CoP–CoM_min). All parameters are defined in Table 1.
https://doi.org/10.1371/journal.pone.0354290.g002
Definitions and clinical significance of gait metrics.
https://doi.org/10.1371/journal.pone.0354290.t001
Data processing involved calculating the positions of key markers, trajectories of the CoP, and coordinates of the CoM in the transverse plane during the stance phase, followed by calculating relative values to compare measurements across participants. Key foot markers, including the markers of the heel, lateral malleolus, medial malleolus, fifth metatarsal head, and first metatarsal head, were analyzed. Coordinates of these markers were adjusted relative to the initial heel position and normalized according to the height of the participants. The CoP and CoM trajectories were adjusted using the same normalization method.
The definitions and clinical significance of all 14 gait metrics are summarized in Table 1, and their derivation from the CoP and CoM trajectories is illustrated in Fig 2. The horizontally projected positions of the CoP relative to the CoM in the global coordinate system in the AP (X-axis) and ML (Y-axis) directions were calculated to understand their dynamic relationship during the stance phase (Fig 2C). AP_CoP–CoM_min and AP_CoP–CoM_max capture the posterior and anterior displacement of the CoP relative to the CoM—constructs conceptually analogous to the trailing and leading limb angles [8,9]—and the ML CoP–CoM parameters characterize mediolateral divergence between the CoP and CoM, related to lateral foot placement and pelvis displacement [11]. These relative values were further used to derive specific metrics including AP_CoP–CoM_min-FO, which quantifies the anterior shift of the CoP relative to the CoM from its most posterior point to foot-off and is a novel parameter not previously reported. Additionally, at the end of the stance phase, the anterior position of the CoP along the X-axis relative to the initial heel position in the global coordinate system (AP_CoP_FO) is a newly defined parameter that characterizes the excessive anterior shift of CoP during preswing, which may reflect toe catching caused by premature swing initiation before sufficient unloading of the limb [20,28]. The direction of foot progression was determined by calculating the angle between the central axis of the foot (from the midpoint between the first and fifth metatarsal heads to the heel) and the X-axis, with angles > 0° indicating a toe-out position, consistent with the role of foot yaw in lateral stability [29].
In contrast, for gait characteristics in the local coordinate system, the foot was divided into the following three regions: hindfoot (heel to malleoli); midfoot (malleoli to metatarsal heads); and forefoot (anterior to metatarsal heads). The percentage of the stance time during which the load was applied to each region was calculated as discrete parameters of CoP progression from hindfoot to forefoot, a feature previously associated with walking ability in hemiparetic gait [15]. To quantify ankle inversion or eversion control at foot-off (ML_CoP_FO), the ML position of the CoP relative to the central axis of the foot (Y-axis) was determined by calculating the perpendicular intersection from the heel to the line connecting the first and fifth metatarsal heads, reflecting a construct related to ankle muscle control of the mediolateral CoP position during stance [17].
To reduce redundancy among gait parameters, we applied a step-by-step variable selection process based on pairwise Pearson correlation coefficients. In each functional category, when two or more variables showed high correlation (|r| > 0.80, more conservative than the r > 0.9 criterion used in a comparable clustering study [30]), the variable with a lower correlation coefficient relative to the others was retained. This process was conducted independently for each feature category (e.g., ML CoP–CoM metrics, AP CoP–CoM metrics, and CoP duration metrics). The final set of variables, which included only those with minimal redundancy, was defined prior to clustering and subsequently used for the clustering analysis (Table 1).
The remaining features were standardized using z-score normalization prior to clustering. Then, a K-means clustering analysis was performed to explore candidate gait subgroups among participants. The optimal number of clusters was primarily determined based on the silhouette coefficient [31] and the elbow method was used as a supplementary tool. To address initialization dependency, k-means clustering was performed using k-means++ seeding with 1,000 independent random initializations, selecting the solution with the lowest inertia as the reference. Cluster stability was further assessed by computing the Adjusted Rand Index (ARI) [32] between the reference solution and 100 repeated runs for k = 3, 4, and 5. The candidate clusters were additionally visualized by projecting the clustering features onto the first two principal components to examine cluster separation in the feature space.
To explore the relative contribution of each metric used in the clustering analysis, we applied the mean decrease accuracy method using the randomForest package in R, which is similar to the method described by Kettlety et al. [31]. The importance scores were normalized to the sum of 100%, thus making it easier to compare the relative influence of each metric on the clustering model. Because LBF is defined as the negative AP GRF phase following the propulsive phase, participants in whom no clear propulsive phase was identified were excluded from all LBF-related analyses. Correlations between GRF data and gait characteristics across all participants were exploratorily analyzed. Data normality was assessed using the Shapiro–Wilk test for all variables included in the correlation analyses. Because GRFs metrics were not normally distributed, Spearman’s rank correlation was used throughout.
Differences in GRF data and gait characteristics among the clusters were determined after the clustering analysis results were analyzed. Because cluster C comprised only four participants, the Kruskal–Wallis rank-sum test was used to determine the presence of significant differences in the variables across groups. To quantify the effect size, the epsilon-squared value was calculated as a measure of association strength. Multiple comparisons among clusters were performed using post hoc analyses with the Steel–Dwass test to determine which clusters significantly differed from each other. All analyses were performed using Python 3.10.12 and R software version 4.4.1.
Seven individuals without PF data (two in cluster B, one in cluster C, and four in cluster D) were excluded from the LBF analysis. A step-by-step variable selection process was conducted to reduce multicollinearity among the 14 predefined gait characteristic indices. As a result, the following variables were excluded because of high correlation (|r| > 0.80) with other metrics: ML_CoP–CoM_min; ML_CoP–CoM_FO (correlated with ML_CoP–CoM_max); AP_CoP–CoM_FO (correlated with AP_CoP–CoM_min); and midfoot_CoP_duration (correlated with forefoot_CoP_duration). Consequently, 10 gait characteristics were used for clustering analyses (Table 1).
Peak and mean EBFs showed the highest correlation with AP_CoP–CoM_max (rs = –0.75 and rs = –0.77). Peak and mean PFs had the strongest correlation with AP_CoP–CoM_min (rs = –0.91 and rs = –0.91) and were significantly correlated with ML_CoP_FO, ML_CoP–CoM_max-FO, and hindfoot_CoP_duration. Peak and mean LBFs were strongly correlated with AP_CoP–CoM_min-FO (rs = 0.79 and rs = 0.78) and AP_CoP_FO (rs = –0.58 and rs = –0.60). Correlations between the other metrics and GRFs are presented in S1 Fig.
The clustering results assigned 31, 17, 4, and 26 participants in clusters A, B, C, and D, respectively (Table 2). The silhouette coefficients were 0.163, 0.178, and 0.173 for k = 3, 4, and 5, respectively, with k = 4 yielding the highest value (S2A Fig). The elbow plot showed a gradual and continuous decline in the sum of squared errors without a distinct inflection point (S2B Fig). Mean ARI values were 0.842 (SD = 0.108), 0.723 (SD = 0.142), and 0.778 (SD = 0.129) for k = 3, 4, and 5, respectively, indicating moderate-to-good stability across initializations (S2C Fig). PCA-based visualization suggested meaningful separation in both the k = 4 and k = 5 solutions (Fig 3), supporting the selection of k = 4 based on the highest silhouette coefficient (S2A Fig). The k = 5 solution subdivided Cluster A into two subgroups, whereas Cluster C showed complete membership overlap across k = 4 and k = 5, further supporting its consistency as a candidate gait subgroup (S2D Fig). The distribution of 10 clustering parameters, including their mean decrease accuracy, and GRF components across clusters are shown in Fig 4, with full numerical summaries provided in the S1 Table for k = 4, S2 Table for k = 5.