Upper Extremity Muscle Activation Pattern Prediction Through Synergy Extrapolation and Electromyography-Driven Modeling

The complex network of the central nervous system (CNS) conveys neural commands to upper and lower extremities. Being activated by neural signals, musculotendon units spanning a specific joint produce the necessary torque to actuate that joint. Actuations in upper and lower extremity joints together create human movement. Brain injuries like stroke can damage this system, preventing the CNS from sending the necessary commands to generate muscle forces. Loss of muscle strength leads to loss in torque generation for upper and lower extremity joints and hinders patients' ability to do daily tasks. It is estimated that one billion people around the world suffer from neurological disorders [1]. Complete recovery of functionality of the upper extremity in patients with neurological disorders is only 25% [2]. A recovering patient requires rehabilitation therapies to strengthen the injured muscles. Patients with neuromuscular disease fail to produce necessary muscle force, thus failing to maintain joint moment required to perform activities of daily living. It is important to evaluate the muscle force values to have an improved understanding of the neuromuscular behavior in human motion. While it is incredibly challenging to measure muscle-tendon forces in vivo, researchers are developing computational approaches to resolve this issue [2–7]. To track a patient's muscle recovery, measuring both individual muscle forces and joint torques is important. Researchers use neuromusculoskeletal models to study muscle-tendon functions, human movements, and joint torques but accurately determining individual muscle forces has been a challenge. Joint torques are commonly estimated from movement data through inverse dynamics computationally, as discussed in Ref. [4]. However, multiple muscle-tendon units about a single joint create redundancy that leads to indeterminacy in muscle force calculation from joint torques. To solve this redundancy issue, an optimization approach based on an objective function is required to calculate the individual muscle forces [5].

Another approach to estimate the muscle forces is through electromyography (EMG) signals. A nerve's stimulation of muscle can be measured as an electrical activity directly from the muscle. Surface EMGs measure this neural activity from surface muscles during human movement [6]. Though Joint torques can be measured using experimental setups such as dynamometers, with measured EMG signals, joint torques can also be computed with [8–11] or without [12–16] a musculoskeletal model. Musculoskeletal models along with experimentally derived EMG data were used by researchers to determine lower [12–15] and upper [17] extremity muscle forces. These studies show muscle activations obtained from EMG signals are crucial to determining joint moments or muscle forces in the absence of an experimental setup. It is important to know the activations of all the muscles involved in a joint movement. Due to difficulties obtaining muscle activations from deep muscles, determining joint moments from superficial muscles only potentially risks undervaluing the net joint torque [17].

Through EMG-driven modeling, determining muscle force requires the knowledge of muscle activation [18]. Neural commands from the CNS induce small electrical activity in muscle during muscle contraction. This activity is directly proportional to the force generated by the muscle-tendon unit in the Hill-type muscle model. During isometric contraction, even though muscle length is constant, muscle force depends on muscle activation and initial muscle fiber length, i.e., $(f(aj, L˜jM))$⁠. Since muscle length is not changing, during isometric contraction, muscle forces do not depend on muscle fiber velocity. The force-velocity and force-length relationships come into play when the muscle is contracting against a load and changing its length. In the Hill-type model, can be dependent on muscle activation, muscle fiber length, and velocity $f(aj, L˜jM, V˜jM)$⁠, where $aj$ is muscle activation, $L˜jM$ and $V˜jM$ are normalized muscle fiber length and velocity of a jth muscle. Either way, it is necessary to obtain the value of activation $aj$ [8,17].

Noninvasive surface electrodes can measure the nerves' stimulation to muscle directly. Through EMG-driven modeling [12–16], obtained activation can be converted to muscle force, as well as joint torque. Even though electromyography is an incredibly powerful method, there are some challenges that limit its full utilization. Surface electrodes cannot measure EMG signals from deep muscles. To measure activations from deep muscles, fine-wire electrodes can be used. For this, a fine needle electrode is inserted into the target muscle to record its electrical activity. The procedure requires the application of local anesthesia to minimize discomfort associated with needle insertion. Once inserted, the electrode can record muscle activation during muscle movement. Determining joint moment through EMG-driven modeling with only surface electrodes may not provide accurate net joint torque [17,19].

