An EMG-marker tracking optimisation method for estimating muscle forces

Abstract

Existing algorithms for estimating muscle forces mainly use least-activation criteria, which do not necessarily lead to physiologically consistent results. Our objective was to assess an innovative forward dynamics-based optimisation, assisted by both electromyography (EMG) and marker tracking, for estimating the upper-limb muscle forces. A reference movement was generated, and EMG was simulated to reproduce the desired joint kinematics. Random noise was added to both simulated EMG and marker trajectories in order to create 30 trials. Then, muscle forces were estimated using (1) the innovative EMG-marker tracking forward optimisation, (2) a marker tracking forward optimisation with a least-excitation criterion, and (3) static optimisation with a least-activation criterion. Approaches (1) and (2) were solved using a direct multiple shooting algorithm. Finally, reference and estimated joint angles and muscle forces for the three optimisations were statistically compared using root-mean-square errors (RMSEs), biases, and statistical parametric mapping. The joint angles RMSEs were qualitatively similar across the three optimisations: (1) $1.63 \pm 0.51$°; (2) $2.02 \pm 0.64$°; (3) $0.79 \pm 0.38$°. However, the muscle forces RMSE for the EMG-marker tracking optimisation ($20.39 \pm 13.24$ N) was about seven times smaller than those resulting from the marker tracking ($124.22 \pm 118.22$ N) and static ($148.15 \pm 94.01$ N) optimisations. The originality of this novel approach is close tracking of both simulated EMG and marker trajectories in the same objective function, using forward dynamics. Therefore, the presented EMG-marker tracking optimisation led to accurate muscle forces estimations.

Appendices

Appendix A

1.1 A.1 The MUSCOD-II software

MUSCOD-II [23] solves optimal control problems based on the direct multiple shooting algorithm [22, 53]. The latter consists in dividing the integration interval into $N$ shorter subintervals, which facilitates and speeds up the convergence of the solution [24]. Additional matching constraints guarantee the continuity of the overall solution over the whole time interval. Inequality constraints are also applied, as, for instance, the ranges of joint angles ($\mathbf{q}$), velocities ($\dot{\mathbf{q}}$), muscle activations ($\mathbf{a}$) and excitations ($\mathbf{e}$):

$$ \begin{gathered} \mathbf{q}_{\mathrm{min}} \leq \mathbf{q} \leq \mathbf{q}_{\mathrm{max}},\\ \dot{\mathbf{q}}_{\mathrm{min}} \leq \dot{\mathbf{q}} \leq \dot{\mathbf{q}}_{\mathrm{max}},\\ 0 \leq \mathbf{a} \leq 1,\\ 0 \leq \mathbf{e} \leq 1. \end{gathered} $$

(9)

In the present study, MUSCOD-II was used with the 4th/5th ODE/DAE Runge–Kutta–Fehlberg solver module, which has a good accuracy level for a given time step [54].

1.2 A.2 Generation of reference muscle excitations

From an anatomical position, the simulated noise-free reference movement mainly consisted of an elbow flexion, hand palm facing upward. The desired joint angles and velocities were defined using the Yeadon quintic spline functions [55]. MUSCOD-II [23] was then used to obtain the reference muscle excitations that produced the desired joint kinematics. Control variables were the muscle excitations ($\mathbf{e}$), and state variables were the joint angles, velocities ($\mathbf{q}, \dot{\mathbf{q}}$) and muscle activations ($\mathbf{a}$). Controls and states variables were jointly optimised with respect to each optimisation objective function and the equation of dynamics, Eq. (4). No objective function was given while generating the optimal noise-free reference excitations with MUSCOD-II. The movement duration was fixed at 1 s. All aforementioned inequality constraints (Eq. (9)) were specified. Specifically, joint angles were forced to respect the desired kinematic values, given as an initial solution at each node of the problem.

Appendix B: Results for the low co-contraction movement

2.1 B.1 Marker and kinematics tracking

The EMG-marker tracking and marker tracking optimisations using MUSCOD-II converged in $25.5 \pm 5.3$ and $73.9 \pm 49.0$ min (mean ± standard deviation of $n = 30$ trials), respectively, for an average of 3.6 million calls of the forward-dynamic function (Intel^® Core™ i5-3570 CPU @3.4 GHz). Comparatively, static optimisation on MATLAB converged in $2.5 \pm 1.0$ min. The average residual actuator in static optimisation was $- 0.17 \pm 0.49~\mbox{N}{\cdot}\mbox{m}$, which is good.

