Clustering analysis across different speeds reveals two distinct running techniques with no differences in running economy

ABSTRACT

Establishing the links between running technique and economy remains elusive due to high inter-individual variability. Clustering runners by technique may enable tailored training recommendations, yet it is unclear if different techniques are equally economical and whether clusters are speed-dependent. This study aimed to identify clusters of runners based on technique and to compare cluster kinematics and running economy. Additionally, we examined the agreement of clustering partitions of the same runners at different speeds. Trunk and lower-body kinematics were captured from 84 trained runners at different speeds on a treadmill. We used Principal Component Analysis for dimensionality reduction and agglomerative hierarchical clustering to identify groups of runners with a similar technique, and we evaluated cluster agreement across speeds. Clustering runners at different speeds independently produced different partitions, suggesting single speed clustering can fail to capture the full speed profile of a runner. The two clusters identified using data from the whole range of speeds showed differences in pelvis tilt and duty factor. In agreement with self-optimisation theories, there were no differences in running economy, and no differences in participants’ characteristics between clusters. Considering inter-individual technique variability may enhance the efficacy of training designs as opposed to ‘one size fits all’ approaches.

KEYWORDS:

• Running performance
• energetics
• kinematics
• unsupervised learning
• machine learning

Introduction

Long-distance running is a widely practiced sport with a growing number of runners participating in competitive events (García-Pinillos et al., 2017). As a result, there is increasing interest in the science underpinning running to enhance performance while reducing running related injuries. For biomechanists and coaches, understanding the relationship between running technique and economy, and how running technique can be improved through training have been central research questions (Folland et al., 2017; Moore, 2016; Tartaruga et al., 2012; van Hooren et al., 2020; van Hooren, Jukic, et al., 2024). These issues are also receiving considerable attention from running consumer technology developers who have found in runners a community keen to include wearable sensors within their sports kit, to track and improve their technique and training (Aroganam et al., 2019; Scheerder et al., 2015).

At first sight, the relationship between technique, running related injuries and running performance seems straightforward from a mechanical perspective: a poor running technique can place excessive strain on the musculoskeletal tissue that may result in injury (Bramah et al., 2018); a more economical running technique (i.e., lower metabolic cost) can allow an athlete to run faster while spending the same amount of energy (Kipp et al., 2019). However, current evidence-based relationships between running technique and injury and economy, respectively, are mostly inconsistent, contrasting or limited at best (Ceyssens et al., 2019; Moore, 2016; van Hooren et al., 2020; Vannatta et al., 2020). Indeed, a recent extensive review with meta-analysis on running technique and economy (van Hooren, Jukic, et al., 2024) identified only consistent trends for a weak association between higher step frequency and better economy, and moderate associations between lower vertical displacement of the centre of mass, higher vertical and leg stiffness and better economy. Inter-individual variability in running technique emerges as a possible contributor to these discrepancies, since runners are thought to self-optimise their movements with practice (Cavanagh & Williams, 1982; Moore et al., 2019; Williams & Cavanagh, 1987), making running technique highly individual. Building upon the self-optimisation principle, the preferred movement path paradigm was proposed in the study of running footwear (Nigg, 2010). According to this theory, runners exhibit a preferred intersegmental coordination pattern that is assumed to be associated with minimal energetic demand and most comfort, and that athletes are reluctant to depart from despite changes in e.g., footwear (Nigg et al., 2017). In related work, Hoerzer et al. (2015) showed how clustering methods could be used to identify functional groups or clusters of runners with a similar technique or style, and how such groups were able to better explain athletes’ response to different footwear compared to generic common groupings such as by age or sex.

The definition of functional groups of runners who share a similar technique has traditionally been based on expert knowledge. For instance, the dual-axis framework has been recently introduced to characterise the full spectrum of running techniques based on stride frequency and duty factor, proposing five main styles: stick (long contact), bounce (long flight), hop (small steps), push (long steps) and sit (somewhere in the middle) (van Oeveren et al., 2021). The definition of these running techniques was based on an extensive review of centre of mass displacement and spatiotemporal running parameters. However, the popularisation of machine learning methods in biomechanics is leading researchers to increasingly adopt data-driven methods such as clustering, which can minimise the introduction of researcher bias in the definition of running functional groups (Phinyomark et al., 2015; Watari et al., 2018). Clustering is a family of unsupervised learning algorithms whose objective is to group datapoints with similar behaviour within a dataset, uncovering potentially a priori unknown patterns within the data. Clustering methods have been embraced in the study of running, with a special focus on injury biomechanics (see Table 1). For instance, Phinyomark et al. (2015) clustered healthy runners based on their running technique and then compared each healthy cluster with a group of runners with patellofemoral pain. Each comparison yielded different conclusions about which aspects of running technique were associated with healthy or injured status, highlighting the importance of considering the existence of different running techniques when comparing healthy and injured individuals. Jauhiainen et al. (2020) clustered runners with different injuries expecting the runners within each cluster would present the same injury. Contrary to their initial hypothesis, runners with the same injury did not necessarily cluster together, stressing again the importance of considering different running techniques when planning injury prevention or rehabilitation strategies.

Table 1. Running clustering studies.

