Sixteen high-performing endurance runners (mean ± SD; age: 21.00 ± 1.31 years; height: 172.00 ± 7.13 cm; body mass: 61.64 ± 8.72 kg) participated in the present study, comprising 9 male and 7 female runners. Sample size was determined using G*Power (version 3.1.9.2, University of Kiel, Germany) for a repeated-measures ANOVA, with the following parameters: significance level (α) = 0.05, statistical power (1 − β) = 0.75, and effect size (Cohen’s f = 0.25) (39, 40). The primary outcome variables were selected based on a prior study that reported statistically significant differences in Achilles tendon loading during the running stance phase (41). Detailed participant information is presented in Table 1. Recruitment was conducted via social media platforms and posters distributed at the university and local running clubs. All participants met the following inclusion criteria: aged 18–35 years; running at least three times per week; completion of at least two half-marathon competitions in the past 6 months; at least 40 km running per week (37); a habitual rearfoot strike pattern; having achieved a self-reported 10-km personal best time of less than 32 min (males) or 36 min (females) (34); and no lower extremity injuries in the past 6 months. This study protocol received approval from the Ningbo University Human Ethics Committee (TY2024034), and all participants provided informed consent.

Prior to testing, all participants performed a 10-min warm-up on the treadmill at 10 km/h and familiarized themselves with the experimental protocol. Participants wore their appropriate running shoes and spandex tights during the data collection. A total of 38 reflective markers (14 mm diameter) were affixed bilaterally on the lower limbs and trunk (Figure 1A). The lower-body marker set included five anatomical segments: pelvis, bilateral thighs, shanks, feet, and torso (42). All markers were placed individually, except for four T-shaped frames on the thighs and calves, each consisting of three tracking points. Furthermore, all reflective markers were attached with skin-medical tape that prevented them from dropping off during the whole running protocol.

Moreover, we used the wireless 32-channel surface electromyogram (EMG) system (Delsys, Boston, MA, USA) electrodes to collect the EMG signals at a sampling of 1,000 Hz during the running stance phases. Skin preparation was performed in accordance with SENIAM guidelines. We placed the electrodes on the rectus femoris (RF), vastus medialis (VM) and lateralis (VL), gastrocnemius medialis (GM) and lateralis (GL), and tibialis anterior (TA) (Figure 1B) to collect the muscle activation during the half marathon running (43). Before the running trials, maximal voluntary contractions (MVC) were measured in all participants (43). A motion capture system consisting of 10 cameras (Oxford Metrics Ltd., Oxford, UK) and two force plates (Kistler Type 9,281 B, Kistler Instrument AG, Winterthur, Switzerland) was synchronized to collect marker trajectories and ground reaction force (GRF) data at 200 Hz and 1,000 Hz, respectively (26, 43) (Figure 1D).

The data collection started with a static calibration trial for model simulation. Participants then ran along a 15-meter overground runway at 4.0 ± 0.5 m/s to obtain pre-running kinematic and GRF data during the stance phase. Thereafter, participants completed the half-marathon running task at a constant speed of 14 km/h on the Zebris FDM-T treadmill system (Zebris Medical GmbH, Isny, Germany). Spatiotemporal parameters were obtained directly from the treadmill, which continuously records right-foot contact and toe-off phases (41). For each running distance (pre-running, 10 km, and 20 km), 30 s of treadmill data were collected for spatiotemporal analysis. After completing 10 km and 20 km on the treadmill, participants returned to the 15-meter overground runway, where synchronized kinematic and GRF data were collected again. Each participant completed four successful running trials through the embedded force plate. The successful trials are defined as the participants running with the entire right foot placed on the force plate at 4.0 ± 0.5 m/s for each participant. Additionally, smart speed devices (Smart Speed, Fusion Sport Inc., Burbank, CA, USA) were positioned on either side of the track to monitor running velocity (43). The running biomechanics trials were conducted 3 times in total: pre-running, post-10 km running, and post-20 km running.

Raw marker trajectories were labelled in Vicon Nexus 1.8.5 (Vicon, Metrics Ltd., Oxford, UK), and the C3D files of the running stance phase for the lower limbs were subsequently exported. All the collection data were then imported into Visual3D (V6.0, C-Motion, Germantown, MD, USA) to calculate joint kinematics and kinetics variables (42). Kinematic data were low-pass filtered at 20 Hz, and ground reaction forces were low-pass filtered at 50 Hz, in accordance with previous studies (41). Raw EMG signals were processed using a zero-lag, fourth-order Butterworth band-pass filter (20–450 Hz) to remove motion artifacts and high-frequency noise. The linear envelopes of the EMG signals were then normalized to the peak value obtained during the maximal voluntary contraction (MVC) for each muscle (43). All EMG signals were then time-normalized to 100 points (1–100%) to allow comparison with simulated muscle activation profiles (43). All biomechanical variables for each stance phase were time-normalized to 101 points utilizing MATLAB software (R2024a; MathWorks, Natick, MA, USA).

