Experimental estimation of energy absorption during heel strike in human barefoot walking

Metabolic energy expenditure during human gait is poorly understood. Mechanical energy loss during heel strike contributes to this energy expenditure. Previous work has estimated the energy absorption during heel strike as 0.8 J using an effective foot mass model. The aim of our study is to investigate the possibility of determining the energy absorption by more directly estimating the work done by the ground reaction force, the force-integral method. Concurrently another aim is to compare this method of direct determination of work to the method of an effective foot mass model. Participants of our experimental study were asked to walk barefoot at preferred speed. Ground reaction force and lower leg kinematics were collected at high sampling frequency (3000 Hz; 1295 Hz), with tight synchronization. The work done by the ground reaction force is 3.8 J, estimated by integrating this force over the foot-ankle deformation. The effective mass model is improved by dropping the assumption that foot-ankle deformation is maximal at the instant of the impact force peak. On theoretical grounds it is clear that in the presence of substantial damping that peak force and peak deformation do not occur simultaneously. The energy absorption results, due the vertical force only, corresponding to the force-integral method is similar to the results of the improved application of the effective mass model (2.7 J; 2.5 J). However the total work done by the ground reaction force calculated by the force-integral method is significantly higher than that of the vertical component alone. We conclude that direct estimation of the work done by the ground reaction force is possible and preferable over the use of the effective foot mass model. Assuming that energy absorbed is lost, the mechanical energy loss of heel strike is around 3.8 J for preferred walking speeds (≈ 1.3 m/s), which contributes to about 15–20% of the overall metabolic cost of transport.


Introduction
The metabolic cost of human walking is substantial. Not only has it been reported to account for about 25% of the daytime energy expenditure of an office clerk [1], it is also well established that walking distance of many people with a locomotor impairment, for instance following a stroke, is limited due to the associated metabolic cost [2][3][4]. Thus, it is important for the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 evaluation and treatment of gait disorders to have a thorough understanding of the processes underlying the metabolic energy expenditure during human locomotion, and to be able to measure the energy associated with these processes [5]. It is accepted that positive and (to a lesser extent) negative muscle fiber mechanical work contribute to the metabolic cost of walking [6]. While it is debated in which phase positive muscle fiber mechanical work is primarily done [6][7][8][9][10][11][12], it is accepted that, during steady-motion level walking, the positive muscle fiber mechanical work serves to compensate for negative muscle fiber mechanical work, for negative work associated with air friction (generally assumed to be negligible), and for mechanical energy lost in the foot-ground contact. There has been recent literature that indicates that soft tissue deformation could play a larger role in cost of transport that was previously thought [13,14]. For running, [14] estimate about 19 J per stance phase during running at 3 m / s associated with soft tissue energy loss. [13] estimated a total 13 J of energy expenditure in the collision phase of walking, where around 7 J could not be accounted for by modeled joint or segment work and thought to be attributed by soft tissue deformation. We want to build on these findings and further investigate a major contributor to this unaccounted energy loss during the collision phase of walking.
In this study, we focus on quantification of the mechanical energy absorbed in the footground contact during steady motion horizontal walking, and in particular on the energy exchange between heel and ground that occurs during heel strike, i.e. during the first 10-50 ms of the double support phase. [15,16]. We define the heel strike to start at heel contact and to end when the vertical foot deformation is maximal or equivalently when the vertical foot deformation speed is zero. Where others, primarily motivated by experimental limitations, defined heel strike to end at 90% or 100% of the impact transient peak seen in the vertical ground reaction force [15,16]. [15] argued that direct experimental estimation of the energy absorbed during heel strike is difficult because it is hard to achieve the required tight synchronization between data on ground reaction force (GRF) and kinematics. As an alternative, these authors analyzed an effective foot mass model similar to the one proposed by [17], and suggested that during heel strike, the behavior is governed by the effective foot mass M eff that is coupled visco-elastically to the ground. Assuming that the velocity of the effective foot mass is zero at the instant of the impact peak in the vertical GRF, they estimated the mechanical energy absorption during heel strike to be 0.8 J. Assuming a muscle mechanical efficiency of 20% for positive muscle fiber mechanical work [18] and assuming all of the energy absorbed by the heel pad is lost, this value would explain about 4% of the metabolic cost per step during steady-motion level walking (defined according to [6]). This suggests that energy loss during heel strike is negligible.
We will show that direct experimental estimation of the mechanical energy absorption during heel strike is possible, using a force-integral method. We will compare the results of the force-integral method to the results obtained using both the original and an improved application of the effective mass model.