Study	Sample	Health status	Variables	Dimensionality reduction	Clustering algorithm	Clusters (n)
Aljohani and Kipp (2020)	32 (17 F) trained	Healthy	0D kinematics, kinetics, EMG and coordination variability	-	SOM	4
Phinyomark et al. (2015)	121 (81 F) trained	Healthy	1D Lower limb joint kinematics	Variable by variable PCA (80%)	AHC	2
Jauhiainen et al. (2020)	291 (146 F) trained and highly trained	Healthy (25) Injured	0D spatio-temporal Segment and joint kinematics	PCA (80%)	AHC	5
Senevirathna et al. (2023)	134 (79 F)	Healthy (58) Injured	Spatio-temporal kinematics and 0D kinetics	PCA (80%)	AHC	2
Hoerzer et al. (2015)	88 (43 F) unspecified	Healthy	Lower limb marker trajectories	-	SOM-k means	8
Watari et al. (2018)	110 (46 F) recreational	Injured	Position derived pelvis acceleration	-	AHC	2
Watari et al. (2021)	85 (42 F) trained	Healthy	1D Ankle/foot kinematics	PCA (99%)	AHC	2
Martin et al. (2022)	53 (32 F) highly trained	Healthy	1D joint kinematics and kinetics, and GRF	PCA (90%)	AHC	5
Forrester et al. (2015)	102 (35 F) recreational	Healthy	Foot strike angle		Regression mixture model	3
Floría et al. (2024)	92 M recreational	Healthy	Stride frequency Duty factor	-	K means	3

F: females; M: males; PCA: Principal Component Analysis; AHC: Agglomerative Hierarchical Clustering; SOM: Self-Organising Maps.

The increasing number of clustering studies (Table 1) show potential to further understand the relation of biomechanics with running economy and performance. The current literature appears to be drifting away from the idea that a single most economical technical model that every runner should pursue exists (Moore, 2016; van Hooren et al., 2020). According to the self-optimisation and preferred movement path theories, clusters of runners with different but similarly economical techniques may exist, but this hypothesis has not been experimentally tested. Also, the association of running technique clusters with physiological factors and anthropometrics has not been explored, despite their importance to running economy and performance (Barnes & Kilding, 2015).

It is well established that changes to technique and training can result in greater economy in some runners while it may not work for others (Barnes & Kilding, 2015; Moran & Wager, 2020). Indeed, Watari et al. (2021) showed how clusters of runners undergoing the same foot core exercise programme reacted differently to the intervention. Cluster-specific training designs may better target individual technical and strength and conditioning needs that can help enhance the effectiveness of the interventions. However, a critical gap persists in the clustering literature prior to developing cluster-specific training: studies have predominantly focused on singular running speeds, disregarding the diverse pacing strategies and racing distances that runners encounter. How well clusters found at one speed generalise to the full spectrum of speeds an athlete runs at remains unanswered. Addressing this issue appears essential before proceeding to develop and test cluster-specific training interventions.

Therefore, the main purpose of this study was to identify clusters of runners who share a similar running technique and compare the clusters for their running economy. It was hypothesised that clusters would not exhibit differences in their economy, in line with the self-optimisation and preferred movement path theories. Additionally, a comprehensive list of demographic, anthropometric, physiological capacity, training and performance characteristics of the clusters were compared. A secondary aim of this study was to compare the clustering partitions of the same runners at three different speeds as well as at the three speeds combined to assess how representative single speed clustering may be of the full spectrum of speeds a runner uses.

Materials and methods

One-hundred runners (50 females) were recruited for this study. Participants were aged 18–50 years, free from neuromusculoskeletal injury in the 6 months preceding the test and with a most recent self-reported race (or equivalent) 10K time under 57:20 (minutes:seconds) for females and 50:00  for males to ensure matched age-grading (USATF, 2020). These inclusion criteria were selected to prevent sex, age, injury status and athletic level from dominating and driving the clustering, since these factors are known to influence running technique (Boyer et al., 2017; Clermont et al., 2017; Willy & Paquette, 2019; Xie et al., 2022). All participants provided written informed consent before taking part in the study and ethical approval was granted by the Research Ethics Approval Committee for Health (REACH) of the University of Bath.

Participants completed a maximal incremental running test on a treadmill (J100, Powerjog Inc., North Charleston, South Carolina, US) at a 1% gradient (Jones & Doust, 1996). The test included four-minute stages of constant speed running to ensure valid blood lactate concentration and maximal respiratory variable measurements (Bentley et al., 2007). The test started at 10 km/h and speed was increased by 1 km/h after every completed stage until volitional exhaustion. Participants were provided with one-minute rest between stages to allow lactate concentration collection and were given frequent verbal encouragement throughout the test to obtain participants’ maximum effort (Andreacci et al., 2010; Midgley et al., 2017; van Hooren, van Der Lee, et al., 2024). Verbal encouragement was systematically provided after every completed minute and 30 seconds before the end of the initial stages using standardised cues (e.g., ‘That is 2 minutes gone, 2 more to go, you are doing really well, keep going’) and a 5-second countdown was provided at the end of each stage. Frequency was increased to every 30 seconds based on higher rate of perceived exertion reported by the runner and became nearly constant (every 5–10 seconds) when runners started exhibiting symptoms of being close to maximal exertion (e.g., ‘Come on’, ‘keep pushing’, ‘you are doing really well’, ‘keep going for as long as you can’, ‘you can do it’, ‘you can finish this one’). To maximise performance and standardise the testing conditions, participants were asked to refrain from strenuous exercise, alcohol and caffeine the 48 hours preceding the test and to attend the test at least 2 hours postprandial. Participants completed the test using their preferred running shoes to prevent influencing their habitual running technique. All the testing sessions took place in the Applied Biomechanics Suite at the University of Bath, where air flow, temperature and humidity were kept consistent through mechanical ventilation.

