The running pattern was defined to be multifactorial and as being the product of overall action of human body as early as 1985 (Subotnick, 1985). Indeed, foot placement, arm swing, body angle, rear leg lift, and stride length were suggested to be considered together (Subotnick, 1985). The running pattern was also described as a global and dynamic system (Gindre et al., 2016). For this reason, the running pattern has been analyzed globally by some researchers. For instance, McMahon et al. (1987) defined running with increased knee flexion and long ground contact time (tc) as Groucho running while Arendse et al. (2004) defined running with aligned acromion, greater trochanter, and lateral malleolus as well as short tc as Pose running. As another example, running with either a midfoot or forefoot strike pattern, short stride length, and with the body slightly leaning forward has been named Chi running (Dreyer and Dreyer, 2009).

More recently, the synthetic review of van Oeveren et al. (2021) proposed that the full spectrum of running patterns could be described by combining two temporal variables: step frequency (SF) and duty factor (DF). The DF variable represents the product of tc and stride frequency, where stride frequency is equivalent to approximately half of SF, with less than 4% differences in step times between right and left sides seen in competitive, recreational, and novice runners between 8 and 12 km/h (Mo et al., 2020). Hence, DF reflects the relative contribution of tc to the running stride (Minetti, 1998; Folland et al., 2017). According to van Oeveren et al. (2021), knowing DF and SF allows to categorize running patterns in one of five distinct categories, namely “stick”, “bounce”, “push”, “hop”, and “sit”, but keeping in mind that running patterns operate along a continuum.

The importance of DF and SF in determining running patterns corroborates previous findings. For instance, Beck et al. (2020) showed that DF is functionally representative of global biomechanical behavior, considering the duration of force production (which takes place during tc) and its cycle frequency (stride frequency). Moreover, DF was used to categorize runners with distinct running patterns (Lussiana et al., 2019; Patoz et al., 2019; Patoz et al., 2020). High and low DF runners were shown to use different running strategies (Lussiana et al., 2019; Patoz et al., 2020). Indeed, low DF runners exhibited a more symmetrical running step, anterior (midfoot and forefoot) strike pattern, and extended lower limb during tc than high DF runners, whereas the latter exhibited greater lower limb flexion during tc, more rearfoot strike pattern, and less vertical oscillation of the whole-body center of mass (COM) to promote forward propulsion than low DF runners (Lussiana et al., 2019; Patoz et al., 2020). Despite these spatiotemporal and kinematic differences, the two DF groups demonstrated similar running economy, indicating the two strategies are energetically equivalent at endurance running speeds (Lussiana et al., 2019). This would suggest that the two DF groups may optimize differently their running pattern, i.e., high DF runners promotes forward propulsion (pulley system) whereas low DF runners optimized the spring-mass model (Lussiana et al., 2019). This statement was further explored by investigating the relationships between DF and force-length relationship and leg stiffness (kleg).

In relation to SF, this variable can reveal individual strategies to increase running speed (Dorn et al., 2012) or achieve top-end running speeds (Salo et al., 2011). Indeed, the consistency in SF was shown to decrease as speed differences increased (tested running speeds: 10–18 km/h) (Patoz et al., 2022) and each runner was shown to self-optimize his step length over SF ratio (Hunter et al., 2017; van Oeveren et al., 2021). Even in subgroups of individuals with similar sprint velocities, a range of SF and step length combinations are present (Hunter et al., 2004). In addition, SF was shown to be more variable in novice than expert runners, independently of the running speed (10 and 15 km/h) (Fadillioglu et al., 2022). Furthermore, Bonnaerens et al. (2021) demonstrated that external forces were lower in recreational runners that run with higher DF and SF values (although non-significant for SF).

These previous studies investigated the association of DF or SF on running biomechanics using summary metrics, i.e., specific temporal focus like foot-strike, mid-stance, or toe-off, of signals such as the whole-body COM trajectory or the lower limb angles during tc (Lussiana et al., 2019; Patoz et al., 2020). This reduction to summary-metric space is not strictly necessary because statistical hypothesis testing can also be conducted in a continuous manner (Pataky, 2012). Indeed, one-dimensional biomechanical curves such as the ground reaction force signals are registrable and their fluctuations can be described and, then, compared expressing them as a function of the normalized stance phase duration (Cavanagh and Lafortune, 1980; Sadeghi et al., 2003). In this case, statistical analysis can be conducted on the original registered curves using statistical parametric mapping (SPM) (Friston et al., 2007), which was recently applied to the field of biomechanics (Pataky, 2010). SPM has the advantages to consider the signal as a whole and presents the results directly in the original sampling space. For this reason, the spatiotemporal biomechanical context is immediately apparent, and allows direct visualization of where do significant differences occur during tc (Pataky, 2012).

