EMG-Assisted Algorithm to Account for Shoulder Muscles Co-Contraction in Overhead Manual Handling

Glenohumeral stability is essential for a healthy function of the shoulder. It is ensured partly by the scapulohumeral muscular balance. Accordingly, modelling muscle interactions is a key factor in the understanding of occupational pathologies, and the development of ergonomic interventions. While static optimization is commonly used to estimate muscle activations, it tends to underestimate the role of shoulder’s antagonist muscles. The purpose of this study was to implement experimental electromyographic (EMG) data to predict muscle activations that could account for the stabilizing role of the shoulder muscles. Kinematics and EMG were recorded from 36 participants while lifting a box from hip to eye level. Muscle activations and glenohumeral joint reactions were estimated using an EMG-assisted algorithm and compared to those obtained using static optimization with a generic and calibrated model. Muscle activations predicted with the EMG-assisted method were generally larger. Additionally, more interactions between the different rotator cuff muscles, as well as between primer actuators and stabilizers, were predicted with the EMG-assisted method. Finally, glenohumeral forces calculated from a calibrated model remained within the boundaries of the glenoid stability cone. These findings suggest that EMG-assisted methods could account for scapulohumeral muscle co-contraction, and thus their contribution to the glenohumeral stability.


Introduction
The glenohumeral (GH) joint is the most mobile joint of the human body. This mobility is enabled by the minimal congruence between the glenoid cavity and the humeral head, which compromises the GH stability. This stability requirement is largely met by the scapulohumeral muscular balance [1]. This balance is ensured when the net joint reaction force is directed towards the glenoid fossa, which requires a fine-tuned coordination between shoulder muscles [2]. While rotator cuff muscles have been reported to act as the GH joint main stabilizers, all shoulder muscles contribute to the GH stability [3]. This intricate stability design puts the GH joint at a higher risk of dislocation, and its muscles at a higher injury hazard. As shoulder pathologies have been linked to numerous occupational settings (e.g., repetitive, high force demand, computer set up), and their impact on quality of life has been reported, it is critical to understand these dysfunctions mechanisms [4]. Accordingly, the prediction of shoulder muscle interactions is a key factor in the development of ergonomic interventions.

Experimental Protocol
The data previously collected [21] from 36 participants (age: 23.1 ± 2.9 years, mass: 68.8 ± 13.2 kg, height: 173.6 ± 9.9 cm, BMI: 22.66 ± 2.9, 17 females and 19 males) with no history of upper limb diseases was used for this study. The research protocol was approved by the University of Montreal Ethics Committee (n • 15-016-CERES-P) and all participants provided their informed consent prior to the experimentation. The experimental protocol fully described by Bouffard et al. [21] and Martinez et al. [22] is briefly summarized.
The participants performed an overhead box lifting task. We limited our study to the right upper limb, assuming that both upper limbs moved symmetrically with regard to the sagittal plane [23]. In line with the kinematic model developed by Jackson et al. [24], 43 reflective skin markers were placed on the thorax and the right upper limb. The kinematics were recorded using an 18 camera Vicon TM motion analysis system (Oxford Metrics Ltd., Oxford, UK). Muscle activity was recorded for ten muscles, with bipolar surface EMG electrodes positioned over the anterior, median and posterior deltoids, biceps, triceps, pectoralis major, latissimus dorsi, upper and lower trapezius and serratus anterior. Indwelling electrodes were inserted into the supraspinatus, infraspinatus and subscapularis. EMG signals were collected at 2000 Hz. They were processed using a second-order bandpass filter, with cut-off frequencies fixed at 20 and 425 Hz, then full-wave rectified. Envelopes were extracted using a 5 Hz low-pass filter. A series of 10 submaximal, and maximal, voluntary contractions (MVC) were performed to validate the electrodes position, and normalize the EMG envelopes, respectively [25]. The forces and moments between the hand and the box (8 × 35 × 50 cm) were measured using a six-dimensional handle force sensor described in previous studies [26]. The box (6 or 12 kg) was lifted from a shelf at hip level and positioned on another shelf at eye level. Each participant performed three trials per load in random order with 30 s rest periods.