Physiological data collection and processing

A lactate analyser (Lactate Plus, Nova Biomedical, USA) was used to record the lactate concentration in blood (~0.7 μl samples) after every completed stage. The device was calibrated weekly with known concentrations following the manufacturer's standards. These data were used to calculate individuals’ lactate threshold by fitting an exponential curve to the data (Santos-Concejero et al., 2014) and using the D_max method (Cheng et al., 1992). Stages under lactate threshold can be assumed to be completed in metabolically steady state conditions (Faude et al., 2009), ensuring valid estimations of metabolic energy expenditure using respiratory variables. During the treadmill test, a portable metabolic analyser (K5, Cosmed, Italy) was used to capture breath-by-breath gas exchange. This device in breath-by-breath mode has been found to obtain excellent reliability against a metabolic simulator (Guidetti et al., 2018). The device was calibrated as per manufacturer’s recommendation prior to every test. The mask was fitted to the participant 5 minutes before commencing the test to allow participants to habituate to breathing wearing the mask. Breath-by-breath data were inspected in 30-second windows to discard aberrant breaths (±2 standard deviations away from the local mean) due to e.g., talking, coughing, swallowing. For every stage completed under lactate threshold, the measured $\dot{V}{O}_2$ and $\dot{V}\mathit{C}{O}_2$ were used to estimate aerobic energy expenditure using the Jeukendrup and Wallis equation for substrate oxidation during exercise (Jeukendrup & Wallis, 2005). The average value over the last minute of each sub-lactate threshold stage was calculated, normalised by body mass and expressed as the aerobic energy expenditure per minute per kg (kcal/min/kg) i.e., running economy (Beck et al., 2018). Additionally, $\dot{V}{O}_2$ data were smoothed using a 15-breath rolling average and the highest smoothed $\dot{V}{O}_2$ during the test was identified and normalised to body mass ($\dot{V}{O}_2\mathit{peak}$) (Robergs et al., 2010).

Biomechanical data collection and processing

A 16-camera motion capture system (Oqus, Qualisys, Sweden) was used to record the three-dimensional trajectories of retroreflective spherical markers attached to the participant (see Appendix) at a sampling frequency of 200 Hz. Marker placement was performed by the same researcher for every participant to prevent systematic differences in segment and joint kinematics due to inter-tester marker placement variability (Pohl et al., 2010). Marker trajectories were low-pass filtered (Butterworth, fourth order, zero-lag) with a cut-off frequency of 10 Hz (Fukuchi et al., 2017) and used to estimate full-body kinematics using a six degrees-of-freedom modelling approach (Visual 3D, C-Motion, C-Motion, Inc, Rockville, MD). The first and last 10 seconds of the trial were discarded to avoid capturing unnatural gait when getting on and off the treadmill. Foot-strike and toe-off were detected based on segment and joint kinematics (Rivadulla et al., 2021). Trials were segmented into strides from foot-strike to foot-strike and strides were time-registered to 201 data points (average stride duration was 141 ± 7 frames). Only right-leg strides were considered. Running technique was characterised by the average pattern of six time-series (1D) variables: the vertical displacement of the centre of mass normalised to leg length (vCOM), trunk-to-pelvis, pelvis, hip, knee, and ankle angles in the sagittal plane; and two discrete (0D) variables: stride frequency normalised to leg length and duty factor. The choice of these variables to characterise running technique was motivated by the fact that during running, most metabolic energy is spent to support body weight and generate forward propulsion (80%) and leg swing (7%) (Arellano & Kram, 2014). Movements in the frontal and transverse plane may be less relevant for running technique and economy (van Hooren, Lennartz, et al., 2024) and are more sensitive to marker placement errors and soft tissue artefact (Osis et al., 2016; Reinschmidt et al., 1997). Therefore, they were not considered in this study.

Data analysis

To maximise sample size included in the analysis we assumed the effects of fatigue on running technique were minimal up to the first stage completed past lactate threshold speed. Out of the 100 participants completing the test, 84 runners (Table 2) were able to complete speeds 10, 11, 12 and 13 km/h in these conditions and were therefore considered. The first stage (10 km/h) was excluded from the analysis as participants were familiarising themselves with running with the testing kit on. Therefore, only data from 84 runners at 11, 12 and 13 km/h were included in the clustering analysis.

Table 2. Runners’ descriptive statistics.

Variable	All		C0 Neutral pelvis	C1 Tilted pelvis	C0-C1 comparison
	F	M
