Humans employ a high degree of redundancy in joint actuation, with different combinations of muscle and tendon action providing the same net joint torque. Both the resolution of these redundancies and the energetics of such systems depend on the dynamic properties of muscles and tendons, particularly their force-length relations. Current walking models that use stock parameters when simulating muscle-tendon dynamics tend to significantly overestimate metabolic consumption, perhaps because they do not adequately consider the role of elasticity. As an alternative, we posit that the muscle-tendon morphology of the human leg has evolved to maximize the metabolic efficiency of walking at self-selected speed. We use a data-driven approach to evaluate this hypothesis, utilizing kinematic, kinetic, electromyographic (EMG), and metabolic data taken from five participants walking at self-selected speed. The kinematic and kinetic data are used to estimate muscle-tendon lengths, muscle moment arms, and joint moments while the EMG data are used to estimate muscle activations. For each subject we perform an optimization using prescribed skeletal kinematics, varying the parameters that govern the force-length curve of each tendon as well as the strength and optimal fiber length of each muscle while seeking to simultaneously minimize metabolic cost and maximize agreement with the estimated joint moments. We find that the metabolic cost of transport (MCOT) values of our participants may be correctly matched (on average 0. 36±0. 02 predicted, 0. 35±0. 02 measured) with acceptable joint torque fidelity through application of a single constraint to the muscle metabolic budget. The associated optimal muscle-tendon parameter sets allow us to estimate the forces and states of individual muscles, resolving redundancies in joint actuation and lending insight into the potential roles and control objectives of the muscles of the leg throughout the gait cycle.
Human walking relies on a complex interplay of several physiological systems, with each exhibiting some degree of redundancy. The nervous system directs muscle contraction while receiving input from many different neural pathways. Muscles work together to produce motion, but different combinations of muscle action can produce the same net torque at a given joint. Tendons provide the interface between muscle and bone, but the energy transferred to the skeleton can come from either the active muscle or the compliance of the tendon. Understanding how humans resolve these redundancies has been a long-standing problem in the fields of neuroscience and biomechanics [1,2]. Knowledge of how the neuromuscular system allocates load during a given task would provide insight into the control objectives that govern its actions. Potential objectives (reviewed in [3]) include joint trajectory planning or minimization of metabolic energy consumption, active muscle volume, or muscle fatigue. However without an adequate understanding of the roles of each component of the system, such control hypotheses remain mere speculation. The roles of individual muscles and tendons in producing motion depend on the neural drive to the muscles and on the force generation properties of both the muscles and the tendons. Several modes of experimental observations provide glimpses of these elements. Electromyography (EMG) can be used to quantify the neural drive to individual muscles, revealing which muscles contribute to a given movement and giving some measure of intensity [4–6]. However it is limited by signal variability, measurement artifacts, and an inexact mapping to physiology and muscle force. Ultrasound probes have recently been used to image the individual motions of some distal leg muscles and tendons in vivo [7–13], but are practical only for small muscles and limited tasks. Motion capture can be combined with a knowledge of anatomy to infer the net movement of muscle-tendon units under much more general circumstances [14]; however breaking the resulting movement profiles into individual muscle and tendon contributions requires knowledge of often unavailable force generation parameters of muscles and tendons. These parameters are typically estimated through cadaver studies, but the scaling of the relevant quantities among different muscles and subjects (not to mention the differences with living specimens) is not well understood [14–16] and can have a significant impact on the resulting modeled dynamics [15,16]. Given the incomplete view afforded by current experimental measures, a unifying theoretical framework that combines the available data modes into a model of neuromuscular function is desirable. Two primary approaches have been taken to address this issue in walking: optimal control and optimal design. Human walking studies based on optimal control [17,18] model the morphology of leg muscle tendon units (MTUs) using literature-based estimates. They infer muscle activation through optimization, choosing control objectives such as metabolic energy minimization and/or motion tracking. They have been successful in predicting joint moments, joint torques, and ground reaction forces but often significantly overestimate the metabolic cost required for locomotion [17,19]. While it is clear that the neural control of the biological system is optimized in some way, it may be infeasible to determine the true objective function for this approach. Many different muscle activation combinations can produce similar muscle torque values, and several different unknown control objectives and neurological factors may contribute at once. Further the underlying uncertainty in the muscle-tendon morphology may result in the excess metabolic cost observed, as improper leveraging of tendon compliance would affect muscle force and state and therefore metabolic estimates. Human walking models based on optimal design utilize the efficiency gains that can be made through MTU parameter tuning. As Lichtwark and Wilson showed [20], experimentally observed muscle-tendon strains may be predicted by maximizing the efficiency of isolated MTUs. This result likely stems from the well-documented ability of tendon to enable muscle to operate economically [21,22]. Krishnaswamy and Herr [23] further explored the potential of optimal design, estimating the torque breakdown of the muscles spanning the ankle during the stance phase of walking at self-selected speed. This work used EMG signals to estimate muscle activation during walking and developed an optimization framework based on the assumption that the morphology of the muscle-tendon units spanning the ankle has evolved to minimize the metabolic cost required for walking at self-selected speed. Its results indicated that one solution set is able to match both human metabolics and kinetics, demonstrating an efficient load-sharing amongst the plantar flexors that qualitatively matched available experimental data. In this work we further the study of [23], modifying and extending it to permit investigation of the full leg. We collected kinematic, kinetic, electromyographic, and metabolic data from five subjects walking at self-selected speed and used them to perform optimizations with prescribed skeletal kinematics. We varied the parameters that determine where each muscle operates on its force-length curve as well as those that shape the force-length relation of each tendon, seeking to simultaneously minimize metabolic cost and maximize agreement with the observed joint moments. Each muscle in the model was modeled as Hill-type [19,24] and driven by activation estimates produced from EMG data. Overviews of the system model and procedure are shown in Figs 1 and 2, respectively. We found that the correct metabolic consumption can be matched with reasonable fidelity in the modeled joint torque through application of a single constraint on the per-muscle metabolic budget. The resulting optimal parameter sets were used to compute muscle force and state, with fascicle length profiles being compared to available experimental measurements. Our results are organized as follows. First we provide our estimated muscle activation profiles, discussing their quality and implications. Second we display the results of our dual objective optimization problem and summarize the methodology used for choosing one optimal solution for each subject. Third we evaluate the optimal solutions, producing estimates of energetic variables and muscle state. We compare these results with available experimental measures, finding quantitative agreement with metabolic data and qualitative agreement with muscle fascicle length data.
Muscle activation provides a scaling factor for the active force generation capability of a given muscle at a given time. As described in Methods, we applied a hybrid approach similar to that of [23] to estimate muscle activation from surface EMG measurements of five participants during walking. A Bayesian algorithm first proposed by Sanger [25] was tuned (as described in Methods) and used to perform a hidden state estimation that effectively determined the neural excitation of each muscle. This method was chosen over more conventional bandpass filtering methods [4–6] because the slightly delayed timing of the profiles it produced more easily allows the production of the observed joint torques. We elaborate on this point further in the Discussion section. The estimated neural excitation was then passed to a shaping filter [26,27] that represents muscle activation dynamics. The resulting average profiles are plotted in Fig 3. As can be seen from this plot, significant variations occurred (both within and among subjects) in the activation estimates obtained from the muscles spanning the hip. This lack of consistency was likely due to some combination of the relatively large depth of these muscles beneath the skin, motion artifacts, and the relative inaccessibility of the area. To prevent this from compromising the ensuing analysis, neural excitation profiles from the wire electrode experiments of [28] were used for the monoarticular muscles spanning the hip. These neural excitation profiles were passed to the activation dynamics from [26,27] for temporal consistency and subsequently used as model input along with the data-based activation estimates from all other muscles. Further details about this procedure are given in Methods and S1 Text. While there is no available ground truth to compare our activation estimates to, we note that the time dependence of the results (when normalized to percent gait cycle) were relatively invariant across trial and subject. The profiles follow the expected build up and decay time scales of muscle activation and, as can be seen in the following, allow the model to produce realistic kinetic and metabolic results. The estimated muscle activations a (t) as well as joint kinematics θjoint (t) derived from motion capture data were used to actuate the full leg model shown in Fig 1 according to the scheme shown in Fig 2. The leg muscle-tendon model M was specified by a set m → of morphological parameters that describe the force generation characteristics of each MTU as well as two parameters that enable passive force generation by the iliofemoral, ischiofemoral, and pubofemoral ligaments as well as other connective tissue at the hip. These ligaments are known to prevent hip overextension and here allow for the recovery of elastic energy in the joint [18,29]; specifically they reduce the load on the iliacus muscle around toe-off. The contributions of each MTU to m → were its maximum isometric force Fmax, an overall scaling factor for tendon slack length lsl and optimal muscle length lopt, its tendon reference strain λref, and its tendon shape factor Ksh. The lumped hip flexor ligament (HFL) acted as a simple rotary spring, parameterized by spring constant KHFL and engagement angle θ0, HFL. Each parameterization of the model generated kinetic (τmod (t) ) and metabolic (C) output costs: M m →, a (t), θ j o i n t (t) → τ m o d (t), C. (1) The parameter vector m → was varied using a stochastic dual objective optimization scheme [30] that simultaneously minimized metabolic cost and the difference between the computed joint moments of the data and those produced by the model. The bounds were specified as described in the Materials and Methods section; they were chosen wide enough that they were not approached by the chosen optimal parameter sets. The solution spaces for this optimization are shown in Fig 4. For each participant, the set of Pareto optimal solutions (i. e. the set of solutions where one would have to compromise on one objective to improve on the other) forms a rounded corner in the objective space. In the ideal case, this corner would be sharp and consist of one solution that optimizes both objectives. However this ideal does not typically occur in noisy (realistic) systems and does not here. To generate predictions based on our model, we chose one optimal parameter set along the Pareto Front for each subject. This was accomplished by evaluating the per-muscle metabolic consumption among all Pareto optimal solutions (Fig 5). Within the set of solutions, those where the metabolic cost was low and the kinetic fit was poor were seen to have uniformly low expenditure per muscle. Those that had the very best kinetic fits but relatively large metabolic costs were seen to be sinking large amounts of metabolic energy into a small number of muscles to produce incremental improvements in the kinetic fit. Such a phenomenon was likely enabled by noise in the data (particularly in the EMG signals) and the imperfect ability of our lumped, partial muscle set to match the force produced by the full set of the human body. The ramp up in metabolic energy seen in this subset of muscles is not physical as it would lead to either rapid fatigue or to the muscle being modeled as much larger than it actually is (since muscle mass scales with Fmax). Hence we excluded solutions that displayed this behavior, choosing as optimal the remaining Pareto optimal solution with best kinetic fit. Mathematically this was achieved by setting a cutoff on the fractional expenditure of the vastus, as the metabolic cost of this muscle group was the largest and increased significantly as kinetic fit improved. The chosen fractional cutoff allowed us to apply one criterion to every participant to match the experimental metabolic cost, as shown in Figs 4 and 5. Table 1 shows how our choice of optimal solution is able to quantitatively match the experimentally-observed metabolic cost of transport in four out of five subjects (and on average) while maintaining acceptable joint moment agreement. The lone subject where quantitative metabolic agreement was not reached displayed only a 6% error. Further details of the cutoff are included in Materials and Methods and in S3 Text while the optimal parameter set for each participant is given at the end of this document. To further evaluate the quality of our joint moment estimates we computed the fractional mean absolute error (FMAE) F M A E = 1 N 1 range (τ e x p) ∑ i = 1 N τ e x p, i - τ m o d, i (2) between the modeled joint moments τmod and the experimentally observed joint moments τobs. These quantities are computed only over the stance phase to facilitate comparison with those generated from another current EMG-driven analysis, [31]. The quantities quoted from [31] represent an average over four different activities (walking, running, side stepping, and crossover) but are generally close to those they generate for walking only (except for the hip, which had lower errors for walking). Their analysis included two treatments; one that considered only one degree of freedom and one that considered multiple degrees of freedom (as in our model). In general our moment fits compare favorably (Table 2). The metabolic expenditure of our model also aligns with previously published results. Across subjects, we found the average efficiency of positive muscle work to be 0. 26 ± 0. 02, consistent with [23,32,33]. We also estimated the metabolic expenditure of the model during different portions of the gait cycle by cross referencing the simulation with the input force plate data. The results are compiled in Table 3 and show a breakdown that is very similar to that simulated in [34]. The optimal muscle-tendon parameter sets estimated by our optimization procedure provide a means to resolve the redundancy in joint actuation for each subject. When applied in conjunction with the computed kinematics and estimated muscle activations, individual muscle force and state may be estimated. Fig 6 shows the torque breakdown for each joint in terms of percentage of body weight times height and averaged over all subjects. Fig 7 shows the trajectories of the muscle fascicle length normalized by lopt and averaged over all subjects. Note that the spread of some of these profiles is not due to a difference in general shape but rather to an overall offset in length, as evidenced by their nearly constant standard deviations. The soleus, hamstring, and vastus fascicle lengths in this plot all agree quite well in both shape and offset with the predictions in [35]. The soleus force and length also agree qualitatively with the projections of [11]. Fig 8 shows the velocity trajectories for each muscle fascicle normalized by its maximal value vmax and again averaged over all participants. As can be seen from the small variation in the shapes of the profiles in these plots, muscle force and state followed similar trajectories across participants. These predictions were also seen to vary little along the Pareto front in the experimentally measured metabolic band of a given subject. In general, muscle fascicle state is extremely difficult to measure. The only currently available means to obtain these profiles is ultrasonography, which is only practical for the relatively short distal muscles of the leg. In Fig 9 the modeled fascicle trajectories of muscles from one subject are compared with the experimental profiles available from published ultrasound studies. In each case the muscle fascicle lengths lm are normalized by their length at heel strike, lmHS. To generate these plots we took the modeled muscle from the subject who most closely matched the average height and weight of the experimental study. Soleus and gastrocnemius profiles came from Ishikawa et al [8], a gastrocnemius profile came from Fukunaga et al [7], and the vastus lateralis profile came from Chleboun [10]. In the plantar flexors, long stretches of nearly isometric operation are observed in mid-stance in both the model and in vivo profiles. In the vastus, the fascicle trajectory is seen to somewhat track the flexion of the knee in both the model and the published data. However while the qualitative trends of each muscle are consistent, quantitative agreement is not observed. We believe that observed differences come from (i) the difference in walking speed between our study and the literature, (ii) natural variation in the kinematics of early stance (which affects the initial muscle length for normalization), and (iii) uncertainty in the breakdown of what constitutes muscle and what constitutes tendon in ultrasound studies. We hope that future experimental methodologies, perhaps employing implantable sensors, will be able to further test our predicted fascicle trajectories.
Several observations may be made about both the methodology employed and the results obtained in this study. On the methodology side, we first address the steps taken to estimate the activations of the muscles in our model. As mentioned above, the conventional approach for estimating muscle excitation based on EMG data involves a bandpass filter [4–6]. We tried this approach with our optimization scheme (on all subjects) but found that it provided an inferior ability to match the observed joint moments compared to the implemented method. Digging further into this we found that two main differences existed between the approaches; (i) the chosen Bayesian method produces a profile that turns on and off more sharply than the signal produced by the bandpass filter and (ii) the signal produced by the Bayesian method consistently lags the bandpassed signal by about 50 ms. The sharpness of the profile does not affect the results significantly as it is mostly washed out by the ensuing activation dynamics and averaging. However the time dependence does matter; this lag enables the build up of muscle force in a manner consistent with the observed joint torques. Interestingly we found that [36] introduced a 40 ms lag to their EMG signal to “account for electromechanical delay between surface EMG and force production. ” This lag was employed by other studies [37,38] and fell within the 10−100 ms range given for these processes in the literature [39,40]. We found that adding a 40 ms lag as in [36] to the excitations produced by bandpassed filtered EMG signals gave performance nearly as good as those provided by the Sanger algorithm, with less variation among gait cycles. However this lagged bandpass method does not lend itself to a biophysical interpretation as clearly as the Bayesian model does. Despite our best efforts at EMG data collection and activation estimation, deficiencies in this part of our data set clearly exist. Surface EMG in general is prone to noise and artifacts, and Sanger’s algorithm (tested only isometrically in his publication [25]) does not remove them. While the effects of inconsistent artifacts was minimized by discarding the EMG data from clearly compromised gait cycles and the use of average trajectories in our model, they were likely not removed entirely. Better results may be obtained in future work through fine wire EMG measurements, which while more invasive are known to produce more reliable signals. Here generic fine wire EMG profiles reported in [28] were used to replace the noisy surface EMG measurements of the muscles spanning the hip. Interestingly our model actually displayed slightly better agreement with the observed hip torque profiles than with those of the other joints, but this is misleading as most of the modeled hip moment came from the hamstrings and the hip flexor ligament (which were unaffected by the generic profiles). A more minor deficiency in the EMG pipeline was normalization by the maximal voluntary contraction (MVC) values. While we do not know of a better alternative, this approach did lead to the normalized excitations of our muscles occasionally exceeding one in fast walking trials. When that occurred we renormalized by the value in the fast walking trial and reprocessed, but the normalization constant could still have been too small in other cases. Fortunately the impact of this scaling is extremely minimal because the estimated activation directly multiplies the maximum isometric force Fmax, which is optimized. A small effect remains because the normalization occurs before the excitation undergoes the activation dynamics (5), but that effect is largely irrelevant to our results. The system identification component of this study produced a methodology for estimating muscle-tendon parameters capable of matching the measured metabolic consumption while producing joint torque profiles that tracked observations reasonably well. One criterion based on the metabolic consumption of the vastus muscle group (the largest in the model) enabled the metabolic match for all participants. While the accuracies of our modeled joint torque profiles compare favorably with other current EMG driven modeling procedures [31], they do tend to underestimate the required joint torques. This characteristic is likely due to the exclusion of some muscles as well as the lumping of some muscle groups. In particular the deficiency exhibited in ankle moment during late stance is likely due to the exclusion of smaller plantar flexor muscles which are known to engage during that time [28]. Similarly the lumping of the three hamstring muscles (semimembranosus, semitendinosus, and biceps femoris long head) may be responsible for some of the deficiencies observed in the knee flexion moments throughout the gait cycle. Further resolution of these and other muscle groups where the muscle activations are not quite concurrent and the muscle-tendon lines of action are not quite aligned would likely improve the predictive power of this approach, allowing more accurate determination of the roles of each muscle during a given task. One notable aspect of the model was the importance of the hip flexor ligament, which produced nearly all of the required hip flexion moment near toe off at no metabolic cost. Its linear form was chosen as in [29] for maximal simplicity, but did not agree with the nonlinear damped form used in previous work [18,41]. Given that it produced torque for free it may have suppressed the required action of the other hip flexors (notably the iliacus) near toe off and produced a better hip moment fit than would have otherwise been possible for the same metabolic cost. However the suppression of other hip flexors could not have been large as these muscles are not strong enough to produce the required torque alone and are not significantly stretched in this time frame. It is also known that the hip ligaments produce the flexion torque necessary to balance the upper body against gravity in a standing position, where the line of gravity passes posterior to the hip joint [42]. They enable people to stand erect and even carry extra weight without significant muscle work at the hip, allowing a low metabolic cost to be maintained [43]. Since the engagement angles we used are consistent with standing and damping in human connective tissue is believed to yield only a slight drift over walking time scales [16], we believe that the contribution of the hip flexor ligament in our model is physiologically reasonable. The optimal muscle-tendon parameters found in this study play different roles in facilitating efficient locomotion. The maximal muscle isometric force (Fmax) must be large enough to meet the torque requirements at each joint, but small enough to keep muscle size and metabolic cost reasonable. The tendon slack lengths lsl and muscle optimal lengths lopt govern the timing of force production, acting in concert with muscle activation. The tendon shape parameters λref and Ksh define the elastic properties of the tendon and are tuned to ensure correct muscle operation and joint actuation. The resulting muscles and tendons together produce an interconnected system capable of producing the joint torque necessary for locomotion at a minimal metabolic cost. Evaluating the velocities at which the muscles in our model contract while activated can lend insight to the goals of their control. As was emphasized in [29,44], muscles minimize metabolic consumption at low speeds (i. e. when operating approximately isometrically). Further, A. V. Hill [45] demonstrated that skeletal muscle maximizes its efficiency while shortening at vCE ≈ −0. 17vmax and its power output at vCE ≈ −0. 30vmax. These three speeds (vCE = 0, −0. 17, −0. 30vmax) are indicated by the horizontal dashed lines in Fig 8. Combining this with Fig 3 allows us to contextualize the modeled velocities of each muscle when activated. At the ankle, the tibialis anterior is seen to operate at low speeds when engaged. This minimizes metabolic cost and is consistent with isometric contractions, as modeled in [29,44]. Both the soleus and gastrocnemius are seen to operate approximately isometrically through their activation in mid-stance before rapidly increasing their contraction velocity in a power stroke toward the end of stance (≈ 60% GC). The required positive work of the plantar flexors is consistent with the need seen in [29,44], but their rapid contraction at the end of stance does not strictly agree with the efficiency goal of the soleus and power goal of the gastrocnemius noted in [23]. This inconsistency amounts to a phase difference; the observation in [23] was based on muscle speed at toe off but our plantar flexor velocities reach approximately the same levels around 53% GC. At the knee, most muscles are seen to operate approximately isometrically. This agrees with [29,44], which note that these muscles primarily serve to modulate the stiffness of the joint in an optimally economical fashion. It is worth noting that the vastus group is seen to re-engage at the end of swing near the the optimal efficiency regime (extending the knee for heel strike) before returning to the low speed regime. At the hip, most muscles are observed to contract at low speed when activated. Two exceptions are the iliacus and the adductor longus, both of which gravitate toward maximal efficiency as they flex the hip around toe off. The adductor magnus may also approach the maximal efficiency regime as it extends the hip around heel strike, but this is less clear. Endo et al [29,44] found that the muscles spanning the hip could not be modeled strictly isometrically, consistent with these observations. Several avenues exist for extending this work in future studies. One route would be to expand the model past the sagittal plane; as shown in [31] matching moments in three dimensions could have an impact on the optimized parameter sets. Another route would be to validate the model by testing under different walking conditions. If more reliable EMG measurements could be obtained, we could train the model based on level-ground walking at self-selected speed and then evaluate the ability of the optimal parameter set to match the observed joint torque profiles and metabolic costs under different conditions (speeds, inclines, etc.) In this case new experimental observations of the change in plantar flexor function across speed [12,13] could be used to evaluate the model. Another possibility is that our overarching hypothesis is incomplete- maybe the body, instead of adapting its morphology to maximize its efficiency at (the presumably most common task of) walking at self-selected speed, actually optimizes with multiple tasks in mind. One could imagine additionally collecting data from a subject running a self-selected speed and adding extra dimensions to the cost function wherein that task is optimized concurrently with walking at self-selected speed. The procedure followed here for self-selected speed walking could thus become part of a larger optimization scheme by including other tasks. While such experiments could lend considerable insight, they will undoubtedly prove challenging due to the difficulties associated with collecting EMG under the stated conditions. One alternative would be to simulate the neural control of the model. Recent forward dynamic models driven by reflexive feedback have shown the ability to walk stably across different terrains [19,46], at different speeds [47,48], and in three dimensions [48,49]. These models are built for a subject with “average” dimensions and could therefore be improved through data-driven customization for individual subjects and the inclusion of more realistic muscle-tendon geometries and morphologies. This would include scaling segment inertias, refining muscle-tendon lines of action and moment arms [14], and optimization of muscle-tendon morphologies to simultaneously minimize metabolic cost and maximize agreement to experimentally observed kinematic and/or kinetic data. In this way a more realistic representation of an individual subject could be obtained, leading to further insights about the roles of individual muscles and tendons during gait and their variation amongst subjects. One could also imagine altering the neural control in this paradigm and using the framework to study movement disorders such as cerebral palsy.
The experiments of this study were conducted in compliance with the principles of the Declaration of Helsinki. The study was approved by the MIT Committee on the Use of Humans as Experimental Subjects (Protocol 1101004266). Prior to the experiments all participants provided written consent for data collection, analysis, and publication. Kinematic, kinetic, electromyographic, and metabolic data were collected at the Harvard University Skeletal Biology Lab. Five healthy adult males participated in the study, with their average height, mass, and age being 1. 77 ± 5 m, 70. 4 ± 6. 3 kg, and 26 ± 2 years, respectively. The required data sets were collected in two phases. First the subjects were outfitted with a portable oxygen consumption mask attached to a Cosmed K4B2 VO2 system. This system employs a standard open-circuit gas analysis technique to estimate metabolic energy consumption based on measurements of oxygen inspired and expired [50]. Three of the subjects were asked to stand still for seven minutes while a basal measurement was recorded, with the rates of the other two participants taken to be the mean of these three closely grouped rates (1. 38 ± 0. 06 W/kg, 1. 60 ± 0. 07 W/kg, 1. 63 ± 0. 06 W/kg). Participants then walked barefoot on an instrumented treadmill for seven minutes at each of six speeds (0. 75 m/s, 1. 00 m/s, 1. 25 m/s, 1. 50 m/s, 1. 75 m/s, and 2. 00 m/s), allowing the variation of metabolic energy expenditure to be measured across speed. The results were quickly tabulated and used to estimate the walking speed where the metabolic cost of transport (MCOT) was minimal. Once the metabolic cost measurements were completed, the oxygen consumption mask and Cosmed system were removed and each participant was outfitted for the second phase. In this phase kinematic, kinetic, and electromyographic data were collected for two minutes of barefoot walking at each of seven speeds; the six listed above and the speed where the subject’s MCOT was found to be minimal. An infrared camera system (8 cameras, Qualisys Motion Capture Systems, Gothenburg, Sweden) was used to track the motion of subjects as they walked in the capture volume. Reflective markers were placed at 43 (bilateral) locations on the participant’s body and their three dimensional trajectories were recorded at 500 Hz. The marker locations were chosen specifically to track joint motion, as prescribed by the Helen Hayes marker model. The ground reaction forces and contact centers of pressure were measured using a split-belt instrumented force plate treadmill (Bertec Corporation, Columbus, OH). Electromyographic signals were collected using a surface system from Motion Lab Systems (Baton Rouge, LA) and electrode placements as dictated in [51]. Fourteen muscles (tibialis anterior, soleus, medial gastrocnemius, vastus lateralis, biceps femoris shorthead, rectus femoris, semimembranosus, biceps femoris long head, illiacus, gluteus maximus (lower), gluteus maximus (upper), gluteus medius, adductor longus, and adductor magnus) on one leg of each subject were recorded, with symmetry being assumed for the other leg. The signals were recorded at the surface using pre-gelled bipolar electrodes (Electrode Store Model BS-24SAF, part number DDN-20), sampled at 1000 Hz, and amplified 20 times by pre-amplifiers (Motion Lab Systems, part number MA-411). Prior to walking trials a maximum voluntary contraction (MVC) trial was conducted for each muscle group wherein the participant was asked to work that particular group as hard as possible. These MVC trials were used for normalization purposes in muscle excitation estimates. The processed data were used to build a dynamic model of the leg during walking (Fig 1). The model includes all muscles that make significant contributions to the components of ankle, knee, and hip torque perpendicular to the sagittal plane during walking. Several muscle groups were lumped together for simplicity; this was deemed appropriate if all muscles within a group had similar lines of action and were activated simultaneously during walking. These included the GAS group (medial and lateral gastrocnemius), the VAS group (vastus lateralis, medialis, and intermedius), the HAM group (semimembranosus, semitendinosus, and biceps femoris long head), and the ILL group (iliacus and psoas). Each muscle group had one effective tendon to represent the net compliance of the muscle’s interaction with the skeleton. In addition to the muscle-tendon units a passive ligament spanning the hip was included to allow for the recovery of elastic energy in that joint. This element represents the contributions of the iliofemoral, ischiofemoral, and pubofemoral ligaments as well as those of other connective tissue at the hip and was modeled using a nonlinear equation in [18]; for simplicity we took it to be a linear spring that only engages around the time the hip goes vertical and begins to extend [29]: τ H F L = - K H F L (θ h i p - θ 0, H F L). (6) Note that both the moment and the angle in this equation were defined with flexion being positive. Below we describe the muscle, tendon, and joint dynamics of the system. An overview of the applied optimization procedure is shown in Fig 2. There were two categories of inputs for each muscle in our model: (i) a (t), lmtc (t), and ri (t) - all estimated from the data and (ii) muscle-tendon morphological parameters m → i which were identified via optimization. Morphological parameters which affect the Hill-type contraction dynamics of the modeled muscles include muscle maximum isometric force Fmax, length where active muscle force is maximal lopt, fiber composition FT, tendon slack length lsl, tendon shape factor Ksh, and tendon reference strain λref. We evaluated the consequences of varying all of these variables, finding differing levels of sensitivity for each. The variables to which the model was insensitive were fixed, leaving Fmax, Ksh, λref, and an overall scaling factor for lsl and lopt as the parameters to optimize. The last factor was chosen to ensure that the muscle operated in reasonable regions of its force-length space while preventing overfitting. The ratio lsl/lopt is known to vary significantly among muscles but not significantly in the same muscle among subjects [14], justifying this choice. As mentioned above, our optimization is built on the hypothesis that the muscle-tendon morphologies of the muscles comprising the leg have evolved to minimize the metabolic cost of walking at the speed where the MCOT is minimal. This assumption is supported by the fundamental importance of bipedal walking as a means of transport and the tendency of humans to walk at or near this speed. For each subject we therefore sought a set of morphological parameters m → i that match the measured joint torques at the ankle, knee, and hip while consuming a minimal amount of metabolic energy. This amounts to a dual objective optimization problem, with metabolic and kinetic cost functions being simultaneously minimized. The specific forms of these two costs are summarized below. As mentioned above, optimal muscle-tendon parameters were identified using a dual objective optimization. Fifty morphological parameters were optimized for each subject; four morphological parameters for each of the twelve modeled muscles and two additional parameters describing the hip flexor ligament. The chosen muscle-tendon parameter bounds are shown in Table 4. The maximum isometric force (Fmax) values for each muscle were constrained to fall in a window surrounding the scaled value from SIMM. The large width of this window helped to compensate for uncertainties in the EMG normalization (via MVC values, as discussed above). The scaling factor for lslack and lopt was chosen to ensure that the muscle fascicle length stayed within reasonable physiological operating ranges. The bounds for Ksh and λref were taken from [16]. The bounds for the spring constant of the hip flexor ligament was chosen so that the ligament could provide anywhere from none to all of the required hip flexion moment near toe off. The engagement angle was chosen so that the tendon could turn on with the hip no more than 10° flexed, as its physiological role is to prevent overextension. The model was constructed using MATLAB and Simulink and integrated using the ode15s (stiff/NDF) solver. Computations were parallelized and carried out using the Mathworks Cloud Center, a computer cluster operated through Amazon Web Services. This resource proved useful as our runs were configured to test many solutions and took about 10 hours to run on a single CPU. The optimization algorithm employed was MATLAB’s gamultiobj, a controlled elitist genetic algorithm that is a variant of the NSGA-II algorithm [30]. This method was chosen because of the likely presence of local minima in both objectives and tuned as shown in Table 5. Optimization tunings were chosen to ensure generations sufficiently large to force the optimizer to thoroughly search the space, eliminating the need for population seeding. Each trial solution was allowed to run for two gait cycles, with only the results of the second gait cycle being considered. This removed initial transient effects and allowed evaluation of the system over the full gait cycle.
Neuromuscular systems often employ redundancy in joint actuation, with different combinations of muscle and tendon action producing the same net joint torque. Both the resolution of these redundancies and the energetics of such systems depend strongly on the force-length relations of muscles and tendons. Many human walking models fail to properly account for elasticity by failing to scale muscle-tendon parameters for different individuals, relying instead on stock values taken from cadaver studies. This can result in inaccurate estimates of metabolic consumption as well as of the forces and states of individual muscles. Instead we estimate muscle-tendon parameters using a data-driven optimization procedure, testing the hypothesis that the human leg has evolved to maximize the metabolic efficiency of walking at self-selected speed. We find that the experimentally observed metabolic consumption can be matched with reasonable joint torque fidelity through the addition of a single constraint on the per-muscle metabolic budget. The associated muscle-tendon parameter sets were used to compute muscle force and state estimates, lending insight into potential roles and control objectives of the major muscles of the leg throughout the gait cycle.
medicine and health sciences
classical mechanics
biological locomotion
biomechanics
optimization
mathematics
pelvis
bioassays and physiological analysis
muscle electrophysiology
research and analysis methods
musculoskeletal mechanics
muscle physiology
musculoskeletal system
hip
joints (anatomy)
connective tissue
biological tissue
electrophysiological techniques
walking
torque
physics
anatomy
electromyography
tendons
physiology
biology and life sciences
physical sciences
motion
Human Leg Model Predicts Muscle Forces, States, and Energetics during Walking