Modeling the Foot-Strike Event in Running Fatigue via Mechanical Impedances

1. Introduction

The human motor system benefits from remarkable muscular redundancies: A motor task is normally performed with the simultaneous involvement of more muscles than strictly necessary. Furthermore, this same task may be performed in multiple ways, with different muscle combinations. From the mechanical viewpoint the musculoskeletal system is indeterminate, whereby the number of unknown muscle forces exceeds the number of available equations. We address in this Chapter the biomechanics of the lower limbs in long-distance running under conditions of developing fatigue. In long-distance running the running speed may result in more than 300 foot-strikes per leg per kilometer. Each such foot-strike evokes an impact loading that results in a vertical shock impulse transmitted upwards through the body and carries with it the potential for injuries in the bone and joint tissues.

Fatigue, or stress, fractures occur in bones in response to repetitive stresses over multiple cycles, when the body’s ability to adapt is exceeded [1,2]. An important factor which affects the incidence of bone stress injury, is exposure to abrupt changes in the bone loading [1], and consequent alteration in the strain distribution [3] with insufficient recovery periods [4]. Other factors include footwear, terrain and intensity of activity or training [1].

Two of the major factors responsible for impulse attenuation at foot- or heel-strike are the shock absorption capacity of the active muscle in the lower limbs, and the cushioning effect of the foot heel-pad tissue. In previous reports we have shown that in long distance running the impact shock load on the lower limbs increases with progressing fatigue [5-8]. One additional question is whether, as a result of fatigue, an imbalance between the activities of the plantar and dorsi flexor muscles of the ankle develops. Such an imbalance would compromise the protective action provided by the muscles to the shank [9-11].

The goal of this research was to characterize the heel-strike shock propagation and attenuation in running by means of a biomechanical model, and to examine changes taking place as a result of running fatigue.

2. Biomechanical modeling of the lower limb

This section deals with the modeling of the heel-strike event. With the development of biomechanical models of human body motion, it has become possible to simulate vertical landing, such as occurring during running, in order to gain insight into intermuscular coordination and to elucidate control strategies of the musculoskeletal system. A common method to deal with this type of problems is to lump together elements of the human body e.g., muscles, tendons, ligaments, bones and joints so that the overall musculoskeletal system is represented as a damped elastic mechanism. Several models describing vertical landing can be found in the literature [12-18]. These models are usually characterized by the presence of elastic springs and viscous dampers, with constant properties and provide a reasonable prediction of the maximal vertical foot/ground reaction force.

Indeterminacy of the locomotor system can be addressed by adopting the lumping method, whereby the material elements of the human body e.g., muscles, tendons, ligaments, bones and joints are lumped together so that the overall musculoskeletal system is represented as a multi-degree-of freedom damped elastic mechanism, interconnecting the masses of the body segments.

The foot- or heel-strike period during landing from fall, during hopping or during the stance phase of running has been generally modeled using one-dimensional models along the vertical direction. [13-15,19-20].

In this study we represent the body segments during heel-strike by a four degree-of-freedom elastically-damped uni-axial biomechanical model. The model thus includes 4 masses connected by elastic stiffnesses with parallel damping elements, as shown in Figure 1. In more details, the masses m_j (j = 1..4) represent, respectively, the foot + shoe, the shank, the thigh and the rest of the body (including the non-supporting leg). Each of the stiffness k_j and damping c_j (j = 1..4) represent, respectively, the lumped effects of the heel-pad + sole, the ankle joint, the knee joint and the hip joint.

3. Model equations

For the above model, the force diagram, as presented in Figure 2, yields the following model equations.