Methods
Participants were asked to walk barefoot at their preferred walking speed on level ground in a large hall where a force platform was mounted in the floor. Ten trials were recorded for each participant. The experiment was focused on heel contact of the right foot. GRF and lower leg kinematics were collected at high sampling frequency, with tight synchronization. This data was used to estimate the work done by the GRF and the energy absorption corresponding to the original and an improved application of the effective mass model. All of the raw data, the processed data, a high speed video made during the pilot experiments together with Matlab scripts used for processing and generating figures can be found in our online repository [19].

Participants
In total 12 (4 male, 8 female) healthy participants (age: 27 ± 5 years; mass: 76.7 ± 15 kg; BMI: 24.5 ± 5.0 kg / m 2) were involved in the experiment. Prior to the experiment, the participants were asked to sign an informed consent form. The experiment was approved by the ethical board of the Human Movements Sciences Faculty at the Free University of Amsterdam.

Instrumentation
GRF was measured using a piezoelectric Kistler force plate (Type 9282B11) mounted permanently in the floor, with a sampling frequency of 3000 Hz. Deformation due to heel strike is primarily attributed to the soft tissue in the heel pad, however [20] has shown by comparing in vivo and in vitro tests of the heel pad and lower leg that deformation of the rest of the lower leg cannot simply be neglected in lower extremity impact tests. Therefore we chose to take the entire foot-ankle deformation into account. The kinematics of 2 LED's were measured using an Optotrak Certus motion capture system, with a sampling frequency of 1295 Hz. The number of Optotrak LEDs was reduced to two in order to maximize Optotrak sample frequency. As the ankle joint angle has been reported to remain more or less constant during heel strike [15], the two LEDs were placed on the skin over the medial surface of the tibia, with an inter-LED distance of around 0.16 m (see Fig 1); skin movement artifact at this location has been reported to be very small [21]. The movement of these LEDs was used to describe the footankle deformation, which is defined as the deformation between the heel contact point and the tibia.
During pilot experiments, an impact test was used to determine the error in synchronization between force measurements and measurements of kinematics. The maximum synchronization error was found to be less than 0.001 s. During pilot experiments the lower leg motion during heel strike was also captured using a high speed camera (sampling frequency of 1200 Hz) in order to determine if the constant ankle angle assumption during heel strike holds (for video see [19]). Visual inspection of the images during heel strike supported the assumption that the change in ankle angle during heel strike is negligible.