In the literature, there are two primary methods for determining activation values for model-based approach, which are inverse dynamics [9,10] and forward dynamic simulation [4,11]. Static optimization is one of the inverse approaches where net joint moments are determined by inverse dynamics approach first. Then, this joint moment is distributed between the muscles crossing that joint. This creates a redundancy problem where multiple solutions may exist for the same net joint moment. This redundancy issue is solved by minimizing sum of squares of the muscle activations [9,10]. As for forward dynamic simulation, a time-varying activation profile can be estimated by tracking the desired kinetics and kinematics (e.g., computed muscle control tool in opensim) [4,11,20,21]. Some researchers implemented a combination of EMG-driven modeling and static optimization to predict missing muscle activations. Zonnino et al. [17] used a combination of EMG-driven forward dynamic estimator and static optimization-based neural model to estimate individual muscle forces. Sartori et al. [22] used a similar approach where missing muscle excitations were predicted using static optimization, but the model was also informed with EMG signals from known muscles. With this hybrid neuromusculoskeletal model, the prediction was done by tracking joint moments across different individuals and motor tasks. However, this hybrid model did not include lower dimensional EMG information, i.e., muscle synergy analysis, which could potentially reduce computation time.

Muscle synergy is another concept that has been introduced to activation prediction recently [23–29]. Rather than neural command being sent to individual muscles, it comes in a batch named synergies. Each synergy consists of two types of data: a time-dependent synergy excitation and a corresponding time-independent synergy vectors for each muscle [30]. The synergies are determined from the EMG data in such a way that the reconstruction of the activation signals from synergies reaches a certain accuracy. Thus, the synergy concept reduces the dimensionality of the EMG dataset. Ajiboye and Weir [26] showed that synergies can be predicted for a specific posture from a subset of postures of hand. Bianco et al. [28] explored the feasibility of estimating muscle excitation of 16 muscles from EMG information available from only eight muscles. Their study proved the effectiveness of muscle synergy analysis to approximate muscle excitation with more than 90% VAF (variance accounted for). However, this study did not include musculoskeletal modeling for the computation. Though muscle synergy analysis is a powerful concept to estimate muscle excitations and activations, most of these studies do not implement biomechanical behavior (EMG-driven musculoskeletal model) of a muscle. Ao et al. [23,24] developed a method (named synergy extrapolation “SynX”) to optimize synergy vector weights by tracking the joint moments obtained from inverse dynamics for lower extremity. This method finds missing muscle excitations one at a time with a combination of EMG-driven musculoskeletal model for gait. Even though some research implemented muscle activation prediction for the upper extremity, either they are missing a wide range of patient-specific tasks or missing inclusion of the musculoskeletal modeling [25]. In some cases, the methodology is missing the concept of muscle synergy analysis. This research aims to further check the feasibility of the approach developed by Ao et al. [23,24] for upper extremity shoulder and elbow joint using a combination of EMG-driven musculoskeletal modeling and SynX to predict muscle activation for the unmeasured muscle. To the knowledge of the authors, no previous research was done that implemented SynX for a wide range of upper extremity rehabilitation tasks. This method tracks inverse dynamics joint moment to determine synergy vector weights and predict muscle activation of selected shoulder and elbow muscles for a healthy subject. This robust methodology can be translated to determine muscle activations of unmeasured deep muscles or patient-specific population. In addition, muscle-tendon parameter values for optimal fiber length, tendon slack length, and maximum isometric force have been personalized to the experimental subject by using the proposed EMG-driven modeling approach.

Figure 1 shows a flow diagram of the summary steps for SynX and EMG-driven modeling. During experimental data collection, EMG muscle excitation and motion capture were recorded for three different tasks using the upper extremity. The experimental protocol was approved by the Internal Review Board at Texas Tech University. Four healthy male subjects (height: 180.2 (±6.81) cm, weight: 75.05 (±11.49) kg, age: 29.25(±2.21) years) were recruited to perform the upper extremity tasks. The first two tasks were single degree-of-freedom (DoF) tasks for shoulder and elbow each. The third task involves both shoulder and elbow joint with five DoF. Recorded EMG excitations went through activation dynamics to generate muscle activations. From n muscle-tendon units, one was marked missing from this step. Therefore, from here, it is assumed to have (n–1) known muscle tendon units. The known muscle activations are directly fed to the muscle synergy analysis. During muscle synergy analysis, time-dependent synergy excitations (⁠$Wm$⁠), and time-independent synergy vectors (⁠$Hm$⁠) were calculated.

