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).

Fig. 1.

Anatomy of the longitudinal arch of the human foot. (A) Radiograph of the longitudinal arch (LA). (B) The plantar aponeurosis attaches to the calcaneus posteriorly, and has two bands with dense, longitudinally oriented fibers. The lateral band attaches to the joint capsule of the fifth tarsometatarsal joint, and the central band divides into five slips that attach distally to the proximal phalanges and nearby subcutaneous tissues. (C) The intrinsic foot muscles are deep to the plantar aponeurosis, and include the flexor digitorum brevis (FDB) and abductor hallucis (AH), which have been linked to LA stiffness. These muscles originate on the calcaneus, and insert onto the proximal phalanges.

Fig. 1.

Anatomy of the longitudinal arch of the human foot. (A) Radiograph of the longitudinal arch (LA). (B) The plantar aponeurosis attaches to the calcaneus posteriorly, and has two bands with dense, longitudinally oriented fibers. The lateral band attaches to the joint capsule of the fifth tarsometatarsal joint, and the central band divides into five slips that attach distally to the proximal phalanges and nearby subcutaneous tissues. (C) The intrinsic foot muscles are deep to the plantar aponeurosis, and include the flexor digitorum brevis (FDB) and abductor hallucis (AH), which have been linked to LA stiffness. These muscles originate on the calcaneus, and insert onto the proximal phalanges.

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.

Participants

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.

Data collection

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.

Fig. 2.

Models of forces on theLAduring running used in this study. (A) Foot marker set, with footwear (Luna MONO sandals) used by participants in the study. (B) Marker set used to create a plane representing the sole of the foot, based on projections of markers on the first metatarsal head, fifth metatarsal head and the posterior calcaneus, with vertical height of the navicular tuberosity marker measured relative to this plane. (C) Marker set used to calculate the projected sagittal plane midfoot angle between the first metatarsal head, the navicular tuberosity and the posterior calcaneus. (D) Beam1 model: ground reaction force (FGRF) and triceps surae force (FTS) load the LA in three-point bending, causing a compressive force (Fbeam1) on the arch and a change in LA height (Δzarch; see B). (E) Beam2 model: compressive force on the arch (Fbeam2) is based on two-point bending from FGRF alone. (F) Midfoot model: FGRF creates a plantarflexion moment about a midfoot ‘joint’ (Mmid) located at the navicular tuberosity, causing angular motion at the joint (Δθmid; see C).

Fig. 2.

Models of forces on theLAduring running used in this study. (A) Foot marker set, with footwear (Luna MONO sandals) used by participants in the study. (B) Marker set used to create a plane representing the sole of the foot, based on projections of markers on the first metatarsal head, fifth metatarsal head and the posterior calcaneus, with vertical height of the navicular tuberosity marker measured relative to this plane. (C) Marker set used to calculate the projected sagittal plane midfoot angle between the first metatarsal head, the navicular tuberosity and the posterior calcaneus. (D) Beam1 model: ground reaction force (FGRF) and triceps surae force (FTS) load the LA in three-point bending, causing a compressive force (Fbeam1) on the arch and a change in LA height (Δzarch; see B). (E) Beam2 model: compressive force on the arch (Fbeam2) is based on two-point bending from FGRF alone. (F) Midfoot model: FGRF creates a plantarflexion moment about a midfoot ‘joint’ (Mmid) located at the navicular tuberosity, causing angular motion at the joint (Δθmid; see C).

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.

Beam models

