Modeling and simulating the neuromuscular mechanisms regulating ankle and knee joint stiffness during human locomotion

Abstract

This work presents an electrophysiologically and dynamically consistent musculoskeletal model to predict stiffness in the human ankle and knee joints as derived from the joints constituent biological tissues (i.e., the spanning musculotendon units). The modeling method we propose uses electromyography (EMG) recordings from 13 muscle groups to drive forward dynamic simulations of the human leg in five healthy subjects during overground walking and running. The EMG-driven musculoskeletal model estimates musculotendon and resulting joint stiffness that is consistent with experimental EMG data as well as with the experimental joint moments. This provides a framework that allows for the first time observing 1) the elastic interplay between the knee and ankle joints, 2) the individual muscle contribution to joint stiffness, and 3) the underlying co-contraction strategies. It provides a theoretical description of how stiffness modulates as a function of muscle activation, fiber contraction, and interacting tendon dynamics. Furthermore, it describes how this differs from currently available stiffness definitions, including quasi-stiffness and short-range stiffness. This work offers a theoretical and computational basis for describing and investigating the neuromuscular mechanisms underlying human locomotion.

Keywords: compliance; electromyography; human leg; neuromusculoskeletal modeling; stiffness.

Copyright © 2015 the American Physiological Society.

Figures

Fig. 1.

The schematic structure of the electromyography (EMG)-driven musculoskeletal model. It comprises five components: musculotendon activation, A; musculotendon kinematics B; musculotendon dynamics, C; joint dynamics, D; and model calibration, E. The EMG-driven model is initially calibrated (E) using experimental EMG-excitations, experimental three-dimensional (3-D) joint angles, and matching experimental joint moments from a set of calibration trials. After calibration, the model is operated in open-loop (i.e., it does not require experimental joint moments to track). The calibrated model is validated on its ability to blindly match experimental joint moments from a set of novel trials that were not used for calibration. The musculotendon activation component (A) maps initial EMG-excitations recorded from 13 muscle groups to activations for 18 musculotendon units (MTUs, Table 1). Subsequently, MTU force and stiffness are determined (C) as a function of muscle activation (A) and experimental 3-D musculotendon kinematics derived from experimental 3-D joint kinematics (B). The resulting MTU forces and stiffness are simultaneously projected (D) on two degrees of freedom including the knee flexion-extension and the ankle plantar-dorsi flexion.

Fig. 2.

A: nominal static properties of muscle tissues. These are determined by the active and passive isometric force-length curves. Fiber force and length values are normalized by the maximal isometric force, Fmax, and optimal fiber length, lom, respectively. B: nominal dynamic properties of muscle tissues. These are determined by the active force-velocity curve. Fiber force and contraction velocity values are normalized by the Fmax and maximal contraction velocity, vmax, respectively. Both active force-length and force-velocity curves are depicted for different levels of muscle activation. C: solution space of the fiber force. This results from combining passive and active force-length and force-velocity curves. Note how muscle force is modulated simultaneously by activation, fiber length, and contraction velocity.

Fig. 3.

A: Profile of the fiber force-length of the soleus muscle over the stance phase of walking (i.e., see continuous colored workloop) with 0% being heel-strike and 100% toe-off events. Two cross-sections of the force-length-velocity-activation solution space along the fiber length axis are depicted in the correspondence of 57.6% stance (red dot) and 72% (black dot) respectively (i.e., see dotted lines). In these two stance points, the respective cross-sections differ with that of each other because they underlie different muscle activations, fiber lengths, and contraction velocities. Furthermore, note how in a specific stance point, the slope derived from the fiber force-length profile differs substantially from that derived from the force-length-velocity-activation cross section. B: force solution surfaces of the soleus muscles in the correspondence of 57.6% stance (red 3-D surface and red dot) and 72% (black 3-D surface and black dot). The fiber force-length profile of soleus muscles (i.e., green work loop) throughout the stance phase of walking (as in A) is also projected on the muscle force solution spaces. C: instantaneous stiffness of soleus (Km, Equation 4) throughout the stance phase of walking. This is predicted from the slope of the fiber force-length-velocity-activation surface in B as the slope along the normalized fiber length directional axis in the correspondence of the instant normalized fiber length predicted by our EMG-driven method (also see materials and methods). D: instantaneous quasi-stiffness (QS, Equation 6) of soleus throughout the stance phase of walking. This is predicted from the tangent of the fiber force-length profile in A. Note how stiffness and QS can predict substantially different slopes in the correspondence of the same stance point. QS relates changes in fiber force with changes in fiber length without accounting for the fact that force may be instead varying because of the activation or contraction velocity components, hence the peaks in (B) relative to the steep tangents in (A). The reported data are relative to subject 1 and are calculated as a function of normalized fiber length and contraction velocity values. Also see Supplemental Video S1.

