Using computed muscle control to generate forward dynamic simulations of human walking from experimental data

The objective of this study was to develop an efficient methodology for generating muscle-actuated simulations of human walking that closely reproduce experimental measures of kinematics and ground reaction forces. We first introduce a residual elimination algorithm (REA) to compute pelvis and low back kinematic trajectories that ensure consistency between whole-body dynamics and measured ground reactions. We then use a computed muscle control (CMC) algorithm to vary muscle excitations to track experimental joint kinematics within a forward dynamic simulation. CMC explicitly accounts for delays in muscle force production resulting from activation and contraction dynamics while using a general static optimization framework to resolve muscle redundancy. CMC was used to compute muscle excitation patterns that drove a 21-degrees-of-freedom, 92 muscle model to track experimental gait data of 10 healthy young adults. Simulated joint kinematics closely tracked experimental quantities (mean root-mean-squared errors generally less than 1 degrees), and the time histories of muscle activations were similar to electromyographic recordings. A simulation of a half-cycle of gait could be generated using approximately 30 min of computer processing time. The speed and accuracy of REA and CMC make it practical to generate subject-specific simulations of gait.     

Sloped muscle excitation waveforms improve the accuracy of forward dynamic simulations

Mathematical models of the muscle excitation are useful in forward dynamic simulations of human movement tasks. One objective was to demonstrate that sloped as opposed to rectangular excitation waveforms improve the accuracy of forward dynamic simulations. A second objective was to demonstrate the differences in simulated muscle forces using sloped versus rectangular waveforms. To fulfill these objectives, surface EMG signals from the triceps brachii and elbow joint angle were recorded and the intersegmental moment of the elbow joint was computed from 14 subjects who performed two cyclic elbow extension experiments at 200 and 300 deg/s. Additionally, the surface EMG signals from the leg musculature, joint angles, and pedal forces were recorded and joint intersegmental moments were computed during a more complex pedaling task (90 rpm at 250 W). Using forward dynamic simulations, four optimizations were performed in which the experimental intersegmental moment was tracked for the elbow extension tasks and four optimizations were performed in which the experimental pedal angle, pedal forces, and joint intersegmental moments were tracked for the pedaling task. In these optimizations the three parameters (onset and offset time, and peak excitation) defining the sloped (triangular, quadratic, and Hanning) and rectangular excitation waveforms were varied to minimize the difference between the simulated and experimentally tracked quantities. For the elbow extension task, the intersegmental elbow moment root mean squared error, onset timing error, and offset timing error were less from simulations using a sloped excitation waveform compared to a rectangular excitation waveform (p<0.001). The average and peak muscle forces were from 7% to 16% larger and 20-28% larger, respectively, when using a rectangular excitation waveform. The tracking error for pedaling also decreased when using a sloped excitation waveform, with the quadratic waveform generating the smallest tracking errors for both tasks. These results support the use of sloped over rectangular excitation waveforms to establish greater confidence in the results of forward dynamic simulations.     

A general-purpose framework to simulate musculoskeletal system of human body: using a motion tracking approach*

Computation of muscle force patterns that produce specified movements of muscle-actuated dynamic models is an important and challenging problem. This problem is an undetermined one, and then a proper optimization is required to calculate muscle forces. The purpose of this paper is to develop a general model for calculating all muscle activation and force patterns in an arbitrary human body movement. For this aim, the equations of a multibody system forward dynamics, which is considered for skeletal system of the human body model, is derived using Lagrange–Euler formulation. Next, muscle contraction dynamics is added to this model and forward dynamics of an arbitrary musculoskeletal system is obtained. For optimization purpose, the obtained model is used in computed muscle control algorithm, and a closed-loop system for tracking desired motions is derived. Finally, a popular sport exercise, biceps curl, is simulated by using this algorithm and the validity of the obtained results is evaluated via EMG signals.     

Evaluation of performance criteria for simulation of submaximal steady-state cycling using a forward dynamic model.