Musculoskeletal Model
A modified OpenSim [13] generic model [19] of the shoulder was used for this study. The model had 10 degrees of freedom (DOF) with 29 musculo-tendon units (MTU) of which 17 actuated the GH joint. Some MTUs were represented by multiple pathways, e.g., anterior, middle, and posterior deltoids. First, an uncalibrated model was created by scaling the generic model to each participant's anthropometry using marker positions during a static pose. The joint generalized coordinates and moments were calculated using inverse kinematics and inverse dynamics, respectively, while the MTUs lengths and moment arms were estimated using muscle analysis in Opensim. The trial data were batch processed using the Pyosim and Pyomeca libraries in Python [27,28].
A personalized model was then created by calibrating the MTU parameters of the scaled uncalibrated model. The calibration step aimed to adjust MTU parameters to minimize the sum of the inverse dynamic joint moment tracking quadratic error and the ratio of GH joint shear to compressive contact forces using a simulated annealing algorithm (See Appendix A). To this avail, the experimental EMG was used as an input of a forward dynamics analysis to predict joint moments for each given set of MTU parameters. For this step, a calibration dataset was used. This dataset consisted of one trial per subject using a 6-kg load. A penalty was applied whenever an MTU's normalized fiber length reached non-physiological values outside the range of [0. 5, 1.5]. For this step, lines of action were grouped based on their function and anatomical classification. Each group received the same EMG signal. Accordingly, lines of action of the same muscle were calibrated using the same EMG signal. Teres minor was grouped with the infraspinatus, as they were expected to have similar functions [29]. Since the experimental data on muscles actuating the sternoclavicular and acromioclavicular joints was limited, we chose to focus this study on the GH joint. The DOFs of the calibrated model included only those of the GH joint, and its muscles were limited to the GH 17 lines of action. Out of these lines, only teres major and coracobrachialis were not calibrated. However, as they are not the main actuators of the GH joint, there influence on the calibration step results was expected to be minimal.
The muscle activations and forces were estimated using three different methods: For all methods, no constraints regarding GH dislocation were implemented while predicting muscle activations. Nevertheless, the calibration process for SOcal and EMGA implicitly tuned the MTU parameters to favor a smaller shear to compressive GH joint reaction ratio during the calibration step. SO and SOcal cost function was the sum of squared activations using the algorithm implemented in Opensim. For these two methods, residual actuators were implemented at all DOFs to ensure finding a solution [30]. The last method (EMGA) sought to adjust existing experimental excitations (i.e., EMG linear envelopes) when available and to synthetize the excitations for the remaining lines of action. Activations were predicted from the excitations using the EMG-to-activation model described by Lloyd and Besier [31]. The adjusted and synthetized excitations simultaneously minimized three terms [32]: • Inverse dynamics joint moments tracking error, • The sum of experimental excitations tracking error for all experimental excitations, • The sum of squared excitations for the 17 lines of actions.
These three terms were combined in a single cost function. The weighting coefficients for each term were chosen to reach a good balance between tracking joint moments and excitations, following the method explained by Sartori et al. [32].

Analysis
The predicted muscle forces from the three models were then used to calculate joint reaction forces at the GH joint, as well as the ratio of shear to compressive joint contact force. The root mean square difference (RMSd) between the inverse dynamics joint moments and those predicted by the EMG-assisted method as well as the Pearson correlation coefficient (R) were calculated to validate its application. The residual actuator moments obtained from SO and SOcal were used for comparison. The activations predicted from SO, SOcal and EMGA were compared with experimental EMGs using RMSd and R. Their respective empirical cumulative distribution functions [33] were also compared using a Kolmogorov-Smirnov dissimilarity index [34]. This index evaluates the probability that the compared samples are drawn from the same distribution. The larger its statistic, the less likely the samples have similar distributions. Coactivation patterns of the rotator cuff muscles between both loads were compared with a Jennrich test on the correlation matrices at each load [35]. Furthermore, the co-contraction was evaluated for muscle pairs known to have synergies: latissimus dorsi and teres major, anterior deltoid and coracobrachialis, median deltoid and supraspinatus, posterior deltoid and infraspinatus [36], using a co-contraction coefficient defined in Equation (1), where E 1 and E 2 are the activation of the less and more activated muscle at each time frame for each muscle pair, respectively [37].
Finally, the GH joint reaction force was expressed in the glenoid local frame. Within this frame, the GH joint reaction ratio, defined as the ratio of the shear to compressive contact forces in the glenoid, was assessed in regard to the reported non-dislocation thresholds [1].