Sample (n)	84		49	35
Demographics
Sex (n, %)	34 (40.5%)	50 (59.5%)	21 F (42.9%)	13 F (37.1%)	Χ^2 = 0.09, p = 0.76
Age (years, $\bar{x}(\mathit{s})$)	33 (10)	35 (10)	35 (10)	33 (9)	U = 973.5, p = 0.29
Anthropometrics
Height (m, $\bar{x}(\mathit{s})$)	1.66 (0.05)	1.79 (0.06)	1.73 (0.08)	1.74 (0.09)	t = −0.37, p = 0.71
Body mass (kg, $\bar{x}(\mathit{s})$)	57 (6.2)	73 (7)	65.7 (10.2)	67.6 (10.7)	t = −0.81, p = 0.43
Trunk length (m, $\bar{x}(\mathit{s})$)	0.73 (0.04)	0.81 (0.04)	0.78 (0.06)	0.78 (0.06)	t = −0.24, p = 0.81
Leg length (m,$\bar{x}(\mathit{s})$)	0.93 (0.04)	0.98 (0.04)	0.96 (0.05)	0.96 (0.05)	t = −0.38, p = 0.71
Pelvis width (m, $\bar{x}(\mathit{s})$)	0.16 (0.01)	0.17 (0.01)	0.17 (0.01)	0.17 (0.01)	t = 0.36, p = 0.72
Thigh length (m, $\bar{x}(\mathit{s})$)	0.42 (0.02)	0.44 (0.02)	0.43 (0.02)	0.43 (0.03)	t = 0.13, p = 0.89
Shank length (m, $\bar{x}(\mathit{s})$)	0.38 (0.02)	0.42 (0.02)	0.41 (0.03)	0.40 (0.03)	t = 0.41, p = 0.68
Foot length (m, $\bar{x}(\mathit{s})$)	0.26 (0.01)	0.29 (0.01)	0.27 (0.02)	0.28 (0.02)	t = −0.26, p = 0.80
Physiology
$\dot{V}{O}_2\mathit{peak}$ (ml/min/kg, $\bar{x}(\mathit{s})$)	51.1 (5.2)	52.3 (7)	51.9 (6.4)	51.8 (6.4)	t = 0.06, p = 0.95
LT (km/h, $\bar{x}(\mathit{s})$)	13.7 (0.9)	14.4 (1)	14.1 (0.9)	14.1 (1.1)	t = 0.06, p = 0.95
Running economy
EE 11 km/h (kcal/min/kg, $\bar{x}(\mathit{s})$)	0.19 (0.02)	0.18 (0.02)	0.19 (0.02)	0.19 (0.02)	C: F(1, 69) = 0.09, p = 0.77 S: F(2, 138) = 628.55, p < 0.001 CXS: F(2, 138) = 0.69, p = 0.50
EE 12 km/h (kcal/min/kg, $\bar{x}(\mathit{s})$)	0.20 (0.02)	0.20 (0.02)	0.20 (0.02)	0.20 (0.02)
EE 13 km/h (kcal/min/kg, $\bar{x}(\mathit{s})$)	0.22 (0.02)	0.21 (0.02)	0.22 (0.02)	0.21 (0.02)
Training
Days-a-week (n, $\bar{x}(\mathit{s})$)	5 (1)	4 (2)	5 (1)	4 (2)	U = 963, p = 0.33
Weekly volume (km, $\bar{x}(\mathit{s})$)	58.7 (23.9)	42.5 (24.7)	50.9 (25.3)	46.4 (26)	U = 980, p = 0.27
Performance
10K time (mm:ss, $\bar{x}(\mathit{s})$)	43:33 (04:49)	39:55 (03:58)	40:53 (04:11)	41:08 (04:57)	t = −0.25, p = 0.80

n: count; $\bar{x}$: mean; $(\mathit{s})$: standard deviation.

The All column is divided in females (F) and males (M), EE: energy expenditureIn the C0-C1 comparison, Χ2: Chi-squared test; U: Mann-Whitney U test for independent groups; t: t-test for independent groups; F: 2-way ANOVA with repeated measures for the main effects of cluster (C), and speed (S) and the interaction effect of these two (CXS).

We used Principal Component Analysis (PCA) to reduce the dimensionality of the data used to characterise participants’ running technique. PCA linearly combines the original data into a new set of uncorrelated variables known as principal components (PCs) which can capture different aspects of the variance within the dataset (Jollife & Cadima, 2016). A small number of PCs can typically capture most of the variance in the original space, providing dimensionality reduction. We chose to keep the number of PCs which explained 99% of the variance in the dataset. This approach was found to be a good compromise between data compression and data representation, as well as cluster separability in a previous study comparing different dimensionality reduction approaches and their effects on clustering of a subset of the present data [Rivadulla et al., (under review)]. To assess how well this PCs were capturing the original data, we looked at the reconstruction of the original data with the retained PCs compared to the original variables as per [Rivadulla et al., (under review)]. The 1D and 0D variables forming each dataset were standardised and combined prior to applying PCA.

Agglomerative hierarchical clustering was then applied using the Euclidean distance as a measure of similarity and the Ward linkage criterion (Ward, 1963) to merge clusters. Initially, each data point constitutes its own cluster, and clusters are iteratively merged while minimising within-cluster variance until all the data points belong to the same cluster, thereby creating a hierarchical structure. To determine the optimal number of clusters, we evaluated various cluster solutions using three internal validity metrics: Silhouette score, Calinski-Harabasz index, and Davies-Bouldin index (Akhanli & Hennig, 2020; Arbelaitz et al., 2013). The Silhouette score measures the cohesion and separability between clusters, ranging from −1 to 1, with higher scores indicating better clustering. The Calinski-Harabasz index assesses clustering quality based on the distance from each point to the global and cluster centroids, while the Davies-Bouldin index compares within-cluster and between-cluster distances. The final number of clusters was determined by majority voting among the three metrics, with priority given to the Silhouette score in the event of a tie, due to its robustness (Arbelaitz et al., 2013).

The dimensionality reduction and clustering process was repeated for each considered speed independently and for the overall dataset, in which we included (i.e., concatenated) data from every speed. We shall refer to these two approaches as ‘single speed’ and ‘multispeed’, respectively. To assess the agreement between the clustering partitions resulting from each single speed and the multispeed approaches, the adjusted mutual information index (AMI) (Vinh et al., 2010) was calculated. The AMI adjusts for the agreement between two different partitions due to chance and ranges between 0 (random agreement) and 1 (perfect agreement).

