Like other animals, humans use their legs like springs to save energy during running. One potential contributor to leg stiffness in humans is the longitudinal arch (LA) of the foot. Studies of cadaveric feet have demonstrated that the LA can function like a spring, but it is unknown whether humans can adjust LA stiffness in coordination with more proximal joints to help control leg stiffness during running. Here, we used 3D motion capture to record 27 adult participants running on a forceplate-instrumented treadmill, and calculated LA stiffness using beam bending and midfoot kinematics models of the foot. Because changing stride frequency causes humans to adjust overall leg stiffness, we had participants run at their preferred frequency and frequencies 35% above and 20% below preferred frequency to test for similar adjustments in the LA. Regardless of which foot model we used, we found that participants increased LA quasi-stiffness significantly between low and high frequency runs, mirroring changes at the ankle, knee and leg overall. However, among foot models, we found that the model incorporating triceps surae force into bending force on the foot produced unrealistically high LA work estimates, leading us to discourage this modeling approach. Additionally, we found that there was not a consistent correlation between LA height and quasi-stiffness values among the participants, indicating that static LA height measurements are not good predictors of dynamic function. Overall, our findings support the hypothesis that humans dynamically adjust LA stiffness during running in concert with other structures of the leg.
The legs of terrestrial animals can be simply modeled as springs of a given stiffness that store and release elastic energy during running, reducing metabolic cost (Blickhan, 1989). Overall leg stiffness, which is determined by the integrated behavior of muscles, tendons and ligaments (Farley et al., 1991), influences several important aspects of gait, including foot contact time and stride frequency (Farley and González, 1996). Animals maintain relatively constant leg stiffness across running speeds (Farley et al., 1993; He et al., 1991), but humans have been shown to adjust leg stiffness in response to compliance or irregularities of the underlying surface to maintain constant running mechanics (Ferris et al., 1998; Grimmer et al., 2008; Kerdok et al., 2002). Leg stiffness is likely modulated by the coordinated action of muscles that control the rotational stiffness of the joints within the leg while also actively absorbing and producing power (Farley and Morgenroth, 1999). Previous studies of running leg stiffness in humans have focused on the contributions of the knee and ankle joints (Arampatzis et al., 1999; Günther and Blickhan, 2002), but so far, the effects of joints distal to the ankle have not been studied. Humans have plantigrade feet composed of multiple joints, but due to its anatomical complexity, the human foot has often been modeled as a single rigid segment. However, recent studies have revealed surprisingly high joint mobility within the human foot (Arndt et al., 2007; Kessler et al., 2019; Lundgren et al., 2008), as well as important mechanical functions for the intrinsic foot joints in gait (Kern et al., 2019; Riddick et al., 2019; Takahashi et al., 2017). It is therefore reasonable to hypothesize that humans may be capable of modulating foot stiffness in coordination with other lower limb joints to control overall leg stiffness during running.
A unique characteristic of the human foot is the longitudinal arch (LA), in which the conformation of the tarsal and metatarsal bones elevates the midfoot, particularly on the foot's medial side (Fig. 1A) (Holowka and Lieberman, 2018). The LA is partly maintained by ligamentous structures including the plantar aponeurosis, a broad sheet of fibrous connective tissue that spans the plantar surface of the human foot (Fig. 1B) (McKeon et al., 2014; Sichting et al., 2020). To investigate the function of the passive ligaments that help maintain the LA, Ker et al. (1987) compressed cadaveric human feet under simulated running loads that were derived from a model of the foot as a beam undergoing three-point bending from ground reaction force anteriorly, and triceps surae force posteriorly. This study found that the LA can function like a spring by storing and releasing elastic energy in ligaments, and that these ligaments cumulatively contribute to overall passive foot stiffness. Recently, Stearne et al. (2016) applied this beam model to humans running in vivo, and found reasonably close agreement between predicted mechanical energy savings in the LA and actual metabolic energy savings during running. These and other studies support the spring-like function of the foot, and suggest that ligamentous structures help determine the passive stiffness of the LA (Alexander, 1991; Huang et al., 1993; Venkadesan et al., 2020; Welte et al., 2018).
The intrinsic and extrinsic foot muscles in humans may be capable of adjusting passive LA stiffness. Although most of these muscles insert onto the toes, their bellies and/or tendons cross multiple joints comprising the LA, potentially making them capable of resisting bending moments at these joints (McKeon et al., 2014). So far, most in vivo experimental investigations have focused on the intrinsic foot muscles, and the two most strongly linked to LA function are the abductor hallucis and flexor digitorum brevis, both of which originate on the calcaneus and insert distally on the phalanges (Fig. 1C). Electromyography (EMG) studies have demonstrated that these muscles are active during walking and running (Kelly et al., 2015; Mann and Inman, 1964; Reeser et al., 1983), can elevate the LA during static compressive loading (Kelly et al., 2014), and can modulate the energetic function of the LA during static loading and stair stepping (Kelly et al., 2019; Riddick et al., 2019). However, Farris et al. (2019) found that anesthetizing these muscles had little effect on the stiffness of the midfoot region in human participants during running, and Kessler et al. (2020) found that relative midfoot stiffness did not correspond to the degree of activation of these muscles during foot contact in a hopping task. Rather than using a beam model, these two studies followed others (e.g. Bruening and Takahashi, 2018) in modeling the human foot as a multi-segment structure with a midfoot ‘joint’ located roughly midway between the calcaneus and metatarsals, and calculated LA stiffness from the moments and motion at this joint using inverse dynamics. While their results suggest that the intrinsic foot muscles are not used to dynamically stiffen the LA, Kern et al. (2019) found that humans can increase LA stiffness when walking with added external loads, although they were unable to identify the mechanism responsible. Additionally, it is possible that extrinsic foot muscles like the tibialis posterior, that insert onto bones of the midfoot but not the toes, could actively stiffen the LA (McKeon et al., 2014). Thus, whether humans are capable of dynamically adjusting LA stiffness in conjunction with overall leg stiffness during running remains unknown.
In addition to potential active contributions, the basic anatomical determinants of LA stiffness are also poorly understood. One such variable of widespread interest is LA height: very high arches are commonly associated with overly rigid feet (‘pes cavus’), and very low arches are associated with highly compliant feet (‘flexible pes planus’) (Franco, 1987; Neumann, 2002). While the mechanistic basis for these associations is not completely clear, low arches are often ascribed to ligamentous laxity, in which the ligaments supporting the LA do not provide the necessary tension to prevent arch compression during weight bearing. Both extremes of LA height have been associated with increased risk of foot injury, with both overly stiff and overly compliant arches thought to be poor shock absorbers (Kaufman et al., 1999; Simkin et al., 1989; Williams et al., 2001). However, there is scant quantitative evidence linking LA height to dynamic stiffness. Zifchock et al. (2006) found that higher-arched individuals have slightly stiffer LAs than lower-arched individuals under static loading conditions, but the actual correlation between arch height and stiffness in their study was weak. Furthermore, Holowka et al. (2018) found that static measurements of LA stiffness were not correlated with dynamic stiffness estimates during walking, revealing considerable uncertainty in the exact relationship between LA height and function during gait.
Here, we experimentally tested whether humans can adjust LA stiffness during running in coordination with other joints of the leg, and whether the height of the LA is a determinant of its dynamic stiffness. To address these questions, we also compared different models of the foot, including beam bending models (Ker et al., 1987; Stearne et al., 2016; Venkadesan et al., 2020) and a model based on midfoot motion and moments. Previous research has demonstrated that humans increase leg stiffness when running at higher stride frequencies in order to reduce vertical displacement of the body's center of mass and thereby enable shorter ground contact times (Farley and González, 1996). This modulation can be accomplished by increasing the stiffness of the lower limb joints, especially the ankle (Arampatzis et al., 2001). A stiffer foot could also contribute in two ways: first, by reducing LA compression during loading, a stiffer foot can help reduce the change in effective limb length, and by extension center of mass displacement; second, given the foot's role in transmitting power from the ankle to the ground, a stiffer foot could allow for more rapid force production during push-off when using shorter contact times. In support of these ideas, Kessler et al. (2020) recently demonstrated that humans increase LA stiffness in correspondence with ankle stiffness during single-leg hopping, when hopping frequency is increased. Thus, we analyzed human participants running at different stride frequencies, and predicted that they would increase LA stiffness with increasing stride frequency in a manner similar to changes in overall leg stiffness. We also tested whether the measured LA heights of participants were associated with dynamic LA stiffness during running, predicting a modest correlation between height and stiffness. Finally, we compared results from the different modeling approaches applied in this study to assess their relative ability to quantify foot kinetics.
MATERIALS AND METHODS
Twenty-seven adult participants (10 female, 22.3±3.4 years, 63.5±7.1 kg; 17 male, 21.9±3.5 years, 81.9±13.4 kg) were included in the study. Inclusion criteria were no musculoskeletal injury that caused pain during running or otherwise affected gait, and the ability to run with stride frequencies 35% above and 20% below preferred frequency for at least 1 min. All experimental procedures were approved by the Harvard University Committee on the Use of Human Subjects, and all participants provided written informed consent.
Prior to experiments, we measured the arch height index (AHI) of the participant's right foot during standing using a custom-machined device following Butler et al. (2008). AHI is calculated as medial midfoot dorsum height divided by medial foot length excluding the hallux, and has been shown to be a robust and repeatable measurement of relative longitudinal arch height (Butler et al., 2008).
We recorded participants running on a split belt forceplate-instrumented treadmill (Bertec, Columbus, OH, USA) using an eight-camera Oqus motion-capture system (Qualisys, Gothenburg, Sweden) in the Skeletal Biology and Biomechanics Laboratory at Harvard University. For all conditions, participants ran at 0.77 Froude (Fr; 2.64±0.09 m s−1) to ensure dynamic similarity regardless of leg length (Alexander and Jayes, 1983), which we measured using right greater trochanter height during standing. We selected 0.77 Fr to correspond with the speeds used in Farley and González's (1996) study of the effects of stride frequency on human leg stiffness. In preliminary tests, we found that individuals could not run comfortably on the treadmill at the necessary stride frequencies while barefoot. Therefore, we had participants run in LUNA Mono sandals (LUNA Sandals, Seattle, WA, USA), which enabled them to run comfortably with their preferred foot strike patterns while allowing us to track markers placed directly on their feet. These sandals consist of a flat and relatively thin layer of EVA foam without any features that would restrict the natural motion of the arch (Fig. 2A). To determine preferred stride frequency, participants ran on the treadmill at 0.77 Fr for several minutes until comfortable, and then we recorded the number of strides taken in 1 min. With the help of a metronome, participants practiced running at +35% and −20% of their preferred stride frequency, until they could maintain these frequencies comfortably. We chose these frequencies because they represented the approximate limits of what most individuals were able to achieve in our preliminary tests, with most finding it considerably more difficult to run at lower as opposed to higher percentages of their preferred frequency.
After this practice session, we placed markers on the pelvis and lower limbs following the marker set described in Cappozzo et al. (1995). We placed an additional marker on the navicular tuberosity to quantify LA stiffness (see below). Unfortunately, the straps of the sandals precluded the application of a more detailed foot marker set, such as those used in multi-segment foot models (e.g. Leardini et al., 2007) (Fig. 2A). We recorded participants standing on the treadmill to determine joint neutral positions for subsequent model building, then had participants run on the treadmill at the high (+35%), low (−20%) and preferred stride frequencies using a metronome, in random order. For each condition, we recorded a 1 min trial once the participant had achieved the desired stride frequency.
We selected 10 strides from each recorded trial for analysis, and used Visual3D v.6 (C-Motion Inc., Germantown, MD, USA) for 3D kinematics and inverse dynamics calculations. We filtered both marker and forceplate data using a fourth order low-pass Butterworth filter with a 15 Hz cutoff frequency, which we deemed appropriate based on residual analysis of the raw data (Winter, 2005). We used the same cutoff frequency for both data types to avoid computational artifacts in inverse dynamics calculations (Derrick et al., 2020; Kristianslund et al., 2012). Kinetic and kinematic variables were imported into custom-written MATLAB scripts (The MathWorks, Natick, MA, USA) for further analysis. All custom-written MATLAB and R code used to process and analyze data is available from the corresponding author upon reasonable request.
LA stiffness models
We calculated LA stiffness using three different models of the foot (Fig. 2). Two were based on the beam bending model originally implemented by Ker et al. (1987), and will be referred to as Beam1 and Beam2 models. The third model we devised is based on a simplified 2D representation of the foot as a two-segment structure with a midfoot joint, which we will refer to as the Midfoot model.
The extent to which triceps surae force actually contributes to bending the foot beam during loading is unclear. Although previous studies have demonstrated that triceps surae activation can increase tension in the plantar aponeurosis (Carlson et al., 2000; Erdemir et al., 2004), the muscle's total contribution to LA bending during running may be somewhat less than that estimated in the original Ker et al. (1987) model, which included all force exerted by triceps surae during running into its calculation. In contrast, Venkadesan et al. (2020) recently modeled the foot as a beam during running while assuming that bending force would be roughly equal to the magnitude of ground reaction force (3 times body weight), resulting in bending forces less than half those in Ker et al. (1987). Welte et al. (2018) used a similar approach to model elastic energy storage in the human LA in vivo at the loads used during walking. To address the uncertainty surrounding triceps surae force contributions to foot bending forces, we applied a two-point bending Beam2 model of the foot where Fbeam2 was equivalent to just the vertical component of ankle joint reaction force (FGRF from Eqn 2 above) (Fig. 2E). Thus, these two beam models provide boundary conditions for the relative contribution of triceps surae force to bending force on the foot.
We calculated Δzarch as the change in vertical height of the marker on the navicular tuberosity relative to the sole of the foot following Stearne et al. (2016) (Fig. 2B). To define the sole of the foot, we used Visual3D to create virtual markers in the standing neutral trial by projecting the positions of the calcaneus, first metatarsal head and fifth metatarsal head markers into an x–y plane parallel with the surface of the treadmill. The relative positions of these virtual markers were subsequently reconstructed in locomotion trials, and used to define a segment coordinate system for the foot sole originating at the virtual calcaneus marker. We transformed the coordinates of the navicular tuberosity marker into this coordinate system, and calculated Δzarch as the change in vertical height of this marker relative to its height at the moment of foot strike (first frame of stance phase).
Stiffness in linearly elastic bodies is calculated as the slope of the relationship between force and displacement. As biological structures like limbs and joints are not true elastic bodies, particularly when acted upon by internal forces from muscles (Latash and Zatsiorsky, 1993; Rouse et al., 2013), all ‘stiffness’ measurements were actually measurements of ‘quasi-stiffness’, in which the structure being measured is not at internal equilibrium, as defined by Latash and Zatsiorsky (1993). As observed in previous in vitro studies (Huang et al., 1993; Ker et al., 1987; Venkadesan et al., 2020), we found that the relationships between Δzarch and Fbeam1,2 were not completely linear, but instead tended to exhibit patterns approximating narrow hysteresis loops with self-intersection. Thus, we calculated LA quasi-stiffness (kbeam) as the slope of the linear least squares regression line fitted to the corresponding values of Δzarch versus Fbeam. This approach is similar to that carried out in previous investigations of limb and joint quasi-stiffness in humans and goats (Clites et al., 2019; Günther and Blickhan, 2002; Kern et al., 2019; Kessler et al., 2020; Shamaei et al., 2013). Recognizing that the quasi-stiffness of a joint or joint complex may change at different points in the gait cycle, we calculated kbeam,load as the slope of the least squares regression line fitted to all values of Δzarch versus Fbeam up to the point of maximum Fbeam. We calculated kbeam,unload using the same approach, but with all values after the point of maximum Fbeam. Similar approaches have been used to assess changing quasi-stiffness values across stance phase in the LA and ankle in previous studies (e.g. Bruening et al., 2018; Kern et al., 2019; Shamaei et al., 2013). Finally, we also calculated the total range of Δzarch as maximum Δzarch minus minimum Δzarch, as well as the maximum Fbeam values during stance phase, to assess their relative effects on kbeam among the different stride frequency conditions.
The LA was also modeled as consisting of two linear segments: the rearfoot, defined by the line between the markers placed on the posterior calcaneus and the navicular tuberosity, and the forefoot, defined by the line connecting the navicular tuberosity to the first metatarsal head (Fig. 2C,F). Midfoot motion, θmid, was calculated as the angle between these segments, projected onto the sagittal plane of the foot segment. This is a highly simplified version of standard multi-segment foot models that permit the calculation of 3D kinetics between rearfoot and forefoot segments (e.g. Bruening and Takahashi, 2018; Kelly et al., 2019), but which require detailed marker sets that we were unable to use. However, one previous study by Kelly et al. (2018) also used a simplified 2D sagittal plane model of the midfoot to estimate LA kinetics during running, and achieved work results that were comparable to those from other studies (Bruening et al., 2018; Stearne et al., 2016). Additionally, Kessler et al. (2019) recently demonstrated that a projected 2D sagittal plane angle used to represent LA motion during running agreed reasonably closely with bone motion captured using biplanar videoradiography, suggesting our approach may provide a reasonable approximation of LA kinematics.
To estimate the external moment acting upon the midfoot ‘joint’, Mmid, we followed some previous 3D multi-segment foot models in assuming that only forces acting upon the forefoot segment create midfoot moments, and that these forces are minimal when the center of pressure is posterior to the midfoot ‘joint’ (Bruening et al., 2018; Farris et al., 2019; Kelly et al., 2018; Kern et al., 2019). When the center of pressure passed anterior to the navicular tuberosity marker, we calculated Mmid as the external moment about the navicular tuberosity marker caused by the ground reaction force vector in the sagittal plane. To do so, we multiplied the vertical and horizontal force components with their respective moment arms and summed their products.
As with the beam models described above, we calculated LA quasi-stiffness values, kmid, as the slopes of the linear least squares regression lines fitted to the corresponding values of the change in midfoot motion (Δθmid) versus Mmid for loading (up to the peak Mmid), and unloading (after the peak Mmid), periods of stance. We also calculated the total range of Δθmid as maximum Δθmid minus minimum Δθmid, and the maximum Mmid value during stance phase, to assess their relative effects on kmid among the different stride frequency conditions.
Leg and joint stiffness models
We calculated leg quasi-stiffness (kleg) in the sagittal plane following a similar approach to previous studies (e.g. Liew et al., 2017). In brief, we measured leg length as the distance from a marker placed on the greater trochanter to the center of pressure under the foot in the sagittal plane, and used this to measure the changes in leg length (Δlleg) during stance phase. We computed compressive force on the leg (Fleg) as the dot product of the sagittal plane ground reaction force and leg length vectors. We then calculated kleg as the absolute value of the slope of the least squares regression line fitted to all values of Δlleg versus Fleg during loading (up to peak Fleg) and unloading (after peak Fleg) periods of stance. This model of leg stiffness assumes that medio-lateral force and displacement values will have negligible effects on quasi-stiffness calculations. Additionally, in the present study, the sole of the sandal worn by participants is necessarily incorporated into the overall leg quasi-stiffness calculations, which likely has some minor effect on kleg values.
For the ankle and knee joints, we calculated flexion–extension excursion angles in the sagittal plane (θankle and θknee), with joint positions measured during the standing neutral trial set as the 0 deg joint angles. We used a standard inverse dynamics approach to calculate net joint moments (Mankle and Mknee), assuming the foot is a rigid segment. Following Günther and Blickhan (2002), we calculated ankle (kankle) and knee (kknee) quasi-stiffness values as the slopes of the least squares regression lines fitted to the respective angular excursion versus net joint moment values during loading and unloading periods of these joints (up to, and after, their respective peak joint moments).
Foot strike angle
As previous studies have demonstrated the effects of different foot strike patterns on midfoot kinematics and stiffness (Bruening et al., 2018; Kelly et al., 2018; McDonald et al., 2016; Perl et al., 2012), we decided to control for this factor in our analysis. To do so, we calculated foot strike angle (FSA) as the projected sagittal plane angle between the long axis of the foot (the vector between the markers placed on the posterior calcaneus and the fifth metatarsal head) and the surface of the treadmill. We computed FSA by deducting the value of this angle at the first frame of stance from the value in the standing neutral position trial, such that positive values indicate initial rearfoot contact, and negative values indicate initial forefoot contact. Although FSAs are often used to distinguish among rearfoot, midfoot and forefoot strikes (Altman and Davis, 2012; Lieberman et al., 2010), we treated FSA as a continuous variable to statistically compare stride frequency conditions.
To compare energy use in the LA during running among the different foot models, as well as with energy use in the leg overall, and in the knee and ankle joints, we calculated work values for the preferred frequency running condition. We calculated compression/angular velocity as the time derivatives of the respective motion variables for each structure (e.g. Δzarch, Δθmid), and calculated net power by multiplying the velocity values by their respective compressive force or moment values (e.g. Fbeam, Mmid). For each structure, we calculated total positive and negative work (W) by integrating the positive and negative power values with respect to time.
After removing steps that had motion artifacts from gaps in marker tracking, we analyzed a minimum of 5 and a maximum of 10 steps per condition per participant (mean±s.d. 9.6±0.9). We averaged values for each variable across all steps within participants for statistical analyses. All statistical tests were carried out using R software (http://www.R-project.org/). We performed Shapiro–Wilk normality tests on each variable and visually checked variable distributions for normality and similarity in variance across conditions. All quasi-stiffness values were log-transformed to achieve normal distributions. We also removed participants with outlier values 3 s.d. above or below sample means from analysis of a given variable. This step was taken because quasi-stiffness values could be affected dramatically by extreme motion or force variables at the beginning or end of stance as a result of things like high impact forces or rapid joint flexion during toe-off.
To test for differences in quasi-stiffness, range of motion/compression, and peak force/moment values for each structure among stride frequency conditions, we used the ‘lme4’ package (Bates et al., 2015) to construct linear mixed effects models, with participant identity set as a random effect. We also included FSA as a covariate in all quasi-stiffness models to account for possible effects of foot strike posture on these variables. We inspected residual plots and q–q plots to assess homoscedasticity and normality, respectively, in model residuals. In each case, these criteria were satisfied, so we performed Type 3 ANOVA tests to test for differences among conditions. When differences were detected, we used the ‘lsmeans’ package (Lenth, 2016) to conduct post hoc pairwise contrasts between frequency conditions with a Holm–Bonferroni P-value correction. FSA was not normally distributed, and could not be made normal by log-transformation. Thus, we carried out a non-parametric Friedman's test to test for differences in FSA among stride frequency conditions, and conducted Wilcoxon signed-ranks tests for post hoc paired comparisons between conditions, adjusting the alpha value for significance following a Holm–Bonferroni correction. For all statistical tests, an alpha value of 0.05 was used as an indicator of significance.
Lastly, we tested for associations between our different measures of LA quasi-stiffness (kbeam, kmid) and static measures of LA height during the preferred frequency running condition. We created ordinary least squares regression models with AHI as the independent variable and quasi-stiffness as the response variable, and carried out Pearson's product-moment correlation coefficient tests to test for associations between variables. We chose not to scale quasi-stiffness variables relative to body size for these analyses, because doing so would entail calculating LA compression as a proportion of static LA height, which is a component of both AHI calculations, thus resulting in auto-correlation.
General features of running frequency conditions
Average values for spatiotemporal, kinematic and kinetic variables for each running frequency are presented in Table 1. All statistical tests are presented in Table S1. On average, participants ran at 81% and 133% of their preferred stride frequencies in the low and high frequency runs, respectively. These differences corresponded to an average of 13% longer contact times in low frequency runs, and 17% shorter contact times in high frequency runs. Maximum vertical ground reaction forces were on average 11–12% lower in high frequency runs than in other conditions, but similar between low and preferred frequency runs. Participants landed with slightly but significantly lower FSAs (P=0.001) when running at high versus preferred frequencies, meaning they tended more towards forefoot strikes in the former and heel strikes in the latter (Fig. S3). There were no other differences in FSA among frequencies.
LA kinematics and kinetics
During stance phase, LA motion (Δzarch, Δθmid) and loading (Fbeam1, Fbeam2, Mmid) variables increased sinusoidally with apices near midstance (Fig. 3). Δzarch and Δθmid typically became negative near the end of stance, indicating that the LA is higher when the foot leaves the ground than at foot strike. Total Δzarch and Δθmid values were significantly different in all comparisons between running frequencies, and were greatest in low frequency runs, and least in high frequency runs (Fig. 3D,J). The magnitude of these differences was large, with low frequency runs averaging 52% greater total Δzarch and 49% total Δθmid than high frequency runs. Peak Fbeam values were significantly greater in preferred and low frequency runs compared with high frequency runs, and were significantly different between preferred and low frequency runs for Fbeam1 but not Fbeam2 (Fig. 3B,F). The relative magnitude of these differences was less than that for total Δzarch, with low frequency runs averaging 23% greater peak Fbeam1 and 12% greater peak Fbeam2 values than high frequency runs. Peak Mmid values were significantly different in all comparisons, and were greatest in low frequency runs but least in high frequency runs (Fig. 3H). The relative magnitude of these differences was less than that of total Δθmid, with low frequency runs averaging 17% higher Mmid values than high frequency runs.
The different foot models yielded generally similar results for LA quasi-stiffness with slight variation (Fig. 4, Table 2). We had to remove three participants from the analysis of kbeam2,load because they exhibited extremely low outlier values due the effects of high impact forces on the force–deformation relationship. In terms of absolute magnitude, kbeam1 values were 3–4 times greater than kbeam2 values, and loading quasi-stiffness values were greater than unloading quasi-stiffness values. During loading, LA quasi-stiffness values were significantly greater in high frequency runs than in low frequency runs for all models. These differences ranged from 24% on average for kbeam2,load to just 5% on average for kmidload. LA loading quasi-stiffness values were also significantly greater in preferred versus low frequency runs for kbeam2,load (11% average difference), and were significantly greater in high versus preferred frequency runs for kmid,load (6% average difference). During unloading, LA quasi-stiffness values were significantly different in all comparisons between running frequencies for each foot model. Generally, these differences were of greater relative magnitude than during the loading period, and ranged from 34% to 26% greater values on average in high versus low frequency runs for kbeam2,unload and kmid,unload, respectively. In most comparisons, for both loading and unloading, the magnitude of the difference between running conditions was greatest for kbeam2 and least for kmid. FSA was positively associated with kbeam1,load and kmid,load, indicating greater stiffness during loading with more heel-first-type foot strike postures.
Leg and joint quasi-stiffness
As expected, patterns in LA quasi-stiffness among running frequencies mirrored the patterns in kleg (Fig. 5A–C, Table 2), which exhibited significant differences in all comparisons between running frequencies for both loading and unloading periods. High frequency runs had the greatest kleg values and low frequency runs had the lowest kleg values during both loading (57% average difference) and unloading (55% average difference) periods. We had to remove four participants from the analysis of kknee,unload because they exhibited extremely low outlier values as a result of the effects of rapid knee flexion at the end of stance during high frequency runs. kankle (Fig. 5E,F) and kknee (Fig. 5H,I) exhibited similar patterns to kleg for loading and unloading, with high frequency runs yielding the greatest values and low frequency runs yielding the lowest. Most differences between running frequencies were significantly different for both joints, except that kknee was not significantly different between preferred and high frequency runs during the unloading phase. FSA was negatively associated with kleg,unload, kankle,unload and kknee,load, but positively associated with kankle,load.
Foot model work comparisons
The different foot models produced markedly different estimates of LA work (Fig. 6). The three-point bending Beam1 model produced negative work values that were on average 2.8 and 3.1 times those calculated from the two-point bending Beam2 and Midfoot models, respectively. Beam1 negative work was slightly greater than average values for the knee and ankle, and more than half the average values calculated for the leg overall. Beam1 also produced positive work values that were on average 3.5 and 1.7 times greater than those of the Beam 2 and Midfoot models, respectively. Beam1 positive work was also on average 39% and 52% of the positive work calculated from the leg and ankle, respectively. In contrast, the Beam2 and Midfoot models produced similar negative work values that were roughly 40% of the values of negative work at the ankle, on average. Positive work values calculated from the Midfoot model were double those calculated from the Beam2 model on average, but both of these models produced positive work values that were considerably lower than those of the knee, ankle or leg overall.
LA height versus stiffness tests
Scatterplots depicting linear regression estimates of the relationships between the LA height index (AHI) and LA quasi-stiffness from the different foot models during loading and unloading periods, along with the results of Pearson's product-moment correlation coefficient tests, are presented in Fig. 7. The Beam1 model showed significant positive relationships between AHI and both kbeam1,load (R=0.46, P=0.016) and kbeam1,unload (R=0.49, P=0.009). However, we found no significant relationships for either the Beam2 or Midfoot models.
In this study, we tested the hypothesis that the LA of the human foot functions like a spring with adjustable stiffness during running, similar to the leg overall. To do so, we experimentally manipulated the stride frequencies of participants during running and used three different models of the foot to estimate LA quasi-stiffness. As predicted, participants adjusted LA quasi-stiffness in a manner similar to overall leg quasi-stiffness, with higher frequency runs tending to yield significantly greater LA quasi-stiffness values during both loading and unloading phases of stance. These results were generally consistent for all three foot models, although none of the models showed significant differences in LA quasi-stiffness during loading between all pairs of conditions. The leg quasi-stiffness results fit with Farley and González's (1996) finding that humans tend to increase stiffness of the leg spring when running at higher stride frequencies, and show that this adjustment holds true for both the loading and unloading portions of stance phase. Additionally, we found that these adjustments are mirrored in changes to the quasi-stiffness at the knee and ankle, which had previously been demonstrated for hopping (Farley and Morgenroth, 1999). Similarly, Kessler et al. (2020) recently demonstrated that when increasing hopping frequency, humans increase midfoot quasi-stiffness in conjunction with the ankle. Their findings lend support to the conclusion that humans coordinate foot, ankle and knee stiffness to adjust the stiffness of the leg spring during running.
Inspection of the LA motion (Δzarch, Δθmid) and loading (Fbeam, Mmid) variables used to calculate LA quasi-stiffness reveals an interesting phenomenon: at higher stride frequencies, the LA tends to exhibit significantly reduced motion and loading. Because lower forces/moments should yield lower stiffness values, these results indicate that LA motion is being restricted disproportionately with loading as stride frequency increases. These results raise the possibility that increased muscle activation could be responsible for increases in dynamic LA stiffness during running, although EMG is needed to confirm this interpretation. This mechanism would be consistent with some previous EMG studies that have suggested that intrinsic foot muscles may be capable of resisting LA compression during running and static loading (Kelly et al., 2014, 2016). However, Farris et al. (2019) recently showed that anesthetizing the intrinsic foot muscles had no effect on midfoot quasi-stiffness during running, suggesting these muscles are not a major determinant of dynamic LA stiffness. Interestingly, Kessler et al. (2020) found that increased midfoot quasi-stiffness during hopping was accompanied by increased activation of one intrinsic foot muscle, abductor hallucis, but only during the flight (non-contact) phase of the hop. They suggested that abductor hallucis pre-activation could explain midfoot quasi-stiffness increases, and thus it is plausible that participants in the present study used a similar mechanism to stiffen the foot. Additionally, extrinsic foot muscles, such as the tibialis posterior, whose tendons cross multiple joints comprising the LA, could also potentially regulate LA stiffness. Thus far, these muscles have received little attention in in vivo studies of LA function in human walking and running, warranting further EMG research.
Outside of the foot muscles, greater engagement of the windlass mechanism during higher frequency runs could passively increase LA stiffness as well. According to the windlass mechanism, greater passive dorsiflexion of the metatarsophalangeal joints increases tension in the plantar aponeurosis, thereby stiffening the LA (Holowka and Lieberman, 2018), although we did not measure metatarsophalangeal joint motion to test this possibility. However, there are several reasons to doubt this mechanism explains our results. First, LA quasi-stiffness increased with stride frequency during the loading phase of stance, when the toes are not being dorsiflexed, and therefore are not engaging the windlass mechanism. Second, Welte et al. (2018) found that toe dorsiflexion actually reduces LA stiffness during static loading, complicating the possible role of the windlass mechanism in dynamic foot stiffness. Thus, further research on metatarsophalangeal joint motion is needed to fully evaluate the role of the windlass mechanism in LA stiffness during running.
Work and quasi-stiffness estimates differed markedly among the foot models used in this study. In particular, the three-point bending Beam1 model produced quasi-stiffness values that were roughly 3–4 times greater than those of the two-point bending Beam2 model, and work values that were up to 3.5 times greater than those of the Beam2 or Midfoot models. These differences were caused by the inclusion of triceps surae force on LA bending, resulting in high Fbeam1 values. According to the Beam1 model, the arch performs greater negative work than the ankle joint, and more than 50% of the negative work in the leg overall, during stance phase of running. This finding is a virtual impossibility, as our calculation of ankle work was based on a rigid foot assumption, which means that LA work estimates were essentially subsumed within overall ankle work calculations (see Kessler et al., 2020; Zelik and Honert, 2018). The Beam1 model produced positive LA work values that were roughly 50% of the ankle positive work, similar to those predicted by Kelly et al. (2018) using a 2D foot model, but far above the 10–15% predicted by Bruening et al. (2018), who used a 3D multi-segment foot model. In contrast, work estimates from Beam2 and Midfoot models were much lower, and produced similar negative work values to one another. The negative work estimated from these models (average of ∼0.2 J kg−1) was somewhat greater than that estimated using the three-point-bending model in Stearne et al. (2016) (∼0.15 J kg−1), and multi-segment midfoot models in Kelly et al. (2018) and Bruening et al. (2018) (∼0.05–0.15 J kg−1). The positive work estimates from the Beam2 and Midfoot models fell between those from the Bruening et al. (2018) and Kelly et al. (2018) studies.
Based on these comparisons, we conclude that the Beam1 model produces work values that are unrealistically high relative to those of the leg and ankle, and that the Beam2 and Midfoot models produce values that are more similar to those of previous studies. Our results suggests that inclusion of all triceps surae force into estimations of LA loading (e.g. Ker et al., 1987; Huang et al., 1993) may lead to considerable overestimates of LA stiffness and work. Paradoxically, Stearne et al. (2016), used a similar model to Beam1, yet estimated much lower amounts of negative work during running. However, whereas our model involved directly calculating work from measured force and deformation of the LA, they introduced an additional computational step requiring a predicted power function based on the static loading tests of Ker et al. (1987). While their model produced estimates of LA energy use that were in rough agreement with measurements of metabolic energy expenditure during running, it entailed an assumed relationship between loading and deformation of plantar tissues that was not entirely derived from their data. Thus, when using a beam model to estimate bending forces on the LA without the computational assumptions of the Stearne et al. (2016) model, we advocate excluding triceps surae force from calculations of work and stiffness (see Venkadesan et al., 2020; Welte et al., 2018). However, we recognize that our simplified model cannot capture the internal dynamics of the foot, and thus acknowledge the possibility that some component of triceps surae force may act to compress the LA during gait. More detailed human foot modeling is necessary to investigate this possibility.
Finally, our results indicate that LA height is not closely correlated with dynamic LA stiffness during running. The Beam1 model did show significant associations between AHI and quasi-stiffness, albeit with relatively low correlation coefficients (R<0.5). However, as discussed above, the Beam1 model may produce inaccurate estimates of LA stiffness, potentially making the Beam2 and Midfoot models more reliable. The correlation coefficient between AHI and kbeam2,unload showed a near-significant trend (P=0.051), suggesting the possibility that a weak relationship might have been detected with a larger sample size, but the Midfoot model showed no such trends. These findings do not support our prediction, the common assumption that high arches indicate stiff feet, and low arches indicate compliant feet (Franco, 1987; Neumann, 2002). However, all but one participant in our study possessed AHI values that fell within 2.5 s.d. of the mean from a previously reported large sample of humans (Butler et al., 2008), and thus it is possible that more extreme outliers in LA height conform more closely to typical expectations concerning LA stiffness. However, the lone outlier in our study had a very low LA but a relatively stiff foot during running, bucking this expectation. Thus, we conclude that within the normal range of anatomical variation, LA height does not serve as a good indicator of dynamic LA stiffness during running, and that other aspects of foot anatomy and potentially muscle activity are more important.
We acknowledge several limitations of our study, particularly regarding the foot models we used to estimate LA stiffness. As previously described, participants wore sandals designed for running, with straps that prevented the use of multi-segment foot model marker sets that are often employed to estimate intra-foot kinetics (e.g. Bruening et al., 2018; Farris et al., 2019; Kern et al., 2019). Consequently, we used a marker set and developed foot models (Beam1 and Beam2) similar to those used in other recent studies of LA energetics in which participants were shod (McDonald et al., 2016; Stearne et al., 2016). This marker set imposed several limitations on these models. First, we had to model the foot as a rigid segment, meaning that we calculated ankle moments and force between the shank and a foot segment encompassing both rearfoot and forefoot segments, rather than just between the shank and rearfoot. Previous research has shown that this leads to substantial overestimates of ankle motion, and hence ankle moments during walking, running and hopping (Kelly et al., 2018; Kessler et al., 2020; Zelik and Honert, 2018), which would result in overestimation of triceps surae force in our Beam1 model calculations. Our Beam2 model does not include triceps surae force, but does assume that the foot is a rigid segment in estimates of ankle joint reaction force, which may have had minor effects on stiffness and work calculations.
An additional source of error in our Beam1 model is estimation of triceps surae moment arm length. We used a kinematics-based approach that produced average moment arm lengths that were comparable to those reported in previous investigations of human triceps surae moment arms (Fig. S1; Maganaris, 2004; McCullough et al., 2011; Sheehan, 2012), but we acknowledge that our approach is susceptible to error due to marker placement. In our sensitivity analysis (Fig. S2), we found that changing moment arm length from −10% to +10% (∼1 cm) changes LA stiffness and work estimates by ∼15%. As the primary focus of this study involved within-participant effects, any errors in moment arm length estimation are unlikely to have affected our findings concerning differences among running frequencies, but they could undermine the accuracy of the absolute values calculated using the Beam1 model.
Our Midfoot model is a highly simplified proxy for 3D multi-segment foot models that calculate six degrees of freedom motion between forefoot and rearfoot segments. In using 2D representations of these segments, as well as the midfoot joint center, our model provides a coarse estimate of LA sagittal plane motion that is restricted to the medial side of the foot. These factors may explain why we calculated much higher ranges of LA motion during running than a recent study by Bruening et al. (2018) that used a multi-segment foot model to estimate LA kinetics. Another limitation of our model is that we were unable to precisely assign ground reaction force to regions posterior and anterior to our designated midfoot joint, and therefore resorted to an approach previously utilized where we calculated the midfoot moment only after the center of pressure under the foot had passed anterior to the joint center (Bruening and Takahashi, 2018; Farris et al., 2019; Kelly et al., 2018; Kern et al., 2019). We also neglected forefoot mass and acceleration in our estimates of midfoot moment, although these would have likely only had trivial effects on our calculations. Despite these limitations, our Midfoot model produced similar average ranges of motion and moments, as well as comparable estimates of LA stiffness to those of the simple 2D foot model used in Kelly et al. (2018), as well as the multi-segment foot model used in Farris et al. (2019). Further, the average work estimates from our Midfoot and Beam2 models were generally in the range of those from previous studies that estimated work within the midfoot (Bruening et al., 2018; Kelly et al., 2018; Stearne et al., 2016). Regardless, the overall goal of this study was to compare foot stiffness estimates across stride frequencies within participants, and across participants with different LA heights. Any effect on the accuracy of our calculations due to our model limitations should have been systematic, and therefore likely would not change the general outcome of this study.
In conclusion, the results presented here show that humans adjust the stiffness of the LA in coordination with the ankle and knee joints to help regulate the overall stiffness of the leg spring during running. While we demonstrated this capability in an artificial manipulation of the stride frequency, it could be examined further by measuring humans running on surfaces with different compliance to determine whether they modulate LA stiffness in a manner similar to the leg overall (Ferris et al., 1998; Kerdok et al., 2002). If so, the LA could serve as a proximate mechanism for optimizing the leg's mechanical response to the underlying surface. For example, when walking or running on soft or pliable surfaces, a stiffer LA might enhance the force transmission necessary for propulsion, while on hard surfaces, a more compliant arch might allow the foot to better dampen potentially damaging ground reaction forces. Recently, Kelly et al. (2016) found some support for this notion, finding less LA deformation and more foot muscle activity when people ran in cushioned shoes compared with running barefoot, which could reflect the use of foot muscles to stiffen the LA when running on a relatively compliant surface. The ability to actively adjust LA stiffness could also be advantageous for other behaviors such as accelerating and decelerating, running with added mass (see also Kern et al., 2019) and moving on uneven surfaces (Grimmer et al., 2008). However, determining the extent to which the adjustments to LA stiffness observed in this study are due to active (i.e. muscle) versus passive (e.g. windlass) mechanisms will require future research with detailed EMG and more sophisticated foot models. Our results do suggest that at least one passive aspect of foot anatomy, LA height, is unlikely to have a major effect on dynamic foot stiffness. Thus, understanding human foot biomechanics, and the adaptive function of the human LA, necessitates continued investigation during dynamic loading conditions common to normal human gait.
We thank Andrew Biewener, Eamon Callison, Tim Kistner, Nicolai Konow, Freddy Sichting, Victoria Tobolsky and Ian Wallace for helpful discussion of this project, Michael Ruiz for assistance with data collection, and Madhusudhan Venkadesan and Andrew Yegian for comments and feedback on earlier drafts of the manuscript.
Conceptualization: N.B.H., D.E.L.; Methodology: N.B.H., A.R.; Formal analysis: N.B.H., A.R.; Investigation: N.B.H., A.R., B.E.S.; Writing - original draft: N.B.H.; Writing - review & editing: N.B.H., A.R., B.E.S., D.E.L.
Support for this research was provided by the American School of Prehistoric Research at Harvard University.
The authors declare no competing or financial interests.