Similarly to the high co-contraction movement, the tracking residuals of the marker trajectories had the same order of magnitude for the three optimisations (EMG-marker tracking: $0.23 \pm 0.10$ cm; marker tracking: $0.24 \pm 0.11$ cm; static optimisation: $0.17 \pm 0.06$ cm). Errors were larger for markers placed on the distal segments of the kinematic chain than for those placed on the proximal segments (Fig. 6).

Fig. 6

Tracking residual of the markers for the three optimisations, averaged across all the markers, across the length of the movement with low co-contraction and across the 30 trials. Note. The EMTO, MTO and SO acronyms stand for the EMG-marker tracking, marker tracking and static optimisations, respectively. (Color figure online)

The bias and RMSE values of the estimated joint angles were similar between the three optimisations (Table 5). The SPM ANOVA thus revealed no significant effect of the Optimisation method on the biases between the reference and the estimated joint angle lasting more than 0.2 s for any DOF (Fig. 7).

Fig. 7

Reference (REF-LC) and mean ± standard deviation of the estimated joint angles obtained with low co-contraction dataset. Grey-shaded zones represent the SPM ANOVA results, that is, the time intervals where there is a significant effect of Optimisation method. (Color figure online)

Table 5 RMSE and bias in degree between the noise-free reference and estimated DOFs of the model with low co-contraction. Note. GH: glenohumeral; SD: standard deviation

2.2 B.2 Muscle activations and forces

The time integral of the squared activations averaged across all the lines of action was 2.9 for the reference, $3.2 \pm 3.2$ for the EMG-marker tracking optimisation, $1.2 \pm 1.1$ for the marker tracking optimisation and $11.4 \pm 8.2$ for static optimisation. Concerning the muscle forces, the EMG-marker tracking RMSE averaged across all the lines of action was $7.61 \pm 4.83$ N with a bias of $2.2 \pm 3.6$ N, meaning a small overestimation (Table 6). RMSE for marker tracking ($34.71 \pm 29.44$ N) and static ($115.51 \pm 75.74$ N) optimisations presented a five- and sixteen-fold increase, respectively, with systematically negative biases for marker tracking optimisation (i.e. forces were underestimated for all muscles) and a positive average bias for static optimisation (Table 6). Muscle forces and activations in static optimisation showed the largest inter-trial variability (see the standard deviations of biases in Tables 6 and 7).

Table 6 RMSE and bias between the noise-free reference and estimated muscle forces with low co-contraction ($n = 30$ trials) for the 20 lines of action of the model. Note. A positive bias corresponds to an overestimation. An asterisk was added next to the muscle lines of action for which a significant effect of Optimisation method was observed on Fig. 8

Table 7 RMSE and bias between the noise-free reference and estimated muscle activations with low co-contraction ($n = 30$ trials) for the 20 lines of action of the model. Note. A positive bias corresponds to an overestimation

For INF, the SPM ANOVA revealed a significant effect of the Optimisation method on the biases between the reference and estimated muscle forces on more than 50% the movement (Fig. 8). For TRI lat., TRI med. and PEC clav., the significant Optimisation method effect was observed on less than 50% of the movement. No significant difference lasting more than 0.2 s was observed for the other muscles. For TRI lat. and TRI med., post hoc comparisons only assessed that the EMG-marker tracking biases were significantly different from the marker tracking ones and from static optimisation ones (i.e. marker tracking and static optimisations biases were never significantly different for these two muscles). For INF, post hoc comparisons only indicated that the marker tracking biases were significantly different from the EMG-marker tracking ones and from static optimisation ones. For PEC clav., post hoc comparisons showed that the EMG-marker tracking biases were significantly different from the marker tracking ones and from static optimisation ones and that the marker tracking and static optimisations biases were significantly different too.

Fig. 8

Reference (REF-LC) and mean ± standard deviation of the estimated muscle forces obtained with the low co-contraction dataset. Grey-shaded zones represent the SPM ANOVA results, that is, the time intervals where there is a significant effect of Optimisation method. (Color figure online)