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

Abstract

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.

1. 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.

While musculoskeletal modelling enables a non-invasive prediction of these forces, it is hindered by the difficulties related to the redundancy of the shoulder musculoskeletal system. The most common method to estimate shoulder muscle activations is the static optimization of various cost functions, with or without a stability constraint for the GH joint [5,6,7,8], such as the sum of activations or the energic cost, which is expected to provide more physiological results [9]. However, the inaccuracies associated with static optimization, mainly due to its assumption of a unique optimal control strategy across individuals [10], led many investigators to consider using electromyographic (EMG) data to inform computations, rather than just for validation [11]. EMG-informed modelling methods have been implemented successfully within the calibrated EMG-informed neuromusculoskeletal modelling toolbox (CEINMS [12]) for OpenSim [13]. These methods seek to track simultaneously the experimental EMG and the various joint moments, in order to predict muscle activation. Additionally, they include a model calibration step, where, the experimental EMG is used as input, and the joint moments are tracked to fine-tune the muscle modelling parameters. CEINMS has been implemented for lower limb joints [12,14], where the co-contraction predicted by this toolbox was consistently higher than the one predicted by static optimization. Additionally, the predicted knee joint reaction forces were similar to experimental ones. This toolbox has also been implemented in its EMG-driven mode for shoulder analysis during isometric tasks [15], where it was concluded that a common cost function underestimates muscle co-contraction in comparison to an EMG-driven approach. However, the EMG-driven mode had already proven its limitations in the analysis of the hip joint due to the insufficient experimental data [12,16]. Additionally, the use of this mode did not account for artefacts related to the envelop extraction [16,17]. Indeed, the implementation of EMG-informed algorithms faces a number of limitations, that are particularly true for the shoulder. The first one being the presence of deep musculature [18] that limits the number of available experimental EMG. The second limitation is related to the quality of the EMG data, particularly because of artefacts introduced with dynamic tasks and envelope extraction. To overcome the first limitation, indwelling electrodes should be used where possible, and unattainable EMG data should be synthetized. As for the second limitation, the input data should be adjusted to account for joint moment tracking errors.

Accordingly, we implemented a shoulder model using OpenSim [13,19] and the EMG-adjusted algorithm of CEINMS [12] for the analysis of a relatively complex dynamic task. Our objective is to predict muscle activation patterns that could account for the stabilizing role of shoulder muscles. It was hypothesized that physiological muscle forces should sum up to a GH joint reaction force that remains within the non-dislocation cone defined by Lippitt and Matsen [1] and Dickerson et al. [20]. The predicted muscle and joint reaction forces were compared with those obtained from a similar model using static optimization. A secondary objective of the study was to evaluate the influence of the calibration on the results by implementing static optimization with a generic and a calibrated model.

2. Materials and Methods 2.1. 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.

2.2. 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:

• Static optimization on the uncalibrated model (SO),
• Static optimization on the calibrated model (SOcal),
• EMG-assisted method on the calibrated model in CEINMS (EMGA).

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].

2.3. 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₁ and E₂ are the activation of the less and more activated muscle at each time frame for each muscle pair, respectively [37].

CC=_E1E2(_E1+_E2)

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].

3. Results 3.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.

3.2. 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 (Figure 3 and Figure 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).

3.3. 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).

4. 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. 4.1. 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].

4.2. 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.

4.3. 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.

4.4. 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.

5. 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.

Figure 1. Left: Net moments of the GH joint calculated by inverse dynamics (ID in blue) and predicted with the EMG-assisted method (EMGA orange), with a schematic of various position of the participant during the box lifting task, right: residual actuators from the generic model (R in blue), calibrated model (R cal in orange), and the difference between the inverse dynamics and the net moments predicted with the EMG-assisted method (ΔEMGAin green). (solid line: mean value, shaded region: 95% confidence interval).

Figure 2. Net Adjusted EMG (EMGA: blue) and experimental EMG (Exp: red). * No experimental data (solid line: mean value, shaded region: 95% confidence interval).