Fig. 4.

A: stiffness predicted by our proposed method (Equation 4) for fibers of the gastrocnemius medialis, tendon, and musculotendon unit. B: short-range stiffness (Equation 7) predicted for fibers of the gastrocnemius medialis, as well as tendon stiffness and resulting musculotendon unit stiffness. Values are averaged across all walking trials performed by all subjects and are reported throughout the stance phase of walking with 0% being heel-strike and 100% being toe-off events.

Fig. 5.

Predicted and experimental joint moments from all validation trials and subjects. The ensemble average (continuous lines) curves are depicted for experimental and predicted joint moments. The standard deviations (shaded area) are also depicted for the experimental joint moments. Experimental and predicted joint moments are reported about two degrees of freedom including knee flexion-extension (KneeFE) and ankle plantar-dorsi flexion (AnkleFE). The reported data are from the stance phase with 0% being heel-strike and 100% being toe-off events.

Fig. 6.

Distribution of the similarity indexes computed between experimental and predicted joint moments from all validation trials and subjects. Similarity indexes include the coefficient of determination (R2) and the root mean square error normalized to the root mean square sum of the experimental joint moments (RMSE). Histograms gather the R2 and RMSE values in intervals with a 0.1-fixed-width in the (zero-to-1) range. In this range, vertical lines highlight a number of variables including the average indexes (i.e., mean) and those within one standard deviation (i.e., SD), the maximum expected RMSE and the minimum expected R2 calculated using the Chebyshev's theorem assuming a 90% confidence interval (i.e., Chebyshev 90%), as well as the most unfavorable index values including R2 = 0.58 (i.e., min R2) and RMSE = 0.704 (i.e., max RMSE).

Fig. 7.

The ensemble average curves (continuous line) and standard deviation (shaded area) for the stiffness predicted about two degrees of freedom including ankle plantar-dorsi flexion (AnkleFE) and knee flexion-extension (KneeFE) during walking across all validation trials and subjects. Data are reported for the stance phase of walking with 0% being heel-strike and 100% being toe-off events.

Fig. 8.

The ensemble average curves (continuous line) and standard deviation (shaded area) for the stiffness predicted about two degrees of freedom including ankle plantar-dorsi flexion (AnkleFE) and knee flexion-extension (KneeFE) during running across all validation trials and subjects. Data are reported for the stance phase of running with 0% being heel-strike and 100% being toe-off events.

Fig. 9.

Variability accounted for (VAF) by all 18 musculotendon units in the model in the net joint stiffness about two degrees of freedom including ankle plantar-dorsi flexion (AnkleFE) and knee flexion-extension (KneeFE). VAF values are averaged from all performed trials and subjects (black bar) and are reported together with the standard deviation (vertical line).

Fig. 10.

Ensemble average joint stiffness (shaded area) and co-contraction ratio (solid line). Data are averaged across all trials performed by all subjects across walking and running trials, respectively. Standard deviation curves are also reported relative to the co-contraction ratio (dotted lines). Data are reported about the ankle plantar-dorsi flexion (AnkleFE) and knee flexion-extension (KneeFE) degrees of freedom and throughout the stance phase with 0% being heel-strike and 100% being toe-off events.

Fig. 11.

Hierarchical structure of how net joint stiffness (A) emerges, in our framework, from the stiffness contribution of the joint constituent muscle-tendon units (MTUs) (B and C) as well as how MTU stiffness (D) results from the stiffness of constituent fibers in series with tendon stiffness (E). Data are averaged across all trials performed by all subjects across walking trials. Data are reported about the knee flexion-extension (KneeFE) degrees of freedom, for the associated MTUs with names as defined in Table 1, and throughout the stance phase with 0% being heel-strike and 100% being toe-off events. Results show that the early stance of walking (i.e., 0% ≤ stance ≤ 25%) underlies a progressive decrease in joint stiffness of knee flexors (B) and a simultaneous increase in net joint stiffness (A) mainly contributed by knee extensors (B). Throughout the midstance (i.e., 25% ≤ stance ≤ 85%), net knee stiffness (A) decreased due to a gradual decrease in joint stiffness of knee extensors (B). In the terminal stance (i.e., stance ≥85%), net knee stiffness (A) reached its lowest value in the correspondence of activity of knee flexors in preparation for the swing phase.