The objectives of this study were twofold. The first was to develop a forward dynamic model of cycling and an optimization framework to simulate pedaling during submaximal steady-state cycling conditions. The second was to use the model and framework to identify the kinetic, kinematic, and muscle timing quantities that should be included in a performance criterion to reproduce natural pedaling mechanics best during these pedaling conditions. To make this identification, kinetic and kinematic data were collected from 6 subjects who pedaled at 90 rpm and 225 W. Intersegmental joint moments were computed using an inverse dynamics technique and the muscle excitation onset and offset were taken from electromyographic (EMG) data collected previously (Neptune et al., 1997). Average cycles and their standard deviations for the various quantities were used to describe normal pedaling mechanics. The model of the bicycle-rider system was driven by 15 muscle actuators per leg. The optimization framework determined both the timing and magnitude of the muscle excitations to simulate pedaling at 90 rpm and 225 W. Using the model and optimization framework, seven performance criteria were evaluated. The criterion that included all of the kinematic and kinetic quantities combined with the EMG timing was the most successful in replicating the experimental data. The close agreement between the simulation results and the experimentally collected kinetic, kinematic, and EMG data gives confidence in the model to investigate individual muscle coordination during submaximal steady-state pedaling conditions from a theoretical perspective, which to date has only been performed experimentally.     

Three-dimensional data-tracking dynamic optimization simulations of human locomotion generated by direct collocation

The aim of this study was to perform full-body three-dimensional (3D) dynamic optimization simulations of human locomotion by driving a neuromusculoskeletal model toward in vivo measurements of body-segmental kinematics and ground reaction forces. Gait data were recorded from 5 healthy participants who walked at their preferred speeds and ran at 2 m/s. Participant-specific data-tracking dynamic optimization solutions were generated for one stride cycle using direct collocation in tandem with an OpenSim-MATLAB interface. The body was represented as a 12-segment, 21-degree-of-freedom skeleton actuated by 66 muscle-tendon units. Foot-ground interaction was simulated using six contact spheres under each foot. The dynamic optimization problem was to find the set of muscle excitations needed to reproduce 3D measurements of body-segmental motions and ground reaction forces while minimizing the time integral of muscle activations squared. Direct collocation took on average 2.7 ± 1.0 h and 2.2 ± 1.6 h of CPU time, respectively, to solve the optimization problems for walking and running. Model-computed kinematics and foot-ground forces were in good agreement with corresponding experimental data while the calculated muscle excitation patterns were consistent with measured EMG activity. The results demonstrate the feasibility of implementing direct collocation on a detailed neuromusculoskeletal model with foot-ground contact to accurately and efficiently generate 3D data-tracking dynamic optimization simulations of human locomotion. The proposed method offers a viable tool for creating feasible initial guesses needed to perform predictive simulations of movement using dynamic optimization theory. The source code for implementing the model and computational algorithm may be downloaded at http://simtk.org/home/datatracking.     

Simulation analysis of muscle activity changes with altered body orientations during pedaling

Testing hypotheses related to the effect of gravitational orientation on neural control mechanisms is difficult for most locomotor tasks, like walking, because body orientation with respect to gravity affects both sensorimotor control and task mechanics. To examine the mechanical effect of body orientation independently from changes in workload and posture, Brown et al. (J. Biomech. 29 p. 1349, 1996) studied pedaling at altered body orientations. They found that subjects pedaling at different orientations changed needlessly their muscle excitations, putatively to preserve body-upright pedaling kinematics. We tested the feasibility of this hypothesis using simulations based on a three biomechanical-function pair organization for control of lower limb muscles (limb extension/flexion pair, extension/flexion transition pair, and foot plantarflexion/dorsiflexion pair), where each pair consists of alternating agonistic/antagonistic muscles. Adjustment of only three parameters, one to scale the muscle excitations of each pair, was sufficient to preserve pedaling kinematics to altered body orientation. Because these adjustments produced changes in muscle excitation and net joint moments similar to those observed in pedaling subjects, the hypothesis is supported. Moreover, the effectiveness of a decoupled gain adjustment procedure where each parameter was adjusted by error in only one aspect of the pedaling trajectory during each iteration (i.e., cadence adjusted the Ext/Flex parameter; peak-to-peak variation in crank velocity over the cycle adjusted the transition parameter; average ankle angle over the cycle adjusted the foot parameter) further supports the distinct function of each muscle pair.     