Joint angles of the ankle, knee, and hip were computed using Cardan angles, with the distal segment relative to the proximal segment in Visual 3D software. The joint moments were calculated using a standard inverse dynamics approach (44), and all kinetic variables were normalized according to each participant’s body mass. Analyses focused on the running stance phase of the right leg, defined as the period from initial ground contact to toe-off, during which the value of ground reaction force exceeded 20 N (43). The range of motion (ROM) for the hip, knee, and ankle joints was also compared during the braking and push-off phases (45). Ankle dorsiflexion ROM was defined as the interval from initial contact to peak dorsiflexion, and ankle plantarflexion ROM as the interval from peak dorsiflexion to toe-off. Knee flexion ROM was defined from initial contact to peak knee flexion, and knee extension ROM from peak knee flexion to toe-off. Hip flexion ROM was defined from initial contact to peak hip flexion, and hip extension ROM from peak hip flexion to toe-off (46). Angular impulse was determined as the time integral of joint moment during the loading and push-off phases of stance (47).

Cumulative impulse ( I C ):

The total number of strides required to complete a 1,000-m distance was denoted as n. For each joint, the cumulative load was calculated as the product of the stance-phase angular impulse in the dominant moment direction and the stride count (n) necessary to cover 1,000 m (48, 49). In addition, we quantified the peak joint moment per stance phase, as well as the kinematic parameters, including peak joint angles and ROM of the ankle, knee, and hip joints during the braking and push-off phases.

The customized musculoskeletal analysis was conducted in OpenSim (Version 4.3, Stanford, CA, USA), which has been widely used to estimate muscle forces and joint contact force during the running stance phases. The 2,392 musculoskeletal model, which comprises 10 rigid body segments, 23 degrees of freedom, and 92 musculotendon units (42), was utilized for all analyses. Marker trajectories exported from Visual3D in “osim” format were imported into OpenSim for musculoskeletal simulation, and the model was scaled for each participant based on anthropometric measurements. Firstly, inverse kinematics (IK) and inverse dynamics (ID) tools were conducted to calculate the joint angle and joint moment. Then, we used the static optimization (SO) tool to compute the lower limb muscle forces during the running stance phases. The analysis tool was applied to compute the patellofemoral contact force during the running stance phases (43). The patellofemoral contact stress was estimated by dividing the compressive component of the joint contact force by the contact area at each knee flexion angle (Figure 1C). Contact area values were obtained from sex-specific MRI-based measurements reported by Besier et al. (50).

The AT force was calculated as the sum of the forces generated by the soleus and the medial and lateral gastrocnemius muscles (Figure 1C) (42). Changes in AT length were estimated by dividing the tendon force by the tendon stiffness, which had previously been measured via ultrasound at 420 N·mm⁻¹ in earlier studies (20, 51). The AT strain was calculated based on previous research using a formula relative to the resting tendon length, which was defined as 250 mm (20).

Simultaneously, we compared the peak stress and strain values per step for the AT and the PFJ, which reflect cumulative loading. Using the number of steps required to cover 1 km, we estimated the CD by applying a weighted impulse for PFJ stress, AT force, and AT strain, as demonstrated in Equation 1. This method, proposed by Firminger et al., incorporates a tissue-dependent weighting factor to estimate kilometer-specific loading for running distances of 0 km, 10 km, and 20 km (8).

In this study, n denotes the number of right-foot contacts, tᵢ represents the start time of the running stance phase, t f the end time of the stance phase, and X_s the measure of internal bone peak impulse loading. The parameter b is a weighting factor that reflects the slope of the power function relating fatigue life to tissue-specific stress or strain. For bone, b was set to 7, and for tendon tissue it was set to 9, based on previous studies (8, 18, 52, 53). The weighted impulse approach emphasizes that load magnitude is a more critical determinant of cumulative tissue strain than the number of loading cycles. The step count per kilometer was calculated by taking the inverse of the stride duration (1/stride duration), which reflects the total number of right-foot contacts included in the cumulative strain analysis of both bone and tendon tissues.

All biomechanical parameters at the three running distances were presented as means and standard deviations. Statistical analyses were performed using SPSS version 25.0 (SPSS Science, Chicago, IL, USA). Data normality was assessed with the Shapiro–Wilk test. One-way repeated-measures ANOVAs with Bonferroni corrections for multiple comparisons were used. The level of significance was set at 0.05.