Ker et al. (1987) and Huang et al. (1993) modeled the LA as a beam undergoing three-point bending from forces applied posteriorly by the triceps surae muscles via their attachment to the calcaneus, and anteriorly by ground reaction force at the foot's point of contact with the ground. To calculate these forces during running in vivo, we slightly modified the methodology from Stearne et al. (2016) (Fig. 2D). Stiffness of the LA can be calculated as bending force on the LA (Fbeam) relative to the vertical compression of the arch (Δzarch). Following previous studies (Ker et al., 1987; Huang et al., 1993; Stearne et al., 2016), the bending force in the Beam1 model, Fbeam1, was calculated as:
formula
(1)
where FTS is the bending force on the LA from the triceps surae, and FGRF is the bending force caused by the ground reaction force vector. In this study, FTS was calculated as:
formula
(2)
where Mankle is the net ankle joint moment in the sagittal plane, which we assume to represent the ankle moment produced by the triceps surae muscles, and rTS is the moment arm of these muscles. We calculated Mankle using inverse dynamics, assuming the foot is a rigid segment, which was a necessary limitation of our analysis given our marker set. We determined rTS as the instantaneous horizontal distance between a marker placed on the heel at the insertion of the Achilles tendon and the virtual ankle joint center (the midpoint between markers placed on the medial and lateral malleoli), when both were projected onto the segment coordinate system for the shank. We multiplied the quotient of Mankle and rTS by the cosine of sagittal plane ankle angle (θankle) to calculate the component of the force produced by the triceps surae that would compress the arch vertically (thereby excluding the shear components of this force). Average moment arm (rTS) calculations for all participants are presented in Fig. S1, and are comparable to those that have been published in previous investigations of human leg muscle moment arms (Maganaris, 2004; McCullough et al., 2011; Sheehan, 2012). We also investigated the sensitivity of this method by adjusting the measured rTS values by ±10%, and assessing the effects on resulting force, quasi-stiffness and work calculations (Fig. S2). Finally, we calculated FGRF using inverse dynamics as the vertical component of the ankle joint reaction force projected into the plane of the sole of the foot. This variable incorporates the effects of acceleration of the foot segment relative to the leg on the forces compressing the arch due to ground reaction force.

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 xy 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.

Midfoot model

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.

Work calculations

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.

Statistical analysis

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.

Table 1.

Average spatiotemporal, kinematic and kinetic variables across running conditions

Average spatiotemporal, kinematic and kinetic variables across running conditions
Average spatiotemporal, kinematic and kinetic variables across running conditions

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.

Fig. 3.

AverageLAkinematics and kinetics during stance phase across running frequencies. (A) Fbeam1 across stance phase and (B) boxplot comparing maximum Fbeam1 among running frequencies. (C) Δzarch across stance phase, with ‘0’ set as the value at initial foot contact, and (D) boxplot comparing total range of Δzarch among running frequencies. (E) Fbeam2 across stance phase and (F) boxplot comparing maximum Fbeam2 among running frequencies. (G) Mmid across stance phase and (H) boxplot comparing maximum Mmid among running frequencies. (I) Δθmid across stance phase and (J) boxplot comparing total range of Δθmid among running frequencies. Lines (A,C,E,G,I) represent the average among participants (N=27) at each percentage of stance phase. Boxplots (B,D,F,H,J) show median, upper and lower quartiles and 1.5× interquartile range. Asterisks indicate significant differences between running condition pairs.

Fig. 3.

AverageLAkinematics and kinetics during stance phase across running frequencies. (A) Fbeam1 across stance phase and (B) boxplot comparing maximum Fbeam1 among running frequencies. (C) Δzarch across stance phase, with ‘0’ set as the value at initial foot contact, and (D) boxplot comparing total range of Δzarch among running frequencies. (E) Fbeam2 across stance phase and (F) boxplot comparing maximum Fbeam2 among running frequencies. (G) Mmid across stance phase and (H) boxplot comparing maximum Mmid among running frequencies. (I) Δθmid across stance phase and (J) boxplot comparing total range of Δθmid among running frequencies. Lines (A,C,E,G,I) represent the average among participants (N=27) at each percentage of stance phase. Boxplots (B,D,F,H,J) show median, upper and lower quartiles and 1.5× interquartile range. Asterisks indicate significant differences between running condition pairs.

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.

Fig. 4.

LAquasi-stiffness calculated from different foot models across running frequencies. Quasi-stiffness values (k) were calculated as slopes of least squares regression lines fitted to motion and force data during loading and unloading periods of stance. (A–C) Beam1 model, (D–F) Beam2 model and (G–I) Midfoot model. Lines (A,D,G) depict average values calculated from all participants (N=27), with arrows indicating loading (upwards) and unloading (downwards) periods. Boxplots (N=27, except E, where N=24) depict stiffness during loading (B,E,H) and unloading (C,F,I). Asterisks indicate significant differences among running condition pairs.

Fig. 4.

LAquasi-stiffness calculated from different foot models across running frequencies. Quasi-stiffness values (k) were calculated as slopes of least squares regression lines fitted to motion and force data during loading and unloading periods of stance. (A–C) Beam1 model, (D–F) Beam2 model and (G–I) Midfoot model. Lines (A,D,G) depict average values calculated from all participants (N=27), with arrows indicating loading (upwards) and unloading (downwards) periods. Boxplots (N=27, except E, where N=24) depict stiffness during loading (B,E,H) and unloading (C,F,I). Asterisks indicate significant differences among running condition pairs.