Therefore, the first purpose of the present study was to investigate the association of DF and SF on the vertical and fore-aft ground reaction force signals for treadmill runs at several endurance running speeds using SPM. The second purpose of this study was to investigate the association of DF and SF on variables derived from the vertical and fore-aft ground reaction force signals, i.e., impact (Fz,impact), active (Fz,max), braking (Fbrake,min), and propulsive (Fprop,max) peaks (Luo et al., 2019). Besides, force related variables that additionally consider the temporal aspect of the running step were also considered because the latter can vary with DF and SF. The third purpose of the present study was to investigate the association of DF and SF on the force-length relationship (Gill et al., 2020) and kleg (Liew et al., 2017).

We hypothesized that 1) a lower DF should be associated to higher vertical and fore-aft ground reaction force fluctuations, and that a lower SF should be associated to higher vertical and fore-aft ground reaction force fluctuations but to a lower extent than for DF (Bonnaerens et al., 2021). Moreover, we hypothesized that 2) both a lower DF and lower SF should be associated to higher peak forces (Fz,impact, Fz,max, Fbrake,min, Fprop,max). Besides, higher DF runners demonstrated a more rearfoot strike pattern (Lussiana et al., 2019; Patoz et al., 2020) but should show lower vertical force than lower DF runners. Hence, we hypothesized that 3) the linearity of the force-length relationship should decrease with increasing DF, due to the higher chance to observe an impact peak when increasing DF, and that a higher DF should be associated to a lower kleg. Furthermore, we hypothesized that 4) a higher SF should correspond to a greater kleg and smaller leg compression, as previously observed (Morin et al., 2007; Coleman et al., 2012; Hobara et al., 2020).

An existing database of 115 recreational runners (Patoz et al., 2021) including 87 males (age: 30 ± 8y, height: 180 ± 6 cm, leg length, measured from motion capture: 86 ± 4 cm, body mass: 70 ± 7 kg, weekly running distance: 38 ± 24 km, and running experience: 10 ± 8y) and 28 females (age: 30 ± 7 years, height: 169 ± 5 cm, leg length: 82 ± 4 cm, body mass: 61 ± 6 kg, weekly running distance: 22 ± 16 km, and running experience: 11 ± 8y) was used in this study. For study inclusion, participants were required to not have current or recent lower-extremity injury (≤1 month), to run at least once a week, and to have an estimated maximal aerobic speed ≥14 km/h (individual estimation). The study protocol was approved by the local Ethics Committee (CER-VD 2020–00334).

After the participants provided written informed consent, retroreflective markers were positioned on the participants (described in Subsec. Data Collection) to record their running biomechanics. For each participant, a 1-s static trial was first recorded while he or she stood in a standard anatomical position on an instrumented treadmill (Arsalis T150–FMT-MED, Louvain-la-Neuve, Belgium) for calibration purposes. Then, a 7-min warm-up run was performed on the same treadmill (9–13 km/h). After a short break (<5min), three 1-min runs (9, 11, and 13 km/h) were performed in a randomized order (1-min recovery between each run). Three-dimensional (3D) kinematic and kinetic data were collected during the static trial and the last 30s of the running trials (83 ± 5 running steps), resulting in more than 20 steps being analyzed (Riazati et al., 2019). All participants were familiar with running on a treadmill, as it was part of their usual training program, and they wore their habitual running shoes (shoe mass: 256 ± 48 g and shoe heel-to-toe drop: 7 ± 3 mm).