A significant main effect of running distance was observed on contact time (F = 43.396, p < 0.001, partial η² = 0.41), step length (F = 10.791, p < 0.001, partial η² = 0.15), step frequency (F = 9.523, p = 0.003, partial η² = 0.13), and stride length (F = 12.250, p < 0.001, partial η² = 0.16) across various running distances (Table 2). Specifically, contact time increased significantly with running distance. At 20 km, the contact time was significantly greater than at 0 km (p < 0.001) and 10 km (p < 0.001). Additionally, after 10 km of running, contact time was significantly longer than at 0 km (p = 0.003). Step frequency also increased significantly with distance, with values at 10 km (p < 0.001) and 20 km (p < 0.001) being higher than at 0 km. Additionally, the step frequency was significantly increased at 20 km (p = 0.002) in comparison to the 10 km running stance phases. In contrast, step length significantly decreased at 20 km (p < 0.001) and 10 km (p < 0.001) in comparison to 0 km. Finally, stride length exhibited a significant decrease with increasing running distance. A significant reduction was observed at 20 km compared to 10 km (p = 0.008).

A main effect of running distance was observed on the ankle push-off phase ROM (F = 7.052, p = 0.004, partial η² = 0.32), knee ROM during the braking phase (F = 16.120, p < 0.001, partial η² = 0.52), and hip ROM during the braking phase (F = 10.842, p < 0.001, partial η² = 0.42) (Table 3; Figure 2). The ankle ROM during the push-off phase was significantly greater at 0 km compared to 10 km (p < 0.001) and 20 km (p = 0.003) (Figure 2H). Conversely, the hip ROM during the braking phase significantly declined at 10 km (p < 0.001) and 20 km (p = 0.006) compared to 0 km (Figure 2F). In contrast, the knee ROM during the braking phase significantly increased at 10 km (p < 0.001) and 20 km (p < 0.001) compared to 0 km (Figure 2E). There were no significant changes observed in the ankle ROM during the braking phase (F = 1.421, p = 0.258, partial η² = 0.09), knee push-off ROM (F = 1.596, p = 0.223, partial η² = 0.10), or hip push-off ROM (F = 0.292, p = 0.669, partial η² = 0.02).

A significant effect of running distance was shown on the peak ankle plantarflexion angle (F = 25.374, p < 0.001, partial η² = 0.29), peak knee flexion angle (F = 11.677, p < 0.001, partial η² = 0.16), peak hip flexion angle (F = 8.563, p < 0.001, partial η² = 0.12), and peak hip extension angle (F = 4.942, p = 0.009, partial η² = 0.07) during the half marathon (Table 3; Figure 2). Specifically, the peak ankle plantarflexion angle at 20 km (p < 0.001) and 10 km (p < 0.001) showed a significant reduction compared to 0 km (Figure 2A). As running distance increased, the peak knee flexion angle significantly increased at 10 km (p < 0.001) compared to 0 km; however, after 20 km, the peak knee flexion angle significantly decreased (p = 0.004) compared to 0 km (Figure 2B). For the hip angles, the peak hip flexion angle significantly increased at 10 km compared to both 0 km (p = 0.002). The peak angle of hip flexion was significantly greater during the 20 km than the 10 km (p < 0.001). On the contrary, the peak hip extension angle was significantly increased during the 10 km rather than the 0 km (p = 0.024) and 20 km (p = 0.003) running stance phases (Figure 2C). No significant changes were observed in peak ankle dorsiflexion (F = 2.074, p = 0.143, partial η² = 0.03) during the half marathon running.

Significant differences were observed in joint moments, moment impulses, and cumulative loading variables during half marathon running (Table 4; Figure 3). A main effect of distance was noted on the peak ankle plantarflexion moment (F = 7.241, p = 0.002, partial η² = 0.10) and the peak hip flexion moment (F = 12.561, p < 0.001, partial η² = 0.17). The results indicated a significant decrease in the peak ankle plantarflexion moment at 10 km (p < 0.001) compared to 0 km (Figure 3A). Furthermore, as the distance increased, the peak ankle plantarflexion moment significantly decreased during the 20 km (p = 0.029) compared to the 10 km distance (Figure 3A). Similarly, the peak hip flexion moment was significantly reduced (p < 0.001) at both 10 km and 20 km compared to 0 km (Figure 3C). Nonetheless, no significant differences were observed in the peak moments of ankle dorsiflexion (F = 1.518, p = 0.223, partial η² = 0.02), knee flexion moment (F = 4.620, p = 0.013, partial η² = 0.07), knee extension moment (F = 0.930, p = 0.435, partial η² = 0.01), or hip extension moment (F = 0.211, p = 0.723, partial η² = 0.01) during the half marathon.