Net Moment Tracking
EMGA tracked moments with an RMSd within the range of the residual actuators used by SO ( Figure 1, Table 1). The calibrated model succeeded in tracking the net moments calculated by inverse dynamics. The two peaks in the residuals actuators when solving SO (at 24% and 42% of the trial) pointed to a possible failure due to the incapacity of MTUs to generate the needed joint moment. These peaks disappeared when using SOcal.

EMG and Muscle Forces
The adjusted excitations mostly had strong correlations (0.30-0.90) with the experimental EMG linear envelopes, with some discrepancies for the posterior deltoid, latissimus dorsi, pectoralis major and the triceps ( Figure 2, Table 2). In comparison, the correlation between the SO predicted activations and the experimental EMG were lower (0.05-0.46). Accordingly, different muscle activation patterns were obtained from SO, SOcal, and the EMGA algorithms ( Figure 3). Compared to SO, SOcal tended to minimize the activation of the rotator cuff muscles, with even larger differences to EMGA results. In general, the activations were smaller and showed less inter-subject variability for SO and SOcal, compared to EMGA (Figures 3 and 4). The dissimilarity index between the experimental EMG distribution and the activation estimated using EMGA and SO was of 0.219 and 0.492 respectively (p < 0.001). The correlation between the different muscle activations followed a similar pattern when using EMGA algorithm, irrespective of the load ( Figure 5). The correlation between the experimental EMG of the infraspinatus, supraspinatus and subscapularis did not differ with the mass (X 2 = 1.94, p > 0.05). Similar results were obtained with the activations predicted with EMGA (X 2 = 3.85, p > 0.05). However, the activation correlations predicted by SO had a significant relationship with the mass (X 2 = 33.98, p < 0.01). The main actuators and the rotator cuff muscles had a higher co-activation for EMGA in comparison with SO ( Figure 6). and the triceps ( Figure 2, Table 2). In comparison, the correlation between the SO predicted activations and the experimental EMG were lower (0.05-0.46). Accordingly, different muscle activation patterns were obtained from SO, SOcal, and the EMGA algorithms ( Figure 3). Compared to SO, SOcal tended to minimize the activation of the rotator cuff muscles, with even larger differences to EMGA results. In general, the activations were smaller and showed less inter-subject variability for SO and SOcal, compared to EMGA (Figures 3 and 4). The dissimilarity index between the experimental EMG distribution and the activation estimated using EMGA and SO was of 0.219 and 0.492 respectively (p < 0.001). The correlation between the different muscle activations followed a similar pattern when using EMGA algorithm, irrespective of the load ( Figure 5). The correlation between the experimental EMG of the infraspinatus, supraspinatus and subscapularis did not differ with the mass (Χ² = 1.94, p > 0.05). Similar results were obtained with the activations predicted with EMGA (Χ² = 3.85, p > 0.05). However, the activation correlations predicted by SO had a significant relationship with the mass (Χ² = 33.98, p < 0.01). The main actuators and the rotator cuff muscles had a higher co-activation for EMGA in comparison with SO ( Figure 6).