Whole-body 3D kinematic data were collected at 200 Hz using motion capture (8 cameras) and Vicon Nexus software v2.9.3 (Vicon, Oxford, United Kingdom). The laboratory coordinate system was oriented such that the x-, y-, and z-axes denoted the mediolateral (pointing towards the right side of the body), posterior-anterior, and inferior-superior axes, respectively. Forty-three and 39 retroreflective markers of 12.5 mm diameter were used for the static and running trials, respectively. They were affixed to the skin and shoes of individuals on anatomical landmarks using double-sided tape following standard guidelines (Tranberg et al., 2011). Synchronized kinetic data (1000 Hz) were also collected using the force plate embedded into the treadmill.

The 3D marker and ground reaction force data (analog signal) were exported in the. c3d format and processed in Visual3D Professional software v6.01.12 (C-Motion Inc., Germantown, MD, United States). The 3D marker data were interpolated using a third-order polynomial least-square fit algorithm, allowing a maximum of 20 frames for gap filling, and were subsequently low-pass filtered at 20 Hz using a fourth-order Butterworth filter. The 3D ground reaction force signal was filtered using the same filter, and down sampled to 200 Hz to match the sampling frequency of the marker data.

A full-body biomechanical model with six degrees of freedom and 15 rigid segments was constructed from the marker set. The segments included the head, upper arms, lower arms, hands, thorax, pelvis, thighs, shanks, and feet. In Visual3D, the segments were treated as geometric objects, assigned inertial properties and COM locations based on their shape (Hanavan, 1964), and attributed relative masses based on standard regression equations (Dempster, 1955). The whole-body COM location was calculated from the parameters of all 15 segments (the whole-body COM was directly provided by Visual3D).

For all biomechanical measures, the values extracted from the 30-s data collection for each participant, including both right and left steps, were averaged for subsequent statistical analyses.

For each running trial, foot-strike, toe-off, and mid-stance events were identified with Visual3D. Foot-strike and toe-off were detected by applying a 20N threshold to the vertical ground reaction force (Smith et al., 2015). Mid-stance was placed at the instant where the fore-aft ground reaction force changed from negative to positive, which permitted to separate the stance phase in a braking and propulsive phase.

Temporal Variables tc, flight time (tf), and swing time (ts) were defined as the time from foot-strike to toe-off of the same foot, from toe-off of one foot to foot-strike of the contralateral foot, and from toe-off to foot-strike of the same foot, respectively.

DF was calculated as DF=tc/tc+ts, where 1/tc+ts represents the stride frequency (Minetti, 1998). SF was defined as the inverse of the sum of tc and tf, i.e., SF=1/tc+tf. Furthermore, SF was normalized by g/L0 (Lieberman et al., 2015; van Oeveren et al., 2021), where g is the gravitational constant and L0 the leg length, calculated as the distance between hip and ankle joint center using the static (calibration) trial.

Braking (tbrake) and propulsive (tprop) times were given as the time from foot-strike to mid-stance and mid-stance to toe-off of the same foot, respectively.

Compression (tcomp) and decompression (tdecomp) times were given as the time from foot-strike to the time where the vertical position of the whole-body COM is at its minimum, i.e., where the vertical ground reaction force is maximum, and from the time where the vertical position of the whole-body COM is at its minimum to toe-off, respectively.

Fz,impact and Fz,max were obtained from the vertical ground reaction force signal (Luo et al., 2019). Noteworthy, an impact peak was not always observed, Fz,impact was quantified in 80% of the running trials. Besides, Fbrake,min and Fprop,max were given by the minimum and maximum values of the fore-aft ground reaction force signal (Luo et al., 2019).

The instantaneous vertical loading rate (LRz) was calculated as the largest slope of the vertical ground reaction force signal between 20 and 80% of the first 15% of the stance phase (Willson et al., 2014). The 15% limit was chosen because an impact peak was not always identified and so that the loading rate of every runner was in the same relative temporal window (Willson et al., 2014). The braking (LRbrake) and propulsive (LRprop) loading rates, because of their relation to running-related injuries (Daoud et al., 2012; Willson et al., 2014; Davis et al., 2016; Johnson et al., 2020), were calculated as the largest slopes of the fore-aft ground reaction force signal between foot-strike and the instant of Fbrake,min and between mid-stance and the instant of Fprop,max, respectively.

The braking (Ibrake) and propulsive (Iprop) impulses were calculated as the integral of the fore-aft ground reaction force signal from foot-strike to mid-stance and from mid-stance to toe-off, respectively (Gottschall and Kram, 2005).

