## 2 Materials and Methods
### Signal Processing and Movement Component Analysis  
All data processing was conducted in MATLAB version 2024a (MathWorks Inc., Natick, MA, USA). The 30 s constant running of each trial of each participant was selected for analysis. An example of original running movement is represented in Supplementary Video S1. Twenty-six markers from each dataset, contributing 78 spatial coordinates (x, y, z), were interpreted as 78-dimensional posture vectors. Each dataset was preprocessed by subtracting the mean posture vector and then normalized to the mean Euclidean distance [13]. The data from all volunteers were concatenated to form an input matrix (100 [sampling rate] \(\times\) 5 [number of trials] \(\times\) 30 [trial duration] \(\times\) 19 [number of participants] \(\times\) 78 [marker coordinates]] for further PCA.  
PCA was calculated using a singular-value decomposition of the covariance matrix through the PManalyzer software [13] to decompose all kinematic marker data retrieved from five running conditions of all participants into a set of orthogonal eigenvectors, i.e., principal components (PCs), where \(k\) denotes the order of movement components. Animated stick figures can be created to characterize each eigenvector's movement pattern, which has been called \"principal movement\" (PM\({}_{\text{k}}\)), visually representing movement synergies that contributed to achieve the running task [13; 26]. Moreover, the actual time evolution (\(t\), time series) of individual PM\({}_{\text{k}}\) is quantified by the PC scores called principal positions (PP\({}_{\text{k}}\)(\(t\))), representing positions in posture space, i.e., the vector space spanned by the PC eigenvectors [13]. Additionally, the acceleration of each PM called principal acceleration (PA\({}_{\text{k}}\)(\(t\))) can also be determined through the second-time differentiations, reflecting the neuromuscular control of each PM, as it is associated with myoelectric activity [7]. The term \"principal\" in the variable names denotes that these variables were obtained through a PCA, and (\(t\)) indicates that these variables are functions of time \(t\)[13].  
In this study, the PM time series were filtered with a third-order zero-phase 6 Hz low-pass Butterworth filter to avoid noise amplification in the differentiation processes. The leave-one-out cross-validation was then used to evaluate the vulnerability of the PM\({}_{\text{k}}\) and the dependent variables to changes in the input data matrix to address validity considerations [13]. The first five PCs, which proved robust and explained 97.7% of the total variance, were selected to test the hypotheses. Each participant-specific relative explained variance (rVAR) of PP\({}_{\text{k}}\) (PP\({}_{\text{k}\_}\)rVAR) and PA\({}_{\text{k}}\) (PA\({}_{\text{k}\_}\)rVAR) calculated from the PP\({}_{\text{k}}\)(\(t\)) and PA\({}_{\text{k}}\)(\(t\)) represent the coordinative structure or composition of running movements and running acceleration, respectively [14]. The PP\({}_{\text{k}\_}\)rVAR and PA\({}_{\text{k}\_}\)rVAR quantify how much (in percent) each PM contributed to the total variance in postural positions and postural accelerations, respectively [14]. It has been suggested in the literature that movement components with limited positional amplitudes, when performed rapidly, significantly affect accelerations and, in turn, the forces acting in the system [27].