Figure 3. Muscle activations (solid line: mean value, shaded region: 95% confidence interval) predicted by EMGA (blue), SO (orange) and SOcal (green) * Teres major activation was synthesized.

Figure 4. Empirical cumulative distribution of all muscle activations (solid line: mean value, shaded region: 95% confidence interval) predicted by EMGA (blue), SO (orange) and SOcal (green): for a given activation (a), the ECDF shows the percentage of data who is at most equal to a.

Figure 5. Activation relationships between the rotator cuff muscles as predicted for the 6 and 12 kg trials: The axes show the activation value of the respective muscles. On the diagonal: the mean distribution of a given muscle activation over the trial. In the upper (lower) triangle: the relationship between two muscles activations throughout the trial. For instance, the EMG algorithm predicted that an increase in the activation of the infraspinatus is correlated with an increase of the supraspinatus's activation. (Mean as predicted by EMGA: right, and SO: left).

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.

Figure 7.(a) GH joint reaction forces (solid line: mean value, shaded region: 95% confidence interval) calculated from the muscle activations predicted from the EMGA (blue), SO (orange) and SOcal (green); (b) Top: GH joint reaction shear/compressive ratio, bottom: ratio of the GH force in regard to the non-dislocation threshold ratios from Lippitt and Matsen [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).

Mass (kg)	GH DOF	RMSd (N.m)	Residual Actuators (N.m)	R
6	Plane of elevation	0.61 ± 0.70	1.29 ± 4.74	0.99 ± 0.03
	Elevation	0.81 ± 1.03	1.23 ± 2.19	0.99 ± 0.02
	Axial rotation	0.74 ± 0.67	1.60 ± 4.42	0.96 ± 0.06
12	Plane of elevation	1.69 ± 7.99	1.00 ± 2.70	0.98 ± 0.09
	Elevation	2.10 ± 10.77	1.54 ± 3.86	0.99 ± 0.06
	Axial rotation	1.24 ± 4.03	1.17 ± 2.26	0.97 ± 0.10

Muscles	RMSd			R
	SO	SOcal	EMGA	SO	SOcal	EMGA
Anterior Deltoid	0.34 ± 0.11	0.27 ± 0.10	0.11 ± 0.07	0.46 ± 0.31	0.53 ± 0.29	0.90 ± 0.09
Median Deltoid	0.25 ± 0.13	0.21 ± 0.07	0.08 ± 0.04	0.37 ± 0.42	0.22 ± 0.39	0.88 ± 0.16
Posterior Deltoid	0.13 ± 0.07	0.13 ± 0.09	0.08 ± 0.05	−0.17 ± 0.26	−0.08 ± 0.30	0.31 ± 0.51
Pectoralis Major	0.09 ± 0.09	0.10 ± 0.09	0.11 ± 0.08	0.17 ± 0.37	0.11 ± 0.36	0.34 ± 0.41
Latissimus Dorsi	0.16 ± 0.16	0.11 ± 0.09	0.12 ± 0.08	−0.13 ± 0.37	−0.19 ± 0.35	0.49 ± 0.44
Infraspinatus	0.32 ± 0.15	0.25 ± 0.13	0.06 ± 0.04	0.32 ± 0.30	0.28 ± 0.29	0.85 ± 0.17
Supraspinatus	0.31 ± 0.15	0.28 ± 0.17	0.06 ± 0.05	−0.05 ± 0.31	−0.13 ± 0.33	0.76 ± 0.30
Subscapularis	0.24 ± 0.20	0.21 ± 0.20	0.06 ± 0.04	0.06 ± 0.34	0.01 ± 0.37	0.56 ± 0.29
Biceps	0.18 ± 0.08	0.14 ± 0.06	0.05 ± 0.03	−0.36 ± 0.30	−0.28 ± 0.33	0.76 ± 0.24
Triceps	0.11 ± 0.07	0.11 ± 0.07	0.07 ± 0.06	−0.13 ± 0.26	−0.04 ± 0.33	0.30 ± 0.57