Fatigue effects on the coordinative pattern during cycling: Kinetics and kinematics evaluation

The aim of the present study was to analyze the net joint moment distribution, joint forces and kinematics during cycling to exhaustion. Right pedal forces and lower limb kinematics of ten cyclists were measured throughout a fatigue cycling test at 100% of PO_MAX. The absolute net joint moments, resultant force and kinematics were calculated for the hip, knee and ankle joint through inverse dynamics. The contribution of each joint to the total net joint moments was computed. Decreased pedaling cadence was observed followed by a decreased ankle moment contribution to the total joint moments in the end of the test. The total absolute joint moment, and the hip and knee moments has also increased with fatigue. Resultant force was increased, while kinematics has changed in the end of the test for hip, knee and ankle joints. Reduced ankle contribution to the total absolute joint moment combined with higher ankle force and changes in kinematics has indicated a different mechanical function for this joint. Kinetics and kinematics changes observed at hip and knee joint was expected due to their function as power sources. Kinematics changes would be explained as an attempt to overcome decreased contractile properties of muscles during fatigue.     

The influence of seat configuration on maximal average crank power during pedaling: a simulation study

Manipulating seat configuration (i.e., seat tube angle, seat height and pelvic orientation) alters the bicycle-rider geometry, which influences lower extremity muscle kinematics and ultimately muscle force and power generation during pedaling. Previous studies have sought to identify the optimal configuration, but isolating the effects of specific variables on rider performance from the confounding effect of rider adaptation makes such studies challenging. Of particular interest is the influence of seat tube angle on rider performance, as seat tube angle varies across riding disciplines (e.g., road racers vs. triathletes). The goals of the current study were to use muscle-actuated forward dynamics simulations of pedaling to 1) identify the overall optimal seat configuration that produces maximum crank power and 2) systematically vary seat tube angle to assess how it influences maximum crank power. The simulations showed that a seat height of 0.76 m (or 102% greater than trochanter height), seat tube angle of 85.1 deg, and pelvic orientation of 20.5 deg placed the major power-producing muscles on more favorable regions of the intrinsic force-length-velocity relationships to generate a maximum average crank power of 981 W. However, seat tube angle had little influence on crank power, with maximal values varying at most by 1% across a wide range of seat tube angles (65 to 110 deg). The similar power values across the wide range of seat tube angles were the result of nearly identical joint kinematics, which occurred using a similar optimal seat height and pelvic orientation while systematically shifting the pedal angle with increasing seat tube angles.     

A theoretical analysis of preferred pedaling rate selection in endurance cycling

One objective of this study was to investigate whether neuromuscular quantities were associated with preferred pedaling rate selection during submaximal steady-state cycling from a theoretical perspective using a musculoskeletal model with an optimal control analysis. Specific neuromuscular quantities of interest were the individual muscle activation, force, stress and endurance. To achieve this objective, a forward dynamic model of cycling and optimization framework were used to simulate pedaling at three different rates of 75, 90 and 105 rpm at 265 W. The pedaling simulations were produced by optimizing the individual muscle excitation timing and magnitude to reproduce experimentally collected data. The results from these pedaling simulations indicated that all neuromuscular quantities were minimized at 90 rpm when summed across muscles. In the context of endurance cycling, these results suggest that minimizing neuromuscular fatigue is an important mechanism in pedaling rate selection. A second objective was to determine whether any of these quantities could be used to predict the preferred pedaling rate. By using the quantities with the strongest quadratic trends as the performance criterion to be minimized in an optimal control analysis, these quantities were analyzed to assess whether they could be further minimized at 90 rpm and produce normal pedaling mechanics. The results showed that both the integrated muscle activation and average endurance summed across all muscles could be further minimized at 90 rpm indicating that these quantities cannot be used individually to predict preferred pedaling rates.     