From the motion capture data, subject specific (parameters optimized) musculoskeletal model calculated the joint kinematics using opensim [11,21] inverse kinematics (IK) tool. $L˜jM$⁠, $V˜jM$⁠, and moment arm (⁠$rj$⁠) are calculated during this step (for jth muscle).

These three time-dependent muscle parameters along with subject optimized time-independent parameter optimal fiber length (⁠$lj0M$⁠), tendon slack length (⁠$ljST$⁠) and maximum isometric force (⁠$Fj0M$⁠) were fed to the EMG-driven model. For a target ith DoF, the EMG-driven model calculated joint moment for known muscles (⁠$M i EMG$⁠). Joint moments were also calculated using opensim inverse dynamics (ID) tool (⁠$M i ID$⁠). The final optimization was done to by tracking this $M i ID$to determine the synergy vector weight (⁠$Hx$⁠) of the unmeasured muscle. Finally, the unknows muscle activation was determined by multiplying this synergy vector (⁠$Hx$⁠) with known muscles' time-dependent synergy weights (⁠$ax=WmHx$⁠). The accuracy of the prediction method was represented in terms of root mean squared error (RMSE) and Pearson's correlation coefficient (R).

Four healthy male participants (Fig. 2) were recruited for the study and were instructed to perform three different tasks (Fig. 3). The first task, referred to as task 1, involved simple elbow flexion–extension with the shoulder held fixed at a 30 deg humerothoracic elevation (abduction) while keeping the hand in a fist. This task primarily activates the flexors of the elbow (biceps, brachialis, and brachioradialis) while keeping the activity of the triceps to a minimum. Task 2 focused solely on the shoulder joint and involved performing a frontal plane abduction–adduction. The participants were asked to keep the elbow joint straight during movement. In addition to visual inspection, elbow joint angles computed from motion capture data were also checked to ensure minimal elbow joint movement. During this task, the palm was facing toward the thigh. The final task, task 3, was a common rehabilitation task where participants were asked to reach for an object on a table. During this task, the participant sat down. To ensure consistency among subjects, the object was placed in such a way that the participants had to fully extend the elbow joint to reach it. Each task was performed three times by each participant.

The motion capture system used in the study consisted of eight Eagle-4 cameras with 10 retroreflective 1-inch diameter markers placed on the arm during each task (Motion Analysis Corporation, Rohnert Park, CA). The cameras had a resolution of four megapixels each and could achieve a maximum frame rate of 500 frames per second. Motion capture was conducted at a frequency of 100 Hz using cortex (Motion Analysis Corporation). A motion capture volume measuring 3.04 × 3.04 square meters and an approximate height of 2.74 meters was chosen. The markers were placed on the anatomical landmarks of the upper extremity, such that they had minimal movement, shown in Fig. 2 based on protocols established in Ref. [18].

The analog EMG signals were recorded using the Trigno Wireless Biofeedback System (DELSYS Incorporated, Natick, MA), which contained 16 surface EMG sensors. Six of the sensors were used to measure the electrical activity of the muscles: three sensors were placed on the elbow muscles (biceps brachii, triceps long head, and triceps lateral head), and the remaining three were placed on the shoulder muscles (anterior deltoid, clavicular pectoralis major, and thoracic latissimus dorsi). Table 1 shows specific sensors and their target muscles tendon units.

Electromyography placement and surface preparation were done according to guidelines established in Ref. [31]. To reduce active noises, any excessive hair from the skin surface was trimmed. An alcohol swab was applied to remove surface oils and other contaminants. Analog EMG was recorded at 1299 Hz during each task using cortex. Before the tasks, a static standing T-pose trial was conducted for scaling purposes in opensim [11,21].