Author Contributions

N.A., C.P., R.M. and D.G.L. developed the programs used for data analysis. N.A. and R.M. conducted data analysis. N.A. wrote the original paper draft. C.P., R.M., D.G.L. and M.B. reviewed the original draft. M.B. supervised the study. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the NSERC Discovery Grant (RGPIN-2014-03912).

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:

_Fcal=_Fτ+|_{FGHshearFGHcomp}|

with_Fτ defined in Pizzolato et al. [12] as:

_Fτ=∑t_Ntrials∑d_NDOFEt,d

where the error_Et,dis defined as:

_Et,d=1_Npoints∑n_Npoints(^{(_Mt,d,n¯-_Mt,d,n)2}var(_Mt,d¯)+∑j_NMTUp(n,j))

p(n,j)={100^{(_l˜mn,j-1)2},if|_l˜mn,j-1|>0.50,otherwise

where_Npointsis the number of data points,_Mt,d,n¯is the inverse dynamics moment, and_Mt,d,n the moment predicted from the EMG-data for the dth DOF for the tth trial at the nth timeframe. Thep(n,j)function is a penalty to prevent non-physiological normalized fiber lengths (_l˜mn,j) that are outside the range of [0.5, 1.5].

Appendix B

Table

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).

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).

Muscles	RMSd	R
Anterior Deltoid	0.11 ± 0.05	0.74 ± 0.17
Median Deltoid	0.08 ± 0.03	0.84 ± 0.10
Posterior Deltoid	0.05 ± 0.03	0.77 ± 0.19
Pectoralis Major	0.02 ± 0.03	0.58 ± 0.20
Latissimus Dorsi	0.02 ± 0.03	0.63 ± 0.25
Infraspinatus	0.05 ± 0.05	0.64 ± 0.38
Supraspinatus	0.05 ± 0.03	0.57 ± 0.41
Subscapularis	0.03 ± 0.03	0.30 ± 0.39
Biceps	0.07 ± 0.03	0.71 ± 0.22
Triceps	0.04 ± 0.04	0.82 ± 0.11

References

1. Lippitt, S.; Matsen, F. Mechanisms of glenohumeral joint stability. Clin. Orthop. Relat. Res. 1993, 291, 20-28.

2. Veeger, H.E.J.; van der Helm, F.C.T. Shoulder function: The perfect compromise between mobility and stability. J. Biomech. 2007, 40, 2119-2129.

3. Ackland, D.C.; Pandy, M.G. Lines of action and stabilizing potential of the shoulder musculature. J. Anat. 2009, 215, 184-197.

4. Pribicevic, M. The Epidemiology of Shoulder Pain: A Narrative Review of the Literature. In Pain in Perspective; InTech: London, UK, 2012.

5. Ambrósio, J.; Quental, C.; Pilarczyk, B.; Folgado, J.; Monteiro, J. Multibody biomechanical models of the upper limb. Procedia IUTAM 2011, 2, 4-17.

6. Engelhardt, C.; Camine, V.M.; Ingram, D.; Müllhaupt, P.; Farron, A.; Pioletti, D.; Terrier, A. Comparison of an EMG-based and a stress-based method to predict shoulder muscle forces. Comput. Methods Biomech. Biomed. Eng. 2015, 18, 1272-1279.

7. Klemt, C.; Prinold, J.A.; Morgans, S.; Smith, S.H.L.; Nolte, D.; Reilly, P.; Bull, A.M.J. Analysis of shoulder compressive and shear forces during functional activities of daily life. Clin. Biomech. 2018, 54, 34-41.

8. Quental, C.; Folgado, J.; Ambrósio, J. A window moving inverse dynamics optimization for biomechanics of motion. Multibody Syst. Dyn. 2016, 38, 157-171.