Table 2.

Average quasi-stiffness values across running conditions

Average quasi-stiffness values across running conditions
Average quasi-stiffness values across running conditions

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.

Fig. 5.

Leg, ankle and knee quasi-stiffness calculated across running frequencies. Quasi-stiffness values were calculated as slopes of least squares regression lines fitted to motion and force data during loading and unloading periods of stance. (A–C) Leg overall, (D–F) ankle joint and (G–I) knee joint. Lines (A,D,G) depict average values calculated from all participants (N=27), with arrows indicating loading (upwards) and unloading (downwards) periods. Boxplots (N=27, except H, where N=23) depict stiffness during loading (B,E,H) and unloading (C,F,I). Asterisks indicate significant differences among running condition pairs.

Fig. 5.

Leg, ankle and knee quasi-stiffness calculated across running frequencies. Quasi-stiffness values were calculated as slopes of least squares regression lines fitted to motion and force data during loading and unloading periods of stance. (A–C) Leg overall, (D–F) ankle joint and (G–I) knee joint. Lines (A,D,G) depict average values calculated from all participants (N=27), with arrows indicating loading (upwards) and unloading (downwards) periods. Boxplots (N=27, except H, where N=23) depict stiffness during loading (B,E,H) and unloading (C,F,I). Asterisks indicate significant differences among running condition pairs.

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.

Fig. 6.

Comparisons of work calculation from foot models, leg and joints (N=27). Average total positive (above the line) and negative (below the line) work per unit body mass across stance phase calculated from foot models (Beam1, Beam2, Midfoot), the leg, ankle joint and knee joint. Error bars represent ±1 s.d. Numbers indicate means±s.d., and units are J kg−1.

Fig. 6.

Comparisons of work calculation from foot models, leg and joints (N=27). Average total positive (above the line) and negative (below the line) work per unit body mass across stance phase calculated from foot models (Beam1, Beam2, Midfoot), the leg, ankle joint and knee joint. Error bars represent ±1 s.d. Numbers indicate means±s.d., and units are J kg−1.

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.

Fig. 7.

Relationship between LA height and quasi-stiffness. Scatterplots depicting relationships between arch height index (AHI) and LA quasi-stiffness calculated from Beam1 (A,D), Beam 2 (B,E) and Midfoot (C,F) models during loading (A–C) and unloading (D–F) periods of stance across all participants (N=27, except B, where N=24). Quasi-stiffness variables were log-transformed to achieve normality. Solid lines indicate least squares linear regression fits, and shaded regions represent 95% confidence intervals. R and P values are from Pearson's product moment correlation association tests.

Fig. 7.

Relationship between LA height and quasi-stiffness. Scatterplots depicting relationships between arch height index (AHI) and LA quasi-stiffness calculated from Beam1 (A,D), Beam 2 (B,E) and Midfoot (C,F) models during loading (A–C) and unloading (D–F) periods of stance across all participants (N=27, except B, where N=24). Quasi-stiffness variables were log-transformed to achieve normality. Solid lines indicate least squares linear regression fits, and shaded regions represent 95% confidence intervals. R and P values are from Pearson's product moment correlation association tests.

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.

Author contributions

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.

Funding

Support for this research was provided by the American School of Prehistoric Research at Harvard University.

