Muscle forces can be selected from a space of muscle recruitment strategies that produce stable motion and variable muscle and joint forces. However, current optimization methods provide only a single muscle recruitment strategy. We modelled the spectrum of muscle recruitment strategies while walking. The equilibrium equations at the joints, muscle constraints, static optimization solutions and 15-channel electromyography (EMG) recordings for seven walking cycles were taken from earlier studies. The spectrum of muscle forces was calculated using Bayesian statistics and Markov chain Monte Carlo (MCMC) methods, whereas EMG-driven muscle forces were calculated using EMG-driven modelling. We calculated the differences between the spectrum and EMG-driven muscle force for 1–15 input EMGs, and we identified the muscle strategy that best matched the recorded EMG pattern. The best-fit strategy, static optimization solution and EMG-driven force data were compared using correlation analysis. Possible and plausible muscle forces were defined as within physiological boundaries and within EMG boundaries. Possible muscle and joint forces were calculated by constraining the muscle forces between zero and the peak muscle force. Plausible muscle forces were constrained within six selected EMG boundaries. The spectrum to EMG-driven force difference increased from 40 to 108 N for 1–15 EMG inputs. The best-fit muscle strategy better described the EMG-driven pattern (R2 = 0.94; RMSE = 19 N) than the static optimization solution (R2 = 0.38; RMSE = 61 N). Possible forces for 27 of 34 muscles varied between zero and the peak muscle force, inducing a peak hip force of 11.3 body-weights. Plausible muscle forces closely matched the selected EMG patterns; no effect of the EMG constraint was observed on the remaining muscle force ranges. The model can be used to study alternative muscle recruitment strategies in both physiological and pathophysiological neuromotor conditions.
Internal forces that physical activity engender on our skeleton through muscles and joints are important for studying human motion  and skeletal mechanics . However, the biomechanical assessment of muscle forces is difficult, because the musculoskeletal system is highly redundant  and the sensorimotor control system is intrinsically variable . Muscle recruitment targets multiple and competing goals, and depends on the task being executed, subjective healthy condition and noise that plagues the sensory inputs and muscle outputs in determining the appropriate motor command . A better understanding of the repertoire of alternative sensorimotor control strategies may be important in studying human motion and skeletal mechanics .
According to the uncontrolled manifold hypothesis, our central nervous system (CNS) uses all the redundant degrees of freedom to ensure flexible and stable motion . Possible muscle synergies can therefore be defined as organizations of muscle forces that stabilize joint torques and motion; or, in other words, alternative solutions to the muscle load-sharing problem. Körding & Wolpert  showed that our CNS probably interprets the problem of optimal performance in a statistical fashion by weighting knowledge gathered from previous experiences and information gathered from multiple sensory modalities. By considering both types of information in the form of prior and likelihood, Bayesian statistics have been shown to properly describe the mechanism behind the generation of movement trajectories , forces  and judgement timing . Likewise, our CNS may solve the muscle load-sharing problem  by recruiting muscles from a space of alternative solutions, ensuring stable motion . However, the large majority of current methods used for calculating muscle forces target, among the infinite possible solutions, the single muscle synergy that minimizes a chosen cost function .
By using different energy- and stress-based cost functions, static optimization methods have been shown to provide muscle force patterns that are in qualitative agreement with recorded electromyography (EMG) . However, static optimization methods are known to underestimate the contribution of balanced agonist–antagonist muscle contractions , the so-called muscle co-contractions, which are essential in a number of circumstances. For example, muscle co-contractions are important (i) in controlling the joint impedance and stability during daily activities  and (ii) for executing motions characterized by rapid changes of joint torque, such as those occurring while landing  and running . Muscle co-contractions have also been found to determine the hip force during activity in terms of magnitude, distribution and timing [13,16]. In cats, the static optimization solution has been shown to be a poor predictor of the soleus and gastrocnemius force pattern .
EMG-driven methods have been developed to calculate muscle force patterns that follow the muscle electrical activity recorded using EMG. EMG-driven muscle forces are calculated by inputting the EMG signal to muscle excitation- and contraction-dynamic models, which are then used to solve the dynamic problem of the motion [18–20]. However, model simplifications and measurement errors cause inconsistencies between the motion being studied and the calculated EMG-driven muscle forces. Inverse EMG-driven models solve the muscle load-sharing problem by forcing the static optimization solution within an arbitrarily defined interval around the calculated EMG-driven muscle force to ensure that a solution to the problem exists . Forward EMG-driven models use optimization-based procedures to tune the model, so that the calculated EMG-driven muscle forces generate the desired motion of the model [18,20]. Lloyd & Besier  used an EMG-driven model of the knee-spanning muscles to estimate the knee torque. Sartori et al.  used a lower-limb model to calculate muscle forces during walking, running, sidestepping and crossover cutting manoeuvres. However, the single representative muscle synergy calculated using both EMG-driven and static optimization methods cannot provide information about the spectrum of muscle synergies driving motion [18–20].
Another possibility is to explore the entire solution space of the muscle load-sharing problem . The space of physiologically possible muscle synergies can be described as the solution space of the inverse, highly indeterminate, linear problem of muscle equilibrium at the joints, which can be geometrically represented by a bounded portion of a hyper-plane in the muscle force domain ; its orientation, offset and boundaries are respectively determined by the muscle lever arms, joint torques and physiological constraints of muscle force. Heino et al.  combined Bayesian statistics and Markov chain Monte Carlo (MCMC) methods for exploring the solution space of highly indeterminate inverse linear problems using software called Metabolica, which uses Bayesian statistics to estimate the posterior probability density function (PDF) for the unknowns and the MCMC algorithm to sample the estimated PDF. By constraining muscle forces between zero and the muscle peak force, this approach has been used to explore the muscle potential to generate force during a single walking frame . No study has investigated the natural unpredictable variability of muscle forces during activity. Combining musculoskeletal models, Bayesian statistics, MCMC sampling methods and EMG-driven muscle force modelling is a viable solution for calculating the spectrum of either potential or physiologically plausible musculoskeletal forces during activity.
The aim of this study was to investigate the repertoire of muscle synergies during walking. The lower-limb joint torques, muscle lever arms, muscle force constraints and the EMG signals for seven walking trials were taken from earlier studies [23,24]. The spectrum of muscle recruitment strategies was calculated using Bayesian statistics and MCMC sampling methods, whereas EMG-driven muscle forces were calculated using the available EMGs and a Hill-type muscle excitation and contraction model. The model was characterized by studying the consistency between EMG-driven forces and the model of the motion, and comparing the proposed method with a commonly used static optimization procedure. Finally, the muscle potential to generate force and the spectrum of physiologically plausible muscle forces were calculated and analysed.
2. Material and methods
2.1. Model development
The model, developed for studying possible muscle recruitment strategies during motion, was based on an earlier lower-limb model of a complete stride . The complete motion capture is available for download at www.physiomespace.com (keyword: LHDL_1stMatchedVolunteer_MOCAP). EMGs were recorded for 15 lower-limb muscles (i.e. gluteus maximus, gluteus medius, rectus femoris, vastus lateralis, vastus medialis, semitendinosus, biceps femoris long head, tibialis anterior, extensor digitorum, extensor hallucis, peroneus longus, soleus, gastrocnemius lateralis, gastrocnemius medialis, flexor digitorum) using a TelEMG system (BTS, Milan, Italy; 2000 Hz). The model, implemented in Matlab (The MathWorks, Natick, MA), is generic in that it can be used to study the muscle load-sharing problem of any musculoskeletal model and task of motion. The analysis is described in four parts: (i) the gait model; (ii) the calculation of the muscle force potential, hereinafter referred to as physiologically possible muscle forces; (iii) the calculation of the spectrum of muscle forces that represent the unpredictable variability of the muscle recruitment process, hereinafter referred to as physiologically plausible muscle forces; and (iv) data analysis.
2.2. The gait model
The musculoskeletal model, including the joint angles and torques while walking, was obtained from earlier studies [23,25,26] (figure 1). In summary, the lower-limb musculoskeletal model was a muscle-actuated articulated system based on the work of Delp et al. , whose anatomy was taken from the computed-tomography images and dissection of an 81-year-old female donor (63 kg weight, 167 cm height). The articulated system was a 13-segment, 15 degree-of-freedom system, actuated by 84 Hill-type muscle–tendon units. The inertial properties were calculated assuming homogeneous bone (1.42 g cm−3) and soft-tissue (1.03 g cm−3) density . The physiological cross-section area (PCSA) was calculated from the muscle volume and length. The peak isometric muscle stress was assumed to be equal to 1.37 MPa, the upper bound of published values . The remaining muscle parameters were based on the work of Delp et al. . The gait simulation used skin-mounted marker trajectories (Vicon Motion Capture, Oxford, UK; 100 Hz) and ground reaction forces at both feet (Kistler Instrument AG, Switzerland; 2000 Hz), recorded following the protocol proposed by Leardini et al. . Joint angles and torques were calculated using the inverse kinematic, dynamic and static optimization algorithms implemented in Opensim . The model yielded joint torques within published values, hip contact forces in agreement with in vivo measurements and muscle force patterns in good qualitative agreement with corresponding EMG recordings [23,25,26].
2.3. Physiologically possible muscle and joint forces
Physiologically possible muscle forces are defined as muscle forces within physiological boundaries, generating the joint torque from inverse dynamics and assuming that muscle activation can range from zero to full activation. Therefore, the muscle's force-generating potential is represented by the boundaries of physiologically possible muscle forces. For each walking frame, the instantaneous equilibrium equation at the joints (2.1), representing the muscle load-sharing problem, was determined by extracting the muscle lever arm, the muscle constraints and the net joint torques from an earlier simulation of walking . The equation takes the form 2.1 where is the matrix of muscle lever arms, is the muscle force vector, is the joint torque vector, and are respectively the lower and the upper muscle force boundaries, m is the number of degrees of freedom of the articulated system and n is the number of muscles in the model. The peak muscle force was calculated using a Hill-type muscle model. The active and passive force–length relationships were taken from the work of Thelen , whereas the force–velocity relationship was taken from the work of Delp et al. . Muscle force vectors within the spectrum were categorized using a single parameter, or muscle co-contraction, defined as the difference between the actual muscle force and the minimal force required to generate a given joint torque. Each muscle force vector, the solution to the muscle recruitment problem, was thus composed by a first minimal co-contraction component, represented by the static optimization solution, and a second component or muscle co-contraction force component. The muscle co-contraction level was assumed to be the fraction between the actual muscle co-contraction force component and the difference between the peak muscle force and the static optimization solution. The lower bound of muscle force, was set to zero, mimicking the muscle inability to sustain compressive forces. The upper bound of muscle force, was defined by studying five uniformly distributed co-contraction levels from zero (i.e. the optimization solution) to full co-contraction (i.e. the peak muscle force vector). Samples of physiologically possible muscle forces, solutions of equation (2.1), were calculated using Metabolica . The software interprets the vector of muscle forces as a multi-variate random variable characterized by its PDF and it samples the calculated PDF using an MCMC algorithm. The vector of muscle force was assumed to be uniformly distributed . Thus, the prior PDF of muscle forces takes the form 2.2 where takes on the value of 1 if all components of the argument are positive and vanishes otherwise. The posterior PDF describing how is distributed is 2.3 meaning that the posterior PDF of muscle forces is proportional to the prior PDF, , and the sensory information about the system states , or likelihood. Stable motion was defined by assuming the joint torque vector , the deterministic value calculated using inverse dynamics. For each walking frame, the MCMC algorithm was used to generate the ensemble of 200 000 samples whose entries are random realizations drawn from equation (2.3). The null space of the matrix , containing the muscle lever arm extracted from the model , is calculated using singular value decomposition. The vector is decomposed into a component lying on the null space and a component orthogonal to . Samples are drawn from the solution space using an MCMC algorithm by separately sampling the component using a hit-and-run algorithm and the component using a Gibbs algorithm .
The hip, knee and ankle reaction forces were calculated using the equation 2.4 where is the joint reaction force vector calculated using inverse dynamics, and is the ith joint-spanning muscle force vector.
2.4. Physiologically plausible muscle forces
Physiologically plausible muscle forces are defined as forces most likely to occur during normal gait and can be seen as a subgroup of the physiologically possible muscle forces. Therefore, physiologically plausible muscle forces were calculated by combining physiologically possible muscle forces and the variability of the muscle electrical activity from repeated EMGs.
Muscle forces were calculated using EMGs and Hill-type excitation- and contraction-dynamic models according to the guidelines proposed by Zajac . The raw EMG signal was band-pass filtered (zero-pole-gain design, eighth order, Butterworth filter) with cut-off frequencies of 10 and 400 Hz to minimize noise owing to motion artefacts and the EMG amplifier . The filtered EMG signal was rectified and low-pass filtered (zero-pole-gain design, second order, Butterworth filter) with a cut-off frequency of 6 Hz  and a 22 ms electromechanical delay, representing the muscle time response to stimuli, applied to synchronize the processed signal with the muscle response . Normalization of the processed EMG signal was then necessary to obtain a signal between zero and 1 representing muscle activation . We scaled the processed EMG signal to match the peak muscle activation calculated using static optimization . The EMG-driven muscle force was calculated using the calculated muscle activation, the active and passive force–length relationships from the work of Thelen  and the force–velocity relationship from the work of Delp et al.  for all seven gait repetitions. The force range, that is, the upper to the lower bound of muscle forces of physiologically plausible muscle forces, was assumed at the 0.68 quantile (i.e. mean ± s.d.) of the EMG-driven muscle force distribution projected onto the solution space of equation (2.1). Samples of physiologically plausible muscle forces were generated using Metabolica by constraining muscle forces within the calculated force range for selected muscles and between zero and the peak physiological force for the remaining muscles.
2.5. Data analysis
Simulations were run on a desktop PC (Windows 7, 64 bit, Intel Xenon E5-2630 v. 2, 2.60 GHz, 64 GB of memory). The gait cycle was divided into clusters of time frames and processed by 12 different CPUs using parallel computing. The speed in generating muscle force samples was output by the code.
The gait model was assessed by comparing the donor's PCSAs with corresponding measurements from donors aged 83 ± 9 years . The muscle lever arm and the joint torques in the model were compared with corresponding published values [36–40]. Calculation of muscle forces was verified by comparing the joint torque generated by the muscles with that calculated using inverse dynamics. The consistency between the EMG-driven muscle forces and the motion was assessed by using the distance between the EMG-driven muscle forces and the solution space of equation (2.1) as the metric; the average force distance over gait and muscles, and the muscle-by-muscle average distance during gait were calculated. The nearest muscle force vector to the EMG-driven force vector, henceforth referred to as the best-fit solution, was used for comparing the ability of the present method with that of a commonly used static optimization procedure in describing the recorded EMG pattern. To this purpose, we calculated (i) the linear regression between the EMG-driven muscle force and the best-fit muscle force and (ii) the linear regression between the EMG-driven muscle force and the static optimization solution obtained by minimizing the squared sum of muscle stress [27,31].
Alternative muscle recruitment strategies were studied in terms of physiologically possible and plausible muscle and joint forces. Physiologically possible forces were assessed by plotting the boundaries of muscle and joint forces for a progressive increase of the muscle co-contraction level. Physiologically plausible muscle forces were calculated by inputting to the model a subset of EMGs ; for this study, we used six of the principal lower-limb muscles spanning the hip, the knee and the ankle (gluteus maximus, rectus femoris, vastus lateralis, biceps femoris long head, tibialis anterior and gastrocnemius medialis). The available EMGs not input to the model were compared with the respective force spectrum from the model.
The process resulted in 103 and 20.6 M of different muscle forces, respectively, representing potential and plausible muscle forces driving walking. The algorithm generated 909 muscle force vectors, solution of equation (2.1), per second per processor.
The model anatomy and motion were in agreement with earlier studies. The muscle lever arms were consistent with those reported in earlier theoretical and experimental studies [37–43] (table 1). The average donor's muscle PCSA was 6.15 cm2, which represents the 25th lower percentile of the elderly population reported by Ward et al.  (table 2). The joint torques pattern was in agreement with that reported by Benedetti et al.  (figure 2). The highest unbalance between the joint torque driving walking and the net joint torque produced by the muscles was 3 × 10−11 Nm. The distance between the model and the EMG-driven force, averaged over gait and muscles, was below 40 N, whereas the peak muscle-by-muscle average distance over gait increased up to 108 N for 15 EMGs input to the model (figure 3). The best-fit solution better represented the EMG pattern than it did the static optimization solution (figure 4). The coefficient of determination between the best-fit solution and EMG-driven muscle forces was R2 = 0.94, and the average error was RMS = 19 N (figure 5). The static optimization solution showed major discrepancies in the muscle force pattern during (i) the early stance phase of walking for the rectus femoris, (ii) the stance-to-swing phase for the rectus femoris and the biceps femoris long head, and (iii) the late swing phase for the gluteus maximus and the tibialis anterior. The coefficient of determination between the static optimization solution and EMG-driven muscle forces was R2 = 0.38, and the average error was RMSE = 61 N.
Physiologically possible muscle synergies comprised muscle forces ranging from zero to the peak muscle force for most muscles. Twenty-seven out of the 34 lower-limb muscles ranged from zero to their peak force whereas seven muscles (gluteus maximus, adductor magnus, semimembranosus, vastus medialis, vastus lateralis, vastus medialis and soleus) could not reach their peak force. The resulting upper and lower boundaries of the hip, knee and ankle force spectrums showed typical double-peak patterns; the upper force boundary reached 11.3, 6.2, 7.6 body-weights (BW), and the lower boundary reached 4.4, 2.5 and 3.5 BW at the hip, the knee and the ankle, respectively. Increasing the upper boundary of the muscle force up to 60% muscle co-contraction caused a proportional shift upwards of the upper force boundary of possible muscle and joint forces. Further increasing the upper bound of muscle forces caused a complex nonlinear response of muscle and joint forces. Negligible changes of the lower force boundary were observed by allowing different muscle co-contraction levels (figures 6 and 7). Physiologically plausible muscle forces well represented the pattern of the muscle electrical activity (figures 8 and 9). The gluteus maximus showed a consistent double-peak activity, reaching its peak values during the early stance (5% gait) and mid-swing (65% gait) phases of walking. The biceps femoris long head peaked at 10%, 50% and 90% gait. The rectus femoris peaked at heel strike (7% gait) and prior to toe-off (50% gait). The vastus lateralis peaked at heel-strike (5% gait) and prior to toe-off (43% gait). The medial gastrocnemius peaked at 40% gait and showed a smaller second peak at mid-swing (70% gait). The tibialis anterior showed a double-peak activity, reaching its peak at early stance (5% gait) and mid-swing (75% gait). The EMG-driven force range for the gluteus medius, vastus medialis, semitendinosus, extensor hallucis, extensor digitorum, peroneus longus, soleus, gastrocnemius lateralis and flexor digitorum, which were not input to the model, was smaller than the calculated force range (figure 9).
The aim of this study was to investigate possible muscle synergies during walking. We used a human gait model in conjunction with Bayesian statistics, MCMC sampling method and EMG-driven muscle force modelling to calculate muscle forces in full respect of physiological and dynamical constraints. The gait model provided reliable information about all the relevant musculoskeletal parameters during walking, including muscle lever arm, muscle size and joint torques [36–39,45]. Muscle forces were calculated in an efficient manner, providing information about (i) the potential muscle and joint forces and (ii) the spectrum of muscle forces consistent with the muscle electrical activity input to the model.
The algorithm generated 909 muscle force samples per second per processor. Therefore, the present method can be used to calculate the spectrum of muscle forces during motion on standard desktop machines and can take advantage of parallel computing on multi-processor systems. The model could well represent, on average, EMG-driven muscle forces. However, the muscle-by-muscle distance between the model solutions and the EMG-driven muscle forces increased for an increasing number of EMG signals input into the model. This inconsistency may explain why EMG-driven models may not offer a solution when several EMGs are input to the model . Other authors optimized the model parameters within physiological boundaries, solving the model consistency problem and ensuring that a solution to the motion problem exists [18,20,46]. However, the optimized model probably provides little information about how the calculated solution represents the subject under study because of the typically large variability of physiological parameters. More work is necessary to understand how model simplifications and input errors influence model calculations. To this purpose, the proposed method is well suited to take in the variability of joint torques and EMGs, caused either by the natural unpredictable variability of motion or by uncertainties on measurements. The best-fit solution in the model more closely represented (R2 = 0.94) the muscle electrical activity than a static optimization procedure (R2 = 0.38) largely accepted for simulating normal walking [11,31], without requiring any assumptions about the adopted sensorimotor behaviour. Therefore, the proposed method can be used to calculate muscle forces when the objective of the sensorimotor behaviour is variable or not known, including in the instances of either physiological or pathophysiological neuromotor conditions.
Physiologically possible muscle synergies include muscle forces ranging from zero to the peak muscle force (figure 6) and joint contact force of up to 11.3 BW at the hip (figure 7). Up to 60% muscle co-contraction caused a linear increase of muscle and joint forces, whereas higher muscle co-contraction caused a nonlinear increase of the same quantities (figures 6 and 7). While the probability for these extreme loading conditions to occur has to be determined, these findings may have implications in studying muscle ability to control joint impedance and stability , the as-yet unresolved fracture mechanism for low-energy osteoporotic fractures , and may reveal important information for the development of exercise therapies for bone health . Physiologically plausible muscle force patterns well represented the muscle electrical activity input to the model (figure 9). Therefore, the proposed approach can be used to study deep aspects of human motion. For example, the calculated spectrum can be used for exploring how different muscles can combine their action in response to the same motor demand.
To the best of the authors' knowledge, this is the first numerical study exploring the spectrum of muscle synergies during motion. The model has been shown capable of yielding kinematic, kinetics, hip contact force and muscle firing patterns in agreement with published patterns for multiple activities [23,25], providing confidence in the reliability of the studied muscle load-sharing problem. The large variability of physiologically possible and physiologically plausible muscle forces is consistent with the known ability of the CNS for adopting very different muscle recruitment strategies [16,49]. Although no measurements of the joint contact force under full muscle co-contractions are available, the range of the calculated hip contact force (3.7–11.4 BW) compares well with the hip contact force of 8.7 BW measured by Bergmann et al.  during stumbling, a value that has been largely attributed to muscle co-contraction rather than to motion dynamics. This provides confidence in the calculated spectrum of muscle forces.
The main limitation of this study is that the majority of the muscle forces that were not constrained between EMG-driven muscle force boundaries (e.g. the semitendinosus) showed a much higher variability (figure 9) than those obtained from repeated EMG recordings, indicating that the calculated spectrum of physiologically plausible muscle forces is larger than that observed in vivo. While it is possible that a reduced number of EMGs input to the model  explain the majority of the muscle force variability, the optimal number and type of EMG signals has to be determined. Second, the processed EMG signal was normalized using the peak muscle activation calculated using static optimization, whereas others normalized the processed EMGs to a maximum voluntary contraction task [18,19]. However, a standardized EMG normalization process has yet to be defined . Third, the present results cannot be generalized owing to the single anatomy used. More research is necessary to solve this limitation. Fourth, the joint torque was set to the deterministic values calculated using inverse dynamics, thus neglecting the joint torque uncertainties attributable to model assumption and measurement errors . However, this allowed the isolated effect of alternative muscle recruitment strategies on calculated muscle and joint forces to be studied. Last, the peak isometric muscle stress was the upper boundary of published values (1.37 MPa ), possibly causing an overestimation of calculated forces. However, the upper boundary of muscle forces is a linear function of the peak isometric muscle stress, whereas the lower boundary is almost invariant . Therefore, the boundaries of muscle and joint forces for every intermediate value of the peak isometric stress can be easily extrapolated.
Despite the study limitations, present findings are important for the biomechanics community in that they provide a viable numerical approach for modelling the stochastic nature of the muscle recruitment process. The present results strengthen the notion that muscle co-contraction is important in studying human motion, and provide a viable numerical approach for studying physiological and pathophysiological conditions characterized by complex sensorimotor behaviours. Moreover, because the proposed approach makes no assumptions on the ‘normality’ of neuromotor control, we expect it to be equally effective in subjects affected by severe neuromuscular pathologies.
This study was supported by the Australian Research Council project (DE140101530) awarded to S.M. The work of E.S was partially supported by NSF-DMS grant no. 1312424.
One contribution of 11 to a theme issue ‘Multiscale modelling in biomechanics: theoretical, computational and translational challenges’.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.