^{*}

Edited by: Yuanyuan Mi, Chongqing University, China

Reviewed by: Shijie Zhao, Northwestern Polytechnical University, China; Karthick P. A., National Institute of Technology, India

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

InverseMuscleNET, a machine learning model, is proposed as an alternative to static optimization for resolving the redundancy issue in inverse muscle models. A recurrent neural network (RNN) was optimally configured, trained, and tested to estimate the pattern of muscle activation signals. Five biomechanical variables (joint angle, joint velocity, joint acceleration, joint torque, and activation torque) were used as inputs to the RNN. A set of surface electromyography (EMG) signals, experimentally measured around the shoulder joint for flexion/extension, were used to train and validate the RNN model. The obtained machine learning model yields a normalized regression in the range of 88–91% between experimental data and estimated muscle activation. A sequential backward selection algorithm was used as a sensitivity analysis to discover the less dominant inputs. The order of most essential signals to least dominant ones was as follows: joint angle, activation torque, joint torque, joint velocity, and joint acceleration. The RNN model required 0.06 s of the previous biomechanical input signals and 0.01 s of the predicted feedback EMG signals, demonstrating the dynamic temporal relationships of the muscle activation profiles. The proposed approach permits a fast and direct estimation ability instead of iterative solutions for the inverse muscle model. It raises the possibility of integrating such a model in a real-time device for functional rehabilitation and sports evaluation devices with real-time estimation and tracking. This method provides clinicians with a means of estimating EMG activity without an invasive electrode setup.