The running kinematics of the resulting clusters were statistically compared using SPM1D for the 1D variables and discrete methods for the 0D variables (α = 0.05). Similarly, clusters were also compared for their running economy, demographics (sex and age), anthropometrics (height, body mass, trunk, leg, thigh, shank and foot length, and pelvis width), physiological capacity ($\dot{V}{O}_2\mathit{peak}$ and lactate threshold), training (training days per week and weekly volume) and performance (10K time). Since the optimal number of clusters was found to be two for every condition, t-tests were used when the normality assumption was met and Mann–Whitney U tests otherwise for every continuous and ordinal 0D variable, whereas non-parametric SPM t-tests were used for the 1D variables. Nominal variables were compared using the ${\text{Chi}}^2$ test. For the multispeed condition and running economy, parametric 2-way ANOVA with repeated measures and 1D non-parametric SPM 2-way ANOVA with repeated measures were used to assess the main effects of cluster, speed and the interaction between these two factors on 0D and 1D variables, respectively. Significant cluster and cluster x speed interaction effects were followed up by Bonferroni adjusted pairwise post hoc tests. For the 1D interaction effects, post hoc tests were conducted on the change of the variable from one speed with respect to (wrt) the previous one (i.e., change in variable at speed 12 km/h wrt to 11 km/h, and 13 km/h wrt 12 km/h). Cohen’s d effect sizes (Cohen, 1977) are reported for pairwise 0D comparisons. All the analyses were performed on Python 3.8, using the scikit-learn 1.3.8 (PCA, clustering and related analyses), pingouin 0.5.3 and spm1d 0.4.20 (statistics) packages. Data, environment file and code to replicate the current results can be found in https://github.com/adrianrivadulla/clustering_runners.

Results

For dimensionality reduction, 20 PCs were enough to account for 99% of the variance for the 11, 12 and 13 km/h datasets whereas 24 PCs were needed for the multispeed condition. These models reconstructed the data with mostly acceptable errors (Appendix). The Silhouette, Calinski-Harabasz and Davies-Bouldin scores agreed in indicating the optimal number of clusters was two (C0 and C1) regardless of the speed condition (Figure 1).

Figure 1. Internal validity scores for clustering partitions with different number of clusters (k) for each single speed and the multispeed dataset.

The agreement between clustering partitions at different speeds and the multispeed condition was generally low, implying that members of each cluster changed at 11, 12, 13 km/h and for the multispeed condition (Figure 2). Since the clustering partitions changed for every single speed dataset, we considered the multispeed dataset the most comprehensive description available of runners’ technique. Therefore, the statistical comparisons between clusters presented hereafter used the clusters obtained by the multispeed model regardless of the speed being compared (comparisons of clusters derived from each individual speed dataset may be found in the Appendix).

Figure 2. Dendrogram for each speed condition. The bottom leaves of the tree represent each single runner in the dataset. The vertical length of the branch is proportional to the distance between the two clusters being merged. Vertical lines underneath each leaf are coloured according to the cluster membership for each participant in the previous condition in descending order as indicated by the black arrows on the sides e.g., vertical lines underneath the leaves at the multispeed dendrogram correspond to the cluster membership at 13 km/h. The reported AMI scores correspond to each cluster membership comparison, with respect to the previous condition.

The two-way ANOVA revealed a main effect of speed for every kinematic variable considered (see Figure 3 and 4) and a main effect of cluster for duty factor, trunk to pelvis, pelvis, hip and knee angles. Cluster follow-up post-hoc tests showed the C0 cluster exhibited greater duty factor at every speed (Figure 2), less trunk-to-pelvis extension (Figure 6) a less anteriorly tilted pelvis (Figure 7), and less hip flexion (Figure 8) throughout the entire gait cycle at every speed but no differences in knee (Figure 9) and ankle (Figure 10) angles. Because of the postural behaviour observed and for the sake of readability, we shall refer to C0 as the ‘neutral pelvis cluster’ and to C1 as the ‘tilted pelvis cluster’. During mid/late stance, there were significant interaction effects between cluster and speed for the vCOM, and hip, knee and ankle angles. Post-hoc analysis following significant cluster x speed interaction effects comparing the change in kinematic variables from one speed to another showed that the tilted pelvis cluster showed greater increase in vCOM (Figure 5), alongside a greater increase in hip (Figure 8) and knee extension (Figure 9) and ankle dorsiflexion (Figure 10) towards the end of the contact phase when changing from 11 to 12 km/h compared to the neutral pelvis cluster. The average pattern for each cluster throughout the gait cycle can be visualised in Figure 11.

Figure 3. 0D variables: Violin plots by cluster for stride frequency and duty factor. Main effects for cluster (C), speed (S) and interaction effect (CxS) are reported in the titles. A significant cluster effect on duty factor was found and posthoc tests results are reported underneath each graph. CT and ST in the Duty Factor y label correspond to contact time and stride time, respectively.

Figure 4. Main effect of speed on 1D variables: Average (and standard deviation cloud) stride pattern (0–100 % of time) for the entire sample at 11, 12 and 13 km/h, respectively. Vertical dotted lines represent the average toe-off for each speed. Grey patches represent areas of the gait cycle where a significant main effect of speed was found. Full statistics are reported in the Appendix.