Alexander
,
R. M.
(
1991
).
Energy-saving mechanisms in walking and running
.
J. Exp. Biol.
160
,
55
-
69
.
Alexander
,
R. M.
and
Jayes
,
A. S.
(
1983
).
A dynamic similarity hypothesis for the gaits of quadrupedal mammals
.
J. Zool.
201
,
135
-
152
.
Altman
,
A. R.
and
Davis
,
I. S.
(
2012
).
A kinematic method for footstrike pattern detection in barefoot and shod runners
.
Gait Posture
35
,
298
-
300
.
Arampatzis
,
A.
,
Brüggemann
,
G.-P.
and
Metzler
,
V.
(
1999
).
The effect of speed on leg stiffness and joint kinetics in human running
.
J. Biomech.
32
,
1349
-
1353
.
Arampatzis
,
A.
,
Brüggemann
,
G.-P.
and
Klapsing
,
G. M.
(
2001
).
Leg stiffness and mechanical energetic processes during jumping on a sprung surface
.
Med. Sci. Sports Exerc.
33
,
923
-
931
.
Arndt
,
A.
,
Wolf
,
P.
,
Liu
,
A.
,
Nester
,
C.
,
Stacoff
,
A.
,
Jones
,
R.
,
Lundgren
,
P.
and
Lundberg
,
A.
(
2007
).
Intrinsic foot kinematics measured in vivo during the stance phase of slow running
.
J. Biomech.
40
,
2672
-
2678
.
Bates
,
D.
,
Mächler
,
M.
,
Bolker
,
B.
and
Walker
,
S.
(
2015
).
Fitting linear mixed-effects models using lme4
.
J. Stat. Softw.
67
,
1
-
48
.
Blickhan
,
R.
(
1989
).
The spring-mass model for running and hopping
.
J. Biomech.
22
,
1217
-
1227
.
Bruening
,
D. A.
and
Takahashi
,
K. Z.
(
2018
).
Partitioning ground reaction forces for multi-segment foot joint kinetics
.
Gait Posture
62
,
111
-
116
.
Bruening
,
D. A.
,
Pohl
,
M. B.
,
Takahashi
,
K. Z.
and
Barrios
,
J. A.
(
2018
).
Midtarsal locking, the windlass mechanism, and running strike pattern: a kinematic and kinetic assessment
.
J. Biomech.
73
,
185
-
191
.
Butler
,
R. J.
,
Hillstrom
,
H.
,
Song
,
J.
,
Richards
,
C. J.
and
Davis
,
I. S.
(
2008
).
Arch height index measurement system: establishment of reliability and normative values
.
J. Am. Podiatr. Med. Assoc.
98
,
102
-
106
.
Cappozzo
,
A.
,
Catani
,
F.
,
Della Croce
,
U.
and
Leardini
,
A.
(
1995
).
Position and orientation in space of bones during movement: anatomical frame definition and determination
.
Clin. Biomech.
10
,
171
-
178
.
Carlson
,
R. E.
,
Fleming
,
L. L.
and
Hutton
,
W. C.
(
2000
).
The biomechanical relationship between the tendoachilles, plantar fascia and metatarsophalangeal joint dorsiflexion angle
.
Foot Ankle Int.
21
,
18
-
25
.
Clites
,
T. R.
,
Arnold
,
A. S.
,
Singh
,
N. M.
,
Kline
,
E.
,
Chen
,
H.
,
Tugman
,
C.
,
Billadeau
,
B.
,
Biewener
,
A. A.
and
Herr
,
H. M.
(
2019
).
Goats decrease hindlimb stiffness when walking over compliant surfaces
.
J. Exp. Biol.
222
,
jeb198325
.
Derrick
,
T. R.
,
van den Bogert
,
A. J.
,
Cereatti
,
A.
,
Dumas
,
R.
,
Fantozzi
,
S.
and
Leardini
,
A.
(
2020
).
ISB recommendations on the reporting of intersegmental forces and moments during human motion analysis
.
J. Biomech.
99
,
109533
.
Erdemir
,
A.
,
Hamel
,
A. J.
,
Fauth
,
A. R.
,
Piazza
,
S. J.
and
Sharkey
,
N.
(
2004
).
Dynamic loading of the plantar aponeurosis in walking
.
J. Bone Joint Surg. Am.
86
,
546
-
552
.
Farley
,
C. T.
and
González
,
O.
(
1996
).
Leg stiffness and stride frequency in human running
.
J. Biomech.
29
,
181
-
186
.
Farley
,
C. T.
and
Morgenroth
,
D. C.
(
1999
).
Leg stiffness primarily depends on ankle stiffness during human hopping
.
J. Biomech.
32
,
267
-
273
.
Farley
,
C. T.
,
Blickhan
,
R.
,
Saito
,
J.
and
Taylor
,
C. R.
(
1991
).
Hopping frequency in humans: a test of how springs set stride frequency in bouncing gaits
.
J. Appl. Physiol.
71
,
2127
-
2132
.
Farley
,
C. T.
,
Glasheen
,
J.
and
McMahon
,
T.
(
1993
).
Running springs: speed and animal size
.
J. Exp. Biol.
185
,
71
-
86
.
Farris
,
D. J.
,
Kelly
,
L. A.
,
Cresswell
,
A. G.
and
Lichtwark
,
G. A.
(
2019
).
The functional importance of human foot muscles for bipedal locomotion
.
Proc. Natl. Acad. Sci.
116
,
1645
-
1650
.
Ferris
,
D. P.
,
Louie
,
M.
and
Farley
,
C. T.
(
1998
).
Running in the real world: adjusting leg stiffness for different surfaces
.
Proc. R. Soc. B Biol. Sci.
265
,
989
-
994
.
Franco
,
A. H.
(
1987
).
Pes cavus and pes planus: analyses and treatment
.
Phys. Ther.
67
,
688
-
694
.
Grimmer
,
S.
,
Ernst
,
M.
,
Günther
,
M.
and
Blickhan
,
R.
(
2008
).
Running on uneven ground: leg adjustment to vertical steps and self-stability
.
J. Exp. Biol.
211
,
2989
-
3000
.
Günther
,
M.
and
Blickhan
,
R.
(
2002
).
Joint stiffness of the ankle and the knee in running
.
J. Biomech.
35
,
1459
-
1474
.
He
,
J. P.
,
Kram
,
R.
and
McMahon
,
T. A.
(
1991
).
Mechanics of running under simulated low gravity
.
J. Appl. Physiol.
71
,
863
-
870
.
Holowka
,
N. B.
and
Lieberman
,
D. E.
(
2018
).
Rethinking the evolution of the human foot: insights from experimental research
.
J. Exp. Biol.
221
,
jeb174425
.
Holowka
,
N. B.
,
Wallace
,
I. J.
and
Lieberman
,
D. E.
(
2018
).
Foot strength and stiffness are related to footwear use in a comparison of minimally- vs. conventionally-shod populations
.
Sci. Rep.
8
,
3679
.
Huang
,
C.-K.
,
Kitaoka
,
H. B.
,
An
,
K.-N.
and
Chao
,
E. Y. S.
(
1993
).
Biomechanical evaluation of longitudinal arch stability
.
Foot Ankle
14
,
353
-
357
.
Kaufman
,
K. R.
,
Brodine
,
S. K.
,
Shaffer
,
R. A.
,
Johnson
,
C. W.
and
Cullison
,
T. R.
(
1999
).
The effect of foot structure and range of motion on musculoskeletal overuse injuries
.
Am. J. Sports Med.
27
,
585
-
593
.
Kelly
,
L. A.
,
Cresswell
,
A. G.
,
Racinais
,
S.
,
Whiteley
,
R.
and
Lichtwark
,
G.
(
2014
).
Intrinsic foot muscles have the capacity to control deformation of the longitudinal arch
.
J. R. Soc. Interface
11
,
20131188
.
Kelly
,
L. A.
,
Lichtwark
,
G.
,
Cresswell
,
A. G.
and
Cresswell
,
A. G.
(
2015
).
Active regulation of longitudinal arch compression and recoil during walking and running
.
J. R. Soc. Interface
12
,
20141076
.
Kelly
,
L. A.
,
Lichtwark
,
G. A.
,
Farris
,
D. J.
and
Cresswell
,
A.
(
2016
).
Shoes alter the spring-like function of the human foot during running
.
J. R. Soc. Interface
13
.
Kelly
,
L. A.
,
Farris
,
D. J.
,
Lichtwark
,
G. A.
and
Cresswell
,
A. G.
(
2018
).
The influence of foot-strike technique on the neuromechanical function of the foot
.
Med. Sci. Sports Exerc.
50
,
98
-
108
.
Kelly
,
L. A.
,
Farris
,
D. J.
,
Cresswell
,
A. G.
and
Lichtwark
,
G. A.
(
2019
).
Intrinsic foot muscles contribute to elastic energy storage and return in the human foot
.
J. Appl. Physiol.
126
,
231
-
238
.
Ker
,
R. F.
,
Bennett
,
M. B.
,
Bibby
,
S. R.
,
Kester
,
R. C.
and
Alexander
,
R. M.
(
1987
).
The spring in the arch of the human foot
.
Nature
325
,
147
-
149
.
Kerdok
,
A. E.
,
Biewener
,
A. A.
,
McMahon
,
T. A.
,
Weyand
,
P. G.
and
Herr
,
H. M.
(
2002
).
Energetics and mechanics of human running on surfaces of different stiffnesses
.
J. Appl. Physiol.
92
,
469
-
478
.
Kern
,
A. M.
,
Papachatzis
,
N.
,
Patterson
,
J. M.
,
Bruening
,
D. A.
and
Takahashi
,
K. Z.
(
2019
).
Ankle and midtarsal joint quasi-stiffness during walking with added mass
.
PeerJ
7
,
e7487
.
Kessler
,
S. E.
,
Rainbow
,
M. J.
,
Lichtwark
,
G. A.
,
Cresswell
,
A. G.
,
D'Andrea
,
S. E.
,
Konow
,
N.
,
Kelly
,
L. A.
and
Kelly
,
L. A.
(
2019
).
A direct comparison of biplanar videoradiography and optical motion capture for foot and ankle kinematics
.
Front. Bioeng. Biotech.
7
,
199
.
Kessler
,
S. E.
,
Lichtwark
,
G. A.
,
Welte
,
L. K. M.
,
Rainbow
,
M. J.
and
Kelly
,
L. A.
(
2020
).
Regulation of foot and ankle quasi-stiffness during human hopping across a range of frequencies
.
J. Biomech.
108
,
109853
.
Kristianslund
,
E.
,
Krosshaug
,
T.
and
van den Bogert
,
A. J.
(
2012
).
Effect of low pass filtering on joint moments from inverse dynamics: implications for injury prevention
.
J. Biomech.
45
,
666
-
671
.
Latash
,
M. L.
and
Zatsiorsky
,
V. M.
(
1993
).
Joint stiffness: myth or reality?
Hum. Mov. Sci.
12
,
653
-
692
.
Leardini
,
A.
,
Benedetti
,
M. G.
,
Berti
,
L.
,
Bettinelli
,
D.
,
Nativo
,
R.
and
Giannini
,
S.
(
2007
).
Rear-foot, mid-foot and fore-foot motion during the stance phase of gait
.
Gait Posture
25
,
453
-
462
.
Lenth
,
R. V.
(
2016
).
Least-squares means: the R package lsmeans
.
J. Stat. Softw.
69
,
1
-
33
.
Lieberman
,
D. E.
,
Venkadesan
,
M.
,
Werbel
,
W. A.
,
Daoud
,
A. I.
,
D'Andrea
,
S.
,
Davis
,
I. S.
,
Mang'Eni
,
R. O.
and
Pitsiladis
,
Y.
(
2010
).
Foot strike patterns and collision forces in habitually barefoot versus shod runners
.
Nature
463
,
531
-
535
.
Liew
,
B. X. W.
,
Morris
,
S.
,
Masters
,
A.
and
Netto
,
K.
(
2017
).
A comparison and update of direct kinematic-kinetic models of leg stiffness in human running
.
J. Biomech.
64
,
253
-
257
.
Lundgren
,
P.
,
Nester
,
C.
,
Liu
,
A.
,
Arndt
,
A.
,
Jones
,
R.
,
Stacoff
,
A.
,
Wolf
,
P.
and
Lundberg
,
A.
(
2008
).
Invasive in vivo measurement of rear-, mid- and forefoot motion during walking
.
Gait Posture
28
,
93
-
100
.
Maganaris
,
C. N.
(
2004
).
Imaging-based estimates of moment arm length in intact human muscle-tendons
.
Eur. J. Appl. Physiol.
91
,
130
-
139
.
Mann
,
R. A.
and
Inman
,
V. T.
(
1964
).
Phasic activity of the intrinsic muscles of the foot
.
J. Bone Jt. Surg. Am.
46
,
469
-
481
.
McCullough
,
M. B. A.
,
Ringleb
,
S. I.
,
Arai
,
K.
,
Kitaoka
,
H. B.
and
Kaufman
,
K. R.
(
2011
).
Moment arms of the ankle throughout the range of motion in three planes
.
Foot Ankle Int.
32
,
300
-
306
.
McDonald
,
K. A.
,
Stearne
,
S. M.
,
Alderson
,
J. A.
,
North
,
I.
,
Pires
,
N. J.
and
Rubenson
,
J.
(
2016
).
The role of arch compression and metatarsophalangeal joint dynamics in modulating plantar fascia strain in running
.
PLoS ONE
11
,
e0152602
.
McKeon
,
P. O.
,
Hertel
,
J.
,
Bramble
,
D.
and
Davis
,
I.
(
2014
).
The foot core system: a new paradigm for understanding intrinsic foot muscle function
.
Br. J. Sports Med.
49
,
290
.
Neumann
,
D. A.
(
2002
).
Ankle and foot
. In
Kinesiology of the Musculoskeletal System
(
ed.
D. A.
Neumann
), pp.
477
-
521
.
Philadelphia
:
Mosby
.
Perl
,
D. P.
,
Daoud
,
A. I.
and
Lieberman
,
D. E.
(
2012
).
Effects of footwear and strike type on running economy
.
Med. Sci. Sports Exerc.
44
,
1335
-
1343
.
R Core Team
. (
2016
).
R: A Language and Environment for Statistical Computing
.
Vienna
:
R Foundation for Statistical Computing
.
Reeser
,
L. A.
,
Susman
,
R. L.
and
Stern
,
J. T.
(
1983
).
Electromyographic studies of the human foot: experimental approaches to hominid evolution
.
Foot Ankle
3
,
391
-
407
.
Riddick
,
R.
,
Farris
,
D. J.
and
Kelly
,
L. A.
(
2019
).
The foot is more than a spring: human foot muscles perform work to adapt to the energetic requirements of locomotion
.
J. R. Soc. Interface
16
,
20180680
.
Rouse
,
E. J.
,
Gregg
,
R. D.
,
Hargrove
,
L. J.
and
Sensinger
,
J. W.
(
2013
).
The difference between stiffness and quasi-stiffness in the context of biomechanical modeling
.
IEEE Trans. Biomed. Eng.
60
,
562
-
568
.
Shamaei
,
K.
,
Sawicki
,
G. S.
and
Dollar
,
A. M.
(
2013
).
Estimation of quasi-stiffness and propulsive work of the human ankle in the stance phase of walking
.
PLoS ONE
8
,
e59935
.
Sheehan
,
F. T.
(
2012
).
The 3D in vivo Achilles’ tendon moment arm, quantified during active muscle control and compared across sexes
.
J. Biomech.
45
,
225
-
230
.
Sichting
,
F.
,
Holowka
,
N. B.
,
Ebrecht
,
F.
and
Lieberman
,
D. E.
(
2020
).
Evolutionary anatomy of the plantar aponeurosis in primates, including humans
.
J. Anat.
237
,
85
-
104
.
Simkin
,
A.
,
Leichter
,
I.
,
Giladi
,
M.
,
Stein
,
M.
and
Milgrom
,
C.
(
1989
).
Combined effect of foot arch structure and an orthotic device on stress fractures
.
Foot Ankle
10
,
25
-
29
.
Stearne
,
S. M.
,
McDonald
,
K. A.
,
Alderson
,
J. A.
,
North
,
I.
,
Oxnard
,
C. E.
and
Rubenson
,
J.
(
2016
).
The foot's arch and the energetics of human locomotion
.
Sci. Rep.
6
,
19403
.
Takahashi
,
K. Z.
,
Worster
,
K.
and
Bruening
,
D. A.
(
2017
).
Energy neutral: the human foot and ankle subsections combine to produce near zero net mechanical work during walking
.
Sci. Rep.
7
,
15404
.
Venkadesan
,
M.
,
Yawar
,
A.
,
Eng
,
C. M.
,
Dias
,
M. A.
,
Singh
,
D. K.
,
Tommasini
,
S. M.
,
Haims
,
A. H.
,
Bandi
,
M. M.
and
Mandre
,
S.
(
2020
).
Stiffness of the human foot and evolution of the transverse arch
.
Nature
579
,
97
-
100
.
Welte
,
L.
,
Kelly
,
L. A.
,
Lichtwark
,
G. A.
and
Rainbow
,
M. J.
(
2018
).
Influence of the windlass mechanism on arch-spring mechanics during dynamic foot arch deformation
.
J. R. Soc. Interface
15
,
20180270
.
Williams
,
D. S.
, III
,
McClay
,
I. S.
and
Hamill
,
J.
(
2001
).
Arch structure and injury patterns in runners
.
Clin. Biomech.
16
,
341
-
347
.
Winter
,
D. A.
(
2005
).
Biomechanics and Motor Control of Human Movement
, 3rd edn.
Hoboken
:
Wiley
.
Zelik
,
K. E.
and
Honert
,
E. C.
(
2018
).
Ankle and foot power in gait analysis: implications for science, technology and clinical assessment
.
J. Biomech.
75
,
1
-
12
.
Zifchock
,
R. A.
,
Davis
,
I.
,
Hillstrom
,
H.
and
Song
,
J.
(
2006
).
The effect of gender, age, and lateral dominance on arch height and arch stiffness
.
Foot Ankle Int.
27
,
367
-
372
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information