9. Praagman, M.; Chadwick, E.K.J.; van der Helm, F.C.T.; Veeger, H.E.J. The relationship between two different mechanical cost functions and muscle oxygen consumption. J. Biomech. 2006, 39, 758-765.

10. Buchanan, T.S.; Shreeve, D.A. An Evaluation of Optimization Techniques for the Prediction of Muscle Activation Patterns During Isometric Tasks. J. Biomech. Eng. 1996, 118, 565-574.

11. Erdemir, A.; McLean, S.; Herzog, W.; van den Bogert, A.J. Model-based estimation of muscle forces exerted during movements. Clin. Biomech. 2007, 22, 131-154.

12. Pizzolato, C.; Lloyd, D.G.; Sartori, M.; Ceseracciu, E.; Besier, T.F.; Fregly, B.J.; Reggiani, M. CEINMS: A toolbox to investigate the influence of different neural control solutions on the prediction of muscle excitation and joint moments during dynamic motor tasks. J. Biomech. 2015, 48, 3929-3936.

13. Delp, S.L.; Anderson, F.C.; Arnold, A.S.; Loan, P.; Habib, A.; John, C.T.; Guendelman, E.; Thelen, D.G. OpenSim: Open-Source Software to Create and Analyze Dynamic Simulations of Movement. IEEE Trans. Biomed. Eng. 2007, 54, 1940-1950.

14. Gerus, P.; Sartori, M.; Besier, T.F.; Fregly, B.J.; Delp, S.L.; Banks, S.A.; Pandy, M.G.; D'Lima, D.D.; Lloyd, D.G. Subject-specific knee joint geometry improves predictions of medial tibiofemoral contact forces. J. Biomech. 2013, 46, 2778-2786.

15. Kian, A.; Pizzolato, C.; Halaki, M.; Ginn, K.; Lloyd, D.; Reed, D.; Ackland, D. Static optimization underestimates antagonist muscle activity at the glenohumeral joint: A musculoskeletal modeling study. J. Biomech. 2019, 109348.

16. Sartori, M.; Reggiani, M.; Farina, D.; Lloyd, D.G. EMG-Driven Forward-Dynamic Estimation of Muscle Force and Joint Moment about Multiple Degrees of Freedom in the Human Lower Extremity. PLoS ONE 2012, 7, e52618.

17. Farina, D.; Negro, F. Accessing the Neural Drive to Muscle and Translation to Neurorehabilitation Technologies. IEEE Rev. Biomed. Eng. 2012, 5, 3-14.

18. Nikooyan, A.A.; Veeger, H.E.J.; Westerhoff, P.; Bolsterlee, B.; Graichen, F.; Bergmann, G.; van der Helm, F.C.T. An EMG-driven musculoskeletal model of the shoulder. Hum. Mov. Sci. 2012, 31, 429-447.

19. Wu, W.; Lee, P.V.S.; Bryant, A.L.; Galea, M.; Ackland, D.C. Subject-specific musculoskeletal modeling in the evaluation of shoulder muscle and joint function. J. Biomech. 2016, 49, 3626-3634.

20. Dickerson, C.R.; Chaffin, D.B.; Hughes, R.E. A mathematical musculoskeletal shoulder model for proactive ergonomic analysis. Comput. Methods Biomech. Biomed. Eng. 2007, 10, 389-400.

21. Bouffard, J.; Martinez, R.; Plamondon, A.; Côté, J.N.; Begon, M. Sex differences in glenohumeral muscle activation and coactivation during a box lifting task. Ergonomics 2019, 1-12.

22. Martinez, R.; Bouffard, J.; Michaud, B.; Plamondon, A.; Côté, J.N.; Begon, M. Sex differences in upper limb 3D joint contributions during a lifting task. Ergonomics 2019, 62, 682-693.

23. Nielsen, P.K.; Andersen, L.; Jørgensen, K. The muscular load on the lower back and shoulders due to lifting at different lifting heights and frequencies. Appl. Ergon. 1998, 29, 445-450.