Figure 5. Cluster (C) and cluster x speed interaction effects (CxS) for the vCOM during the stride (foot-strike to foot-strike, 0–100%): Top row: Average vCOM displacement (and standard deviation cloud) for each cluster. Bottom row: Average change in vCOM displacement (and standard deviation cloud) for each cluster. Areas of the gait cycle where significant differences were found in post-hoc analysis following a significant cluster x speed interaction are highlighted in grey. Vertical dotted lines indicate the end of the average contact phase of the stride.

Figure 6. Cluster (C) and cluster x speed interaction effects (CxS) for the trunk-to-pelvis angle during the stride (foot-strike to foot-strike, 0–100%): Top row: Average trunk-to-pelvis angle (and standard deviation cloud) for each cluster. Areas of the gait cycle where significant differences were found in post-hoc analysis following a significant cluster effect are highlighted in grey. Bottom row: Average change in trunk-to-pelvis angle (and standard deviation cloud) for each cluster. Vertical dotted lines indicate the end of the average contact phase of the stride.

Figure 7. Cluster (C) and cluster x speed interaction effects (CxS) for the pelvis tilt angle during the stride (foot-strike to foot-strike, 0–100%): Top row: Average pelvis angle (and standard deviation cloud) for each cluster. Areas of the gait cycle where significant differences were found in pot-hoc analysis following a significant cluster effect are highlighted in grey. Bottom row: average change in pelvis angle (and standard deviation cloud) for each cluster. Vertical dotted lines indicate the end of the average contact phase of the stride.

Figure 8. Cluster (C) and cluster x speed interaction effects (CxS) for the hip angle during the stride (foot-strike to foot-strike, 0–100%): Top row: Average hip angle (and standard deviation) cloud for each cluster. Areas of the gait cycle where significant differences were found in pot-hoc analysis following a significant cluster effect are highlighted in grey. Bottom row: Average change in hip angle (and standard deviation cloud) for each cluster. Areas of the gait cycle where significant differences were found in pot-hoc analysis following a significant cluster x speed interaction are highlighted in grey. Vertical dotted lines indicate the end of the average contact phase of the stride.

Figure 9. Cluster (C) and cluster x speed interaction effects (CxS) for the knee angle during the stride (foot-strike to foot-strike, 0–100%): Top row: Average knee angle (and standard deviation cloud) for each cluster. Bottom row: Average change in knee angle (and standard deviation cloud) for each cluster. Areas of the gait cycle where significant differences were found in pot-hoc analysis following a significant cluster x speed interaction are highlighted in grey. Vertical dotted lines indicate the end of the average contact phase of the stride.

Figure 10. Cluster (C) and cluster x speed interaction effects (CxS) for the ankle angle during the stride (foot-strike to foot-strike, 0–100%): Top row: Average ankle angle (and standard deviation cloud) for each cluster. Bottom row: Average change in ankle angle (and standard deviation cloud) for each cluster. Areas of the gait cycle where significant differences were found in pot-hoc analysis following a significant cluster x speed interaction are highlighted in grey. Vertical dotted lines indicate the end of the average contact phase of the stride.

Figure 11. Time-synchronised visualisation of cluster average stride kinematics (C0 = neutral pelvis, in orange; C1 = tilted pelvis, in green) using a modified version of the 3DGaitModel2392 model (Delp et al., 2007) with just the right leg (only right leg data were considered in this study). The cluster average of each kinematic variable was prescribed to the model locking every other degree of freedom at 11 km/h (top), 12 km/h (middle) and 13 km/h (bottom). The segments of the model have been scaled to the average cluster segment length. The selected frames include foot-strike, mid-stance, toe-off, leg recovery, mid swing and foot strike.

Despite the differences in their running technique, clusters did not exhibit significant differences in their running economy for each of the considered speeds (Table 2 and Figure 12).

Figure 12. Running economy by cluster for the three speeds considered. Main effects of cluster (C), speed (S) and interaction effect (CxS) are reported in the title.

Similarly, clusters did not show differences in their demographics, anthropometrics, physiological capacity and training and performance (Table 2 and Figure 13).

Figure 13. Demographics, anthropometrics, physiological capacity, training and performance variables by cluster. No significant differences were found between clusters. Full statistics are reported in Table 2.

Discussion

This study aimed to cluster runners based on a comprehensive selection of kinematic variables describing running technique and then compare the resulting groups of runners for their kinematics, running economy, demographics, anthropometrics, physiological capacity, training and performance. The novelty of our work resides in the use of clustering to investigate the relationship between running technique and economy, and to provide an extensive pool of anatomical, physiological, training and performance characteristics to discuss their potential association with running technique-based clusters. We identified two clusters of runners who showed distinct kinematics: the ‘neutral pelvis’ cluster, whose runners showed less pelvic tilt, less trunk-to-pelvis extension and less hip flexion throughout the gait cycle and greater duty factor at every speed compared with the ‘tilted pelvis’ cluster. Athletes using these two running techniques did not show differences in their running economy in support of our hypothesis based on the self-optimisation and preferred movement path theories. Clusters did not show differences in their demographics, anthropometrics, physiological metrics, training and performance either suggesting each of these factors alone may not maximise kinematic differences in running technique. We also assessed whether clustering the same runners at different speeds independently would have led to consistent clustering partitions. Our results showed that this was not the case and therefore, clusters of runners found at one single speed may not generalise to a runner’s entire speed profile.