Data analysis
Optotrak data was filtered using a low pass forward-backward filter (cut-off frequency 100 Hz) in order to remove high frequency noise. Two coordinate systems were used, the first being the global XYZ coordinate system with origin O, the second being the local xyz coordinate system with origin L that moves with the rigid part of the foot-lower-leg system (see Fig 1). This local system was used to define a heel point H, the position of which is constant in this local xyz coordinate system (r H is constant). The net GRF, as well as the X and Y-coordinate of the point of application of the vertical component of the GRF ('center of pressure' COP) were calculated using standard methods using the 8 Kistler force plate output channels (4 for the vertical force components and 2 for each horizontal component corresponding to the lines of action through the four sensors embedded in the force plate) [22,23].
The location of the heel point H was defined using this COP position relative to the global XYZ coordinate system as a basis. First the time-variant COP position was transformed from the global XYZ coordinate system to the local xyz coordinate system by establishing and subsequently applying the time-variant rotation and translation between these coordinate systems [24]. The median xyz-position of these local COP positions surrounding heel impact was taken as the constant location of the heel contact point r H relative to the local xyz coordinate system. The first sample included in this calculation was the first sample in which the vertical GRF exceeded 25% of the impact peak force; the last sample included in this calculation was chosen such that the impact peak is in the middle of the samples used. The r H vector was calculated for each trial separately. The time-variant location of the heel point R H relative to the global XYZ coordinate system is found by applying the inverse time-variant rotation and translation readily established between the global and local coordinate system to the constant r H [24].
Force-integral method. Direct estimation of the work done by the GRF during heel strike requires that foot-ankle deformation was determined from the kinematic data. To that end, we have defined a 'heel point' H, the position of which is time invariant relative to the xyz coordinate system (r H ) and time variant relative to the XYZ coordinate system (R H ). We assumed that foot-ankle deformation S equals the displacement of this heel point H relative to the XYZ coordinate system (S(t) = R H (t) − R H (t 0 ), where t 0 is the instant of heel strike at which GRF is zero) (see Figs 1 and 2). Subsequently, the displacement of this heel point relative to the XYZ coordinate system was calculated from the kinematic data. Finally, energy W absorbed by the foot-ankle using the force-integral method was calculated by numerically integrating the GRF The direction X is lateral, Y is anterior and Z is vertical up. A local coordinate system xyz is defined, where the origin L is situated at an LED kinematic marker. This is the lower of two markers placed directly on the skin over the medial surface of the tibia. The local z direction points to the second LED marker, the y direction is perpendicular to the z axis and lies in the plane of progression (defined as a plane parallel to the global YZ-plane), the x direction is defined orthogonal to the yz-plane to complete the local right-handed orthogonal coordinate system. The position of the heel point H is derived from the location of the center of pressure as obtained from the force plate data during heel strike, as measured in the global coordinate system XYZ. These global center of pressure points are transformed to the local coordinate system xyz. The median xyz-position of these transformed local center of pressure points is taken as the position of the heel point r H , which is a time-invariant vector in the local lower-leg coordinate system xyz. The global position of the heel point R H is time variant, due to the movement of the local coordinate system in the global coordinate system. vector F over the foot-ankle deformation vector S: Where t e is the instant when the vertical heel velocity is zero (Ż H (t e ) = 0). We refer to the interval from t 0 to t e as 'heel strike'. The values of these integrals were determined using a trapezoidal method. Improved application of the effective mass model. The energy change was also determined using the effective foot mass M eff model [15,17]. The M eff in the effective foot mass model represents the part of the body mass that participates in the initial phase after footground impact and is only acted upon by the GRF. The effective foot mass model is concerned only with motions and forces in the vertical Z direction. In the original effective mass model, it is also assumed that the period of heel strike ends at the time of the 'impact peak' (t p ) in the GRF and it is further assumed that the vertical velocity Ż H is zero at the impact peak time t p , making the model applicable to kinematic and force data that are not synchronized [15]. Below we will argue that, due to the viscous properties of the heel pad, the impact peak in the force data does not occur at zero vertical heel velocity. Since we have kinematic and force data which is tightly synchronized, we are not forced to make this assumption. Therefore in this improved application of the M eff model, we used the synchronized kinematic data to identify the end of heel strike to be at Ż H (t e ) = 0. The effective foot mass corresponding to this improved application of the model M e was determined from the impulse-momentum Where Ż H (t 0 ) and Ż H (t e ) are the vertical heel velocities at the beginning and end of heel strike respectively, which are ascertained by using a gradient approach on the vertical heel point position Z H (which is the vertical component of R H ), and g is the gravitational acceleration (taken as g = 9.81 m / s 2). The total change in energy ΔE e , consisting of a change in kinetic ΔE kin and potential ΔE pot energy, of this M e were determined in order to indirectly determine work W done on the heel pad.
Where t e is end time of heel strike and (Z H (t e ) − Z H (t 0 )) is the maximal vertical foot-ankle deformation, which is the vertical component of S(t e ).
In order to investigate the influence of this improvement on the application of the effective mass model and be able to compare our results to [15], we evaluated the original model under the assumption that heel strike ends at the instant of the impact force peak t p and that Ż H (t p ) = 0 to estimate the corresponding effective mass M p_CS and energy absorption due to heel strike ΔE p_CS and also evaluated the improved application of the effective mass model until the instant of the impact force peak t p (where Ż H (t p ) 6 ¼ 0) to estimate the corresponding effective mass M p and energy absorption due to heel strike ΔE p .