Table 4 shows the angular impulse and cumulative impulse mechanics of the ankle, knee, and hip joints in the sagittal plane during half-marathon running. The results revealed a significant decrease in the ankle plantarflexion angular impulse at 10 km (p = 0.007) compared to 0 km (F = 6.485, p = 0.004, partial η² = 0.09) during the running stance phases (Figure 3D). However, the angular impulse for ankle plantarflexion significantly increased at 20 km when compared to 10 km (p < 0.001) (Figure 3D). There were no significant changes observed in the knee flexion stance phase angular impulse (F = 1.681, p = 0.193, partial η² = 0.03) (Figure 3E) or the hip flexion stance phase angular impulse (F = 1.455, p = 0.238, partial η² = 0.02) (Figure 3F). Additionally, the cumulative impulse for ankle plantarflexion (F = 7.181, p = 0.003, partial η² = 0.10) significantly increased with distance, showing a significant difference at 20 km compared to 0 km (p < 0.001) (Figure 3G). Conversely, the cumulative impulse for ankle plantarflexion significantly decreased at 10 km (p = 0.018) compared to 0 km as the distances increased (Figure 3G). The knee flexion cumulative impulse also increased significantly with distance (F = 4.242, p = 0.020, partial η² = 0.06), indicating that the knee flexion cumulative impulse was greater at 20 km compared to both 0 km (p = 0.034) and 10 km (p = 0.002) during the running stance phases (Figure 3H). Furthermore, the hip flexion cumulative impulse significantly increased with distance (F = 7.538, p = 0.001, partial η² = 0.11), demonstrating that the hip cumulative impulse was greater at 20 km compared to 0 km (p < 0.001) and 10 km (p = 0.009) during the running stance phases (Figure 3I).

To verify the validity of the musculoskeletal model, we compared the time- and magnitude-normalized EMG data to the modeled muscle activation. Figure 4 illustrates a comparison of muscle activation results acquired from the EMG sensor with those simulated by OpenSim across different running distances. The results indicated that the activation levels of the six selected muscles were predominantly consistent, demonstrating that the musculoskeletal model established in this study exhibits a significant level of reliability.

The main effects of running distance on peak force, peak strain, force impulse, strain impulse, cumulative weight force impulse, and cumulative weight strain impulse for the AT were identified (Table 5; Figure 5). A significant main effect of running distance was observed on the AT peak force (F = 6.711, p = 0.002, partial η² = 0.10), AT peak strain (F = 56.523, p < 0.001, partial η² = 0.47), AT peak force impulse (F = 14.276, p < 0.001, partial η² = 0.19), and AT peak strain impulse (F = 30.005, p < 0.001, partial η² = 0.32) during the half marathon. The peak AT force significantly diminished at 10 km (p < 0.001) in comparison to 0 km. Conversely, the AT force significantly increased after 20 km (p = 0.021) when compared to 10 km (Figure 5A). The AT peak strain exhibited a significant decrease with increasing distances, particularly at 10 km (p < 0.001) and 20 km (p < 0.001) during the running stance phases compared to 0 km (Figure 5B). However, at 20 km, the AT peak strain significantly increased (p < 0.001) compared to 10 km. Regarding the AT force impulse, a significant decrease was noted at 10 km compared to 0 km (p = 0.015) (Figure 5D). The AT force impulse was significantly greater at 20 km compared to both 0 km (p = 0.011) and 10 km (p < 0.001) during the running stance phases (Figure 5D). Additionally, a significant reduction in AT peak strain impulse was found at 10 km compared to both 0 km (p < 0.001) and 20 km (p < 0.001) (Figure 5E). There were significant main effects of distance on the cumulative weight AT force impulse (F = 6.415, p = 0.003, partial η² = 0.09) and the cumulative AT strain impulse (F = 51.638, p < 0.001, partial η² = 0.45) (Figure 5G). The results indicated that with increased running distance, the cumulative weight AT force impulse significantly declined during the 10 km run compared to 0 km (p = 0.001) and 20 km (p = 0.005) (Figure 5G). The cumulative AT strain impulse significantly decreased during the 10 km (p < 0.001) and 20 km (p < 0.001) runs compared to 0 km (Figure 5H). Furthermore, the cumulative AT strain impulse significantly increased at 20 km (p < 0.001) compared to 10 km (Figure 5H).

The main effects of running distance on peak stress, stress impulse, and CD stress impulse were observed in the PFJ during the half marathon (Table 5; Figure 5). Specifically, both the peak stress impulse (F = 6.645, p = 0.003, partial η² = 0.10) and the CD of the PFJ impulse (F = 14.428, p < 0.001, partial η² = 0.19), demonstrated significant differences throughout the half- marathon. As the distance increased, the peak PFJ stress impulse was significantly greater at 10 km (p = 0.037) and 20 km (p = 0.002) compared to 0 km (Figure 5F). Furthermore, the CD of PFJ stress was also significantly greater at 10 km (p < 0.001) and 20 km (p < 0.001) in comparison to 0 km (Figure 5I). However, no significant differences were found in the PFJ peak stress variables across the three checkpoint distances (F = 1.495, p = 0.229, partial η² = 0.02) (Figure 5C).