- Research
- Open access
- Published:

# A fair and EMG-validated comparison of recruitment criteria, musculotendon models and muscle coordination strategies, for the inverse-dynamics based optimization of muscle forces during gait

*Journal of NeuroEngineering and Rehabilitation*
**volume 18**, Article number: 17 (2021)

## Abstract

Experimental studies and EMG collections suggest that a specific strategy of muscle coordination is chosen by the central nervous system to perform a given motor task. A popular mathematical approach for solving the muscle recruitment problem is optimization. Optimization-based methods minimize or maximize some criterion (objective function or cost function) which reflects the mechanism used by the central nervous system to recruit muscles for the movement considered. The proper cost function is not known a priori, so the adequacy of the chosen function must be validated according to the obtained results. In addition of the many criteria proposed, several physiological representations of the musculotendon actuator dynamics (that prescribe constraints for the forces) along with different musculoskeletal models can be found in the literature, which hinders the selection of the best neuromusculotendon model for each application. Seeking to provide a fair base for comparison, this study measures the efficiency and accuracy of: (i) four different criteria within the static optimization approach (where the physiological character of the muscle, which affects the constraints of the forces, is not considered); (ii) three physiological representations of the musculotendon actuator dynamics: activation dynamics with elastic tendon, simplified activation dynamics with rigid tendon and rigid tendon without activation dynamics; (iii) a synergy-based method; all of them within the framework of inverse-dynamics based optimization. Motion/force/EMG gait analyses were performed on ten healthy subjects. A musculoskeletal model of the right leg actuated by 43 Hill-type muscles was scaled to each subject and used to calculate joint moments, musculotendon kinematics and moment arms. Muscle activations were then estimated using the different approaches, and these estimates were compared with EMG measurements. Although no significant differences were obtained with all the methods at statistical level, it must be pointed out that a higher complexity of the method does not guarantee better results, as the best correlations with experimental values were obtained with two simplified approaches: the static optimization and the physiological approach with simplified activation dynamics and rigid tendon, both using the sum of the squares of muscle forces as objective function.

## Introduction

Determination of muscle forces during gait is of great interest to extract the principles of the central nervous system (CNS) control, to facilitate assessment of pathological gait, or to estimate the loads on bones and joints (prevention of injuries in sports, surgical planning to reconstruct diseased joints) [1,2,3]. The invasive character of in vivo experimental measurements, and the uncertain relation between muscle force and EMG, makes computer modeling and simulation a useful substitutive approach. Determination of muscle forces by computer modeling and simulation was extensively treated and numerous approaches can be found in the literature to solve the redundancy problem of the muscle recruitment, as well as to represent the musculotendon actuator dynamics [4,5,6,7,8], where each author highlights the advantages of his own approach. Ambrosio and Kecskemethy [4] suggested that the physiological criteria improve human motion prediction, Hardt [5] defended the linear constraint programming, Anderson and Pandy [6] showed that static and dynamic optimization solutions for gait are practically equivalent, Millard et al. [7] offered a damped equilibrium musculotendon model to improve the physiological optimization, and Shourijeh et al. [8] found that forward static optimization was a suitable method for solving forward dynamic musculoskeletal simulations. However, results do not depend only on the approach used, but also, on the experimental data collection and on the musculoskeletal model used, which makes more difficult for the readers to select objectively which approach to use for a certain application [9].

To objectively compare different approaches, it is necessary to test them under the same conditions. When proposing a new approach, authors generally make a comparison with experimental measurements in order to validate it [10, 11], and, in some cases, they compare their results with those provided by a previous approach which gives confidence to readers [12, 13]. In some applications, a benchmark problem can be found that establishes some defined conditions, so that researchers can get a fair comparison [14]. In the case of the resolution of the muscle force-sharing problem, the authors only found a benchmark where the computational speed and biological accuracy of three musculotendon models was compared during simple muscle-driven simulations [7]. The physiological effect of static and dynamic optimization during gait was compared by Pandy et al. [10] and De Groote et al. [15], but none of the studies offered an experimental validation that allowed to conclude which method provided the most realistic results.

Few years ago, a grand challenge competition to predict in vivo knee loads was organized by some researchers who shared their experimental data collections for the analysis and its evaluation [16]. However, the musculoskeletal modeling could differ between participants so that, by using a different multibody model (with different degrees of freedom) and different muscle geometry (which implies different arm moments), the results could not dissociate the effect of neuromusculotendon models.

In this context, the aim of the present work is to compare the efficiency and accuracy of: i) four different cost functions within the static optimization approach (where the physiological character of the muscle, which affects the constraints of the forces, is not considered); ii) one static and three physiological representations of the musculotendon actuator dynamics: activation dynamics with elastic tendon, simplified activation dynamics with rigid tendon, and rigid tendon without activation dynamics; iii) a synergy-based method; all of them within the framework of inverse-dynamics based optimization. The best cost function (criterion) obtained in (i) was used for (ii) and (iii). Motion/force/EMG gait analyses were performed on ten healthy subjects. A musculoskeletal model of the right leg actuated by 43 Hill-type muscles was scaled to each subject and used to calculate joint moments, musculotendon kinematics and moment arms. Therefore, the muscle force-sharing problem was solved under the same conditions and using the same inputs. Muscle activations were then estimated using the different approaches, and these estimates were compared with EMG measurements which served as experimental reference.