Joint Reaction Forces
The magnitude of the estimated GH force was higher with EMGA compared to SO and SOcal. However, EMGA solutions did not present any sudden peaks in the force amplitude for the duration of the test, unlike SO. These SO peaks were smoothed out when using SOcal (Figure 7a). Peaks also appeared for the joint reaction ratio throughout the test with SO (Figure 7b). A maximum ratio was reached close to the maximum extension of the GH joint with a mean value of 0.65 and 0.27 for the EMGA and SOcal, respectively. In comparison, the joint reaction ratio had an erratic behavior when estimated from SO. The direction of the GH force was directed within the stability cone for 100% of the trials for the SOcal and EMGA models, but not for the SO model (Figure 7b).          Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 16 Figure 6. Coactivation coefficients between various muscle pairs (solid line: mean value, shaded region: 95% confidence interval) as predicted by EMGA (blue), SO (orange) and SOcal (green): the higher the coefficient, the higher is the muscle pair's coactivation.

Joint Reaction Forces
The magnitude of the estimated GH force was higher with EMGA compared to SO and SOcal. However, EMGA solutions did not present any sudden peaks in the force amplitude for the duration of the test, unlike SO. These SO peaks were smoothed out when using SOcal (Figure 7a). Peaks also appeared for the joint reaction ratio throughout the test with SO (Figure 7b). A maximum ratio was reached close to the maximum extension of the GH joint with a mean value of 0.65 and 0.27 for the EMGA and SOcal, respectively. In comparison, the joint reaction ratio had an erratic behavior when estimated from SO. The direction of the GH force was directed within the stability cone for 100% of the trials for the SOcal and EMGA models, but not for the SO model (Figure 7b).

Discussion
This study compared muscle activation patterns and GH joint reaction forces obtained from three different neuromusculoskeletal modelling approaches: SO, SOcal, and EMGA. As hypothesized, shear to compressive GH joint reaction forces predicted with EMGA remained within the stability cone boundaries, unlike those predicted with SO. Similar results were obtained with SOcal which pointed out the importance of musculoskeletal calibration.  [1] in the glenoid frame (mean ± 1 std: solid ± dashed black lines): If the GH joint reaction shear/compressive ratio is within the cone, the joint is considered as stable. The data shown are from the EMGA (blue), SO (orange) and SOcal (green, solid line: mean value, shaded region: 95% confidence interval).

Discussion
This study compared muscle activation patterns and GH joint reaction forces obtained from three different neuromusculoskeletal modelling approaches: SO, SOcal, and EMGA. As hypothesized, shear to compressive GH joint reaction forces predicted with EMGA remained within the stability cone boundaries, unlike those predicted with SO. Similar results were obtained with SOcal which pointed out the importance of musculoskeletal calibration.

Net Moment Tracking and Calibration
Calibration is a key step to tune MTU parameters that vary non-linearly across subjects, and thus is critical to EMG-informed models [38]. In this study, the difference in residual moments between SO and SOcal pointed out the impact of calibration. For SO, the presence of peaks in the residual moment pointed the incapacity of the model to generate the required inverse dynamics moment. This could be the result of the cost function used for static optimization. Indeed, for a similar task, the distribution of muscle activations revealed a weak activation of the majority of the muscles (<20% MVC), while the task was mainly led by the superior trapezius, the anterior and median deltoids, and the infraspinatus [39]. Additionally, the forces were mainly generated by the previously mentioned muscles as well as the latissimus dorsi. Another possibility is that the model required such an activation pattern due to its overall weakness. In other words, the model was unable to generate the moments responsible for the upper limb movement relying solely on its MTUs activations. With the calibration, the maximal isometric force was increased for all muscles. Thus, it is possible that the generic muscle parameters were not representative of our participants, and thus would require further calibration. However, the choice of the calibration cost function is also crucial. Limiting the calibration objective function to tracking the inverse dynamics moments only, might lead to a non-physiological increase of the MTU's maximal isometric force [14,40].