A theoretical analysis of an optimal chainring shape to maximize crank power during isokinetic pedaling

Previous studies have sought to improve cycling performance by altering various aspects of the pedaling motion using novel crank–pedal mechanisms and non-circular chainrings. However, most designs have been based on empirical data and very few have provided significant improvements in cycling performance. The purpose of this study was to use a theoretical framework that included a detailed musculoskeletal model driven by individual muscle actuators, forward dynamic simulations and design optimization to determine if cycling performance (i.e., maximal power output) could be improved by optimizing the chainring shape to maximize average crank power during isokinetic pedaling conditions. The optimization identified a consistent non-circular chainring shape at pedaling rates of 60, 90 and 120 rpm with an average eccentricity of 1.29 that increased crank power by an average of 2.9% compared to a conventional circular chainring. The increase in average crank power was the result of the optimal chainrings slowing down the crank velocity during the downstroke (power phase) to allow muscles to generate power longer and produce more external work. The data also showed that chainrings with higher eccentricity increased negative muscle work following the power phase due to muscle activation–deactivation dynamics. Thus, the chainring shape that maximized average crank power balanced these competing demands by providing enough eccentricity to increase the external work generated by muscles during the power phase while minimizing negative work during the subsequent recovery phase.     

Evaluation of a method to scale muscle strength for gait simulations of children with cerebral palsy

Cerebral palsy (CP) is a neurological disorder that results in life-long mobility impairments. Musculoskeletal models used to investigate mobility deficits for children with CP often lack subject-specific characteristics such as altered muscle strength, despite a high prevalence of muscle weakness in this population. We hypothesized that incorporating subject-specific strength scaling within musculoskeletal models of children with CP would improve accuracy of muscle excitation predictions in walking simulations. Ten children (13.5 ± 3.3 years; GMFCS level II) with spastic CP participated in a gait analysis session where lower-limb kinematics, ground reaction forces, and bilateral electromyography (EMG) of five lower-limb muscles were collected. Isometric strength was measured for each child using handheld dynamometry. Three musculoskeletal models were generated for each child including a ‘Default’ model with the generic musculoskeletal model’s muscle strength, a ‘Uniform’ model with muscle strength scaled allometrically, and a ‘Custom’ model with muscle strength scaled based on handheld dynamometry strength measures. Muscle-driven gait simulations were generated using each model for each child. Simulation accuracy was evaluated by comparing predicted muscle excitations and measured EMG signals, both in the duration of muscle activity and the root-mean-square difference (RMSD) between signals. Improved agreement with EMG were found in both the ‘Custom’ and ‘Uniform’ models compared to the ‘Default’ model indicated by improvement in RMSD summed across all muscles, as well as RMSD and duration of activity for individual muscles. Incorporating strength scaling into musculoskeletal models can improve the accuracy of walking simulations for children with CP.     

Static and dynamic optimization solutions for gait are practically equivalent

The proposition that dynamic optimization provides better estimates of muscle forces during gait than static optimization is examined by comparing a dynamic solution with two static solutions. A 23-degree-of-freedom musculoskeletal model actuated by 54 Hill-type musculotendon units was used to simulate one cycle of normal gait. The dynamic problem was to find the muscle excitations which minimized metabolic energy per unit distance traveled, and which produced a repeatable gait cycle. In the dynamic problem, activation dynamics was described by a first-order differential equation. The joint moments predicted by the dynamic solution were used as input to the static problems. In each static problem, the problem was to find the muscle activations which minimized the sum of muscle activations squared, and which generated the joint moments input from the dynamic solution. In the first static problem, muscles were treated as ideal force generators; in the second, they were constrained by their force-length-velocity properties; and in both, activation dynamics was neglected. In terms of predicted muscle forces and joint contact forces, the dynamic and static solutions were remarkably similar. Also, activation dynamics and the force-length-velocity properties of muscle had little influence on the static solutions. Thus, for normal gait, if one can accurately solve the inverse dynamics problem and if one seeks only to estimate muscle forces, the use of dynamic optimization rather than static optimization is currently not justified. Scenarios in which the use of dynamic optimization is justified are suggested.     