The markers in motion capture data were labeled and filtered at cutoff 6 Hz frequency in cortex [16]. EMG postprocess was done using visual-3d (C-MOTION, Washington, DC). Raw EMG data first were high-pass filtered with cutoff frequency of 50 Hz, then rectified and low-pass filtered with cutoff frequency of 6 Hz [30]. A fourth-order Butterworth filter smoothed out the data. Once processed, EMG data were normalized using the maximum activation value of each muscle.

Recorded EMG excitations went through first-order activation dynamics that are provided as a supplemental file.

Moment calculation implemented in Eq. (5) requires muscle-tendon length (⁠$ljMT$⁠), muscle-tendon velocity (⁠$vjMT$⁠), and moment arm (⁠$rj$⁠). A generic upper extremity opensim musculoskeletal model is selected for this study [20]. The model includes 7 DoF: three at the shoulder, two at the elbow, and two at the wrist. The joints are actuated by 50 Hill-type muscle-tendon units crossing the shoulder, elbow, forearm, and wrist [20].

The first step to use the model is scaling, to make the generic musculoskeletal model geometrically equivalent to the subject. Static marker data along with subject's mass scales the geometry of body segments, muscle actuators, and inertial properties of the model. Inverse Kinematics tool in OpenSim calculated joint kinematics by tracking marker positions from experimental data. As tasks in the experiment did not require wrist movements, those coordinates in the model were locked. Time-dependent variables muscle-tendon length (⁠$ljMT$⁠), muscle-tendon velocities (⁠$vjMT$⁠), and moment arms (⁠$rj$⁠) were computed in this step.

As the above equations generate an over-constrained linear system, it is solved in least square sense to get the values of optimal fiber length (⁠$lj0M$⁠) and tendon slack length (⁠$ljST$⁠) of the jth muscle.

where $C ij(t)= [rj(t). [aj(t).fl(L˜jM(t)).fv( V˜jM(t))+fp (L˜jM(t))]cosαj]$

Muscle synergies are calculated using non-negative matrix factorization with built-in matlab function nnmf. The recorded EMG signals were removed one at a time and treated as unmeasured, while the remaining (n–1) muscle tendon units were treated as measured.

Muscle synergy analysis is a method of reducing the dimensionality of EMG data. This is done by extracting a low-dimensional set of time-varying measured synergy excitations (W_m) and a corresponding time-invariant synergy vector containing weights (H_m) from the known muscle's EMG data. The synergy vectors describe the relation between synergy excitation and muscle activation. The accuracy of the non-negative matrix factorization is determined by variance accounted for (VAF) [34]. A minimum of 90% VAF was ensured to construct the muscle synergy vectors.

Activation for the unmeasured muscle (a_n) is constructed using the known muscles synergy excitations (W_m) and a time-invariant synergy vector containing weights (H_x) of the unmeasured muscle. The unmeasured synergy vector weights (H_x) are calculated iteratively by tracking experimental inverse dynamics joint moments (Eq. (10)) through optimization [23,24]. For this step, H_x is given a lower bound to be between 0 and 1. The unmeasured synergy vector weights (H_x) are being initialized by the optimization using randomly chosen values between 0 and 1. The unmeasured muscle activation is also constrained to be between 0 and 1.

Final optimization is done using global optimization manner with particle swarm algorithm [35,36]. $M i EMG$represents EMG-driven joint moment for ith DoF. Here m represents total number of discretized time points. In the optimization formulation, H_x is the design variable vector. Again, for one unmeasured muscle, tasks 1 and 2, H_x has a dimension of 2 × 1, since only two synergies (S = 2) are enough to reach 90% VAF. Here S represents the number of synergies for a specific task. For task 3, the number of synergies required to be four (S = 4), and so H_x has a dimension of 4 × 1 here. Therefore, there are two design variables for tasks 1 and 2 and for task 3 there are four design variables. With optimized $Hx$⁠, the final activation was calculated using $an=Wm Hx$⁠, where W_m is taken from known muscle set synergies.

The muscle activation calculation was done individually for each participant using the model and then averaged for presenting in the result section in terms of mean and standard deviation. Comparison between the calculated muscle activation and experimental EMG activation was performed using two metrics: RMSE and Pearson's correlation coefficient (R). This analysis is done using matlab.