Statistical analysis
The results of the 10 trials per participant were grouped together and the mean and standard deviation were calculated of these groups. The mean of the subject means gives the overall mean, the mean of the subject standard deviations was used as a measure for the intrasubject variability and the standard deviation of the subject means was used as a measure for the intersubject variability.
A Kolmogorov-Smirnov test on the subject means for the energy absorption during heel strike (W, W Z , ΔE e , ΔE p and ΔE p_CS ) indicated that the underlying distributions are unlikely to be normal (p(W) = 2.8E-10, p(W Z ) = 1.9E-9, p(ΔE e ) = 2.1E-9, p(ΔE p ) = 7.6E-7, p(ΔE p_CS ) = 1.8E-5). Therefore a Friedman test was performed to test the statistical significance of the difference between the methods. Wilcoxon signed rank tests were used to investigate specific differences. For these paired tests a Bonferroni correction was applied to set the level of significance to 0.01.

Parameter sensitivity analysis
A parameter sensitivity analysis was performed on the estimationr H of the heel point H position r H in the local coordinate system. A new heel point Q was defined, which is located approximately 1.7 cm away from point H by moving 1 cm downwards (-z), 1 cm posterior (-y) and 1 cm lateral (x), which is assumed to be further away than any reasonable inaccuracies expected in ascertaining the heel point H. In order to verify the validity of this assumption, we have compared the location of this point Q to the various estimated locations of the COP points during heel strike in the local xyz coordinate system. Point Q was used in order to calculate the energy absorption during heel strike in the same manner as for point H and the results for these two points were compared.

Results
During heel strike, the variation in inter-marker distance was in the same order of magnitude as the accuracy of the Optotrak system (0.1 mm). The preferred walking speed of the participants was consistent over trials (Table 1), and the average value was consistent with previous research [25].
A typical GRF measurement is shown in Fig 3. In this study we were interested only in the early part (% 0.03 s) of the double support phase (% 0.1 s), where the heel strike occurs. In Fig  4 the GRF is shown together with the foot-ankle deformation for the same measurement as Fig  3. The zero vertical heel velocity occurred substantially later than the impact peak; this substantial difference between t p and t e was seen in all the measurements ( Table 1). As can be seen from Fig 4, the vertical foot-ankle deformation indicates that the vertical heel velocity remained more or less constant during the initial part of heel strike (Table 1: on average Ż H (t 0 ) = −0.57 m / s and Ż H (t p ) = −0.52 m / s ) and decreased quickly after the impact force peak. The maximum vertical foot-ankle deformation S Z (t e ) was around −13 ± 2 mm, which is a compressive deformation. The corresponding shear deformations S X (t e ) and S Y (t e ) were −7 ± 4 mm and 14 ± 7 mm respectively, together bringing the total foot-ankle deformation at the end of heel strike to |S(t e )| = 21 ± 6 mm. Fig 5 shows the progression of the COP in the local xyz coordinate system together with the defined heel point H and point Q used for a parameter sensitivity analysis (see Section Parameter sensitivity analysis). This figure shows that the point H was located within the area of the COP samples and point Q was located at a considerable distance from it.
Using the force-integral method (Eq 1) the total work done by the GRF vector was found to be W = −3.8 ± 1.7 J. Herein the largest component was the vertical component W Z = −2.7 ± 1.1 J, then the lateral component W X = −0.8 ± 0.7 J and the anterior component W Y = −0.4 ± 0.3 J. For some trials there were components of this calculated work which had a small positive value (6% for W X , 13% for W Y , 0% for W Z of all trials), which results in a small positive mean value for a single subject for a single component (subject 3: W X = 0.016 J), which must be due to measurement and/or modeling errors.
The original effective foot mass model predicted an energy absorption of ΔE p_CS = −0.8 ± 0.4 J, with corresponding foot mass M p_CS = 4.1 ± 0.6% of bodyweight. The improved application of the effective foot mass model resulted in a substantially higher energy absorption of ΔE p = −1.4 ± 0.1 J, with corresponding foot mass M p = 18 ± 5% of bodyweight and ΔE e = −2.5 ± 0.9 J, with corresponding foot mass M e = 11 ± 3% of bodyweight.
In this study the overall difference between the methods of determining energy absorption due to heel strike (W, W Z , ΔE e , ΔE p , ΔE p_CS ) was found to be statistically significant (χ 2 (59) = 45.9, p = 2.6E-9). Where the result of the vertical component of the force-integral method (W Z ) and the improved application of the effective foot mass model (ΔE e ) did not differ significantly (Z = −1.2, p = 0.24), however all the other method combinations did show statistically significant differences (Z = −3, p = 2.2E-3).
In order to determine the sensitivity of the results for the choice of heel point, results were also calculated for point Q (see Section Parameter sensitivity analysis); energy absorption for point Q was found to be 2-5% smaller than for the default heel point and results of statistical tests were similar.
In Fig 6 the relation between maximum vertical foot-ankle compression |S Z (t e )| and energy absorption per heel strike W Z is shown together with experimental results collected by [26] (from [27][28][29][30]). In these studies, in vivo tests were performed on human constrained lower legs with bare feet to determine compressive mechanical properties of the heel pad (and lower leg [20]). Typically, energy input was varied and resulting maximal lower leg compression was