Estimation of accuracy of patient-specific musculoskeletal modelling: case study on a post polio residual paralysis subject

For patients with patterns ranging out of anthropometric standard values, patient-specific musculoskeletal modelling becomes crucial for clinical diagnosis and follow-up. However, patient-specific modelling using imaging techniques and motion capture systems is mainly subject to experimental errors. The aim of this study was to quantify these experimental errors when performing a patient-specific musculoskeletal model. CT scan data were used to personalise the geometrical model and its inertial properties for a post polio residual paralysis subject. After having performed a gait-based experimental protocol, kinematics data were measured using a VICON motion capture system with six infrared cameras. The musculoskeletal model was computed using a direct/inverse algorithm (LifeMod software). A first source of errors was identified in the segmentation procedure in relation to the calculation of personalised inertial parameters. The second source of errors was subject related, as it depended on the reproducibility of performing the same type of gait. The impact of kinematics, kinetics and muscle forces resulting from the musculoskeletal modelling was quantified using relative errors and the absolute root mean square error. Concerning the segmentation procedure, we found that the kinematics results were not sensitive to the errors (relative error < 1%). However, a strong influence was noted on the kinetics results (deviation up to 71%). Furthermore, the reproducibility error showed a significant influence (relative mean error varying from 5 to 30%). The present paper demonstrates that in patient-specific musculoskeletal modelling variations due to experimental errors derived from imaging techniques and motion capture need to be both identified and quantified. Therefore, the paper can be used as a guideline.     

Estimation of accuracy of patient-specific musculoskeletal modelling: case study on a post polio residual paralysis subject

For patients with patterns ranging out of anthropometric standard values, patient-specific musculoskeletal modelling becomes crucial for clinical diagnosis and follow-up. However, patient-specific modelling using imaging techniques and motion capture systems is mainly subject to experimental errors. The aim of this study was to quantify these experimental errors when performing a patient-specific musculoskeletal model. CT scan data were used to personalise the geometrical model and its inertial properties for a post polio residual paralysis subject. After having performed a gait-based experimental protocol, kinematics data were measured using a VICON motion capture system with six infrared cameras. The musculoskeletal model was computed using a direct/inverse algorithm (LifeMod software). A first source of errors was identified in the segmentation procedure in relation to the calculation of personalised inertial parameters. The second source of errors was subject related, as it depended on the reproducibility of performing the same type of gait. The impact of kinematics, kinetics and muscle forces resulting from the musculoskeletal modelling was quantified using relative errors and the absolute root mean square error. Concerning the segmentation procedure, we found that the kinematics results were not sensitive to the errors (relative error<1%). However, a strong influence was noted on the kinetics results (deviation up to 71%). Furthermore, the reproducibility error showed a significant influence (relative mean error varying from 5 to 30%). The present paper demonstrates that in patient-specific musculoskeletal modelling variations due to experimental errors derived from imaging techniques and motion capture need to be both identified and quantified. Therefore, the paper can be used as a guideline.     

The effect of perturbing body segment parameters on calculated joint moments and muscle forces during gait

This study examined the effect of body segment parameter (BSP) perturbations on joint moments calculated using an inverse dynamics procedure and muscle forces calculated using computed muscle control (CMC) during gait. BSP (i.e. segment mass, center of mass location (com) and inertia tensor) of the left thigh, shank and foot of a scaled musculoskeletal model were perturbed. These perturbations started from their nominal value and were adjusted to ±40% in steps of 10%, for both individual as well as combined perturbations in BSP. For all perturbations, an inverse dynamics procedure calculated the ankle, knee and hip moments based on an identical inverse kinematics solution. Furthermore, the effect of applying a residual reduction algorithm (RRA) was investigated. Muscle excitations and resulting muscle forces were calculated using CMC. The results show only a limited effect of an individual parameter perturbation on the calculated moments, where the largest effect is found when perturbing the shank com (MS_com,shank, the ratio of absolute difference in torque and relative parameter perturbation, is maximally −7.81 N m for hip flexion moment). The additional influence of perturbing two parameters simultaneously is small (MS_{mass+com,thigh} is maximally 15.2 N m for hip flexion moment). RRA made small changes to the model to increase the dynamic consistency of the simulation (after RRA MS_com,shank is maximally 5.01 N m). CMC results show large differences in muscle forces when BSP are perturbed. These result from the underlying forward integration of the dynamic equations.     