24. Jackson, M.; Michaud, B.; Tétreault, P.; Begon, M. Improvements in measuring shoulder joint kinematics. J. Biomech. 2012, 45, 2180-2183.

25. Dal Maso, F.; Marion, P.; Begon, M. Optimal Combinations of Isometric Normalization Tests for the Production of Maximum Voluntary Activation of the Shoulder Muscles. Arch. Phys. Med. Rehabil. 2016, 97, 1542-1551.e2.

26. Blache, Y.; Begon, M.; Michaud, B.; Desmoulins, L.; Allard, P.; Dal Maso, F. Muscle function in glenohumeral joint stability during lifting task. PLoS ONE 2017, 12, e0189406.

27. Martinez, R. Pyomeca/Pyosim: Interface between OpenSim and the Pyomeca Library. Available online: https://github.com/pyomeca/pyosim (accessed on 4 December 2018).

28. Martinez, R.; Michaud, B.; Begon, M. Pyomeca: Biomechanics toolbox for python. J. Open-Source Softw. 2020. (submitted for publication).

29. Inman, V.T.; Saunders, J.B.; Abbott, L.C. Observations of the Function of the Shoulder Joint. Clin. Orthop. Relat. Res. 1996, 330, 3-12.

30. Hicks, J.L.; Uchida, T.K.; Seth, A.; Rajagopal, A.; Delp, S.L. Is My Model Good Enough? Best Practices for Verification and Validation of Musculoskeletal Models and Simulations of Movement. J. Biomech. Eng. 2015, 137, 020905.

31. Lloyd, D.G.; Besier, T.F. An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo. J. Biomech. 2003, 36, 765-776.

32. Sartori, M.; Farina, D.; Lloyd, D.G. Hybrid neuromusculoskeletal modeling to best track joint moments using a balance between muscle excitations derived from electromyograms and optimization. J. Biomech. 2014, 47, 3613-3621.

33. Kvam, P.H.; Vidakovic, B. Nonparametric Statistics with Applications to Science and Engineering; Wiley Series in Probability and Statistics; Wiley-Interscience: Hoboken, NJ, USA, 2007; ISBN 978-0-470-08147-1.

34. Breakthroughs in Statistics; Kotz, S.; Johnson, N.L. (Eds.) Springer Series in Statistics; Springer: New York, NY, USA, 1992; ISBN 978-0-387-97566-5.

35. Jennrich, R.I. An Asymptotic χ2 Test for the Equality of Two Correlation Matrices. J. Am. Stat. Assoc. 1970, 65, 904-912.

36. Kapandji, A.I.; Tubiana, R. Physiologie Articulaire: Schémas Commentés de Mécanique Humaine: 805 Dessins Originaux de l'auteur; Maloine: Paris, France, 2005; ISBN 978-2-224-02647-9.

37. Rudolph, K.S.; Axe, M.J.; Snyder-Mackler, L. Dynamic stability after ACL injury: Who can hop? Knee Surg. Sports Traumatol. Art. 2000, 8, 262-269.

38. Winby, C.R.; Lloyd, D.G.; Kirk, T.B. Evaluation of different analytical methods for subject-specific scaling of musculotendon parameters. J. Biomech. 2008, 41, 1682-1688.

39. Martinez, R.; Assila, N.; Goubault, E.; Begon, M. Sex differences in upper limb musculoskeletal biomechanics during a lifting task. Appl. Ergon. 2020, 86, 103106.

40. Hoang, H.X.; Pizzolato, C.; Diamond, L.E.; Lloyd, D.G. Subject-specific calibration of neuromuscular parameters enables neuromusculoskeletal models to estimate physiologically plausible hip joint contact forces in healthy adults. J. Biomech. 2018, 80, 111-120.

41. Gagnon, D.; Larivière, C.; Loisel, P. Comparative ability of EMG, optimization, and hybrid modelling approaches to predict trunk muscle forces and lumbar spine loading during dynamic sagittal plane lifting. Clin. Biomech. 2001, 16, 359-372.