For this study, we aimed to recruit a sample of runners with a relatively consistent level of running capability, while acknowledging the diversity in sex and age typically present in the running population (Nesburg et al., 2023). Biomechanical differences in runners from different sexes (Xie et al., 2022) and ages (DeVita et al., 2016; Paquette et al., 2018) and interactions between age and sex (Boyer et al., 2017; Phinyomark et al., 2014), as well as athletic level (Clermont et al., 2017) are indeed typically used to justify runners’ subgrouping in biomechanics research. Despite the cluster differences in running kinematics and in agreement with previous clustering studies (Phinyomark et al., 2015; Senevirathna et al., 2023), there was a lack of evidence for any association between the two clusters and their underlaying demographics. It should be noted that although some studies have found differences in knee and hip flexion between males and females, sex differences may be primarily related to frontal and transverse plane kinematics (Phinyomark et al., 2014) which were not included in our analysis. Similarly, biomechanical changes associated with age within the considered age range (18–50) can be less apparent when training volume is matched (Paquette et al., 2018). Overall, it appears that sex and age may not be key determinants of maximal kinematic variance in running technique.

The observed differences in pelvic tilt, trunk-to-pelvis extension and hip angle between the neutral and the tilted pelvis cluster could be described as a near-constant offset throughout the entire gait cycle as opposed to gait phase specific. These could be two different strategies to keep the trunk upright depending on pelvic tilt: The more anterior pelvic tilt, the more trunk-to-pelvis extension is needed to keep the trunk up and the more flexed the hip is throughout the stride. Radiographic studies have shown pelvic tilt to be considerably variable in healthy individuals (standing pelvic tilt = −13 ± 6°, range [4.5, −27]) (Vialle et al., 2005). Indeed, in a retrospective analysis of the pelvic tilt during standing in the static calibration trial, we found that runners in the neutral pelvis cluster exhibited a more neutral pelvis compared to the tilted pelvis cluster (−12 ± 4° and −16 ± 3°, respectively; t = 4.75, p < 0.001, d = 1.06), which may contribute to explaining the observed differences during running. No other anthropometrics or anatomical differences were found between clusters. Insights into more complex anatomical aspects such as bony morphology, the multiple strong muscles actuating the pelvis and neuromuscular control strategies may be required to better understand the observed differences in pelvis orientation during standing and running. The large variability in features studied and methodological approaches taken in previous running clustering studies (Table 1) prevent direct comparison with our current results. Considering the dual-axis framework (van Oeveren et al., 2021), the two found running techniques would fall somewhere in between the bounce (long flight time, closer to the tilted pelvis cluster) and stick (long contact time, closer to the neutral pelvis cluster) styles, also described as aerial and terrestrial runners in previous studies (Lussiana et al., 2017), since differences between clusters were found in duty factor but not in stride frequency. Nevertheless, the five running styles proposed in the dual-axis framework did not emerge in our dataset, in agreement with a recent running clustering study (Floría et al., 2024). Further research is warranted to elucidate how well the dual-axis framework can capture running styles and to integrate expert knowledge and data-driven approaches.

There were no significant differences in running economy at any speed between the two clusters suggesting the participants from this study were able to find two distinct running techniques at a similar energy cost, confirming our initial hypothesis. This finding aligns with the self-optimisation and preferred movement path theories, and reinforces the importance of running technique variability in research and training. Duty factor and pelvic tilt played a key role in defining the two kinematic patterns. It is believed that smaller duty factor can improve running economy (Nijs et al., 2022) and that pelvic tilt needs to be minimised to conserve energy and maintain efficiency (Novacheck, 1998; Schache et al., 1999; van Hooren, Lennartz, et al., 2024). Our results do not appear to support those statements like Lussiana et al. (2017) and van Hooren, Jukic, et al. (2024), who found no differences in running economy between aerial (small duty factor) and terrestrial (large duty factor) runners. It is also worth noting that the neutral pelvis cluster showed greater duty factor than the tilted pelvis cluster, which may be described as ‘the bitter with the sweet’. Running requires the coordinated action of the full body, and studies looking at one variable in isolation without considering or controlling the other components of the system may also contribute to the contrary findings about running technique and economy. Similarly, although greater pelvic tilt has been found in recreational runners compared to competitive runners (Clermont et al., 2017), the athletic level of the present sample can be defined as recreationally active to trained runners (McKay et al., 2022) with no differences in the physiological capacity ($\dot{V}{O}_2\mathit{peak}$ and lactate threshold), training (running days a week, weekly distance run), and performance (10K times) between clusters. Ultimately, these results seem to build upon the idea that the search for one single most economical running technique model that every runner should strive for may be ill-posed, and that greater personalisation might be more appropriate (Moore, 2016; van Hooren et al., 2020).

The current results advocate for embracing inter-individual running technique variability in training too. ‘One size fits all’ training approaches may fail (or be suboptimal) to promote improvements to athletes’ running technique and economy. Coaches may choose to enhance a runner’s strengths while working on the weaknesses. For instance, despite the apparent lack of association of the tilted pelvis technique with worse economy and lower athletic-level in the present study, it is unclear whether this pattern may be detrimental compared to the neutral pelvis technique at faster speeds where greater hip extension may be required or as technique adapts to fatigue in a longer run. Under the assumption that a reduced duty factor and less anterior pelvic tilt can improve economy (Nijs et al., 2022; Schache et al., 1999; van Hooren, Lennartz, et al., 2024), neutral pelvis runners may benefit from maintenance strength work on the muscles stabilising the pelvis while trying to reduce duty factor since this may help at faster speeds (van Hooren, Jukic, et al., 2024). Runners in the tilted pelvis group on the other hand could increase core strength and hip mobility exercises in their training so their more anteriorly tilted pelvis position does not become a potential limitation at faster speeds or as fatigue develops. Nevertheless, it is important to stress that external validation and larger samples of runners are key to further establish the practical significance of these clusters. Working closely with coaches and running researchers is key in this sense to make the most of blending data-driven analysis and expert knowledge.