EMG and Muscle Forces
The difference between adjusted and experimental EMGs were comparable to the differences between the EMG envelopes obtained from two repetitions of the 6 kg trials (See Appendix B). Accordingly, the adjusted EMGs were considered representative of the experimental EMGs. The muscles that exhibited the largest differences were the posterior deltoid, pectoralis major, latissimus dorsi, and triceps. In the second half of the lift, the triceps was mostly involved in the elbow flexion, thus the increased length of the MTU was balanced by a reduction of its EMG, to reduce its contribution to the net moments of the GH joint. Pectoralis major and latissimus dorsi generated torques at the GH, sternoclavicular and acromioclavicular joints. Since the calibration was performed by matching only GH joint torque, the error in calibration for these muscles was probably larger. Similar results have been reported for the lower limb. The activations predicted with a multi-joint (hip, knee, ankle) model were more accurate than using a single joint model [16]. Overall, the agreement of the EMGs offered a preliminary validation of the activations predicted with EMGA.
As shown in previous studies on the upper limb, lower limb, or the spine, activation patterns predicted by EMG-informed algorithms differ from SO based predictions [15,41]. SO, based on least squares activation, tends to favor muscles with large moment arms [42] and underestimates co-contraction [43]. While this could be easily observed for simple humeral motions such as abduction or flexion, it was less straightforward for our lifting task due to its kinematic complexity. Overall, SO reduced activation for most muscles, in comparison to activations predicted by the EMGA method. This tendency was even more evident on the overhead phase, where the supraspinatus and subscapularis activations were minimized with SO and SOcal, even though such positions are known to be at higher risk for injury, and thus would require more activation of the rotator cuff muscles [44][45][46]. Similarly, during the concentric phase of the task (0%-30%), teres major was expected to assist the latissimus dorsi in the extension of the humerus [47], which was not predicted when using SO or SOcal. Accordingly, the difference in the empiric activation distribution was probably due to the lack of a co-contraction component in SO and SOcal. This limitation of SO could also be inferred from the random correlations between the rotator cuff muscles that were not consistent between trials. Whereas, the activations predicted by the EMGA method showcased correlations between the rotator cuff muscles irrespective of the task and load, as well as between agonist and antagonist muscles. Additionally, the observed loop-shape co-activation was reported experimentally as a typical force sharing pattern in cat walking [48]. Muscle activations were obtained without the implementation of a non-dislocation constraint, which reduced the effect of mathematical penalties and restrictions that might not be physiological. Indeed, different muscle recruitment trends could be predicted depending on the implementation of a stability constraint [20,26]. The implications of such a constraint could be detrimental when studying unstable shoulders, that might physiologically fail to ensure stability, but would succeed to do so in a constrained model.

Joint Reaction Forces
The pattern of the GH joint reaction forces obtained from EMGA were similar to that reported in-vivo using an instrumented prosthesis [49,50] or via modelling [7], with higher amplitudes close to the maximum elevation and a higher contribution of the compressive component of the force. Our study however reported a higher force magnitude (mean peak value: 2907 N vs. about 1000 N). This could partly be explained by the difference in the lifted weights (6 vs. 2 kg), the age and fit of the studied population, and the kinematic differences between the two studies [51]. Additionally, the use of an implant, for the in-vivo study might have shifted the GH center of rotation and muscles moment arms [50]. The GH force calculated by SO had random peaks throughout the trials. However, these peaks disappeared with the calibration (SOcal). They occurred when GH main actuators (i.e., deltoid and pectoralis major) had activation peaks. Since these main actuators have a large contribution to the net joint reaction force, their activation patterns would reflect on the GH joint reaction amplitude. The model calibration, as discussed in Section 1, enabled the MTUs to generate the inverse dynamics joint moment without saturating these main actuators, and thus eliminated these peaks from the GH joint reaction force with SOcal.
The muscle activations obtained from SO led to an unstable glenohumeral joint, as evidenced by the excursions outside the non-dislocation stability cone. The GH stability is probably attained, in part, by the co-contraction of rotator cuff muscles [52]. However, these muscles tend to have a smaller moment generating capacity. Thus, SO solutions would not favor their activation in comparison to that of the main GH actuators. Thereby, SO might not generate an activation that would prevent the humerus from sliding under the deltoid elevating action. From a stability perspective, the SO shoulder model was the most vulnerable at the end of the concentric phase along the superior-inferior direction, which was not found in our previous study based on SO [26]. The MTUs definition and the analyzed task are relatively similar between both studies. However, Blache et al. [26] optimized the maximal isometric forces to track the experimental EMG (genetic algorithm combined with SO), which partly resembles our calibration step. Accordingly, implementing SOcal, the GH instability decreased, despite the use of static optimization to predict activations. This emphasizes the importance of a subject-specific MTU model for studying complex systems such as the shoulder.
The GH joint reaction forces predicted with the EMGA method were on the boundary of the non-dislocation cone obtained from the ratios reported on cadaveric data when a compressive force of 50 N was applied [1]. The ratios are expected to decrease as the compressive force is increasing [53]. However, as far as we know, no data are reported in the literature for our complete range of compressive forces (366-2813 N). Additionally, these ratios were obtained when the humerus was positioned in the center of the glenoid, with the arm in a neutral position only, whereas these ratios have been reported to be a function of the joint's 3-DOF [53]. Finally, our model did not include the biceps long tendon. This tendon was reported to oppose the luxation of the humerus as the biceps is contracting [36], thus this contribution might account for the increased risk of dislocation observed in our study.