Table 2 shows the number of synergies obtained for each task and their corresponding VAF. For the first two tasks, a minimum of two muscle synergies are found to be adequate to explain over 90% of the muscle activations, whereas, for the third task, at least four muscle synergies are required to achieve the same level of explanation.

Table 3 shows a summary of prediction accuracy for all tasks and target muscles in terms of RMSE and Pearson's correlation coefficient (R). The most accurate prediction was found for task 1 biceps long head (R = 0.9980) and the lowest to be thoracic latissimus dorsi for task 2 (R = 0.5188).

A comparison of the predicted and actual muscle activation for task 1 was illustrated in Fig. 4, and the results were examined for four subjects. The biceps long head muscle's activation pattern was predicted with high accuracy using the proposed method, as demonstrated by a low RMSE of 0.0306 (±0.0208) and a high Pearson's correlation coefficient of 0.9980 (±0.001), making it the most precise prediction among all tasks.

After successfully predicting the elbow muscle activation, the method was applied to predict activations of the shoulder muscles. Specifically, task 2 involved shoulder abduction in the frontal plane, which was deemed to be a single DoF task. During this task, two muscles were intentionally excluded from the calculation, one at a time. Meaning, when anterior deltoid was unknown, the thoracic latissimus dorsi was known, and vice versa. Figure 5 demonstrates that the proposed method could accurately predict the activation patterns of the anterior deltoid (with an RMSE of 0.1334 (±0.0817) and a Pearson's correlation coefficient of 0.9240 (±0.0619)) and thoracic latissimus dorsi muscles (with an RMSE of 0.2624 (±0.0487) and a Pearson's correlation coefficient of 0.5188 (±0.3494)).

To verify the method's predictability, a final task was conducted that required the coordination of both the shoulder and elbow joints. The reaching task (task 3) involved a total of 5 DoF, with three at the shoulder and two at the elbow. As shown in Fig. 6, the proposed method demonstrated better predictability for the anterior deltoid compared to the biceps long head. Again, the computation was done for one muscle at a time. The anterior deltoid showed Pearson's correlation coefficient of 0.7134 (±0.2507) with RMSE of 0.2964(±0.0957), whereas the biceps long head exhibited a Pearson's correlation coefficient of 0.5413 (±0.2997) with RMSE of 0.3067 (±0.1352).

This study evaluates SynX [23] method and EMG-driven modeling to predict unmeasured muscle activation by minimizing the difference between moments obtained by EMG-driven modeling and inverse dynamics for elbow and shoulder joint. Three tasks related to upper extremity rehabilitation were selected to evaluate the methodology. Primary results show that this method can predict the activations successfully for a variety of muscles as long as the muscle has a certain level of excitation. To evaluate and compare the accuracy of the predictions, the methodology marks one of the surface muscles as unmeasured and compares the predicted activation with the real EMG data of that muscle.

Tasks 1 and 2 had a lower number of muscles compared to task 3. For example, task 1 had seven muscle-tendon units. Due to the low number of muscles, only two synergies were enough to reach a VAF of 90%. In the case of task 3, with 13 muscle-tendon units, four synergies were required to reach a VAF past 90%. Authors believe, with higher number of synergies, the activation prediction will get better. But increased number of synergy vectors compromises the power of dimensionality reduction of muscle activations.

The authors used normalization of the EMG signals by each muscle's maximum activation value instead of using maximal voluntary isometric contraction (MVC) normalization. Normalization by maximum activation value is a muscle-specific normalization method that takes into account the individual differences in muscle size and contraction ability. This approach eliminates the issues associated with MVC normalization, such as the hardness obtaining true maximum contraction ability of a muscle for an individual [37]. However, it should be noted that normalization by maximum activation value may still have some limitations, such as the possible variation in maximum activation value due to different contraction types or fatigue levels [37]. The focus of this study is mainly to determine pattern of muscle activation, rather than actual value. That is why the EMG is normalized in max per muscle sense. The decision is taken due to variability and arguments in different EMG normalization methods. We believe a different normalization method such as MVC normalization or max overall normalization will still be effective in current methodology.