Discussion
In this study, we estimated the mechanical energy absorption during heel strike, by directly calculating the work done by the GRF on the foot by using the force-integral method. Our results on the relation between foot-ankle compression and compressive energy absorption are similar to previous mechanical lower leg property results [20] collected by [26]. Similar values for energy absorption results, due the vertical force only, were obtained from an improved application of the effective mass model of [15]; we found no significant difference in the results using the improved application of the effective foot mass model and the vertical component of the force-integral method. We did however find significant differences between these methods and the (total) force-integral method and the original effective foot mass model. We were able to reproduce very similar results using the original effective mass model to those reported by [15], however these are substantially lower than the values obtained using the other methods in this study.
Without synchronized kinematic and force data, [15] had to make assumptions about the relation between force data and kinematic data during heel strike. A key assumption in their approach was that the vertical heel velocity is zero at the instant of the impact peak and chose this instant as the final instant in their analysis. However this assumption is not consistent with our measurements and should be expected to be incorrect on theoretical grounds  whenever the relative damping is non-negligible [31]. The observation that application of the original effective foot mass model in this study yielded values for energy absorption that are similar to those reported by [15] suggests that our heel point adequately represents the heel center studied by [15]; the observation that the original effective mass model yielded values for energy absorption that are substantially lower than those obtained with both the improved application of the effective mass model and the force-integral method proposed here, indicates that the erroneous assumption in the original effective mass model had substantial consequences for the estimated energy absorption.  1), where the subscript X corresponds to mediolateral, Y to anterior-posterior and Z to vertical component (same measurement as in Fig 3). The start of heel strike is shown by a vertical dotted line indicated by t 0 , the instant of the vertical force peak is indicated by t p and the end of heel strike, which is defined as zero vertical heel velocity, by t e . https://doi.org/10.1371/journal.pone.0197428.g004 Experimental estimation of energy absorption during heel strike A possible source of overestimation of the foot-ankle deformation is skin motion artifact, since the markers could both have moved downward relative to the bone after impact. However we saw no evidence of damped oscillation after impact in the raw marker data (raw data can be found at [19]) and the marker distance variation was negligible. It is therefore assumed that skin movement artifact is not the dominant error source. Another limitation of our research is that we assume that the contact area during heel strike is a constant point, whereas in reality rolling of the contact point, albeit with small amplitude, is likely to occur during heel strike. This means that a small part of the estimated foot-ankle deformation, should be in fact attributed to the movement of the contact point in the local coordinate system due to rolling of the contact point. Furthermore, analysis of the sensitivity of energy absorption for the position of the heel point has indicated that the methods for determining energy absorption proposed here do not depend strongly on the choice of the heel contact point. We therefore expect that these limitations will only have a relatively small contribution to the methodical error.
Our experimental approach depends on the assumption that the lower leg plus foot constitute a single body with constant configuration (no joint rotation) during heel strike; based on Experimental estimation of energy absorption during heel strike this assumption, the time-invariant location of the heel point in the local xyz coordinate system was determined. This assumption was previously made by [15], and was supported by inspection of high speed video data of the foot during heel strike, made during pilot experiments [19]. However we would reason that the error in our estimates is dominated by the neglected  Experimental estimation of energy absorption during heel strike ankle rotation during initial contact, which is expected to be smaller than 5˚ [32]. Estimating the distance of the heel point to the ankle at 60 mm, this implies that in the worst case the heel point would move 5 mm in the local coordinate system due to the neglected ankle rotation. Considering that the vector from ankle joint axis to heel point is oriented nearly vertically during heel strike, we expect the effect of the ankle joint rotation on vertical displacement of the heel point to be very small. This would result in a small overestimation of the foot-ankle deformation and therefore also a small overestimation on the energy absorption estimation. We would therefore suggest to investigate other marker configurations, with at least one marker on the foot, in particular when studying abnormal gait patterns. We expect that the limitations discussed above result in a small bias in the results of both methods, therefore the comparison between methods is not compromised.
The compressive mechanical characteristics of the heel and lower leg (energy absorption as a function of maximum lower leg compression [20]) as estimated in this study are quite similar to those reported in previous research [26]. For some subjects the foot-ankle compression is slightly higher in the current study. This could be due to the fact that body mass and BMI of the participants in this study varied widely, whereas half of the data points reported in [26] were obtained from a single male recreational runner. Furthermore, activities such as running have been shown to influence the mechanical heel pad properties [33]. Differences could also be related to methodological issues, both in our study and in the studies evaluated by [26].
Given that the vertical component of the force-integral method proposed in this study and the improved application of the effective mass model yield similar values for energy absorption, one may ask which of these two methods is to be preferred in future studies. It is clear that the force-integral method entails less assumptions than the effective mass method, and on that ground alone the force-integral method is to be preferred. Furthermore, both methods require tightly synchronized high-frequency force and position data. The improved application of the effective mass model could possibly be useful in answering other research questions related to foot strike mechanics. However, we note that the movement of the heel point (which is assumed to represent the point mass in the effective mass model) during heel strike is not correctly predicted by the improved application of the effective mass model; in particular, the acceleration of the heel point is by no means proportional to the GRF, as they should be according to Newton's second law. The force-integral method has as added strength that the total energy absorption due to heel strike, including the approximate shear losses corresponding to the lateral and anterior components of the GRF, can be estimated. All together, we suggest to use the force-integral method outlined in this study in future attempts to quantify the energy absorbed by the foot or heel pad during heel strike.
It should be noted that in this study we have only quantified the amount of energy absorbed by the heel pad, foot and ankle during its compression in the impact stage, the initial footground contact. In later stages of the stance phase other parts of the heel and the rest of the foot will deform and absorb energy. This study sheds no light on the question to what extent this energy is buffered in elastic structures, allowing it to be returned to the body at a later stage, and to what extent it is converted into heat by viscous structures. However, this issue has been addressed by [26], who suggest that 75-95% of the energy absorbed by the heel pad and lower leg is dissipated; only a small fraction of the energy is stored in elastic structures. It is an open question if the elastic energy stored in the foot can be effectively returned to the system at a later stage of the stance phase. Although it has been shown that energy stored in a rotational spring at the ankle can be beneficial to the metabolic cost of walking [34,35], it remains to be investigated to what extent this occurs for elastic energy stored in the heel pad, foot and ankle. At this point, we estimate the mechanical energy loss per step due to heel strike to be about 3.8 J.
Results of this study have implications for the importance of energy loss during heel strike in the energetics of human walking. According to this study, energy loss is around 3.8 J of mechanical energy per step during barefoot steady-motion level walking at preferred velocity. Assuming a mechanical efficiency for positive muscle fiber work of 20% [18], and considering that the total metabolic energy expenditure per step is in the order of 120 J [6], this corresponds to 15-20% of the total metabolic energy expended per step. In contrast to previous suggestions, mechanical energy lost during foot-ground contact contributes substantially to the metabolic energy cost of walking.

Conclusion
We conclude that direct estimation of the work done by the ground reaction force is possible and preferable over the use of the effective foot mass model. The mechanical energy loss of heel strike is around 3.8 J for preferred walking speeds (% 1.3 m / s ), which contributes to about 15-20% of the overall metabolic cost of transport.