Determining residual reduction algorithm kinematic tracking weights for a sidestep cut via numerical optimization

Musculoskeletal modeling allows for the determination of various parameters during dynamic maneuvers by using in vivo kinematic and ground reaction force (GRF) data as inputs. Differences between experimental and model marker data and inconsistencies in the GRFs applied to these musculoskeletal models may not produce accurate simulations. Therefore, residual forces and moments are applied to these models in order to reduce these differences. Numerical optimization techniques can be used to determine optimal tracking weights of each degree of freedom of a musculoskeletal model in order to reduce differences between the experimental and model marker data as well as residual forces and moments. In this study, the particle swarm optimization (PSO) and simplex simulated annealing (SIMPSA) algorithms were used to determine optimal tracking weights for the simulation of a sidestep cut. The PSO and SIMPSA algorithms were able to produce model kinematics that were within 1.4° of experimental kinematics with residual forces and moments of less than 10 N and 18 Nm, respectively. The PSO algorithm was able to replicate the experimental kinematic data more closely and produce more dynamically consistent kinematic data for a sidestep cut compared to the SIMPSA algorithm. Future studies should use external optimization routines to determine dynamically consistent kinematic data and report the differences between experimental and model data for these musculoskeletal simulations.     

Dynamic optimization of human walking

A three-dimensional, neuromusculoskeletal model of the body was combined with dynamic optimization theory to simulate normal walking on level ground. The body was modeled as a 23 degree-of-freedom mechanical linkage, actuated by 54 muscles. The dynamic optimization problem was to calculate the muscle excitation histories, muscle forces, and limb motions subject to minimum metabolic energy expenditure per unit distance traveled. Muscle metabolic energy was calculated by slimming five terms: the basal or resting heat, activation heat, maintenance heat, shortening heat, and the mechanical work done by all the muscles in the model. The gait cycle was assumed to be symmetric; that is, the muscle excitations for the right and left legs and the initial and terminal states in the model were assumed to be equal. Importantly, a tracking problem was not solved. Rather only a set of terminal constraints was placed on the states of the model to enforce repeatability of the gait cycle. Quantitative comparisons of the model predictions with patterns of body-segmental displacements, ground-reaction forces, and muscle activations obtained from experiment show that the simulation reproduces the salient features of normal gait. The simulation results suggest that minimum metabolic energy per unit distance traveled is a valid measure of walking performance.     

Diurnal variations in cycling kinematics

Physiological and biomechanical constraints as well as their fluctuations throughout the day must be considered when studying determinant factors in the preferred pedaling rate of elite cyclists. The aim of this study was to monitor the diurnal variation of spontaneous pedaling rate and movement kinematics over the crank cycle. Twelve male competitive cyclists performed a submaximal exercise on a cycle ergometer for 15 min at 50% of their W(max). Two test sessions were performed at 06:00 and 18:00 h on two separate days to assess diurnal variation in the study variables. For each test session, the exercise bout was divided into three equivalent 5-min periods during which subjects were requested to use different pedal rates (spontaneous cadence, 70 and 90 rev min(-1)). Pedal rate and kinematics data (instantaneous pedal velocity and angle of the ankle) were collected. The results show a higher spontaneous pedal rate in the late afternoon than in the early morning (p < 0.001). For a given pedal rate condition, there was a less variation in pedal velocity during a crank cycle in the morning than in the late afternoon. Moreover, diurnal variations were observed in ankle mobility across the crank cycle, the mean plantar flexion observed throughout the crank cycle being greater in the 18:00 h test session (p < 0.001). These results suggest that muscular activation patterns during a cyclical movement could be under the influence of circadian fluctuations.