Force variables were all normalized by BW.

The spring-mass characteristics of the lower limb were assessed by computing the force-length relationship (Gill et al., 2020), i.e., the force vector projected along the leg as function of the leg compression/decompression during stance, and kleg (Liew et al., 2017), calculated using both the compression and decompression of the human body (Gill et al., 2020) and adapted from Liew et al. (2017). More explicitly, compressive (kleg, comp) and decompressive (kleg, decomp) leg stiffnesses were given by the maximum of the force vector projected along the leg (Fleg,max) divided by the maximum leg compression (∆Lcomp) and decompression (∆Ldecomp) during stance, respectively. Following the definition of the spring-mass model, i.e., a massless spring attached to a point mass located at the whole-body COM (Blickhan, 1989), the leg length was represented by the magnitude of a 3D leg vector defined from the whole-body COM to the center of pressure of the foot. The center of pressure being subject to large fluctuations for low vertical force values, a 200N vertical threshold was used for foot-strike and toe-off events in this specific case (see supplementary materials). ∆Lcomp and ∆Ldecomp were given by the difference between the leg length at foot-strike and the minimum value of the leg length and by the difference between the leg length at toe-off and the minimum value of the leg length, respectively. A leg angle (θleg) was calculated as the angle between the leg vector and anterior-posterior axis, and evaluated at foot-strike (θleg,FS) and toe-off (θleg,TO) (Coleman et al., 2012).

Fleg,max was normalized by BW, ∆Lcomp and ∆Ldecomp were expressed in absolute and relative (as a percentage of participant’s height) units and similarly for kleg, comp and kleg, decomp.

