Foot strike pattern has long been a central topic in running biomechanics. Historically, humans running barefoot predominantly adopted a non-rearfoot strike (NRFS; midfoot or forefoot strike), while the widespread use of cushioned running shoes has led most runners to adopt a rearfoot strike (RFS) pattern characterized by initial heel contact (Bramble and Lieberman, 2004; Lieberman et al., 2010). At present, only about 25% of runners use NRFS, while the majority adopt RFS, which is typically associated with a pronounced impact peak in the vertical ground reaction force (vGRF) and has been linked to impact-related injuries (Patoz et al., 2019). By contrast, NRFS runners exhibit a smoother landing pattern without a distinct impact peak. Although NRFS may increase loading on the Achilles tendon and triceps surae due to greater ankle plantarflexion (Kulmala et al., 2013), it has been suggested to reduce the risk of patellofemoral pain and tibial stress syndrome and may improve running economy under certain conditions; however, evidence remains mixed (Bovalino and Kingsley, 2021).

Gait-retraining interventions conducted over 4–12 weeks, particularly under barefoot or minimalist footwear conditions, have enabled habitual rearfoot strikers (HRFS) to adopt a non-rearfoot strike (NRFS) pattern (Deng et al., 2024; Warne et al., 2014; Yang et al., 2020; Zhang et al., 2021). These interventions have been shown to alter foot-strike angle, vertical loading rates, plantar pressure, joint kinetics, and muscle–tendon mechanics. HRFS can also acutely adopt an NRFS pattern and reproduce some of the joint kinematic and running-mechanics characteristics observed in habitual non-rearfoot strikers (HNRFS) (Valenzuela et al., 2015; Williams et al., 2000). However, these adaptations appear to be outcome-specific, and successful modification of foot-strike pattern does not necessarily produce uniform changes across biomechanical variables. Thus, although gait retraining and acute foot-strike modification may enable HRFS to reproduce some observable biomechanical characteristics of HNRFS, it remains unclear whether the underlying temporal organization and dynamic stability of running also become comparable.

Running is a nonlinear dynamical system characterized by small but persistent stride-to-stride fluctuations, whose temporal organization may provide information beyond mean kinematics or the magnitude of variability alone (Hunter et al., 2023; Stergiou and Decker, 2011). Sample entropy (SampEn) quantifies the regularity and predictability of a time series, with higher values indicating less repetitive and less predictable movement patterns; however, its direction should not be interpreted as inherently beneficial or detrimental (Yentes and Raffalt, 2021). In contrast, the largest Lyapunov exponent (LyE) quantifies the rate at which initially neighboring trajectories diverge in reconstructed state space, with higher values indicating lower local dynamic stability (Dingwell et al., 2001; Rosenstein et al., 1993). These measures characterize different dynamical properties of a time series: SampEn is based on the likelihood that similar sequential patterns remain similar when their length increases, whereas LyE measures the short-term divergence of initially neighboring trajectories. Consequently, greater temporal regularity does not necessarily imply greater local dynamic stability, and the two measures need not change in parallel. Assessing both measures therefore provides complementary information about the temporal organization and local dynamic stability of movement (Crowley et al., 2021; Strongman and Morrison, 2021).

Evidence concerning nonlinear movement dynamics following foot-strike modification remains limited. Ekizos et al. (2017) reported that an acute transition from shod to barefoot running, accompanied by a more anterior foot-strike pattern, was associated with reduced local dynamic stability. Similarly, Hollander et al. (2021) found lower local dynamic stability in novice barefoot runners than in shod runners, although stability did not change over the 8-week intervention period. However, because the transition from shod to barefoot running was accompanied by a more anterior foot strike, the independent effects of footwear and foot-strike modification could not be distinguished. Furthermore, existing studies have focused predominantly on the largest Lyapunov exponent (LyE) as the primary nonlinear metric, without jointly evaluating movement regularity and local dynamic stability to characterize dynamic behavior comprehensively. Another critical gap is that prior work has not directly compared movement dynamic profiles across three distinct conditions: acute NRFS adoption, retrained NRFS, and habitual NRFS. From a motor-learning perspective, adopting an unfamiliar foot-strike pattern may involve exploratory adjustments and progressive reorganization of movement coordination (Krakauer et al., 2019; Wood et al., 2024). It therefore remains unclear whether acute and retrained NRFS adoption represent similar movement dynamic states, and how each state compares with habitual NRFS in terms of temporal organization and local dynamic stability.