## Methods

### Experimental data collection

Ten subjects (seven males, three females, age 42 ± 16 years, height 173 ± 16 cm, body mass 73 ± 26 kg) were recruited for this study. All subjects gave written informed consent for their participation. Subjects walked at their self-selected speed (1.1 ± 0.2 m/s) along a walkway with two embedded force plates (AccuGait, sampling at 100 Hz; AMTI, Watertown, MA, USA). The motion was captured using 12 optical infrared cameras (OptiTrack FLEX 3, also sampling at 100 Hz; Natural Point, Corvallis, OR, USA) that computed the position of 36 optical markers (Fig. 1). An extended Kalman filter (EKF) was used to filter the marker trajectories and reconstruct the motion with a process noise variance of 1 m/s2 and a cutoff frequency of 20 Hz [17]. Additionally, 9 surface EMG signals were recorded from the right leg at 1 kHz (BTS, FREEEMG, Quincy, MA, USA). The participants' right leg was shaved, the skin was cleaned with alcohol, and the electrodes were placed according to the guideline presented in [18]. EMG signals were acquired from the following muscles: tibialis anterior, vastus medialis, vastus lateralis, gastrocnemius medialis, gastrocnemius lateralis, semitendinosus, biceps femoris, gluteus maximus and gluteus medius. Each EMG signal was rectified, filtered by singular spectrum analysis (SSA) with a window length of 250 [19] (equivalent to the common forward and reverse low-pass 5th order Butterworth filter with a cut-off frequency of 6 Hz), and, then, normalized with respect to its maximal value. This cut-off frequency value is consistent with the ranges reported in previous lower extremity studies using EMG data [20,21,22,23,24].

### Musculoskeletal model

The human body was modeled as a three-dimensional multibody system formed by rigid bodies (Fig. 1, left and center). The model consisted of 18 anatomical segments [25]: two hindfeet, two forefeet, two shanks, two thighs, a pelvis, a torso, a neck, a head, two arms, two forearms, and two hands. The segments were linked by ideal spherical joints, thus defining a model with 57 degrees of freedom (DOFs). The axes of the global reference frame were defined as follows: *x*-axis in the anterior–posterior direction, *y*-axis in the medial–lateral direction, and *z*-axis in the vertical direction. The computational model was defined with 228 mixed (natural + angular) coordinates. The subset of natural coordinates comprised the three Cartesian coordinates of 22 points and the three Cartesian components of 36 unit vectors, thus yielding a total of 174 variables.

Matrix-R formulation [26] was used to perform an inverse-dynamics analysis to obtain the joint torques along the motion by means of the in-house developed MBSLIM library [27] programmed in FORTRAN, as described in [28]. Once the joint torques were computed, it was assumed that 43 right leg muscles contributed to the following six right-leg inverse-dynamics moments: the three rotational DOFs at the hip, the flexion/extension DOF at the knee, and the plantar/dorsi flexion and inversion/eversion at the ankle. Muscles were modeled as one or more straight-line segments with via points. These points corresponded to the attachments of muscle and tendon to bone and were defined as the origin (i.e., proximal attachment) and insertion (i.e., distal attachment). Muscle properties and local coordinates for these points were obtained from OpenSim (model Gait2392) [29] and scaled to each subject from the generic reference OpenSim model, as commented further in 2.4.

### Optimization problem

As introduced before, the fundamental problem is that there are more muscles serving each degree of freedom of the system than those strictly necessary from the mechanical point of view. In this case, there are 43 muscles at the leg to actuate 6 degrees of freedom (other degrees of freedom of the leg are controlled by joint structures as bones and ligaments, yielding a reaction moment instead a drive torque). Consequently, there is an infinite number of solutions for this problem and, in order to reproduce the specific strategy of muscle coordination adopted by the CNS, optimization is used.

The inverse-dynamics based optimization problem that serves to determine the muscle forces at each time-point can be formulated in general form as:

where \(C\) is the cost function, \({\mathbf{Q}}^{ID}\) is the vector of inverse-dynamics joint moments at the right leg (where the force-sharing problem is addressed), \({\mathbf{F}}^{MT}\) is the vector of muscle forces, \({\mathbf{J}}\) is the Jacobian whose transpose projects the muscle forces into the joint drive torques space, \(F_{i}^{Min}\) and \(F_{i}^{Max}\) are the instantaneous minimum and maximum allowed forces in muscle *i*, respectively, and *m* is the number of muscles. Expression of the objective function *C* depends on the muscle recruitment criterion used. In the literature, several muscle recruitment criteria have been suggested to represent the CNS behavior. In this work, four of them have been considered in the context of static optimization.

#### Static optimization (SO)

##### Nonlinear polynomial criteria

The polynomial criterion can be written as

where \(k_{i}\) denotes a positive weighting factor and \(w\) is the power of the polynomial. According to [30], the muscle force prediction obtained by minimizing the sum of muscle stresses raised to a power *w* whose value ranges between 1.4 and 5.1 is physiologically analogous to minimizing muscle fatigue. As Anderson and Pandy did in their study [10], a power of 2 was chosen.

##### Criterion I—minimization of the sum of the squares of muscle forces

##### Criterion II—minimization of the sum of the squares of relative muscle forces

with \(F_{0}^{M}\) the maximum isometric force from [29].