All data are presented as the mean ± standard deviation. Pearson correlation coefficient (r) between DF and SF together with corresponding 95% confidence interval (lower, upper) were computed at the three running speeds separately. Correlations were considered very high, high, moderate, low, and negligible when absolute r values were between 0.90–1.00, 0.70–0.89, 0.50–0.69, 0.30–0.49, and 0.00–0.29, respectively (Hinkle et al., 2002). In this study, collinearity between DF and SF was prevented because r was smaller than 0.7 (Table 1) (Van Oeveren et al., 2019). The association of DF and SF on the vertical and fore-aft ground reaction force signals (along the entire stance phase) was examined using SPM and linear regression for each tested speed. Bonferroni correction was employed to consider that three running speeds were tested. To compare participants, the stance phase was normalized and therefore expressed in percentage. Besides, residual plots were inspected and no obvious deviations from homoscedasticity or normality were observed. Hence, the association of DF and SF (covariates) and running speed on temporal, ground reaction force, and stiffness related variables was evaluated using a linear mixed effects model fitted by restricted maximum likelihood. The within-subject nature was controlled for by including random effects for participants (individual differences in the intercept of the model). The fixed effects included running speed (categorical variable) and DF and SF (continuous variables). The linearity of the force-length relationship was quantified using the coefficient of determination (R2) during both leg compression (Rcomp2) and decompression (Rdecomp2). However, the calculation of R2 was modified so that Rcomp2 and Rdecomp2 values were computed by comparing the compression and decompression force-length relationships to the perfectly elastic compression and decompression lines, i.e., linear relations obtained using slopes equal to kleg, comp and kleg, decomp, respectively (Gill et al., 2020). In other words, R2 evaluates how far the force-length relationship is from a linear model obtained using kleg. Statistical analysis was performed using spm1D (v0.4.6, https://spm1d.org) (Pataky, 2012), Python (v3.7.4, http://www.python.org), and Jamovi (v1.6.23, https://www.jamovi.org) with a level of significance set at p ≤ 0.05.

TABLE 1. Duty factor (DF) and step frequency (SF), as well as their Pearson’s correlation coefficient (r) together with their 95% confidence interval (lower, upper) and statistical significance (P ≤ 0.05), indicated in bold, for three tested speeds.

The increase of the running speed from 9 to 13 km/h was accompanied with a decrease of DF of 13.3 ± 3.8% and an increase of SF of 5.9 ± 3.1%. The correlation between DF and SF was low but significant at all tested speeds (r ≤ 0.32; p < 0.001; Table 1).

tc, tbrake, tprop, tcomp, and tdecomp significantly increased with increasing DF while tf decreased (p < 0.001; Table 2). These six variables significantly decreased with increasing SF (p < 0.001; Table 2). Besides, tc and tdecomp decreased with increasing speed while tf increased with increasing speed (p ≤ 0.02; Table 2).

TABLE 2. Temporal variables for runners at endurance running speeds. Significant differences (P ≤ 0.05) identified by linear mixed effects modeling are indicated in bold. Note: values are presented as mean ± standard deviation. DF: duty factor, SF: step frequency, tc: contact time, tbrake: brake time, tprop: propulsion time, tcomp: compression time, tdecomp: decompression time, and tf: flight time. SF covariate was normalized by g/L0, where g is the gravitational constant and L0 the leg length. Up (↑) and down (↓) arrows indicate positive and negative effects of the covariate, respectively. † and ‡ Significantly different from the value at 13 km/h.

The vertical ground reaction force signal was significantly negatively related to DF at all tested speeds (stance range: 0 and 15–100% at 9 and 11 km/h, and 0 and 14–100% at 13 km/h; Figure 1). Similar findings were obtained for SF but to a lower extent (stance range: 60–99% at 9 km/h, 59–99% at 11 km/h, and 67–83% at 13 km/h; Figure 2).

FIGURE 1. Statistical parametric mapping (SPM) analysis, i.e., t-statistics (SPM), of the linear relationship between the vertical ground reaction force (Fz) and the duty factor (DF) along the running stance phase at (A) 9 km/h, (B) 11 km/h, and (C) 13 km/h. In the upper panels, Fz, expressed in body weight (BW), is depicted for each participant (the color depends on the DF value) and for the mean (black line) ± standard deviation (dashed black line) over all participants. In the lower panels, the black dashed horizontal lines represent the critical (parametric) threshold while the portion of the running stance phase which is statistically significant (p ≤ 0.017; Bonferroni correction was applied to take into the three tested speeds) is given by the gray shaded area.

FIGURE 2. Statistical parametric mapping (SPM) analysis, i.e., t-statistics (SPM), of the linear relationship between the vertical ground reaction force (Fz) and the step frequency (SF) normalized by g/L0, where g is the gravitational constant and L0 the leg length, along the running stance phase at (A) 9 km/h, (B) 11 km/h, and (C) 13 km/h. In the upper panels, Fz, expressed in body weight (BW), is depicted for each participant (the color depends on the DF value) and for the mean (black line) ± standard deviation (dashed black line) over all participants. In the lower panels, the black dashed horizontal lines represent the critical (parametric) threshold while the portion of the running stance phase which is statistically significant (p ≤ 0.017; Bonferroni correction was applied to take into the three tested speeds) is given by the gray shaded area.

The fore-aft ground reaction force signal was significantly positively related to both DF and SF in the first 50% of the stance (negative fore-aft force) and negatively related to both DF and SF in the last 50% of the stance at all tested speeds (stance range for DF: 5–11, 27–34, and 69–100% at 9 km/h, 7–12, 29–35, and 71–100% at 11 km/h, and 6–13 and 68–100% at 13 km/h; Figure 3; stance range for SF: 15–33 and 68–95% at 9 km/h, 14, 19–35, 47–52, and 70–98% at 11 km/h, and 14–28 and 71–89% at 13 km/h; Figure 4).

FIGURE 3. Statistical parametric mapping (SPM) analysis, i.e., t-statistics (SPM), of the linear relationship between the fore-aft ground reaction force (Fy) and the duty factor (DF) along the running stance phase at (A) 9 km/h, (B) 11 km/h, and (C) 13 km/h. In the upper panels, Fy, expressed in body weight (BW), is depicted for each participant (the color depends on the DF value) and for the mean (black line) ± standard deviation (dashed black line) over all participants. In the lower panels, the black dashed horizontal lines represent the critical (parametric) threshold while the portion of the running stance phase which is statistically significant (p ≤ 0.017; Bonferroni correction was applied to take into the three tested speeds) is given by the gray shaded area.

FIGURE 4. Statistical parametric mapping (SPM) analysis, i.e., t-statistics (SPM), of the linear relationship between the fore-aft ground reaction force (Fy) and the step frequency (SF) normalized by g/L0, where g is the gravitational constant and L0 the leg length, along the running stance phase at (A) 9 km/h, (B) 11 km/h, and (C) 13 km/h. In the upper panels, Fy, expressed in body weight (BW), is depicted for each participant (the color depends on the DF value) and for the mean (black line) ± standard deviation (dashed black line) over all participants. In the lower panels, the black dashed horizontal lines represent the critical (parametric) threshold while the portion of the running stance phase which is statistically significant (p ≤ 0.017; Bonferroni correction was applied to take into the three tested speeds) is given by the gray shaded area.

Fz,max, Fz,impact, Fprop,max, LRprop, and Ibrake significantly decreased with increasing DF while Fbrake,min, LRbrake, and Iprop significantly increased (p ≤ 0.01; Table 3). Fz,impact, Fbrake,min, LRprop, and Ibrake significantly increased with increasing SF while, Fprop,max and Iprop significantly decreased (p ≤ 0.02; Table 3). Considering absolute values, all the ground reaction force variables significantly increased with increasing speed (p ≤ 0.005; Table 3).

TABLE 3. Ground reaction force variables for runners at endurance running speeds. Significant differences (P ≤ 0.05) identified by linear mixed effects modeling are indicated in bold.

The force-length relationships of all participants, colored according to their DF and SF, are depicted in Figures 5, 6, respectively, for each tested speeds and separately for the compression and decompression phases. Rcomp2 significantly decreased with increasing DF or running speed, and increased with increasing SF (p ≤ 0.007; Table 4), while there was no change of Rdecomp2 with DF, SF, and speed.

FIGURE 5. Force-length relationship, i.e., ground reaction force projected along the leg (Fleg) as function of the leg compression/decompression, for each participant [the color depends on the duty factor (DF) value] and for the mean (black line) over all participants during the running stance phase, at three running speeds, and expressed using (A) SI units and (B) normalized units, i.e., body weight (BW) for Fleg and percentage of runners’ height for leg compression.

FIGURE 6. Force-length relationship, i.e., ground reaction force projected along the leg (Fleg) as function of the leg compression/decompression, for each participant [the color depends on the step frequency (SF) value; SF was normalized by g/L0, where g is the gravitational constant and L0 the leg length] and for the mean (black line) over all participants during the running stance phase, at three running speeds, and expressed using (A) SI units and (B) normalized units, i.e., body weight (BW) for Fleg and percentage of runners’ height for leg compression.

TABLE 4. Linearity of the force-length relationship during leg compression (Rcomp2) and decompression (Rdecomp2). Significant differences (p ≤ 0.05) identified by linear mixed effects modeling are indicated in bold.

Fz,max, ∆Lcomp, kleg, comp, and kleg, decomp significantly decreased with increasing DF while ∆Ldecomp and θleg,TO significantly increased (p ≤ 0.03; Table 5). ∆Lcomp, ∆Ldecomp, and θleg,TO significantly decreased with increasing SF while kleg, comp, and kleg, decomp significantly increased (p < 0.001; Table 5). Fleg,max, ∆Ldecomp, kleg, decomp, θleg,FS |, and θleg,TO significantly increased with increasing speed (p ≤ 0.005; Table 5).

TABLE 5. Stiffness variables for runners at endurance running speeds. Significant differences (P ≤ 0.05) identified by linear mixed effects modeling are indicated in bold.

DF was previously analytically shown to be inversely proportional to the maximum of an approximated, based on a sine-wave model (Beck et al., 2020), vertical ground reaction force signal (Morin et al., 2005). This previous knowledge is further expanded by the present results which showed that Fz,max is significantly negatively related to DF (Table 3), and corroborates previous findings which showed that DF was negatively correlated to Fz,max (Bonnaerens et al., 2021). This suggests that DF should be inversely related to Fz,max without using a sine-wave model to approximate the vertical ground reaction force. Moreover, the present study extends to the fact that a lower DF results in a larger vertical ground reaction force during most of the stance (∼15–100%; Figure 1) but after the 15% temporal window representative of the “impact” phase (Willson et al., 2014). Therefore, the SPM analysis additionally revealed that the association between DF and the vertical ground reaction force signal is not only given at Fz,max but through almost the entire stance (after the impact phase; ≥15%). This result suggests that the shape of the vertical ground reaction force during the impact phase is not affected by the DF. This result might be attributed to the fact that the vertical ground reaction force signal is given by the force contributions of two discrete body mass components, i.e., a distal mass composed of the foot and shank and the remaining mass (Clark et al., 2017; Udofa et al., 2019). Hence, the impact phase, represented by the distal mass in this model, might not be affected by the DF. This study also showed that for the runners having a visible Fz,impact, this Fz,impact was significantly larger for lower than higher DF values (Table 3). This discrepancy could be explained by fact that the impact peak might be happening at a different instant of the running stance phase (within the first 15%) depending on individuals.

Similarly, a lower SF resulted in a larger vertical ground reaction force, but only at the end of the stance (∼65–95%; Figure 2). This was likely related to the longer step length for running at the same speed. In fact, it has previously been shown that a larger vertical ground reaction force (i.e., support force) produces a larger step length (Weyand et al., 2000; Dorn et al., 2012). Our SPM analysis demonstrated that this larger support force was located only at the end of the stance. Indeed, Fz,max was not related to SF (Table 3). Non-etheless, the reason why this larger support force was located at the end of the stance could not readily be explained. Besides, Fz,impact significantly increased with increasing SF (Table 3). This result contradicts previous findings which observed a decrease of the impact peak with increasing SF (Lieberman et al., 2015). However, these findings were obtained when asking individuals to voluntarily increase their SF. Hence, this could lead to a different running pattern than the spontaneous running pattern of runners with a naturally high SF.

The present study reported no association of DF and SF on LRz (Table 3). This could partly follow from the fact that the SPM analysis did not report any significant association between DF and SF and the vertical ground reaction force during the impact phase (the first 15% of the stance). This result corroborates the absence of correlation between LRz and both DF and SF at slow running speeds, as reported by Bonnaerens et al. (2021). However, assuming that DF is partly related to foot-strike pattern, i.e., the higher the DF, the more likely that this runner is a rearfoot striker (Lussiana et al., 2019; Patoz et al., 2020), this result contradicts the result of a meta-analysis which reported higher LRz for rearfoot than non-rearfoot strikers (Almeida et al., 2015).

The fore-aft ground reaction force signal was positively related to DF around ∼5–10% of the stance and to both DF and SF around ∼25–35% (positively) and ∼70–90% (negatively; Figures 3, 4). The positive association of DF on the fore-aft ground reaction force signal reported by the SPM analysis around ∼5–10% of the stance can be explained by the foot-strike pattern. Indeed, fore-foot strikers were shown to have a negative spike on the fore-aft ground reaction force signal around ∼5–10% of the stance (Nordin et al., 2017) and DF was related to the footstrike pattern (Lussiana et al., 2019; Patoz et al., 2020). However, the association of DF on the fore-aft force signal around ∼5–10% of the stance was not accompanied by an association of DF on the vertical force signal at the same percentage of the stance. This suggests that the effect of DF during the impact phase was more important in the fore-aft than vertical force signal. The other two significant regions are around the braking and propulsive peaks (Fbrake,min and Fprop,max), which were also significantly related to DF and SF (Table 3). These results partly corroborate previous observations, which showed that the peak braking force was correlated to DF but not to SF (Bonnaerens et al., 2021). Moreover, they confirm that larger ground reaction forces during propulsion are needed to lift and accelerate the body during stance to generate longer step lengths (Schache et al., 2014). As previously suggested (van Oeveren et al., 2021), combining vertical and horizontal ground reaction forces into a single vector could be useful to properly characterize their orientations and actions and carefully describe the relationship of this single vector with DF and SF, especially at the end of the stance.

The linearity of the force-length relationship was higher for lower DF and SF than for higher DF and SF runners during the leg compression but there was no difference during the leg decompression (Table 4). This means that higher DF and SF values were associated to more variations of the instantaneous compressive stiffness, i.e., the slope for each pair of point during the leg compression. However, the decompressive stiffness during the leg decompression was independent of DF and SF (Table 4). This result corroborates the choice made by several authors to use the decompression phase instead of the compression one to calculate the vertical stiffness (Cavagna et al., 1988; Schepens et al., 1998). Deviation from linearity of the force-length relationship among individuals was also reported by Gill et al. (2020). Indeed, these authors reported that the linearity of the force-length curve was foot-strike index (foot-strike pattern) dependent and that this curve should be investigated before using the spring-mass model. Furthermore, these authors suggested that for R2<0.95, it may be more appropriate to segment the stance phase and to individually investigate the different subphases. Hence, the deviation from linearity observed herein during the leg compression for higher than lower DF and SF runners suggest that the stiffness should be split into several phases during the leg compression and thus invalidate the usage of kleg, comp for these runners. However, the linearity observed during leg decompression for all participants suggest that kleg, decomp could be used.

This study reported that kleg, decomp and Fleg,max significantly increased with decreasing DF while ∆Ldecomp decreased (Table 5). Hence, the elastic energy (Eel), which could be calculated as Eel=Fleg,max2kleg,decomp2 using the definition of the spring-mass model, increased with decreasing DF. Furthermore, the compression was more vertical and tc was shorter for lower than higher DF runners. High DF runners could be characterized by a slow stretch-shortening cycle (runners with tc longer than 250 m) while low DF runners by a fast one (Vogt and Hoppeler, 2014). These results, together with the higher linearity of the force-length curve during the compression phase observed for lower than higher DF runners suggest that a lower DF runner better optimizes the spring-mass model than a higher DF runner. On the contrary, θleg,TO significantly increased with increasing DF (Table 5). These results suggest that the higher the DF, the higher the promotion of forward propulsion of the body. This compensates for the lower utilization of the spring-mass model of higher than lower DF runners and corroborates previous findings (Lussiana et al., 2019; Patoz et al., 2020). These findings bring further evidence and reinforce previous statements that low DF runners rely more on the optimization of the spring-mass model whereas high DF runners promotes forward propulsion (pulley system) (Lussiana et al., 2019; Patoz et al., 2020).

This study further revealed that the higher the SF, the larger kleg, decomp and the smaller ∆Ldecomp (Table 5), which corroborates previous findings (Coleman et al., 2012; Hobara et al., 2020). Moreover, θleg,TO significantly decreased with increasing SF (Table 5). These results suggest that higher SF runners better optimize the spring-mass model than lower SF runners, confirming that SF seems to be an indirect factor influencing kleg through its effect on tc (Morin et al., 2007).

Most of the variables studied herein reported an opposite association of DF and SF covariates (Tables 3-5). In other words, for most of the variables, if DF had a positive association on a given variable, then SF had a negative association on the same variable, and vice versa. This observation sounds counter-intuitive because SF is analytically associated to DF, i.e., DF=0.5 tc SF. However, though significant, correlations between DF and SF were low at all tested speeds (Table 1). Hence, the direct association of SF covariate on a given variable is more important than the indirect association caused by the relationship between SF and DF. Besides, the low correlations between SF and DF tend to reduce the direct association of a covariate on a given variable. Noteworthy, correlations between DF and tc were high (r ≥ 0.78) and significant (p < 0.001) and correlations between DF and tf were very high (r ≥ 0.95) and significant (p < 0.001). Hence, these results corroborate that DF and SF can be viewed as two variables that complement each other and that should be used together to describe the full spectrum of running patterns (van Oeveren et al., 2021).

A few limitations to the present study exist. Few findings of this study were obtained using the spring-mass model, which include many assumptions and limitations (Blickhan, 1989; McMahon and Cheng, 1990; Farley and González, 1996) that may restrict our conclusion on the underlying mechanisms. However, due to the methodological challenges associated to in vivo measurements under dynamic conditions to understand the role of muscle-tendon unit during running, the use of spring-loaded inverted pendulum model seems rational and relevant. In addition, the running speeds were limited to endurance speeds representative of the running speeds employed by recreational runners during endurance running training (Selinger et al., 2022) and experimental trials were performed on a treadmill. Similar results might also be obtained using overground running trials because spatiotemporal parameters between motorized treadmill and overground running are largely comparable (Van Hooren et al., 2020). However, it was also concluded that participants behaved differently when attempting to achieve faster speeds overground than on a treadmill (Bailey et al., 2017). Therefore, further studies should investigate the association of DF and SF on running kinetic using additional conditions, i.e., faster speeds, positive and negative slopes, and different types of ground.

To conclude, this study revealed that the lower the DF and the higher the SF, the more the runner relies on the optimizatio