42. Robertson, D.G.E.; Caldwell, G.E.; Hamill, J.; Kamen, G.; Whittlesey, S.N. Research Methods in Biomechanics, 2nd ed.; Human Kinetics: Champaign, IL, USA, 2014; ISBN 978-0-7360-9340-8.

43. Cholewicki, J.; McGill, S.M.; Norman, R.W. Comparison of muscle forces and joint load from an optimization and EMG assisted lumbar spine model: Towards development of a hybrid approach. J. Biomech. 1995, 28, 321-331.

44. Palmerud, G.; Forsman, M.; Sporrong, H.; Herberts, P.; Kadefors, R. Intramuscular pressure of the infra- and supraspinatus muscles in relation to hand load and arm posture. Eur. J. Appl. Physiol. 2000, 83, 223-230.

45. Ward, S.; Hentzen, E.; Smallwood, L.; Eastlack, R.; Burns, K.; Fithian, D.; Friden, J.; Lieber, R. Rotator Cuff Muscle Architecture: Implications for Glenohumeral Stability. Clin. Orthop. Relat. Res. 2006, 448, 157-163.

46. Wattanaprakornkul, D.; Halaki, M.; Boettcher, C.; Cathers, I.; Ginn, K.A. A comprehensive analysis of muscle recruitment patterns during shoulder flexion: An electromyographic study. Clin. Anat. 2011, 24, 619-626.

47. Jeno, S.H.; Varacallo, M. Anatomy, Back, Latissimus Dorsi. In StatPearls; StatPearls Publishing: Treasure Island, FL, USA, 2019.

48. Arslan, Y.Z.; Jinha, A.; Kaya, M.; Herzog, W. Prediction of muscle forces using static optimization for different contractile conditions. J. Mech. Med. Biol. 2012, 13, 1350022.

49. Bergmann, G.; Graichen, F.; Bender, A.; Kääb, M.; Rohlmann, A.; Westerhoff, P. In vivo glenohumeral contact forces-Measurements in the first patient 7 months postoperatively. J. Biomech. 2007, 40, 2139-2149.

50. Westerhoff, P.; Graichen, F.; Bender, A.; Halder, A.; Beier, A.; Rohlmann, A.; Bergmann, G. In vivo measurement of shoulder joint loads during activities of daily living. J. Biomech. 2009, 42, 1840-1849.

51. Blache, Y.; Begon, M. Influence of Shoulder Kinematic Estimate on Joint and Muscle Mechanics Predicted by Musculoskeletal Model. IEEE Trans. Biomed. Eng. 2018, 65, 715-722.

52. Yanagawa, T.; Goodwin, C.J.; Shelburne, K.B.; Giphart, J.E.; Torry, M.R.; Pandy, M.G. Contributions of the Individual Muscles of the Shoulder to Glenohumeral Joint Stability During Abduction. J. Biomech. Eng. 2008, 130.

53. Klemt, C.; Nolte, D.; Grigoriadis, G.; Federico, E.D.; Reilly, P.; Bull, A.M.J. The contribution of the glenoid labrum to glenohumeral stability under physiological joint loading using finite element analysis. Comput. Methods Biomech. Biomed. Eng. 2017, 20, 1613-1622.

54. Ludewig, P.M.; Reynolds, J.F. The Association of Scapular Kinematics and Glenohumeral Joint Pathologies. J. Orthop. Sports. Phys. Ther. 2009, 39, 90-104.

55. Soslowsky, L.J.; Malicky, D.M.; Blasier, R.B. Active and passive factors in inferior glenohuimeral stabilization: A biomechanical model. J. Shoulder Elbow Surg. 1997, 6, 371-379.

56. Halder, A.; Kuhl, S.; Zobitz, M.; Larson, D.; An, K. Effects of the Glenoid Labrum and Glenohumeral Abduction on Stability of the Shoulder Joint Through Concavity-Compression: An in Vitro Study. J. Bone Joint Surg. Am. 2001, 83, 1062-1069.