$$\{\begin{array}{l}m1\ddot{x}1=-k1x1-c1\dot{x}1-k2(x1-x2)-c2(\dot{x}1-\dot{x}2)-m1g\\ m2\ddot{x}2=k2(x1-x2)+c2(\dot{x}1-\dot{x}2)-k3(x2-x3)-c3(\dot{x}2-\dot{x}3)-m2g\\ m3\ddot{x}3=k3(x2-x3)+c3(\dot{x}2-\dot{x}3)-k4(x3-x4)-c4(\dot{x}3-\dot{x}4)-m3g\\ m4\ddot{x}4=k4(x3-x4)+c4(\dot{x}3-\dot{x}4)-m4g\end{array}$$E1

with initial conditions:

$$xi(0)=0\text{ ; }\dot{x}i(0)=v0=-1m/s\text{, }i=\mathrm{1,2,3,4}$$E2

and gravitational acceleration

$$g=-9.81m/s^2$$

These values rely on reported landing velocities between -0.8 m/s to -1.2 m/s for running speeds of 3.5 m/s (comparable to the speeds of this study), while wearing various types of running shoes [21-23].

The above masses are expressible in terms of the total body mass from anthropometric data [24].

$$\{\begin{array}{l}{m}_1=0.0145\cdot m\\ m_2=0.0465\cdot m\\ m_3=0.100\cdot m\\ m_4=1-(m_1+m_2+m_3)=0.839\cdot m\end{array}$$E3

From the simultaneous recording of the foot ground reaction forces and accelerations on the masses m_j, information about the rise time of the peak acceleration can be obtained.

Figure 1.

Lumped model including 4 masses connected by elastic stiffnesses with parallel dampings. The masses m_j (j = 1..4) represent, respectively, the foot + shoe, the shank, the thigh and the rest of the body (including the non-supporting leg). Each of the stiffness k_j and damping c_j (j = 1..4) represent, respectively, the lumped effects of the heel-pad + sole, the ankle joint, the knee joint and the hip joint.

Figure 2.

The force diagram for the model elements presented in Figure 1.

4. Experimental setup

Information about the impulsive loading along the skeletal elements in long-distance running can be non-invasively obtained from the foot-ground reactive forces [25] and, more directly, by measuring the transient accelerations on the shank caused by impact.

4.1. Impact accelerations

Non-invasive in vivo measurements of acceleration and impact transmission along the human body were made by externally attaching light-weight, high-sensitivity accelerometers at strategic points including bony prominences, such as the tibial tuberosity below the knee area, the greater trochanter near hip level and the sacrum area at the lower back [13,15,22,26-29].

In this study, each subject was instrumented with two light-weight (4.2 grams) uniaxial (Kistler PiezoBeam, type 8634B50, Kistler, Winterthur, Switzerland), skin-mounted accelerometers connected to a coupler (Kistler Piezotron, type 5122). One was attached on the tibial tuberosity, and the second - on the sacrum. To achieve good reliability of the measurements by means of bone-mounted accelerometers, the accelerometers were pressed onto the skin in closest position to the bony prominences of the tibial tuberosity and the sacrum, by means of two elastic belts passed in a horizontal plane around the shank and the waist, respectively. The tensions of the belts were well above the level in which the acceleration trace for a given impact force became insensitive to the accelerometer attachment force, thus ensuring stability of the accelerometer as well as consistency of the readings and reproducibility of the data [13,30].

The shank accelerometer was aligned with the axis of the tibia to provide the axial component of the tibial acceleration and the accelerometer on the sacrum was oriented along the spine. These accelerometers allowed us to acquire the shock accelerations propagated in the longitudinal directions of the tibia and the spine. As earlier reported, such attachment is suitable for faithfully measuring the amplitude of shock acceleration [5-8].

Force platforms, type Kistler Z-4035, were used for the simultaneous recordings of the foot-ground reaction forces and acceleration.

4.2. Running fatigue tests

An overview of the experimental setup is shown in Figure 3. For examining the effect of global fatigue due to running, the subjects were asked to run on a Quinton Q55 treadmill.

Figure 3.

Description of the experimental apparatus with a subject running on a treadmill. Two accelerometers are attached, one on the tibial tuberosity (Tib. Tub.) and the other on the sacrum (Sac.) ; the accelerometers’ data are sampled through an amplifier and A/D converter to a PC. Likewise, the respiratory data are sampled and collected on-line.

Global, or metabolic fatigue is associated with the development of metabolic acidosis following an endurance exercise and is accompanied by a decrease in the end tidal carbon dioxide pressure (PETCO₂) [31]. In long distance running metabolic fatigue is reached when the running speed exceeds the anaerobic threshold [31].

Running was thus for a duration of 30 min and at a speed exceeding the anaerobic threshold speed of each subject by 5%. Before the test a 15 min warming up running on the treadmill was performed. In this study, the average running speed for all 14 subjects was 3.53 m/s (SD, 0.19). It should be noted, however, that in addition to global fatigue, local fatigue in a muscle may also take place as a result of an intensive activity of this muscle. This type of fatigue is reflected by certain changes in its electromyogram (EMG) signal in the time and/or frequency domains [32]. Local fatigue was not considered in this study.

Respiratory data were collected from a Sensor-Medics 4400 device and included $${\text{<mover accent="true">V <mo>˙</mo></mover>}}_{{\text{O}}_{\text{2}}}$$- Oxygen consumption, $${\text{<mover accent="true">V <mo>˙</mo></mover>}}_{{\text{CO}}_{\text{2}}}$$- CO₂ production, $${\text{<mover accent="true">V <mo>˙</mo></mover>}}_{\text{E}}$$- minute ventilation, PETCO₂ - end-tidal CO₂ pressure, $$\dot{V}E/\dot{V}{O}_2$$- ventilatory equivalent for oxygen, $$\dot{V}E/\dot{V}C{O}_2$$- ventilatory equivalent for CO₂. The anaerobic threshold was determined as the point of initial increase of $$\dot{V}E/\dot{V}{O}_2$$ and$$\dot{V}E/\dot{V}C{O}_2$$, which just precedes the initial decline of PETCO₂ [31]. In a previous study, we have shown that 30 min running at a speed exceeding the anaerobic threshold is a sufficient time to induce general fatigue [8].

The respiratory data were evaluated at each of the 1^st, 5^th, 10^th, 15^th, 20^th, 25^th and 30^th min of running and the accelerometer and force platform data were online sampled at 1667 Hz sampling rate. The model parameters were estimated, however, at the 1^st, 15^th and 30^th min of running.

The dynamics of acceleration build-up at heel-strike is shown in Figure 4 where the simultaneous recordings of the tibial tuberosity acceleration and the ground reaction force (GRF) are shown in two time scales: complete running cycle (panel a) and zooming-in on the heel-strike event (panel b). In this case the build-up time to the tibial tuberosity peak acceleration was ~ 30 ms. It is also noted that the ground reaction force exhibits two peaks: a smaller one shortly after heel-strike and a larger one (~ 2.5 body weights), towards the middle of the running cycle.

Figure 4.

Typical foot-ground reaction force (GRF) (transient curve) measured simultaneously with tibial tuberosity acceleration (spiky curve). a. complete stance phase; b. heel-strike period only.

Typical accelerometer traces for the tibial tuberosity (right leg) and sacrum are shown in Figure 5, for a complete running cycle, i.e., from heel-strike of the right foot till the next heel-strike of the same foot. Two major differences should be noted between the traces: (a) intensity of ~ 8 g in the tibial tuberosity versus less than 3 g in the sacrum; (b) while the tibial tuberosity exhibits one major peak within the first 50 ms of the running cycle from heel-strike and originating from the heel-strike of the right foot, the sacrum acceleration, due to its central location, exhibits two comparable positive peaks within the running cycle, reflecting each of the right and left heel-strikes.

Figure 5.

Typical accelerometer traces for the tibial tuberosity (T.T.) (right leg) and sacrum, for a complete running cycle, i.e., from heel-strike of the right foot till the next heel-strike of the same foot.

5. Parameter estimation of the model

The mechanical properties of biological material are, in general, multiple variable-dependent. Specifically stiffness, in addition to its being non-linear e.g. strain dependent, often depends on the deformation rate. This is also the case with bones [33], tendons and ligaments [34], cartilage [35] and muscle [36]. Damping too may be position-dependent.

Due to nonlinearity of the stiffness/damping properties of the joints of the leg [e.g. 20,37], we were not generally able to estimate the model parameters while assuming that they remain constant over the heel-strike period. Thus, the heel-strike period was divided into two equal periods (22 ms each) and the parameters were estimated separately for each of these periods, using the Gauss-Marquardt [38-39] method of non-linear estimation. For the first period the initial conditions were as prescribed in equation (2), and for the second period the initial conditions used were the values reached at the end of the first period.

Figure 6 shows the model prediction of the shank mass (m₂ in Figure 1) acceleration (continuous line) versus the experimentally measured tibial tuberosity acceleration (denoted as +). The results of two different subjects are demonstrated (vertical partition of the Figure) for two time points of running (horizontal partition): 1^st min, fatigue-free, and 15^th min of the running test. The subtle, brief, discontinuity in the traces at mid-cycle represents the transition between the two parameter estimation periods.

Figure 6.

The model prediction of the m₂ acceleration (continuous line) versus the experimentally measured tibial tuberosity acceleration (denoted as +). The upper panels represent the solution for the 1^st min of running and the lower panels represent the solution for the 15^th min of running. The subtle, short, discontinuity in the traces at mid-cycle represents the transition between the two parameter estimation periods, with piece-wise constant stiffness each. Each of the left and right panels represents a different subject.

6. Model results and model sensitivity analysis

Table 1 shows the sample results of the stiffness and damping parameters for the two time zones. A sensitivity analysis can provide an indication to the quality of the estimation of parameters. Thus, sensitivity of the m₂ and m₄ model results to each of the estimated parameters was performed by two-fold multiplying and/or dividing each of the estimated parameters separately and for each time zone. The two-fold variation demonstrated a strong sensitivity of m₂ to k₁ and k₂ in the first time zone and to parameters k₁, k₂, k₃ and c₃ in the second time zone. Small to medium sensitivity was obtained in the second time zone to parameters c₁ and k₄, repectively. The m₄ acceleration was not sensitive in the first time zone to none of the parameters, while in the second time zone it demonstrated strong sensitivity to parameters k₂ and k₃ and small sensitivity to parameters k₄ and c₃.

Table 1.

Summary of the sensitivity of the m₂ and m₄ accelerations to the model parameters. A twofold change up and down was checked in the two time solution-zones, separately for each parameter. Parameters for which the accelerations were not sensitive were checked also for a change in one order of magnitude (up and down). The results of this latter test are shown in parentheses (* = negligibly low effect; ** = low to medium effect; *** = high effect, i.e. more than 10% change in peak acceleration of m₂ or m₄).

In cases of low- or no-sensitivity the parameters were varied, up and down, by one order of magnitude with the results shown in parentheses. The one order of magnitude variation in the parameter values did not evoke sensitivity beyond that of the twofold variation, except for c₁ on m₂ and k₁ on m₄.

Due to the fact that the parameter estimation of the model coefficients was performed in short time intervals, fractions of the heel strike event, the damping coefficients disclosed a high variability. Better estimations would have probably resulted if the period of estimation was higher. Accordingly, the stiffness parameters k₁ and k₂, in the first time zone, and k₁, k₂ and k₃ in the second time zone, considered more reliably estimated, were reported in what follows.

Difficulties in estimating the damping coefficients are not unusual due to their expected low values. It has been reported that in repetitive physical activity, such as in running, the subject bounces on the ground in a spring-like manner [17, 40-45]. Depending on the range of joint flexion and on the frequency of motion, a considerable amount of elastic energy can be stored and re-used. It has been shown that the dissipated energy in muscles increase when the amplitudes of joint movement are increased [46]. It has also been also commented that the utilization of stored elastic energy depended on the shortness in latency between the stretch and shortening phases of the muscles [47]. Accordingly, during the ground-contact period of running, the leg was modeled as a one-dimensional four-degree-of-freedom piece-wise linear spring, with no damping. During heel-strike, the joints did not have a damping effect, to contribute to energy dissipation.

Summary of the values of these parameters for the two solution time zones is given in the rightmost column of Table 2. The values, averages for 10 subjects (SD), are calculated at each of the 1^st, 15^th and 30^th min of running to evaluate the effect of fatigue. The asterisks indicate a significant change p < 0.05 between the values of zone 1 and zone 2. The differences between k₁ in zone 1 and k₁ in zone 2 were significant with p < 0.0007, for each of the 1^st, 15^th and 30^th min of running. For the k₂ parameter a statistically significant difference was obtained for the 1^st min of running only, with p < 0.009. The differences in parameter values due to fatigue were not statistically different, despite the notable differences in the averages. The reason was the big variability among the tested subjects. On the individual level, though, the differences due to fatigue were statistically significant in most of the cases.

It should be noted that the stiffness k₁, relating to the heel-pad may be alternatively obtained from the ratio between the foot-ground reaction force and the heel-pad deformation. Past reports using this method have reported an increased stiffness with the heel-pad deformation, as occurs during heel-strike [27,48-52]. These studies have reported an approximately tenfold increase in heel-pad stiffness, i.e., from 6.6 – 135 kN/m to 77 – 1430 kN/m, using different methodologies including actual running in conjunction with measurements of the heel pad by means of X-rays; heel-strike simulation by in vivo pendulum impact, or Instron measurements made on cadavers. Despite the differences with the method applied in this study, the stiffness k₁ indicates similar values with those obtained by other methods. Our model provides further the fatigue effect on k₁, which is found opposite in the first time zone, where the stiffness decreases with fatigue, to the second time zone where the stiffness is found to increase.

Table 2.

Parameter means (S.D) of 10 subjects: k₁, k₂ in ZONE 1 (upper part of Table), and k₁, k₂, k₃ in ZONE 2 (lower part of Table). Parameter-values are from 3 time points, i.e. at the 1^st, 15^th and 30^th min of running. Parameter-means change with fatigue was significant with p>0.148. The asterisks in the lower part of Table signify significant changes (p<0.05) from ZONE 1 to ZONE 2.

The m₁ displacement obtained from our results revealed that the foot sinks about one cm in the first time zone of heel-strike, followed by a bouncing back phase. Previous reports also indicated a 1 cm deformation of the heel [e.g. 48]. The elastic properties of the heel-pad of various mammals were studied and it has been found that full up and down oscillations might result before actual settling down of the foot [53]. At cases, the foot may even temporarily lose contact with the ground during these oscillations. In our case, the oscillation phase would succeed the first time zone.

The obtained average knee stiffness k₃ of 15.8 kN/m was within the range of 6.89 – 112.0 kN/m, previously reported [12-14,16,19,22,29,45,54-55].

Further exploration of the effect of fatigue in the course of running was performed by correlating, using linear regression, the tibial tuberosity peak acceleration to the parameter values obtained from the model results. Figure 7 shows the correlation for which the Pearson’s coefficient was statistically different from zero (with 95% significance). The parameters are k₁ in the first and second zones, upper and middle rows, respectively, and k₃ in the second time zone of stance. The 1^st, 15^th and 30^th min of running are shown in the left, middle and right columns, respectively. Table 3 summarizes the regression results for the cases displayed in Figure 7. The results show the following: (a) For k₁ in time zone 1, the Pearson correlation r_p starts off by a low value of 0.26 in the 1^st min of running but increases in the 15^th (0.49) and the 30^th (0.44) min of running, indicating that a higher peak in the tibial tuberosity acceleration is associated with a lower k₂ in time zone 1; (b) there is a high correlation (0.86) between the peak tibial tuberosity acceleration and the k₁ parameter in time zone 2 at the initial stage of running (1^st min). This correlation, however, was lower in the 15^th min (0.27) and in the 30^th min (0.49) of running, thus suggesting that high peak acceleration at the tibial tuberosity is associated with a higher k₁ value in time zone 2; (c) a medium correlation (0.73) was found between k₃ and peak acceleration at the tibial tuberosity in the 1^st mi of running. But here too this correlation was decreased with the development of fatigue, suggesting that higher peak acceleration at the tibial tuberosity is associated with a high k₃ parameter value in time zone 2.

Figure 7.

Simple linear-regression made to express the relationship between the tibial tuberosity peak acceleration and parameter values (results of 10 subjects). The regression was performed for the 1^st min (leftmost panels), the 15^th min (middle panels) and for the 30^th min (rightmost panels). The correlation coefficient of the presented regression lines was significantly different from zero, at 95% significance level. The 3 upper panels - for k₁ in ZONE 1 (Z1); 3 middle panels for k₁ in ZONE 2 (Z2); and 3 lower panels - for k₃ in ZONE 2.

It has been shown that in running with shoes, the foot is restricted from bulging sideways, thus limiting the vertical deformation to an average of 5.5 mm, as opposed to 9 mm when running barefoot [48,56]. This explains the higher stiffness during the first time zone compared to bare foot running. It also explains the lower stiffness during the second time zone compared to bare foot running. It has also been shown by that better energy absorption and impact shock attenuation is associated with lower stiffness [51].

Table 3.

Linear regression results (Y= aX+b) of Figure 7, at 95% confidence intervals for coefficients a and b ; r_p = Pearson correlation coefficient ; and r² = coefficient of determination.

The correlation found between low stiffness in the first time zone to the high stiffness in the second time zone is obvious from the anatomy of the heel pad, which consists of nearly closed collagen cells filled with fatty cells [27,48]. The vertical orientations of these cells, together with the high viscosity of the fat tissue are the major factors responsible for the absorption of the impact energy at heel strike. Initially, the fat flows sideways and small loads result in high deformation (low stiffness). In the second time zone, after the heel pad has already considerably deformed, further increase in deformation provokes a high load, thus high stiffness. The effect of fatigue could be explained by means of the heating effect during the course of running. With nearly 80 heel-strikes per min, the whole running duration of 30 min results in some 2500 heel-strikes, each of which causing a rapid deformation of the heel pad and during which the fatty tissue frictions while squeezed out of the collagen cells. The accumulated heat due to friction reduces the viscosity and the vertical displacement is accelerated, causing a reduction in stiffness during the first zone of heel strike. In the second zone, however, the thinner remaining tissue together with the underlying bone evokes an increase in stiffness.

7. Conclusion

Stress fractures in long bones of the lower limbs are believed to originate from repetitive and/or excessive loading, such as may take place in long-distance running at a speed exceeding the anaerobic threshold. In the present study the average running distance per test was 6.30 km (30 min of running at the average speed of 12.6 km/h) in agreement with the definition of ‘long distance’ [57]. We have measured and analyzed the following: respiratory data to monitor global fatigue; and accelerometry, to provide quantitative information on loading of the major segments of the lower limb. While providing accelerometer boundary values for the model system, accelerometry is an advantageous method due to its being non-invasive.

We have addressed a major fatigue-related factor taking part in exposing the shank to stress fractures risk: the decline in end tidal carbon dioxide pressure, the latter expressing metabolic fatigue [31,58]. The mechanical consequence of fatigue in long-distance running is two-fold: enhanced impact acceleration due to global fatigue and muscle activity imbalance due to local fatigue before and during foot contact, resulting in the development of excessive shank-bone bending stresses and higher risk of stress injury [11].

While departing from the stiffness constancy concept, the model revealed that a correct and sufficient variability of the joint stiffness is a two-region piece-wise constant stiffness indicating that a higher order of non-linearity is not necessary. This result should be considered meaningful in those problems where the constant stiffness representation is not sufficient and in cases where the system’s representation has to be improved. Joint stiffness is dominated by muscular activation [59-60] and as the joints stiffen, they undergo smaller angular displacements during the ground-contact phase, also resulting in smaller excursion of the hip and higher leg stiffness. Thus, since stiffness is related to muscle activation, the piece-wise constant stiffness obtained solution also provides, through the obtained stiffness profiles, an insight into the patterns of the muscular activation in the legs’ joints.

The fact that the simple model of a piece-wise constant stiffness can predict major features of the running exercise makes it an effective tool for future designing of artificial legs and robots and also for the development of more accurate control strategies.

Acknowledgement

This Chapter is partly based on an MSc Thesis of second author DD, carried out in JM’s Biomechatronics Laboratory, Department of Biomedical Engineering, Technion – Israel Institute of Technology, under the joint supervision of JM and Prof. Yacov Ben-Haim.