Comparing the multispeed and single speed clustering resulting partitions raised important considerations for the comprehensive representation of runners in clustering studies. The three single-speed clustering models and the multispeed model agreed in highlighting the distinct kinematic patterns between clusters previously discussed. However, the clustering solution at 11 km/h displayed a notable shift in cluster composition compared to the other speeds and members of the tilted pelvis cluster at this speed showed greater vCOM displacement and plantar flexion around foot contact (see Appendix). These results underscore the impact that clustering at a single speed may have on subsequent decision-making in training design and injury research. In this study, we proposed a simple alternative which was including data from several speeds in the clustering dataset. Due to the high correlations between data from different speeds, including three speeds in the multispeed dataset had a minimal impact in terms of the number of PCs required to capture 99% of the variance compared to single speed datasets (e.g., 24 PCs as opposed to 20 PCs) while still providing an appropriate representation of the data based on the reconstruction errors (Appendix). This approach circumvents the cluster agreement issues of single speed clustering, had minimal methodological considerations (e.g., small extra data collection and processing time) and provides a more comprehensive description of running technique. Furthermore, analysis across speeds allowed us to investigate whether runners with different kinematic patterns adapted to speed differently. There were general changes to running technique as speed increased (Figures 3–4) which were aligned with previous literature (Fukuchi et al., 2017; Orendurff et al., 2018). The tilted pelvis cluster also showed a greater increase in vCOM displacement (Figure 5) alongside greater hip (Figure 8) and knee (Figure 9) extension and ankle dorsiflexion (Figure 10) towards the end of the contact phase when changing from 11 to 12 km/h. Although the magnitude of these differences was marginal if we consider the reliability of marker-based motion capture (McGinley et al., 2009; Pohl et al., 2010), the question of whether these different cluster responses will accentuate at faster speeds remains.

It must be noted that the identified clusters may be sample and method dependent and other clustering algorithms may result in different clustering partitions. We chose agglomerative hierarchical clustering with Ward’s linkage, a popular algorithm in the data mining literature and in previous running biomechanics clustering studies (Table 1), but this approach is known to be sensitive to outliers and to favour spheroidal clusters regardless of the natural configuration of the data (Zimmerman et al., 1982). Density-based methods such as HDBSCAN (Campello et al., 2013) may help dealing with outliers and finding shape assumption free clusters, but they typically require larger datasets than what we had available. Larger samples and a thorough exploration of clustering algorithms may help reinforce the present clusters or to find more robust ones. The resulting clusters are variable dependent. We focused on sagittal plane variables because most body movements in running occur in this plane and energy expenditure during running is mostly determined by sagittal plane motion (Arellano & Kram, 2014; van Hooren, Lennartz, et al., 2024). As such, it is believed different running patterns should emerge in these variables (van Oeveren et al., 2021). The treadmill protocol using a 1% gradient was designed to replicate the energetic demands of overground running (Jones & Doust, 1996) but more recent work has indicated that a 0% gradient can be appropriate to replicate the oxygen cost of running at submaximal speeds (Miller et al., 2019). A 0% gradient can be more representative of overground running kinematics and may be considered in future studies investigating running technique and economy. We allowed participants to use their preferred footwear to capture their habitual running pattern. However, footwear can affect lower limb biomechanics (Nigg et al., 2012) and running economy (Hoogkamer et al., 2018), potentially acting as a confounding factor in the clustering process and running economy comparison. Information about the shoe-type used by the runners could have added to the discussion. Including kinematic variability, kinetic, and muscle activity variables could also help both enrich the definition of running technique and the discussion of the found clusters, representing an area to be explored in future studies. Finally, the identified clusters respond to a mathematical optimum (i.e., maximum Silhouette score) and close work with running coaches is key for the external validation of the present results and to maximise the usefulness and impact of the present results in running practice.

Conclusions

We have investigated running technique in healthy trained runners and identified two distinct running patterns characterised by a neutral pelvis and larger duty factor (neutral pelvis cluster) versus a more tilted pelvis and shorter duty factor (tilted pelvis cluster). These differences were close to constant throughout the stride and we suspect trunk-to-pelvis angle and hip flexion are adjustments to compensate for variability in pelvic tilt to keep the trunk upright. These patterns were identified at three different speeds so they may represent individual preferences for such patterns regardless of the speed athletes are running at. These two techniques seemed similarly economical and they were not associated with the demographics, anthropometrics, physiological capacity and athletic performance level of the cluster members. Coaches are advised to consider inter-individual differences in running technique to enhance runners’ strengths while working on their weaknesses when designing their training programmes to achieve greater response from their athletes as opposed to ‘one size fits all’ approaches. From a methodological perspective, clustering at different speeds independently led to different cluster memberships being formed and additional kinematic differences being found between clusters. We suggest researchers should consider including several speeds in their clustering studies to obtain a more comprehensive representation of runners’ technique, since single speed clustering of the same athletes at different speeds can lead to different cluster make-up and affect subsequent analysis and decision-making.

Acknowledgments

The authors would like to thank Mr. Andrew Hammocks for his assistance in data collection for this study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/14763141.2024.2372608

Additional information

Funding

This work was supported by a PhD studentship funded by the University of Bath and NURVV, Ltd.