This methodology possesses some limitations. One of them is to determine one muscle at a time. The total number of muscles is represented as n. Among them, (n–1) are marked as known (activation value known). The one unmeasured muscle's activation is calculated. The methodology is developed to calculate this one unknown muscle. For task 3, for instance, only one of the thirteen muscle-tendon units was identified as missing, and the methodology predicted the activation of that muscle. Thus, this approach has the potential to be expanded to estimate the activation of a greater number of muscles by introducing multiple muscles in the optimization cost function (Eq. (10).

Another limitation of this study is the optimization of subject-specific muscle-tendon parameters. There are three parameters that need to be subject-specific in a musculoskeletal model: optimal fiber length, tendon slack length, and maximum isometric force. These parameters cannot be measured directly for a specific subject, rather they need to be calibrated for a specific subject. In this study, the optimal fiber length and tendon slack length were optimized according to the subject's anthropometry. But maximum isometric force is not an anthropometric property, rather it is closely related to force generation capability of a muscle. The maximum isometric force parameter values were optimized for the specific subjects by tracking the inverse dynamics joint moment. In the literature, it can be found all three parameters may be obtained by tracking the inverse dynamics of joint moments [24]. The authors acknowledge that the model is not fully subject-specific due to the non-normalization of activation dynamics parameters. In this study, these parameters were obtained from the literature [8], but a subject-specific determination of these parameters will improve the joint moment prediction by EMG-driven modeling, thus improving the muscle activation prediction.

Results show a high variation of muscle activation between subjects, collected through EMG sensors. Muscle activation can vary depending on the type of task being performed. Different tasks require different levels of muscle activation and coordination, resulting in changes in muscle activation patterns. For example, when reaching for an object with the upper extremity, the muscle activation patterns can vary between different subjects [38]. The amount of force required to reach an object may differ depending on factors such as the object's weight, distance from the body, and height. Additionally, individuals may have different strategies for completing the task, such as using their shoulder muscles more or relying more on their elbow muscles. These differences in strategy can result in differences in muscle activation patterns, even when completing the same task [39]. Similarly, tasks that involve balance or stabilization may require activation of certain muscles to a greater extent compared to tasks that do not involve these aspects. Thus, the way a task is performed can significantly impact muscle activation patterns. Even though there are constraints in task performance, kinematics variations exist. Here it shows one example for task 1 where the elbow flexion–extension is shown in Fig. 7.

Figure 7 shows the standard deviation and mean for elbow angle among four subjects in task 1. During flexion, the deviation seems lower among subjects. But as it reaches its maximum joint angle at the end of flexion, different subjects reach different peak values. Therefore, high deviation can be noticed during extension, which is normal as range of motion varies among subjects.

Instead of principal component analysis, authors chose to use non-negative matrix factorization for synergy extrapolation. The argument behind this is that previous studies established non-negative matrix factorization as the most appropriate method for synergy extrapolation [40], although Ao et al. [23] extended to principal component analysis as the best method to determine muscle activations of the lower extremity. Further research could test multiple synergy evaluation methods to identify the most accurate approach for the upper extremity. In addition, for the upper extremity, we will investigate machine learning algorithms [41,42].

For the proof of concept, the developed methodology was applied to healthy subjects to calculate unmeasured muscle activations. As discussed in the introduction, the methodology can be translated to clinical populations to measure muscle activations, as well as muscle forces, as a means of quantifying muscle recovery during rehabilitation.

This study evaluated a method to predict unmeasured muscle activation for a variety of upper extremity tasks. The results demonstrated that the proposed method could predict the activation of unmeasured muscle with high accuracy, particularly for the biceps long head muscle in task 1 and the anterior deltoid muscle in tasks 2 and 3. The study also found that a minimum of two muscle synergies were adequate for single DoF tasks, whereas four synergies were required for higher DoF. Even though the model is not fully subject-specific, the proposed method showed promising results for predicting muscle activations in specific upper extremity tasks like elbow flexion–extension, shoulder abduction–adduction, and reaching task. Future studies could focus on further validating the model for more diverse subjects and complicated upper extremity tasks including wrist joints.

• National Science Foundation (Grant Nos. CBET # 1703093 and 2014278; Funder ID: 10.13039/100000001).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.