Muscle contractions generate tension and, as a result, a moment about the human joint. Knowing the muscle forces or activation signals during human movement can assist in understanding the underlying biomechanical systems (Crowninshield and Brand,

Forward simulation of a musculoskeletal system with multiple muscles.

Conversely, the inverse dynamic problem yields either (I) muscle tensions or (II) EMG signals from a predefined kinematic motion, introducing a redundancy resolution problem (Anderson and Pandy,

Inverse simulation of a musculoskeletal system with multiple muscles, using

In conventional dynamic approaches, muscle activations or EMG signals are determined by an optimization algorithm (

Machine learning is an instrument or model for solving complicated mathematical problems without knowing the analytical relationship between the inputs and outputs (Rane et al.,

Recently, machine learning models have been developed to estimate skeletal muscle tensions without explicit modeling of the physical muscle behaviors (Arjmand et al.,

Many previous studies (Arjmand et al.,

A successful study by Gonzalez-Vargas et al. (

the development of an optimized recurrent neural network configuration to estimate EMG signals from kinematic and dynamic data;

the evaluation and acquiring of necessary biomechanical data for the estimation, including joint angle, velocity, acceleration, joint torque, and activation torque;

the assessment of estimation possibility of EMG channels for a given joint;

the evaluation of the general and subject-based model in the estimation model; and

an explanation of how InverseMuscleNET can be used in applications for biomechanical analysis (inverse dynamic musculoskeletal simulation), post-rehabilitation analysis (comparison of subjects' EMG signals against those of a healthy individual via estimated activation signals), sports engineering/optimization, functional electrical stimulation control (application of the estimated healthy limb muscle activation to the functional electrical stimulation probe), and assessment and biofeedback (Gonzalez-Vargas et al.,

In this paper, firstly, data preparation is introduced. Secondly, the machine learning model, configuration optimization, and evaluation algorithms are described. Finally, results are discussed and evaluated.

Experimental data of object manipulation in the sagittal plane was used for training the inverse muscle model. The protocol for data collection and preparation is discussed in the current section. The motor task is visualized in

A depiction of EMG electrodes and the retroreflective position marker placements by Whittaker et al. (

Whittaker et al. (

The 3D upper-limb motion and muscle EMG signals were recorded at the same time. As suggested by Avers and Brown (

The following sections describe how the upper-limb position, external load, and EMG signals are used in the data preparation (

Schematic of data preparation for training InverseMuscleNET. The inputs of InverseMuscleNET are joint angle, joint velocity, joint acceleration, joint torque, and activation torque. The outputs of the InverseMuscleNET are the predicted EMG signals.

A 10 Hz low-pass filter was used to eliminate the recording noise in the 3D joint angle (Euler angles for shoulder joint) data. Joint velocities were calculated using (A) the numerical first derivative of position (Euler angle), (B) the transformation of position (Euler angle) derivatives to angular velocity by Equation (1), and (C) a low-pass filter with a 20 Hz cut-off frequency. Joint accelerations were estimated using another numerical derivative (first) of the joint velocities and a 30 Hz low-pass filter. In this paper, the joint angle θ(

An upper limb skeletal dynamic model was used within inverse dynamic simulations to estimate the required joint moments. The model had 3-DoF of rotation at the shoulder, only flexion/extension at the elbow, axial rotation (pronation/supination) of the forearm, and no wrist joint. The upper-arm and lower arm body segment inertial parameters (BSIP) were estimated from Dumas et al. (_{h} and was used as one of the inputs of the inverse muscle mapping. In Equation (2),

Equation (2) is solved for Q, one component of which is _{h}.

Nasr et al. (

The shoulder elevation joint torque (_{h}) was converted to an activation torque using a Muscle Torque Generator (MTG) model (Inkol et al., _{act}, _{a}, _{v}, and _{p} are the activation torque, the function of position-scaling, the velocity-scaling, and the passive torque respectively. The position-scaling, the velocity-scaling, and the passive functions were adopted from experimental study by McNally and McPhee (_{act} was defined for the shoulder flexor/extensor and was used as one of the five biomechanical input signals for the inverse muscle mapping.

The raw EMG data was filtered using seven steps: (1) a 20–500 Hz band-pass filter (signals lower than 20 Hz were ignored to eliminate motion artifacts, and higher than 500 Hz were ignored as they had minor power spectral density Reaz et al.,

Machine learning is an adaptive system that trains by using interconnected nodes or neurons in a layered composition. RNNs can outperform forward networks for direct muscle modeling (EMG to biomechanical signals) (Nasr et al.,

The nonlinear autoregressive with external input neural (NARX) network was selected as an RNN learning model. NARX is a recurrent network with feedback connections that are commonly used in time-series modeling. It predicts the next output signal value of a nonlinear dynamic system. A schematic of the NARX network is shown in

where _{y}) and former values of the input signal _{u}). _{y} and _{u} are the number of the previous values of the output signal and input signal. As shown in _{y}) and operate on current sequences

Schematic of a two-layer NARX network as an RNN learning model for EMG signal prediction. _{u} and _{y} are the input signal former values and the output signal former values, respectively.

The four primary configuration parameters of the RNN model are the: (1) number of hidden layers, (2) number of nodes in each hidden layer, (3) number of input signal former values, and (4) number of output signal former values. Using many of these parameters would lead to a high variable model that would require a time-consuming training process on a large dataset. One of the goals of this paper is to find the optimum configuration that can map kinematic and dynamic signals to EMG signals with limited training time and a limited dataset.

The datasets consisted of the five input signals and one output channel for 17 subjects. The inputs of the RNN model are biomechanical signals such as (1) joint angle, (2) joint velocity, (3) joint acceleration, (4) joint torque, and (5) activation torque. The first three kinematic signals (joint angle, velocity, and acceleration) are calculated in section 2.1. The last two dynamic signals (joint torque and activation torque) are calculated in sections 2.2 and 2.3, respectively. The outputs of the RNN model are the predicted EMG signals. The filtered EMG prepared in section 2.4 is used for the model output.

There are two methods for assessing and evaluating a machine learning model: (I) subject-based estimation and (II) general-model performance. The subject-based estimation test is used as a full estimation model. It employs a data set from one subject that has not been utilized for training purposes. Thus, 14 subjects' data sets were used for training, two different subjects' data sets were used to validate the model, and one subject's data set was used in the estimation test. This method is used to optimize the RNN configuration in sections 4.1, 4.4 and is used in

Primarily, out-of-sample (OOS) evaluation was used for the specific selection of training and validation sets. Specifically, cross-validation (CV) is a method that comparatively divides the databank into two sets of training and validation. K-fold CV acquires

To solve the nonlinear least-squares problem of the mapping model, the Levenberg-Marquardt backpropagation algorithm was used. The backpropagation algorithm calculates the Jacobian of the training performance concerning the weight and bias variables of the network. The initial adaptive mu, mu decrease factor, mu increase factor, and maximum mu value were 0.001, 0.1, 10, and 1e10, respectively. The training performance was evaluated using the mean squared normalized error (MSE) function. Parallel processing with the MEX 4 workers was used to make the training process faster. The training time was limited to 1-h, and the epoch (the number of passes of the entire training dataset to the model) was set to infinity. The linear regression of targets relative to outputs (R) was used as the mapping performance in this paper. A PC with an Intel^{®} Core™ i7-3370 CPU @ 3.40GHz and 16.0 GB of memory was used for training the RNN model.

The estimation performance relies on the NARX network configuration, consisting of the number of hidden layers _{l}, the nodes in each layer _{n}, the input signal former values _{u}, and the output signal former values _{y}. These configuration parameters depend on the nature of the input signals (biomechanical signals) and the output signals (EMG signals). This research tested different possible configurations and evaluated the estimation accuracy.

For preliminary optimization of the machine learning model's configuration, all available biomechanical signals (joint angle, velocity, acceleration, torque, and activation torque) were used as the input signals. All the EMG signals were used as the output signals. The input signal sensitivity and the estimation possibility of the output signals were analyzed after the preliminary optimization of the RNN configuration.

The training condition was limited to a maximum of 1-h and 200 maximum epochs. Different NARX network configurations were tested to evaluate the mapping or pattern estimation of the RNN model. The minimum and maximum optimization limits, along with steps, are presented in

The optimization variables, limits, and steps of RNN configuration.

Number of input signal former values | _{u} |
1 | 10 | 1 |

Number of output signal former values | _{y} |
1 | 10 | 1 |

Number of hidden layers | _{l} |
1 | 3 | 1 |

Number of nodes in each hidden layer | _{n} |
10 | 50 | 5 |

EMG signals were measured from 11 sites over the right arm muscles. The goal of this research was to estimate EMG signals from biomechanical signals. For the mentioned pick and place motion in the sagittal plane, some of the muscles may not be required for the mentioned motion. For evaluating the EMG signals estimation possibility, each EMG signal was selected as the output of the RNN model, and the estimation error was calculated as an assessment criterion. This assessment has been done using the general-model for comprehensive results and conclusions.

The joint angle, velocity, torque, and activation signal are the variables in the mathematical muscle model (Winters,

The sensitivity analysis to input biomechanical signals was done using a sequential backward selection algorithm for evaluating the sensitivity and dominant signals. The sensitivity analysis algorithm is shown in

Ignore the signal index

If

The signal index j with the lowest mapping error or highest mapping performance is labeled as the least dominant biomechanical signal and is deleted from the biomechanical input pool (new {

If the biomechanical pool is not empty, go to step (1).

Schematic of the sequential backward selection algorithm for evaluating the sensitivity analysis and dominant biomechanical signals for estimating the EMG channels.

The number of input signals was reduced using the information from section 3.6. The number of output signals was decreased by employing the results of section 3.5. Since the input and the output signals were optimized to the dominant and non-negligible ones, the initial configuration was optimized according to the new conditions of the input and output signals. The final RNN configuration optimization was done using the best results of sections 3.4 and 4.4.

The optimum configuration of the RNN model, the analysis of dominant input signals, and the estimation possibility of the EMG signals are presented in the following sections.

The mentioned algorithm tested 1,701 possible configurations, and the top 9 configurations are presented in _{l}) was 1, not 2 or 3. Thus, the mapping of biomechanical signals to EMG signals should be limited to a linear function. Secondly, the optimal number of nodes in each layer (_{n}) was approximately 20 to 30. Fewer or extra nodes in each layer lead to lower performance. Fewer nodes means a fewer number of weights and biases, leading to lower regression accuracy. In contrast, extra nodes in a model require more time and data to train (Lee et al., _{u}) was between 3 and 6 which shows that the history of biomechanical signals is important for pattern regression of EMG signals. Finally, the number of the output signal former values (_{y}) was required to be between 4 and 7 to give the best performance for the mapping. The number of output signal former values tries to model the dynamics of the EMG signals. It shows that the current EMG signal value is related to the previous EMG signal values. We have used the first RNN configuration in

Top performance of RNN model configurations.

_{l} |
_{n} |
_{u} |
_{y} |
||||||
---|---|---|---|---|---|---|---|---|---|

1 | 20 | 3 | 5 | 1,046 | 46 | 91.1 | 85.5 | 80.6 | 89.9 |

1 | 20 | 6 | 6 | 1,466 | 45 | 90.9 | 85.5 | 80.0 | 89.6 |

1 | 25 | 5 | 5 | 1,556 | 42 | 91.5 | 83.6 | 73.9 | 89.5 |

1 | 20 | 6 | 3 | 1,106 | 40 | 92.3 | 76.1 | 75.5 | 89.4 |

1 | 25 | 6 | 5 | 1,681 | 41 | 91.3 | 81.1 | 78.2 | 89.3 |

1 | 25 | 4 | 7 | 1,731 | 40 | 90.2 | 83.4 | 83.7 | 89.0 |

1 | 20 | 4 | 4 | 1,026 | 41 | 90.0 | 84.6 | 80.5 | 88.8 |

1 | 30 | 2 | 3 | 1,056 | 41 | 90.0 | 80.2 | 79.1 | 88.2 |

1 | 20 | 6 | 5 | 1,346 | 51 | 86.7 | 85.5 | 82.4 | 86.3 |

An example of using the RNN model with the preliminary optimized configuration for mapping the biomechanical signals to EMG signals is shown in

The performance of the preliminary configuration of InverseMuscleNET on test data to estimate EMG signals for 11 upper limb muscles.

In

The mapping performance of 11 EMG signals for the general-model. All regression (training and validation) accuracies are shown on the top, and the K-fold cross-validation regression accuracies (average, maximum, and minimum) are shown on the bottom.

The results of the sensitivity analysis of the biomechanical signals are shown in

Sensitivity analysis of 5 input biomechanical signals.

1 | 5 | ✓ | ✓ | ✓ | ✓ | ✓ | 83.9 | 91.0 | |||

2 | 4 | ✓ | ✓ | ✓ | ✓ | 81.6 | 85.8 | 1 | |||

4 | ✓ | ✓ | ✓ | ✓ | 82.3 | 86.6 | 3 | ||||

4 | ✓ | ✓ | ✓ | ✓ | 83.8 | 87.9 | 5 | 5 | Joint acceleration | ||

4 | ✓ | ✓ | ✓ | ✓ | 82.6 | 87.0 | 4 | ||||

4 | ✓ | ✓ | ✓ | ✓ | 82.1 | 86.0 | 2 | ||||

3 | 3 | ✓ | ✓ | ✓ | 79.9 | 83.2 | 2 | ||||

3 | ✓ | ✓ | ✓ | 80.9 | 83.9 | 3 | |||||

3 | ✓ | ✓ | ✓ | 81.0 | 84.4 | 4 | 4 | Joint velocity | |||

3 | ✓ | ✓ | ✓ | 79.2 | 82.5 | 1 | |||||

4 | 2 | ✓ | ✓ | 75.8 | 82.0 | 2 | |||||

2 | ✓ | ✓ | 76.0 | 82.9 | 3 | 3 | Joint torque | ||||

2 | ✓ | ✓ | 69.3 | 79.0 | 1 | ||||||

5 | 1 | ✓ | 65.5 | 77.8 | 2 | 2 | Activation torque | ||||

1 | ✓ | 58.2 | 76.9 | 1 | 1 | Joint angle |

The final RNN configuration optimization utilized all five biomechanical signals and 6 of the 11 EMG signals. The configuration optimization was similar to section 3.4, but the 6 EMG signals consisted of Anterior Deltoid (ADEL), Infraspinatus (INFR), Middle Deltoid (MDEL), Pectoralis Major (PECC), Serratus Anterior (SERR), and Upper Trapezius (UTRA). These were selected according to the observations in section 4.2. Suppose the rest of the signals (those five with lower regression accuracy) are needed for specific clinical purposes/applications for different tasks. In that case, more data with different motions and tasks are required for training. The training properties were similar to the previous step. However, the maximum training time was set to 2 h, and the maximum epoch was set to 200 times. The training could be stopped after 100 epochs if the validation performance had not changed.

In this case, the training was stopped after 134 iterations instead of 200 iterations since the best performance occurred at epoch 34 (

The RNN training (blue), validation (green), and test (red) performance (MSE) for 134 iterations. The exported model is based on epoch 34. The training stops after 100 iterations after the best validation performance.

Error histogram indicating the difference between actual EMG signal and estimated output.

Regression or scatter plot of training, validation, test, and all data to RNN model.

From the optimization, the number of hidden layers (_{l}) was 1. The machine learning model consisted of two layers, one hidden layer, and one output layer. The number of nodes in the hidden layer (_{n}) was 20. The number of input signal former values (_{u}) was 3 and number of the output signal former values (_{y}) was 5. Since the sample rate of all input and output signals was 50 Hz, the input signals and the feedback signals should consist of 0.06 and 0.1 s of the vector input (containing current and previous states).

For example, by feeding the biomechanical signals of the test subject to the model and having the initial value of EMG signal at 0 s, the estimation of 6 dominant EMG signals for 60 s has been presented in

The performance of final optimized InverseMuscleNET configuration on test data for the estimation of EMG signals from 6 upper limb muscles.

The model's accuracy was slightly increased when using the elbow joint angle, velocity, and acceleration in addition to the previous shoulder joint angle, velocity, acceleration, torque, and activation torque. The training, validation, test, and all regression accuracy increased from 91.0, 82.6, 88.9, and 89.2% to 91.3, 83.3, 88.9, and 90.9%, respectively. The improvement is not much but can help the accuracy of the estimation model. We propose using the joint biomechanical variables relevant to the muscles that impact them.

Optimal control and musculoskeletal modeling have been studied as potential tools for functional rehabilitation and sports evaluation. Despite significant progress, musculoskeletal modeling has encountered some significant problems. The muscle within a musculoskeletal system is complex and challenging to model, and the static optimization simulation requires a time-consuming iterative problem solution based on some assumptions (Williams and Constandinou,

This study intended to discover how to predict muscle activations accurately and efficiently using a machine learning model, specifically a recurrent neural network. Machine learning is an instrument or model for solving complicated mathematical problems without knowing the analytical relationship between inputs and outputs (Rane et al.,

Recently, machine learning models have been developed to estimate skeletal muscle tensions without explicit modeling of the physical behaviors of muscles (Arjmand et al.,

Schematic of using InverseMuscleNET within a forward dynamic simulation of the skeletal model for biomechanical analysis and sports engineering/optimization.

The inverse muscle model was a NARX network, selected as an RNN learning model. The RNN includes a feedback connection to provide memory capacity for time-varying data. Successfully, Dao (

Moreover, the RNN model was used to estimate muscle activation with fewer input signals in the sensitivity analysis section. These fewer input signals demonstrated the flexibility of the RNN model in the estimation and muscle model. This model could be trained based on the availability of input data, consisting of joint angle, velocity, acceleration, torque, and activation torque. It is noteworthy to mention that this RNN flexibility could not be done with static optimization approaches.

Developing a machine learning model has some challenging aspects, including (I) determining the optimal configuration and (II) requiring a rich and representative dataset for training. We have tested multiple topologies with the same training conditions to determine the model topology's optimal configuration. The number of hidden layers (1), the number of nodes in each layer (20), the number of input signal former values (3), and the number of the output signal former values (5) have been optimized to have the highest mapping accuracy. Although having a model with more parameters may have more accuracy, the model may have a delay and require more data for training. Having fewer parameters in the trained network allows the network to be saved and straightforwardly used in real-time.

Forming a rich database for model training was established on the data collection for a repetitive manual task from 17 healthy right-handed young individuals (Whittaker et al.,

Another strength of this study relates to the rich number of involved muscles. Having access to 11 muscles provides the opportunity to evaluate the estimation possibility. The dominant and most accurate estimation occurred with the Anterior Deltoid (ADEL), Infraspinatus (INFR), Middle Deltoid (MDEL), Pectoralis Major (PECC), Serratus Anterior (SERR), and Upper Trapezius (UTRA) for the described task and motion.

Analysis of using five biomechanical signals for the estimation provides the possibility of sensitivity analysis. The order of most essential signals to least dominant ones was as follows: joint angle, activation torque, joint torque, joint velocity, and joint acceleration. This result demonstrated that the estimation of activation signals was relevant to joint angle, activation torque, and joint velocity. The results from the sensitivity analysis concur with previous classical muscle mathematical model studies. The muscle pennation angle and muscle-tendon length have been found to be sensitive to the joint angle and velocity (Winters,

Since we found that the input signal and the feedback signal should consist of 0.06 and 0.1 s of the vector input (containing current and previous states) in the optimization of the RNN configuration process, the estimation of EMG signals relies on the history of the previous value of EMG signals and the biomechanical signals.

The estimated and the measured EMG signals in

The main limitation of the mentioned model for inverse muscle stimulation is the amount of training data. The machine learning model depends significantly on training data, like the human brain that tunes the neurons' weights and bias using motion and response. As far as we are concerned, there were 14 sets of input/target pairs from fourteen subjects, two to three subjects' data for validation, and one set for training. The data consisted of shoulder flexion/extension movement (the combination of external load, slow, and a fast motion), each of them in 12 trials of the pick-and-placing task. Having more data for more ranges of motion would facilitate a more general model.

A model estimating the shoulder muscle activation was developed using an optimized recurrent neural network. This estimation was based on biomechanical (kinematic and kinetic) input signals. This study suggests that the classical static optimization of the dynamic model could be replaced with an inverse muscle model using a machine learning model. This approach could increase potential decision support tools for functional rehabilitation with real-time estimation and muscle activation or EMG signal tracking.

The datasets presented in this article are not readily available due to ethical and privacy reasons. Requests to access the datasets should be directed to Clark R. Dickerson,

The studies involving human participants were reviewed and approved by the University of Waterloo Office of Research Ethics. The patients/participants provided their written informed consent to participate in this study.

AN and JM contributed to the conception and design of the study. AN organized the interpretation of data, the model, and the code. AN, KI, SB, and JM performed the statistical analysis. AN, SB, and KI wrote the first draft of the manuscript and wrote sections of the manuscript. KI and JM performed revising it critically for important intellectual content. All authors contributed to manuscript revision, read, and approved the submitted version.

This research is supported by funding from the Canada Research Chairs Program and the Natural Sciences and Engineering Research Council of Canada (NSERC).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

The authors would like to thank Prof. Clark R. Dickerson (University of Waterloo, Kinesiology and Health Sciences) for providing the experimental data.