Accordingly, the present study aimed to compare lower-limb joint kinematics, movement regularity, and local dynamic stability among three observed states: HNRFS, HRFS acutely adopting NRFS (HRFS-AN), and HRFS after a 12-week GR program (HRFS-GR). We hypothesized that: (1) discrete joint kinematics would not differ substantially among the three states; (2) HRFS-AN would exhibit less regular movement patterns and lower local dynamic stability than HNRFS; and (3) HRFS-GR would demonstrate movement regularity and local dynamic stability comparable with those of HNRFS.

A simulation-based sample-size analysis was conducted using the random-intercept linear mixed-effects model specified for the primary analysis (Outcome ~ State + (1 | participant) (Green and MacLeod, 2016). The calculation targeted the primary within-participant contrast between HRFS-GR and HRFS-AN, using a two-sided α of 0.05. The expected effect was based on the sagittal-plane knee SampEn change reported by Hunter et al. (2021) (approximate dz = 0.962). The analysis indicated that 11 HRFS participants with observations under both conditions were required to achieve 80% power. A total of 33 male recreational runners participated in the current study, including 19 HRFS who expressed an interest in transitioning their foot strike pattern (age: 34 ± 8 years; body mass: 68.0 ± 8.5 kg; height: 1.74 ± 0.05 m) and 14 HNRFS (age: 34.0 ± 10.2 years; body mass: 70.0 ± 8.9 kg; height: 1.73 ± 0.04 m). The inclusion criteria were as follows: (1) running at least 20 km per week in the past 6 months without changing their habitual foot strike pattern, (2) no musculoskeletal injury within the past 6 months, and (3) abstinence from alcohol and caffeine within 24 h prior to testing. All participants provided written informed consent, and the study protocol was approved by the institutional ethics committee.

During the first laboratory visit, participants ran naturally on a treadmill while wearing the experimental shoes (Nike Pegasus 34, 285 g, size 43, heel-to-toe drop: 12 mm). Their habitual foot-strike patterns were recorded in the sagittal plane using the slow-motion video mode of an iPhone 12 Pro at 240 frames·s⁻¹ and were classified frame by frame according to the foot position at initial ground contact, following the criteria described by Hasegawa et al. (2007). RFS was defined as heel first contact without simultaneous midfoot or forefoot contact, whereas NRFS was defined as midfoot or forefoot first contact or simultaneous contact of the heel and midfoot or the entire sole. The posterior one-third of the shoe was considered the heel region. Two experimenters jointly reviewed each recording and reached a consensus classification. On the basis of video analysis, runners were categorized into HRFS or HNRFS. An additional assessor separately re-evaluated all recordings using the same criteria (agreement = 90.9%; Cohen’s κ = 0.81); all disagreements involved NRFS classifications, which were subsequently confirmed by three-dimensional motion capture during the formal biomechanical testing.

Before formal testing, HRFS participants completed an individualized familiarization period at the testing speed, during which they practiced the intended NRFS pattern. Familiarization continued until they verbally reported feeling comfortable maintaining the pattern (Valenzuela et al., 2015). HNRFS participants completed the warm-up without altering their habitual foot-strike pattern.

A total of 36 retro-reflective markers (14 mm diameter) were placed bilaterally on the following anatomical landmarks: anterior superior iliac spine, posterior superior iliac spine, iliac crest, and greater trochanter (pelvis); medial and lateral femoral condyles (knee); medial and lateral malleoli (ankle); and 1st and 5th metatarsal heads, distal phalanx of the hallux, and calcaneus (foot). In addition, rigid clusters that consisted of 3 noncollinear markers were affixed onto the lateral thigh and shank of each leg.

During formal testing, the participants ran at 9 km·h⁻¹ on an instrumented treadmill (Bertec FIT, USA) while wearing the experimental shoes. The HRFS participants were instructed to adopt NRFS acutely (HRFS-AN), while HNRFS runners maintained their natural strike pattern. Data collection began only after an investigator visually confirmed a stable running pattern, operationally defined as maintaining the intended foot-strike pattern without observable abrupt changes in step frequency or stride length. Data were then recorded for 30 s. Marker trajectories were captured using an eight-camera Vicon T40 system (Oxford Metrics, UK) at 200 Hz. The HRFS participants who completed the 12-week GR protocol returned for post-tests, during which they ran with their self-selected post-training strike pattern (HRFS-GR). Data collection procedures were identical to those of baseline testing.

HRFS runners voluntarily completed a 12-week GR program designed to promote the adoption of NRFS. Minimalist shoes (Vibram FiveFingers; 3 mm rubber outsole, no midsole, zero heel–toe drop; minimalist index: 92%) were used to simulate barefoot running while providing plantar protection. Training volume was progressively increased based on individual weekly mileage (Lieberman et al., 2010; Zhang et al., 2021), with all the participants completing three minimalist shoe sessions per week. During Weeks 1–2, runners completed 10% of their weekly mileage in minimalist shoes by using NRFS, with the remaining mileage performed while wearing conventional cushioned shoes and using their habitual strike pattern. Minimalist shoe mileage was subsequently increased by 10% per week, reaching 100% by Week 11, and maintained throughout Week 12 (Figure 1). Total weekly mileage was kept constant throughout the retraining period. The participants were additionally instructed to perform foot and ankle strengthening exercises under the supervision of a rehabilitation specialist to facilitate musculoskeletal adaptation. Exercises targeted intrinsic and extrinsic foot and ankle musculature and were performed for approximately 20 min per session. Training adherence was monitored through weekly mileage logs, video submissions of foot strike patterns, and supervised training sessions. Competitive races were not permitted during the retraining period.

Kinematic data and vGRFs were processed using Visual3D (C-Motion, Inc., Rockville, MD, USA). Marker trajectories were low-pass filtered using a fourth-order Butterworth filter at 7 Hz, and vGRF signals were filtered at 50 Hz (Deng et al., 2024; Zhang et al., 2021). Initial contact was identified when the vGRF exceeded 30 N, and a gait cycle was defined as the interval between two successive initial contacts of the dominant foot. The middle 20 s (5–25 s) of each 30-s trial were retained for analysis, from which the first ten consecutive complete gait cycles of the dominant limb, defined as the preferred kicking limb, were selected.

Joint angles were calculated using Euler angles between adjacent segments. Foot strike angle was defined as the angle between the line that connected the first metatarsophalangeal joint and calcaneus markers and the global anterior–posterior axis, offset by the corresponding angle measured during static standing (Altman and Davis, 2012). Angles <8° were classified as NRFS and angles >8° as RFS. To quantify whole-body movement dynamics, the position of the center of mass (COM) was estimated using the centroid of the pelvic markers, including the bilateral anterior superior iliac spine and posterior superior iliac spine (Napier et al., 2020; Saini et al., 1998). The COM position was calculated as the mean three-dimensional coordinates of these markers.

Vertical average loading rate (VALR) and vertical instantaneous loading rate (VILR) were calculated from vGRF during early stance. A point of interest (POI) was defined as the point immediately before the vGRF slope decreased below 15 body weights per second, with the additional requirement that vGRF magnitude exceeded body weight (Futrell et al., 2020). VALR was defined as the average slope between 20% and 80% of the most linear portion of the vGRF curve between initial contact and the POI, while VILR was defined as the maximum slope between 20% and 100% of the vGRF–time curve over the same interval.

Discrete variables were extracted from 10 consecutive gait cycles and averaged for statistical analysis. Within-subject variability of the discrete joint angle measures was quantified using the coefficient of variation (CV) and calculated as the standard deviation (SD) divided by the mean across the 10 gait cycles.

Regularity and complexity of movement were quantified using SampEn, computed from the continuous joint angle time series of the hip, knee, and ankle in the sagittal, coronal, and transverse planes (McCamley et al., 2018; Richman and Moorman, 2000). In addition, whole-body movement variability was quantified using SampEn calculated from the COM position time series in the anterior–posterior, medial–lateral, and vertical directions. Time series were extracted from the 10 s of each trial, yielding 2000 data points per joint and plane at a sampling frequency of 200 Hz. SampEn was calculated as follows:

where N is the number of data points, m is the embedding dimension, and r is the similarity tolerance expressed as a proportion of the signal SD. A represents the number of vector pairs of length m + 1 within tolerance r, and B denotes the number of vector pairs of length m. For primary presentation, m = 2 was selected to reduce the dimensionality of pattern matching and maintain sufficient template matches within the 2000-point time series, while r = 0.25 × SD was selected to reduce the risk of sparse matches across the continuous joint-angle and COM signals. Because SampEn estimates may depend on parameter selection, three additional combinations (m = 2, r = 0.10 × SD; m = 3, r = 0.10 × SD; and m = 3, r = 0.25 × SD) were analysed to evaluate parameter sensitivity (Yentes and Raffalt, 2021). Results obtained using m = 2 and r = 0.25 × SD are presented in the main manuscript, and results for the alternative parameter combinations are provided in Supplementary Figure S1, Supplementary Figure S2, Supplementary Figure S3.

Local dynamic stability was evaluated using the largest Lyapunov exponent (LyE), calculated from the ten cycles of running time series of right hip, knee, and ankle joint angles in the sagittal, coronal, and transverse planes, as well as from the COM position time series (Supplementary Table S1). The LyE was estimated using the Rosenstein algorithm (Rosenstein et al., 1993). Before state-space reconstruction, each time series was standardized to zero mean and unit SD. Time delay was determined as the first local minimum of average mutual information, or the global minimum within the evaluated lag range if no local minimum was found. Embedding dimension was estimated using the false-nearest-neighbours method as the smallest dimension (2–10) yielding ≤1% false neighbours, with relative and absolute thresholds of 10 and 2, respectively. The Theiler window was defined as the mean signal period derived from the dominant non-zero frequency in Welch’s power spectral density estimate, and candidate neighbours were required to be separated by more than this window. Mean logarithmic divergence was calculated over approximately 1.5 gait cycles. LyE was estimated as the least-squares slope over a stride-scaled early divergence window extending from sample 1 to approximately one gait cycle (110–157 samples, corresponding to 0.55–0.785 s at 200 Hz), ensuring consistent and objective estimation across participants and signals:

where d(t) denotes the separation between neighboring trajectories in the reconstructed state space at time t, and λ represents the slope of the mean logarithmic divergence curve over the specified early divergence interval.

All outcomes are reported as mean ± SD. For each outcome, a random-intercept linear mixed-effects model (LMM) was fitted, with participant included as a random effect to account for the repeated observations in the HRFS participants. A three-level group–condition factor, hereafter termed State, was included as a fixed effect, with levels HRFS-AN, HRFS-GR, and HNRFS. State represented the three observed group–condition combinations rather than three independent participant groups. The model was specified as: Outcome ~ State + (1 | participant)

For each outcome, the linear mixed-effects model evaluated the overall effect of State and the following three contrasts: HRFS-GR minus HRFS-AN, HNRFS minus HRFS-AN, and HNRFS minus HRFS-GR. The p values for the three contrasts were adjusted within each outcome using the Holm correction. Contrast results are reported as estimated mean differences with 95% confidence intervals and Holm-adjusted p values. Statistical significance was set at a two-sided α of 0.05. All statistical analyses were performed in Python using statsmodels package version 0.14.6 (Seabold and Perktold, 2010).

A total of 19 RFS runners enrolled in GR, 2 discontinued due to triceps surae soreness, and 1 was excluded for repeated protocol nonadherence. Subsequently, 16 completed the post-intervention test, while 3 failed to convert to NRFS, and thus, were excluded, leaving 13 participants (age: 34.2±8.6 years; body mass: 70.0±8.5 kg; height: 1.73 ± 0.07 m) in the HRFS-GR group for analysis.

Table 1 presents the descriptive statistics and LMM comparisons for the spatiotemporal, foot strike angle, and ground reaction force parameters across the three states.

All three observed states were classified as NRFS, with mean foot strike angles below 8°. A significant overall effect of State was observed for foot strike angle (p < 0.001, Table 1). Compared with HRFS-AN, foot strike angle was higher in HRFS-GR (MD = 4.52°; p < 0.001) and HNRFS (MD = 2.95°; p = 0.011). No significant difference was observed between HRFS-GR and HNRFS (p > 0.05).

Significant overall effects of State were also observed for VILR (p = 0.028) and VALR (p = 0.025). Compared with HRFS-AN, HRFS-GR exhibited higher VILR (MD = 17.32 BW/s; p = 0.012) and higher VALR (MD = 14.93 BW/s; p = 0.016). Neither VILR nor VALR differed significantly between HNRFS and HRFS-AN or between HNRFS and HRFS-GR (all p > 0.05). No significant overall State effects were observed for cadence, step length, contact time, or vGRF (all p > 0.05; Table 1).

Table 2 presents the discrete lower-limb joint angle and COM ROM results across the three states. Significant overall effects of State were observed for the hip sagittal-plane minimum angle (p = 0.046), knee sagittal-plane minimum angle (p = 0.035), knee sagittal-plane ROM (p = 0.013), and ankle coronal-plane maximum angle (p = 0.034). For the hip sagittal-plane minimum angle, HNRFS exhibited a higher value than HRFS-GR (MD = 6.36°, p = 0.040). Knee sagittal-plane ROM was higher in HRFS-GR than in HRFS-AN (MD = 7.79°; p = 0.038) and higher in HNRFS than in HRFS-AN (MD = 14.25°; p = 0.035). For the ankle coronal-plane maximum angle, HNRFS exhibited a higher value than HRFS-GR (MD = 6.51°; p = 0.044). Although a significant overall State effect was observed for the knee sagittal-plane minimum angle, none of its pairwise contrasts remained significant after Holm adjustment. No significant overall State effects were detected for the remaining discrete joint angle parameters or COM ROM variables (all p > 0.05; Table 2).

Table 3 presents the CV results for the lower-limb joint angle and COM variables across the three states. Significant differences in joint-angle CV parameters were observed among the three states. At the hip, the sagittal-plane minimum angle CV was lower in HRFS-GR than in HRFS-AN (MD = −0.33; p = 0.004). Overall State effects were also observed for the hip sagittal-plane ROM CV (p = 0.042) and transverse-plane ROM CV (p = 0.042), although none of their pairwise comparisons were significant.

At the knee, HRFS-GR exhibited lower CV values than HRFS-AN for the sagittal-plane minimum angle (MD = −0.09; p = 0.008), sagittal-plane ROM (MD = −0.14; p = 0.002), transverse-plane maximum angle (MD = −0.19; p = 0.049), and transverse-plane ROM (MD = −0.07; p < 0.001).

At the ankle, the sagittal-plane maximum angle CV was lower in both HRFS-GR (MD = −0.35; p = 0.005) and HNRFS (MD = −0.35; p = 0.005) than in HRFS-AN. HRFS-GR also exhibited lower sagittal-plane ROM CV (MD = −0.13; p = 0.036) and transverse-plane ROM CV (MD = −0.10; p = 0.034) than HRFS-AN. No significant overall State effects were observed for the remaining joint-angle or COM CV parameters (all p > 0.05).

Figure 2 presents the SampEn values and pairwise comparisons for the hip, knee, ankle, and COM across the three states (m = 2, r = 0.25).

At the hip, a significant overall effect of State was observed for transverse-plane SampEn (p = 0.012). HRFS-GR exhibited higher SampEn than HRFS-AN (MD = 0.07; p = 0.012), while HNRFS exhibited lower SampEn than HRFS-GR (MD = −0.05; p = 0.050).

At the knee, a significant overall State effect was also observed for knee sagittal-plane SampEn (p < 0.001). Compared with HRFS-AN, SampEn was lower in HRFS-GR (MD = −0.02; p = 0.004) and HNRFS (MD = −0.05; p < 0.001). HNRFS also exhibited lower SampEn than HRFS-GR (MD = −0.03; p < 0.001).

At the ankle, a significant overall State effect was observed for sagittal-plane SampEn (p < 0.001). SampEn was lower in HRFS-GR (MD = −0.04; p = 0.002) and HNRFS (MD = −0.05; p = 0.002) than in HRFS-AN, whereas HRFS-GR and HNRFS did not differ significantly.

For the COM, a significant overall State effect was observed for vertical SampEn (p = 0.002). Vertical COM SampEn was lower in HRFS-GR (MD = −0.01; p = 0.014) and HNRFS (MD = −0.03; p = 0.014) than in HRFS-AN, with no significant difference between HRFS-GR and HNRFS. No significant overall State effects were observed for the remaining joint or COM SampEn outcomes (all p > 0.05).

Figure 3 presents the LyE values and pairwise comparisons for the hip, knee, ankle, and COM across the three states.

At the hip, a significant overall effect of State was observed for coronal-plane LyE (p = 0.040). HNRFS exhibited lower coronal-plane LyE than HRFS-GR (MD = −0.24; p = 0.038), while HRFS-AN did not differ significantly from either HRFS-GR or HNRFS. No significant overall State effects were observed for hip sagittal- or transverse-plane LyE.

At the knee, a significant overall State effect was observed for sagittal-plane LyE (p < 0.001). Compared with HRFS-AN, sagittal-plane LyE was lower in HRFS-GR (MD = −0.18; p = 0.017) and HNRFS (MD = −0.38; p < 0.001). HNRFS also exhibited lower sagittal-plane LyE than HRFS-GR (MD = −0.20; p = 0.044). No significant overall State effects were observed for knee coronal or transverse-plane LyE.

At the ankle, a significant overall State effect was observed for coronal-plane LyE (p = 0.044); however, none of the pairwise comparisons were significant. No significant overall State effects were observed for ankle sagittal or transverse-plane LyE.

For the COM, a significant overall State effect was observed for vertical LyE (p = 0.003). Vertical COM LyE was lower in HRFS-GR (MD = −0.14; p = 0.019) and HNRFS (MD = −0.34; p = 0.009) than in HRFS-AN, while HRFS-GR and HNRFS did not differ significantly. No significant overall State effects were observed for COM LyE in the anterior–posterior or medial–lateral directions.

This study compared lower-limb kinematics, movement regularity, and local dynamic stability across the HRFS-AN, HRFS-GR, and HNRFS states. The first hypothesis was largely, but not completely, supported: most spatiotemporal and discrete kinematic outcomes did not differ significantly among states, although significant State effects were observed for foot-strike angle, vertical loading rates, selected joint-angle variables, and several CV outcomes. The second hypothesis was supported for selected nonlinear outcomes: compared with HNRFS, HRFS-AN exhibited higher sagittal-plane knee and ankle SampEn, higher vertical COM SampEn, and higher sagittal-plane knee and vertical COM LyE. The third hypothesis was partially supported: HRFS-GR and HNRFS did not differ significantly for most outcomes, but differences remained in selected hip and knee SampEn and LyE measures. All three states had mean foot-strike angles below 8° and were therefore classified as NRFS, indicating that the HRFS participants were able to adopt the intended strike pattern both acutely and after retraining. Cadence, step length, contact time, and peak vGRF did not differ significantly among states, and most discrete joint-angle and COM ROM outcomes also showed no significant differences.

These findings are broadly consistent with previous studies showing that habitual rearfoot strikers can acutely reproduce several sagittal-plane kinematic characteristics of habitual forefoot or non-rearfoot strikers (Valenzuela et al., 2015; Williams et al., 2000). Previous gait-retraining studies have also reported reductions in impact peak and vertical loading rates relative to habitual RFS running (Futrell et al., 2020; Yang et al., 2020). Although the present study did not assess the habitual RFS condition, neither HRFS-AN nor HRFS-GR differed significantly from HNRFS in peak vGRF, VILR, or VALR, providing no evidence of greater impact loading in either modified NRFS state than in habitual NRFS under the present testing conditions. Nevertheless, VILR and VALR were higher in HRFS-GR than in HRFS-AN. This difference may partly reflect the more negative foot-strike angle in HRFS-AN, indicating a more pronounced anterior foot strike, whereas the foot-strike angle of HRFS-GR was closer to that of HNRFS.

The nonlinear measures further differentiated acute NRFS adoption from retrained and habitual NRFS running. Compared with HRFS-AN, both HRFS-GR and HNRFS exhibited lower SampEn in sagittal-plane knee and ankle and vertical COM movement, together with lower LyE in sagittal-plane knee and vertical COM movement. Given the more pronounced anterior foot strike adopted during HRFS-AN, the higher LyE in these selected trajectories may be compatible with the challenge of acutely adopting a more pronounced and unfamiliar strike strategy. This interpretation is indirectly supported by Ekizos et al. (2017; 2018), who reported reduced local dynamic stability when runners acutely shifted toward a more anterior foot strike, whereas this effect was not observed following a longer period of gradual transition. In the present study, the lower LyE observed after 12 weeks of retraining similarly suggests that the greater local divergence associated with the acute condition may diminish with continued practice. Directly comparable evidence regarding SampEn following foot-strike modification remains scarce. However, selected reductions in COM-velocity SampEn have been observed after accommodation to walking with unstable footwear (Buchecker et al., 2018), providing indirect evidence that repeated exposure to a modified locomotor task may alter the temporal organization of movement. Therefore, the lower SampEn observed after retraining may reflect the development of a more regular movement strategy through continued practice.

SampEn, LyE, and CV should therefore be interpreted as complementary rather than interchangeable measures of movement behaviour (Dingwell et al., 2001; Stergiou and Decker, 2011). SampEn characterizes the temporal regularity and complexity of a movement trajectory, with lower values indicating a more regular and predictable temporal structure (Richman and Moorman, 2000; Yentes and Raffalt, 2021), whereas LyE quantifies the rate of local trajectory divergence, with lower values indicating greater local dynamic stability (Dingwell et al., 2001; Rosenstein et al., 1993). In contrast, CV describes the consistency of discrete kinematic extrema across gait cycles without considering the temporal organization of the continuous trajectory. In the present study, the concurrent lower SampEn and LyE in sagittal-plane knee and vertical COM motion indicates that these specific trajectories were both more regular and less locally divergent in HRFS-GR and HNRFS than in HRFS-AN. However, sagittal-plane ankle SampEn differed without a corresponding LyE difference, while transverse-plane hip SampEn was higher in HRFS-GR than in HRFS-AN. Similarly, reductions in CV after retraining were observed in selected sagittal- and transverse-plane outcomes, but the overall CV pattern did not correspond fully with the nonlinear findings. These joint- and plane-specific results demonstrate that changes in movement regularity, local dynamic stability, and the variability of discrete kinematic events do not necessarily occur together. Accordingly, the nonlinear measures provided additional information that was not captured by CV, rather than simply being more sensitive, and their combined use enabled a multidimensional characterization of movement adaptation.

Despite the overall similarity between HRFS-GR and HNRFS, residual differences remained in selected nonlinear outcomes. HRFS-GR exhibited higher SampEn in sagittal-plane knee and transverse-plane hip motion and higher LyE in sagittal-plane knee and coronal-plane hip motion than HNRFS. Because HNRFS was an independent cross-sectional reference group rather than the longitudinal endpoint of the HRFS participants, these differences may reflect distinct adaptation histories and between-group characteristics rather than an insufficient retraining duration alone. Importantly, greater similarity to HNRFS should not automatically be regarded as a superior outcome. Movement variability may support adaptability, and its association with musculoskeletal injury remains inconsistent (Baida et al., 2018; Stergiou and Decker, 2011). Bonacci et al. (2020) also reported increased entropy based joint variability after gait retraining despite clinically meaningful improvements in pain. Thus, higher SampEn could reflect continued exploration or adaptability, whereas higher LyE indicates greater local trajectory divergence; neither finding alone establishes an adverse response or increased injury risk. Because injury incidence, pain, tissue loading, and performance were not prospectively assessed, the clinical significance of these residual differences remains uncertain.

Although SampEn and LyE cannot establish injury risk, adverse symptoms may occur during the transition to an NRFS pattern, particularly when foot-strike transition is combined with minimalist footwear exposure. In the present study, two participants withdrew during the early retraining phase because of calf and foot pain. Because the intervention also included foot and ankle strengthening exercises, these symptoms cannot be attributed to foot-strike modification, minimalist footwear, or any single intervention component. Nevertheless, these observations support gradual, individually monitored progression during gait retraining.

The current study has several limitations. First, nonlinear outcomes may be influenced by methodological choices. SampEn outcomes depend on parameter selection, and although multiple combinations of embedding dimension (m) and tolerance (r) were examined and a representative setting was selected, comparisons are only valid under identical parameter conditions. Likewise, LyE estimates may vary depending on reconstruction settings. However, to reduce methodological bias, LyE was calculated using the same Rosenstein algorithm for all time series, while time delay and embedding dimension were determined for each time series according to established recommendations, the corresponding time delay and embedding dimension values for each parameter have been provided in the Supplementary Material. Second, only male runners were included to improve sample homogeneity and reduce potential confounding effects related to sex-specific physiological fluctuations, which may influence training-load control during the 12-week gait retraining intervention. However, this limits the generalizability of the findings to female runners. Third, although a 12-week GR protocol was adopted, no post-intervention follow-up was conducted. Therefore, the long-term retention of changes in movement regularity and kinematic patterns after removal of external guidance remains unknown. Future studies should incorporate follow-up assessments, and workload monitoring and prospective symptom or injury surveillance, to better clarify the durability and practical significance of GR-induced adaptations during the transition process.

Acute NRFS adoption was characterized by a more anterior foot strike than HRFS-GR and HNRFS and lower VILR and VALR than HRFS-GR, whereas peak vGRF and most spatiotemporal and discrete kinematic outcomes were similar among states. Compared with acute adoption, HRFS-GR showed greater regularity in selected knee, ankle, and COM trajectories and greater local dynamic stability in knee and COM motion, broadly resembling the pattern observed in HNRFS. These within-participant changes suggest progressive adaptation to the newly learned NRFS pattern, although selected knee and hip differences remained between HRFS-GR and HNRFS. Overall, the findings suggest that a gradual 12-week gait-retraining program may be an appropriate approach for transitioning habitual rearfoot strikers to NRFS.