##### Criterion III—minimization of the sum of the squares of muscle stresses

with \(PCSA\) the physiological cross sectional area from [31].

##### Min/max criterion

The min/max criterion distributes the collaborative muscle forces in such a way that the maximum relative muscle force is as small as possible. Therefore, the largest endurance for a task is attained when the maximum relative muscle force [32] or the maximum muscle stress [33] is as small as possible. The min/max criterion takes the form:

For this study, the following criterion is used:

##### Criterion IV—minimization of the largest relative muscle force

For SO, the physiological behavior of the musculotendon actuator dynamics is not considered, so, the limit values of the muscular forces are \(F_{i}^{Min} = 0\) and \(F_{i}^{Max} = F_{0,i}^{M}\).

##### Physiological approach (PHY1)

At physiological level, musculotendon actuator dynamics introduces muscle force constraints. Whereas the static optimization approach disregards these constraints in order to simplify the problem, the so-called physiological approach [34] takes them into consideration. This approach applies optimization techniques at each time-point (Fig. 2, right), and prescribes minimal and maximal constraints for the forces by extrapolating the force values from the previous time-point through feasible muscle dynamics (Fig. 2, left).

The dynamics of musculotendon actuators can be divided into two parts. First, the activation dynamics which corresponds to the transformation of a neural excitation sent by the brain into an activation of the contractile apparatus. Activation dynamics is described by a first-order ordinary differential equation that contains the relationship among the muscle activation \(a\), its derivative \(\dot{a}\), and the neural excitation \(u\) as:

with \(\tau = \tau_{act}\) when \(a(t_{k - 1} ) \le u(t_{k} )\) and \(\tau = \tau_{deact}\) when \(a(t_{k - 1} ) > u(t_{k} )\). The activation and deactivation time constants \(\tau_{act}\) and \(\tau_{deact}\) are set to 15 ms and 50 ms, respectively [35, 36].

Second, this activation is transformed into a muscle force by the second phase, the contraction dynamics. The force generated by a muscle is constrained by its force–length-velocity properties, related to the Hill-type musculotendon model used (Fig. 3), which is defined by this second differential equation:

The musculotendon length \(l^{MT}\) and velocity \(v^{MT}\) depend on the position and velocity of the body segments and, in turn, the generated tendon force \(F^{MT}\) affects the motion of the body segments. Thus, there exists interaction between muscles and body segments.

The complete musculotendon dynamics can be expressed as a system of two differential equations which can be written, in a simplified form, as.

This system is used to define the minimal and maximal constraints for the forces by extrapolating the force values from the previous time-point using feasible muscle dynamics, integrating (10) with \(u(t) = 0\) to estimate \(F^{Min} (t)\) and with \(u(t) = 1\) to estimate \(F^{Max} (t)\) (the muscular excitation is assumed to be constant between the two time frames). Using the physiological approach, the initial activations and muscles forces are needed. The determination of initial activations and muscular forces is based on the static condition which states that the initial fiber velocity of each muscle is set to zero and, \(F^{Min}\) and \(F^{Max}\) correspond to \(a = 0\) and \(a = 1\), respectively.

Integration is carried out with Matlab *Ode23t*. Two integrations per muscle are required at each time step, which make the optimization and integration process heavy and slow. By programming the muscular functions of the Hill-type musculotendon model into a FORTRAN mex file, the computational time is reduced by a factor of 10. However, the computational time is still long and, in addition, the high tendon stiffness makes really difficult to use this approach (a suitable scaling of muscle parameters is needed) [37]. Therefore, in order to simplify the problem while keeping some physiological characteristics, most authors prefer to use a Hill-type musculotendon model with a rigid tendon [13, 15, 38].

### Physiological approach with rigid tendon

In this way, the tendon length is constant and, consequently, the muscle fiber length and velocity depend only on the musculoskeletal geometry as well as on body segment configurations (which affect \(l^{MT}\) and \(v^{MT}\)) and not the musculotendon force. Consequently, the force–length-velocity allowed is expressed as:

Use of the rigid tendon model avoids the two integrations needed to calculate the limits of the muscle force at each instant. Then, in order to further reduce the computational burden, the first-order ordinary differential Eq. (8) used to estimate the muscular activation, \(a\), can be simplified as follows:

### Time response considered (PHY2)

In order to keep the muscular time response relation given by (8), the first-order ordinary differential equation can be converted into:

with \(\tau = \tau_{act}\) when \(a(t_{k - 1} ) \le u(t_{k} )\), \(\tau = \tau_{deact}\) when \(a(t_{k - 1} ) > u(t_{k} )\) and \(\Delta t\) is the time step.

Therefore, the minimal and maximal muscular force constraints of the optimization problem can be obtained through integration with any of the described approaches by extrapolating the force values from the previous time-point using feasible muscle dynamics.

### Time response ignored (PHY3)

However, authors who consider the tendon as a rigid element usually choose to ignore the muscular time response and assume that:

In this work, the minimization of the sum of the squares of muscle forces (Criterion I) was used as objective function for the three physiological models.

#### Synergy optimization (SynO)

The fact that synergies take a high dimensional control space and reduce it to a low dimensional space is potentially useful for reducing the amount of indeterminacy when estimating muscle forces via optimization. For this reason, some authors started to investigate how to include it to solve the muscle force-sharing problem.