Limitations
This study had limitations that should be considered. The model calibration has been done with a task similar to those analyzed. Since this task required low sternoclavicular and acromioclavicular joint net moment, our study was limited to the GH joint. For future studies, a more appropriate calibration trial should be used, to account for the contribution of the trapezius and serratus anterior, as they are critical in rotator cuff pathologies [54]. Furthermore, the need of indwelling electrodes to calibrate the model is a practical limitation of the EMG-informed approaches for the upper limb. Further studies should be done to understand the sensitivity of each parameter of the calibration, and thus refine its process and reduce the number of required EMG signals, whether they require surface or indwelling electrodes. Additionally, while the symmetrical activity is confirmed in the literature for manual handling tasks [23], it would be interesting to compare calibration results as well as activation patterns between the dominant and non-dominant sides. The musculoskeletal model had some limitations. The first limitation was related to the geometry, as it was difficult to validate the muscle lengths and moment arms for all 36 participants due to large variability in their morphological features. Secondly, the GH joint stability was ensured only by MTUs, which neglected the effect of the labrum and ligaments [55]. However, it was expected that, in our range of motion labrum and ligaments, contributions were much smaller than muscles contributions [1,56]. Finally, the GH joint reaction magnitude in this study was relatively larger, it is expected that an improved calibration step would yield a joint reaction magnitude closer to those reported in the literature.

Conclusions
In conclusion, the activations predicted by EMGA were more representative of rotator cuff muscles co-contractions, as well as co-contractions between the main GH actuators and stabilizers. The calibration of the model enabled the prediction of a shear to compressive GH joint reaction force that remained within the boundaries of the non-dislocation cone, without the implementation of a non-dislocation constraint. These findings highlight the importance of subject-specific models and EMG-informed algorithms in the study of the shoulder complex in general, and more particularly in the prediction of the scapulohumeral muscular balance. A deeper understanding of the scapulohumeral muscular pathomecanisms should help the development of adequate ergonomic interventions in the work environment. The EMGA method could also be implemented for other populations that are at high risk of upper limb injury such as athletes and manual wheelchair users.  Acknowledgments: This research was undertaken thanks, in part, to funding from the Canada First Research Excellence Fund through the TransMedTech Institute.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The MTU parameters relative to the activation dynamics (EMG-to-activation) and contraction dynamics (activation-to-force) were personalized using a simulated annealing algorithm to minimize two components. The first one expressed the error between joint moments predicted by inverse dynamics and those calculated by the EMG-driven algorithm. The second component minimized the ratio of shear to compressive joint reaction forces at the GH joint. The objective function was defined as: with F τ defined in Pizzolato et al. [12] as: where the error E t,d is defined as: where N points is the number of data points, M t,d,n is the inverse dynamics moment, and M t,d,n the moment predicted from the EMG-data for the d th DOF for the t th trial at the n th timeframe. The p(n, j) function is a penalty to prevent non-physiological normalized fiber lengths ( l m n,j ) that are outside the range of [0.5, 1.5].
Appendix B Table A1. RMSd and R between the experimental EMG of two different trials with a 6 kg load (mean ± 1 std for n = 36 participants).