The synergy optimization (SynO) approach used in [13] estimates muscle forces during human walking using synergy-constructed muscle activations, similar to the more complex approach proposed in [39]. SynO finds muscle forces that match the inverse-dynamics joint moments as closely as possible through the moment tracking error term in the cost function. In SynO, synergies couple muscle activations across time frames, requiring the optimization to be performed over all the time frames simultaneously as follows:

where \(T_{{f \times n_{S} }} (T_{p} )\) and \(V_{{n_{S} \times m}} \,\) are the time-varying synergy activations defined by B-spline nodes, and the corresponding time-invariant synergy vectors, respectively. Each muscle activation synergy is composed of a single time-varying synergy activation defined by *p* = *(f-1)/5* + *1* (nearest integer, *f* = number of frames) B-spline nodal points along with its corresponding time-invariant synergy vector defined by *m* = 43 weights specifying inter-muscle activation coupling. Thus, for *n*_{S} synergies (*n*_{S} = 3 in this study), the number of design variables is *n*_{S}**(p* + *m)*. Muscle synergy quantities are used as the design variables for synergy optimization. Each optimization problem is theoretically over-determined. However, in practice, the problems remain under-determined since neighboring time frames are not completely independent from one another.

Using these design variables, the SynO cost function is formulated as follows:

where \(a_{ij}^{*}\) are the synergy-based muscle activations, and \(\lambda_{ij,pen} = \left\{ \begin{gathered} 0{ 0} \le a^{*}_{ij} \le 1 \hfill \\ 10^{5} {\text{ otherwise}} \hfill \\ \end{gathered} \right.\) are penalization factors for muscle *i* at the time frame *j* to ensure that muscle activations stay between 0 and 1. While previous approaches enforce the muscle forces to exactly reproduce the inverse-dynamics joint moments through its equality constraints, the SynO approach minimizes the error between \(Q_{{}}^{ID}\) and \(Q_{{}}^{MT}\), being \(Q_{{}}^{MT}\) the joint moments produced by the muscle forces estimated by SynO. A scale factor \(\beta = 100\) is applied to achieve the best compromise between joint moment tracking and activation minimization [40]. In his previous work [40], Michaud found that best correlations with experimental EMG patterns were obtained using three synergies for SynO, and that the resulting mean joint intersegmental moment matching between \(Q_{{}}^{ID}\) and \(Q_{{}}^{MT}\) across subjects was higher than 96%. For this reason, the SynO approach will be used in this study with only three synergies.

The objective function is programmed as a FORTRAN mex file to reduce computation time (16 times faster than the original Matlab function). Linear equality constraints enforce that the sum of weights within each synergy vector is equal to 1, which makes the synergy construction unique, while lower bound constraints enforce the synergy activation B-spline nodes and synergy vector weights to be non-negative. The same musculotendon model used for the PHY3 approach is used here.

#### Subject-specific scaling of musculotendon parameters

Due to the sensitivity of physiological approaches [37], a suitable scaling of musculotendon parameters is needed. In addition to the high tendon stiffness which makes implementation really difficult, some Hill-muscle equations become numerically stiff when numerical singularities are approached [7]. Since these conditions are often encountered during a simulation, to prevent that the solver gets stuck at points that were numerically feasible yet not physiologically sound [41] (which slows the process of numerical integration), a scaling correction was applied. Length parameters were scaled in two steps. First, for each muscle, the tendon slack length (\(l_{S}^{T}\)) and the optimal muscle fiber length (\(l_{0}^{M}\)) were scaled with a scale factor calculated as the relation between the subject’s musculotendon length in standing position and that of the generic model in the same position. As the pennation angle of the reference, \(\alpha_{0}\), is kept, the scaled distance between the aponeuroses of muscle origin and insertion, \(w\) (which remains constant during the muscle contraction), is given by:

Then, because tendons are so stiff that their lengths do not change significantly during movement, the approximated muscle fiber length, \(l^{{M*}}\), can be calculated for each muscle during the complete gait cycle as follows:

with \(\alpha (t) = \arctan \left( {\frac{w}{{l_{{}}^{MT} (t) - l_{S}^{T} }}} \right)\). Finally, in order to keep the normalized muscle lengths (\(\overline{{l_{{}}^{M} }} = \frac{{l^{M} }}{{l_{0}^{M} }}\)) within the physiological optimal conditions (\(0.5 < \overline{{l^{M} }} \le 1.2\)) [42], the final scaled \(l_{0}^{M}\) was set to the maximum approximated muscle fiber length along the motion.

### Optimization protocol and EMG comparison

Optimization seeks to find the best solution from all the feasible ones by minimizing the objective function. Finding the global minima of a function is really difficult because of the many local minima. In order to get the best possible results, the following protocol is used in this work. Although each optimization problem is solved using the Matlab’s *fmincon* nonlinear constrained optimization algorithm, five global optimizations are run using Matlab’s *ga* genetic optimization algorithm with a population size of 50 to provide random initial guesses for *fmincon*. The solution with the lowest objective function value is chosen as initial guess for the initial time point. Thereafter, as muscle activation is normally smooth and continuous during gait, the optimal solution from the previous time frame is used as the initial guess for the current time frame [43].

Matching between estimated muscle activations and EMG was quantified via cross-correlation using the Pearson correlation coefficient *r* (Matlab’s function *corrcoef*) with a maximum time delay of 150 ms [44]. The correlation coefficient *r* was chosen to compare muscle activations and EMG data so as to focus on shape rather than on magnitude discrepancies, as there is no direct relationship between EMG and muscle force amplitude [45, 46].

## Results

The different approaches presented in this study were compared with EMG measurements for the ten healthy subjects. Normalized muscle activations during a gait cycle of one healthy subject estimated by SO using all the criteria are plotted in Fig. 4 along with the corresponding normalized EMG measurements. Comparison of muscle activations estimated with the different criteria are significantly different but show some similarities too.

Then, to observe the physiological effect of the musculotendon model, four approaches using criterion I were compared: normalized muscle activations during a gait cycle of one healthy subject estimated by static optimization (SO-I) and the three physiological approaches (PHY1, PHY2 and PHY3) are represented in Fig. 5. While PHY3 presented distinct results, muscle activations estimated by SO-I, PHY1 and PHY2 were very similar. PHY1 and PHY2 showed almost the same results.

Furthermore, the normalized muscle activations during a gait cycle of one healthy subject estimated by PHY3 and synergy optimization with 3 synergies (SynO3) are compared in Fig. 6 to highlight the effect of the synergy structure. Both used the same musculotendon model (physiological approach with rigid tendon and activation time response ignored). However, results are significantly different.

Mean across subjects Pearson correlation coefficient *r* values between EMG vs. muscle activations of all the approaches of this study for the nine muscles are reported in Table 1. Mean correlations of the many approaches did not present as significant differences as those observed previously for one subject. Mean values of the different approaches are close, between 0.61 (SO-IV and SynO3) and 0.74 (SO-I and PHY2). Approaches SO-I, PHY1 and PHY2 show almost the same correlations (means of 0.73 and 0.74). The paired sampled t-test realized (Table 2) showed no statistical differences (p < 0.05) between SO-I, SO-II, PHY1, PHY2 and PHY3, while differences were observed with the other approaches.

Mean across muscles Pearson correlation coefficient *r* for the ten subjects reported in Table 3 offered similar results as Table 1. Good correlations (*r* > 0.60, in italics) and close results were obtained for most approaches. In Table 4, the paired sampled t-test yielded the same conclusions as Table 2: SO-I, SO-II, PHY1, PHY2 and PHY3 are statistically similar (p < 0.05).

Besides, the computational efficiency of the different approaches studied in this work is compared in Table 1. All calculations were performed on an Intel® Core™ i7-6700 K processor running at 4.00 GHz with 16 GB of RAM. With a mean computational time of 2.5 s to simulate a gait cycle, SO-I is the fastest approach, while PHY1 is the slowest one (225.2 s).

Finally, the resulting joint reaction forces at hip, knee and ankle for the different muscle recruitment approaches are compared (Fig. 7 and Table 4). Similar joint reaction forces are obtained with the non-synergy-based approaches: SO-I, PHY1 and PHY2 offer almost the same results, while SO-II and PHY3 show joint reactions slightly higher than their counterparts at hip and knee levels. However, the joint reaction forces at hip, knee and ankle calculated from the muscle forces estimated with SynO3 are much higher than those obtained with the other approaches (Table 5).

## Discussion

This work offers a comparison of the efficiency and accuracy of: (i) four different criteria; (ii) three different physiological representations of the musculotendon actuator dynamics; (iii) a synergy-based method; all of them in the framework of inverse-dynamics based optimization. All the approaches were used under the same conditions by taking the same inputs from motion/force/EMG gait analyses performed on ten healthy subjects. Results obtained with the different methods do not present large discrepancies. Higher complexity of the method does not guarantee better results, as the best correlations with experimental values were obtained with two simplified approaches.

First, muscles activations obtained from SO and four different criteria exhibit visually different shapes along with some similarities. In addition, mean across subjects Pearson correlation coefficient *r* values between EMG vs. muscle activations of all the criteria do not present significant discrepancies. The best correlations were obtained with the simplest and fastest criterion (SO-I), which yielded a correlation of 74%, while the worst correlations were obtained with the most involved criterion (SO-IV).

Second, it was observed that the physiological representation of the musculotendon actuator dynamics does not affect the estimation of muscle forces during gait. Muscle activation shapes, experimental correlations and joint reaction forces are almost the same as those obtained through the non-physiological method (SO). The paired sampled t-tests demonstrated that SO-I, PHY1, PHY2 and PHY3 are statistically similar (p < 0.05). The same conclusion was drawn by De Groote et al. [15], Anderson and Pandy [6] and Millard et al. [7]. However, for faster, higher-powered tasks, like running or jumping, a compliant tendon model could be preferable. Moreover, despite its disadvantages (harder to implement and higher computational time), the physiological approach served to implement some Hill-based energy expenditure methods [47, 48], since it provides the muscular variables required as inputs. Despite the physiological realism of these approaches, it must be pointed out that the Hill-based muscle dynamics does not include the activation of gamma motor neurons and stretch receptors (i.e., proprioceptive receptors), which can induce minor activation in the stretched muscles during gait [49]. Consequently, some discrepancies can be explained by the EMG-linear envelope extraction procedure that may not properly demodulate neural excitations from the motor neurons action potentials, or because the noise contamination from cross-talk and movement artifacts, two intrinsic limitations in surface EMG measurements in addition of the inability to access deep muscles [23].

Third, as previously observed in [40], the synergy structure imposed within the SynO approach did not improve prediction of muscle activations during gait. This approach showed significant differences with the best approaches offered in this work (SO-I, PHY1, PHY2 and PHY3). The muscle synergy hypothesis has been notoriously difficult to prove or falsify [50], and results of this study do not allow to draw a conclusion in this regard. It can only be said that the SynO approach offers reasonable prediction of muscle activations and that its reduced dimensional control space could be beneficial for applications such as epidural electrical stimulation [51] or motion control and prediction [38].

Finally, all the estimated joint reaction forces at the hip were higher than the direct experimental measurements reported in the literature [52,53,54]. Brand et al. reported that hip contact-force predictions in the literature are higher than force measurements because of modeling assumptions [54]. In this work, and in the literature [9, 54], it has been shown that, paradoxically, physiological representation of the musculotendon actuator dynamics increases rather than reduces the discrepancies between force predictions and measurements, due to its constraints. Same conclusion can be drawn for the SynO approach, which disproportionately increases the joint forces due to its imposed synergy structure and reduced dimensional control space. Shourijeh and Fregly observed that the joint stiffness results were visibly different between the SynO and SO solutions, and that the stiffness decreased as the number of synergies was increased [13].

## Conclusion

In conclusion, this study evaluated several approaches to predict muscle activations during gait by comparing them with EMG measurements obtained experimentally, and found that higher complexity of the method does not guarantee better results. No significant differences among predicted EMG patterns within different physiological representations of the musculotendon were found. However, the simplified physiological approach with rigid tendon and activation time considered, presented the best accuracy and a very competitive computational time.

## Availability of supporting data

The datasets generated for this study are available on request to the corresponding author.

## References

Michaud F, Mouzo F, Lugrís U, Cuadrado J. Energy Expenditure Estimation During Crutch-Orthosis-Assisted Gait of a Spinal-Cord-Injured Subject. Front Neurorobot. 2019;13:1.

M. Daniel, “Mathematical simulation of the hip joint loading,” Czech technical university in Prague, 2004.

Gervais B, Vadean A, Raison M, Brochu M. Failure analysis of a 316L stainless steel femoral orthopedic implant. Case Stud Eng Fail Anal. 2016;5–6:30–8.

Ambrosio JAC, Kecskemethy A. Multibody dynamics of biomechanical models for human motion via optimization.,” in

*In J.C. Garcia Orden, J.M. Goicolea, J. Cuadrado (Eds.) Multibody Dynamics – Computational Methods and Applications*, Dordrecht: Springer, 2007, pp. 245–270.Hardt DE. Determining muscle forces in the leg during normal human walking - an application and evaluation of optimization methods. J Biomech Eng. 1978;100:72–8.

Anderson FC, Pandy MG. Static and dynamic optimization solutions for gait are practically equivalent. J Biomech. 2001;34(2):153–61.

Millard M, Uchida T, Seth A, Delp SL. Flexing computational muscle: Modeling and simulation of musculotendon dynamics. J Biomech Eng. 2013;135:2.

Shourijeh MS, Mehrabi N, McPhee J. Forward static optimization in dynamic simulation of human musculoskeletal systems: a proof-of-concept study. ASME J Comput Nonlinear Dyn. 2017;12(5):051005.

Wesseling M, et al. Muscle optimization techniques impact the magnitude of calculated hip joint contact forces. J Orthop Res. 2015;33:3.

Pandy MG, Anderson FC, Hull DG. A parameter optimization approach for the optimal control of large-scale musculoskeletal systems. J Biomech Eng. 1992;114(4):450–60.

Heintz S, Gutierrez-Farewik EM. Static optimization of muscle forces during gait in comparison to EMG-to-force processing approach. Gait Posture. 2007;26(2):279–88.

Michaud F, Lugris U, Ou Y, Cuadrado J, Kecskemethy A. Comparison of forward-dynamics approaches to estimate muscular forces in human gait. In:

*The 4th Joint International Conference on Multibody System Dynamics*, 2016, pp. 4–5.Shourijeh MS, Fregly BJ. Muscle synergies modify optimization estimates of joint stiffness during walking. J Biomech Eng. 2020;142:1.

González M, González F, Luaces A, Cuadrado J. A collaborative benchmarking framework for multibody system dynamics. Eng Comput. 2010;26(1):1–9.

De Groote F, Kinney AL, Rao AV, Fregly BJ. Evaluation of direct collocation optimal control problem formulations for solving the muscle redundancy problem. Ann Biomed Eng. 2016;44(10):2922–36.

Fregly BJ, et al. Grand challenge competition to predict in vivo knee loads. J Orthop Res. 2012;30(4):503–13.

Cuadrado J, Michaud F, Lugrís U, Pérez Soto M. Using accelerometer data to tune the parameters of an extended kalman filter for optical motion capture: Preliminary application to gait analysis. Sensors. 2021. https://doi.org/10.3390/s21020427.

Criswell E, Cram JR. Cram’s Introduction to Surface Electromyography. 2011; 340–376.

Romero F, Alonso FJ, Gragera C, Lugrís U, Font-Llagunes JM. Estimation of muscular forces from SSA smoothed sEMG signals calibrated by inverse dynamics-based physiological static optimization. J Brazilian Soc Mech Sci Eng. 2016;38(8):2213–23.

Buchanan TS, Lloyd DG, Manal K, Besier TF. Estimation of muscle forces and joint moments using a forward-inverse dynamics model. Med Sci Sports Exerc. 2005;37(11):1911–6.

Shourijeh MS, Flaxman TE, Benoit DL. An approach for improving repeatability and reliability of non-negative matrix factorization for muscle synergy analysis. J Electromyogr Kinesiol. 2016;26:36–43.

Bisi MC, Stagni R, Houdijk H, Gnudi G. An EMG-driven model applied for predicting metabolic energy consumption during movement. J Electromyogr Kinesiol. 2011;21(6):1074–80.

Sartori M, Farina D, Lloyd DG. Hybrid neuromusculoskeletal modeling to best track joint moments using a balance between muscle excitations derived from electromyograms and optimization. J Biomech. 2014;47(15):3613–21.

Lloyd DG, Besier TF. An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo. J Biomech. 2003;36(6):765–76.

Lugrís U, Carlín J, Pàmies-Vilà R, Font-Llagunes JM, Cuadrado J. Solution methods for the double-support indeterminacy in human gait. Multibody Syst Dyn. 2013;30(3):247–63.

García de Jalón J, Bayo E. Kinematic and Dynamic Simulation of Multibody Systems. Berlin, Heidelberg: Springer-Verlag; 1994.

Dopico D. MBSLIM: Multibody Systems en Laboratorio de Ingeniería Mecánica, http://lim.ii.udc.es/MBSLIM. 2016.

Lugris U, Carlin J, Luaces A, Cuadrado J. Gait analysis system for spinal cord injured subjects assited by active orthoses and crutches. J Multi-body Dyn. 2013;227(4):363–74.

Delp SL, et al. OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE Trans Biomed Eng. 2007;54(11):1940–50.

Crowninshield RD, Brand RA. A physiologically based criterion of muscle force prediction in locomotion. J Biomech. 1981;14(11):793–801.

Brand RA, Pedersen DR, Friederich JA. The sensitivity of muscle force predictions to changes in physiologic cross-sectional area. J Biomech. 1986;19(8):589–96.

Rasmussen J, Damsgaard M, Voigt M. Muscle recruitment by the min/max criterion — a comparative numerical study. J Biomech. 2001;34(3):409–15.

An KN, Kwak BM, Chao EY, Morrey BF. Determination of muscle and joint forces: a new technique to solve the indeterminate problem. J Biomech Eng. 1984;106(4):364–7.

Pipeleers G, et al. Dynamic simulation of human motion: Numerically efficient inclusion of muscle physiology by convex optimization. Optim Eng. 2008;9(3):213–38.

Thelen DG, Anderson FC, Delp SL. Generating dynamic simulations of movement using computed muscle control. J Biomech. 2003;36(3):321–8.

Winters JM. An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models. Ann Biomed Eng. 1995;23(4):359–74.

Scovil CY, Ronsky JL. Sensitivity of a Hill-based muscle model to perturbations in model parameters. J Biomech. 2006;39(11):2055–63.

Meyer AJ, Eskinazi I, Jackson JN, Rao AV, Patten C, Fregly BJ. Muscle synergies facilitate computational prediction of subject-specific walking motions. Front Bioeng Biotechnol. 2016;4:77.

Gopalakrishnan A, Modenese L, Phillips ATM. A novel computational framework for deducing muscle synergies from experimental joint moments. Front Computat Neurosci. 2014;8:153.

Michaud F, Shourijeh MS, Fregly BJ, Cuadrado J. Do muscle synergies improve optimization prediction of muscle activations during gait? Front Comput Neurosci. 2020;14(July):1–12.

Van Campen A, Pipeleers G, De Groote F, Jonkers I, De Schutter J. A new method for estimating subject-specific muscle-tendon parameters of the knee joint actuators: a simulation study. Int J Numer Method Biomed Eng. 2014;30(10):969–87.

Garner BA, Pandy MG. Estimation of musculotendon properties in the human upper limb. Ann Biomed Eng. 2003;31(2):207–20.

Shourijeh MS, Mehrabi N, McPhee J. Forward Static optimization in dynamic simulation of human musculoskeletal systems: a proof-of-concept study. J Comput Nonlinear Dyn. 2017;12(5):051005.

Del Vecchio A, Úbeda A, Sartori M, Azorín JM, Felici F, Farina D. Central nervous system modulates the neuromechanical delay in a broad range for the control of muscle force. J Appl Physiol. 2018;125(5):1404–10.

Hof AL. The relationship between electromyogram and muscle force. Sport. 1997;11(3):79–86.

Buchanan TS, Lloyd DG, Manal K, Besier TF. Neuromusculoskeletal modeling: Estimation of muscle forces and joint moments and movements from measurements of neural command. J Appl Biomech. 2004;20(4):367–95.

Umberger BR, Gerritsen KGM, Martin PE. A model of human muscle energy expenditure. Comput Methods Biomech Biomed Engin. 2003;6(2):99–111.

Bhargava LJ, Pandy MG, Anderson FC. A phenomenological model for estimating metabolic energy consumption in muscle contraction. J Biomech. 2004;37(1):81–8.

Duysens J, Beerepoot VP, Veltink PH, Weerdesteyn V, Smits-Engelsman BCM. Proprioceptive perturbations of stability during gait. Neurophysiol Clin. 2008;38(6):399–410.

Tresch MC, Jarc A. The case for and against muscle synergies. Curr Opin Neurobiol. 2009;19(6):601–7.

Slawinska U, et al. Comment on restoring voluntary control of locomotion after paralyzing spinal cord injury. Science. 2012;338(6105):328–328.

Rydell NW. Forces acting on the femoral head-prosthesis: a study on strain gauge supplied prostheses in living persons. Acta Orthop Scand. 1966;37(sup88):1–132.

Bergmann G, Graichen F, Rohlmann A. Hip joint loading during walking and running, measured in two patients. J Biomech. 1993;26(8):969–90.

Brand RA, Pedersen DR, Davy DT, Kotzar GM, Heiple KG, Goldberg VM. Comparison of hip force calculations and measurements in the same patient. J Arthroplasty. 1994;9(1):45–51.

## Acknowledgements

Authors would like to acknowledge the different funding supports. Moreover, they would like to acknowledge the subjects for their voluntary participation in this project.

## Funding

This work was funded by the Spanish MICIU under project PGC2018-095145-B-I00, co-financed by the EU through the EFRD program, and by the Galician Government under grant ED431C2019/29. Moreover, F. Michaud would like to acknowledge the support of the Spanish MINECO by means of the doctoral research contract BES-2016–076901, co-financed by the EU through the ESF program.

## Author information

### Authors and Affiliations

### Contributions

FM designed the experiments with the supervision of UL and performed the experiments with the help of ML, FM derived the models and analyzed the data. FM and JC wrote the manuscript in consultation with UL and ML. All authors read and approved the final manuscript.

### Authors' information

Florian Michaud graduated in Engineering (Mechanical and Ergonomy) at University of Belfort-Montbeliard (France) in April 2014. In June 2014 he joined the Laboratory of Mechanical Engineering of University of La Coruña, granted by the research funds of the Laboratory, to work on the biomechanics of human motion. Since May 2017 until January 2020 he enjoyed a doctorate grant from the Spanish Government. In 2018 he carried out a five-month pre-doctoral stay at Rice University in Houston, USA. In January 2020 he obtained his doctoral degree. Since then he continues working, now as a doctor, on applications of multibody dynamics to the biomechanics of human motion.

Mario Lamas graduated in Biomedical Engineering at University of Barcelona (Spain) in July 2017, and got his master degree in Biomedical Engineering, with an intensification in Biomechanics, at University of Zaragoza (Spain) in December 2018. In March 2019 he joined the Laboratory of Mechanical Engineering of University of La Coruña, granted by the research funds of the Laboratory. His work focuses on the analysis and simulation of human motion.

Urbano Lugris graduated in Mechanical Engineering from the University of La Coruña (Spain) in March 2003. Then, he joined the mechanical department of consulting company Norcontrol S.A., where he carried out technical assistance and quality control tasks for wind power plants and steel structures. In May 2004 he started his doctoral thesis at the Laboratory of Mechanical Engineering, granted by the Spanish Government, and got is Ph.D. in November 2008. His work has focused on real-time formulations for multibody dynamics with flexible bodies, spending three months at Universidad de Sevilla in 2005, and three more at University of Stuttgart in 2006. Since January 2011 he is Associate Professor, and his work is currently focused on biomechanics of human motion. In 2014 he spent three months at University Duisburg-Essen.

Javier Cuadrado graduated and obtained his PhD in Mechanical Engineering from University of Navarra in 1990 and 1993, respectively. In 1994, he joined the Mechanical Engineering department of University of La Coruña as Associate Professor, becoming Full Professor in 2000. In 2002, he founded the Laboratory of Mechanical Engineering (http://lim.ii.udc.es) at the mentioned department. His research is oriented to the computational kinematics and dynamics of multibody systems with rigid and flexible links, and their applications, mainly in the automotive, biomechanical, machinery and maritime industrial sectors, having published more than 50 papers in journals and 150 in conferences. He has participated in projects funded by the European Commission, the European Space Agency, and the US Army Research Office. He has supervised ten doctoral theses. He is or has been Associate Editor of three journals, Mechanism and Machine Theory, Journal of Computational and Nonlinear Dynamics, and Mechanics Based Design of Structures and Machines, and member of the Editorial Board of two other journals, Multibody System Dynamics, and Journal of Multi-body Dynamics. He chaired the Technical Committee for Multibody Dynamics of IFToMM (International Federation for the Promotion of Mechanism and Machine Science) since its creation in 2005 until 2013, and is secretary of IMSD (International Association for Multibody System Dynamics) since its creation in 2010.

### Corresponding author

## Ethics declarations

### Ethics approval and consent to participate

The studies involving human participants were reviewed and approved by the Committee of Ethics of the University of A Coruña. The participants provided their written informed consent to participate in this study.

### Consent for publication

The participants provided their written informed consent for publication.

### Competing interests

No conflicts of interest lie with any of the authors.

## Additional information

### Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

## About this article

### Cite this article

Michaud, F., Lamas, M., Lugrís, U. *et al.* A fair and EMG-validated comparison of recruitment criteria, musculotendon models and muscle coordination strategies, for the inverse-dynamics based optimization of muscle forces during gait.
*J NeuroEngineering Rehabil* **18**, 17 (2021). https://doi.org/10.1186/s12984-021-00806-6

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12984-021-00806-6