Mapping Individual Motor Unit Activity to Continuous Three-DoF Wrist Torques: Perspectives for Myoelectric Control

The surface electromyography (EMG) decomposition techniques provide access to motor neuron activities and have been applied to myoelectric control schemes. However, the current decomposition-based myoelectric control mainly focuses on discrete gestures or single-DoF continuous movements. In this study, we aimed to map the motor unit discharges, which were identified from high-density surface EMG, to the three degrees of freedom (DoFs) wrist movements. The 3-DoF wrist torques and high-density surface EMG signals were recorded concurrently from eight non-disabled subjects. The experimental protocol included single-DoF movements and their various combinations. We decoded the motor unit discharges from the EMG signals using a segment-wise decomposition algorithm. Then the neural features were extracted from motor unit discharges and projected to wrist torques with a multiple linear regression model. We compared the performance of two neural features (twitch model and spike counting) and two training schemes (single-DoF and multi-DoF training). On average, 145 ± 33 motor units were identified from each subject, with a pulse-to-noise ratio of 30.8 ± 4.2 dB. Both neural features exhibited high estimation accuracy of 3-DoF wrist torques, with an average $\text{R}^{{2}}$ of 0.76 ± 0.12 and normalized root mean square error of 11.4 ± 3.1%. These results demonstrated the efficiency of the proposed method in continuous estimation of 3-DoF wrist torques, which has the potential to advance dexterous myoelectric control based on neural information.


I. INTRODUCTION
M YOELECTRIC control techniques have found wide applications in human-machine interfaces such as prosthetics, exoskeletons, and gaming [1], [2], [3], [4]. The bio-signal used in these systems, electromyography (EMG), contains rich information about muscle contraction and thereby could be used to decode the human motion intention. Among current myoelectric control schemes, state-machine and pattern recognition are the two most representative and popular methods. The state-machine was first proposed, realizing a few motions' simple control [5]. Pattern recognition has been applied to myoelectric control since the 1980s [6]. From then on, various classification models have been investigated to improve the control performance [7], [8], which could discriminate tens of motions or gestures with high accuracy (>95%). However, these two control schemes could only identify the discrete motions. Currently, one of the most challenging issues is how to decode the continuous kinematics/kinetics of human movements.
The common strategies for continuous control are based on regression models such as linear regression, Gaussian mixtures, and neural networks [2], [9], [10], [11]. The EMG features, in the time or frequency domain, are usually extracted and projected to the continuous force/angle of multiple degrees of freedom (DoFs) using these regression models [12], [13]. In regression-based continuous control schemes, the EMG signals are usually regarded as stochastic signals, where the movement kinematics/kinetics are reflected by the statistical properties. Therefore, it is difficult to interpret the neural drive transferred into the muscles. Another continuous control scheme is the simultaneous and proportional control (SPC), which is proposed based on muscle synergy [14], [15]. The SPC scheme takes muscle coordination into consideration during human movements. Based on the non-negative matrix factorization (NMF), the muscle synergies are extracted from several EMG channels, allowing for the observation of activation patterns and their activation intensities across multiple muscles. Nevertheless, the implementation of NMF is based on the EMG features, which could merely reflect the partial neural drive [16], [17].
Although SPC-related studies have gained tremendous momentum in the past decade, achieving natural and intuitive EMG-based control is still challenging. The EMG This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ decomposition techniques provide another possible solution. Using the EMG decomposition, the motor unit action potential train (MUAPT) could be identified [18], directly representing the motor neuron discharges derived from the spinal cord. After verifying the decomposition algorithm itself [19], [20], [21], the MUAPT-based control schemes have been rapidly applied to the continuous estimation of movement kinematics, such as the angle and force/torque of fingers, wrists, and knees [22], [23], [24], [25]. On the one hand, the discharge properties of motor units could be projected to the activation level of each DoF linearly [22], [26]. The discharge rate of motor units has been demonstrated to be highly correlated with kinematics/kinetics [27]. On the other hand, the neural drive can be estimated from MUAPTs and sent into a neuromuscular-skeletal model [25], [28]. The EMG decomposition techniques contribute to a better understanding of the control mechanism underlying human movements. Furthermore, the derived MUAPT-based control scheme has shown promising results in intuitive and dexterous myoelectric control.
However, the current MUAPT-based control schemes still have several limitations. The previous studies mainly investigated the feasibility of the MUAPT-based continuous control during single-DoF or two-DoF tasks [22], [29], [30]. In addition, compared with the conventional continuous control, the experimental protocol in MUAPT-based control is comparatively simple, which cannot satisfy the practical application. To overcome these limitations, this study aimed to map the motor unit discharges, identified from high-density surface EMG signals, to the multi-DoF and continuous wrist torques. A segment-wise algorithm [31], which is suitable for decomposing EMG signals of multi-DoF tasks, was used to decode the MUAPTs. Then we proposed a continuous multi-DoF estimation method and validated its performance in individual DoFs' movements and their various combinations. The estimation performance of two training methods and two neural features were tested with a multiple linear regression model.

A. Subjects
Eight healthy subjects with no neurological or psychiatric disorders (5 males and 3 females, aged 25 ± 3 years, all righthanded) participated in the experiment. All participants signed informed consent prior to participation. The experimental protocol as well as the informed consent were approved by the local ethics committee of Shanghai Jiao Tong University (approval number B2020026I) and followed the Declaration of Helsinki.

B. Experimental Setup
In the experiment, the subjects performed the isometric contractions of the wrist without joint movements. The torques and EMG signals were recorded concurrently. We used an external trigger channel to synchronize the two signals and designed a graphical user interface (GUI) to provide the real-time display of the torque signals for subjects (Fig. 1).
The wrist torque curves were displayed to the subject in real time and indicated by the cursor's position in a GUI, as shown in Fig. 1c. The clockwise and counterclockwise rotation of the arrow represented Pro and Sup (DoF1). The left-right and updown movements of the arrow represented the activation of Fle-Ext (DoF2) and Rad-Uln (DoF3), respectively.
2) EMG Signals: High-density surface EMG signals were recorded from the forearm muscles with three electrode grids (ELSCH064NM3, OT Bioelettronica, Italy). Each grid contained 64 channels (8 × 8) with a 10 mm inter-electrode distance in two directions. The electrode grids were mounted around the proximal third of the forearm (Fig. 1d), covering most wrist-related muscles. The EMG signals were acquired and amplified by a multi-channel amplifier (EMG-USB2+, OT Bioelettronica, Italy). The gain was set as 500 with a sampling rate of 2048 Hz.

C. Experimental Protocol
Before the experiment, the subjects had enough time to familiarize themselves with the experimental protocol and the GUI. The maximum contraction torque of each DoF was measured and used to normalize the activation space of the cursor. During the experiment, the subject sat in a chair with his/her forearm placed into the torque measurement device and supported by an armrest. The forearm and palm were secured in the device to ensure the isometric contraction of wrist movements.
The subjects were required to perform wrist motor tasks involving 1-3 DoFs. The experimental protocol included 7 sessions, as described in Tab. I. In each session, we designed different tasks to imitate various circumstances in realistic motions as much as possible. Sessions 1-3, 4-6, and 7 covered the 1-3 DoF wrist tasks, respectively. In each task, the subject performs the corresponding movements cyclically, with a frequency of about 0.5 Hz. The contraction intensity of each movement was medium, which is not strictly requested. For each task, the subject performed two repeated trials. Sometimes one more trial was required depending on the completeness of the task, while always two trials were used for the following analysis. Each trial lasted for 20 seconds. Rest periods of 1-2 minutes between sessions were allowed to avoid muscle fatigue. The duration of the experiment was about 1 hour for each subject.

D. Data Analysis
The flow chart of data analysis is shown in Fig. 2. After preprocessing, we decomposed the EMG signals from 1-DoF tasks into MUAPTs using a segment-wise decomposition method [31]. The obtained separation matrix was used to decompose the multi-DoF EMG signals in a pseudo-online manner. The neural features were extracted and projected to the torques using a linear regression model. 1) Preprocessing: The torque signals were first up-sampled to 2048 Hz by interpolation to match the frequency of EMG signals, and then low-pass filtered by 20 Hz to remove the baseline noise. For the EMG signals, we used a 4th-order The flow chart of data analysis. The EMG signals from 1-DoF tasks were first decomposed offline to obtain the separation matrix. Then we decomposed the EMG signals of multi-DoF tasks in a pseudo-online manner. Two features, twitch force and spike count, were extracted from the training data and projected to the wrist torques using multiple linear regression. The regression model was validated on the testing data.
Butterworth filter with a bandpass frequency of 20-500 Hz and a comb filter with a cutoff frequency of 50 Hz. Several EMG channels (usually < 5) with excessive noise, where the root mean square value of EMG signals was higher than three times the average, were removed. The excessive noise might be caused by poor contact between electrodes and skin.
2) EMG Decomposition: A segment-wise decomposition method was used to decompose the EMG signals [31]. The EMG signals from each DoF were regarded as a segment and decomposed into MUAPTs individually using the convolution kernel compensation (CKC) algorithm [32], [33]. Briefly, the generation of EMG signals can be modeled as a multi-inputmulti-output system, where the EMG signals are described as the convolutive mixture of a series of motor unit discharges and their action potentials. The unknown mixture matrix is compensated, and the spike train of each motor unit is estimated as:ŝ whereȳ is the extended EMG signal with an extending factor of 10, c s jȳ = E(ȳ(t)s T j (t)) is the cross-correlation vector between the jth spike train s j and the EMG signals, Cȳȳ = E(ȳ(t)ȳ T (t)) is the correlation matrix of the extended EMG signals, and t indicates the sample point in the time domain. E(·) denotes the mathematical expectation. In this study, the cross-correlation vector (c s jȳ ) was estimated iteratively with the natural gradient descent algorithm [33]. The crosscorrelation vectors calculated from each segment were grouped to estimate the spike trains from the global EMG that was obtained by cascading the EMG signals from each DoF [31]. The estimation formula was the same as Eq. 1, while theȳ was replaced with the extended global EMG. Then the discharges were extracted fromŝ j (t) using the K-means++ clustering.
In this study, the 1-DoF tasks were decomposed into MUAPTs using the segment-wise method. The first-trial EMG signals of each DoF were cascaded and decomposed, as described above. The decomposition was applied to each grid individually. The MUAPTs identified from each grid were grouped after removing the duplications. Fig. 3 shows an example of the decomposition result of 1-DoF tasks. The separation matrix and other decomposition parameters (e.g., clustering centroids for spike extraction) were kept. For the 2/3-DoF tasks in the training data and all the tasks in the testing data, we directly multiplied the preprocessed EMG signals by the separation matrix to estimate the MUAPTs [23], [34]. The discharges were extracted based on the obtained clustering centroids. Due to the lack of ground truth, the decomposition accuracy was evaluated based on a signal-based metric, pulse-to-noise ratio (PNR) [35].
3) Feature Extraction: The MUAPT was decoded in the form of a series of motor unit firing sequences and the corresponding action potential waveforms. From MUAPTs, we extracted two kinds of neural features. One is based on a muscle force model, where the motor unit mechanically generates a twitch force at every firing instant. The other was extracted by simply counting the firings in a sliding window.
The twitch force of each firing was regarded as an impulse response and modeled as in [36]. A critically damped, secondorder system can reasonably approximate the twitch force of the ith motor unit. where P i is the peak amplitude of the twitch force, T i is the rise time (from the discharge instant) to the peak force of the impulse response. The relationship between twitch force and contraction time is approximated as an inverse power function in the form as: where T L represented the longest contraction time desired for the motor neuron pool, which was set as 90 ms in this study. c is a constant set to 4.2.
In the twitch force model, the gain of each firing varies as a function of the firing rate to simulate the non-linear force behavior. The gain in the jth firing of ith motor unit is formulated as: where I S I j indicates the inter-spike interval between the jth and ( j + 1)th firing. As a result, the twitch force of each firing and the force curve of each motor unit (F i (t)) is modeled as: where t i j indicates the jth firing instant of the ith motor unit. Fig. 4 illustrates the generation of the twitch force curve for each motor unit. In this study, the peak amplitude of each twitch was regarded as the same (P i = 1), while the contribution of each motor unit to the torques was described by the regression coefficient in the multiple linear regression model. A sliding window was applied to feature extraction. For the twitch feature, the twitch force curve in the sliding window was first generated depending on Eq. 6, and then, averaged across all the time instants in the sliding window. For the spike feature, the number of firings in the sliding window was counted for each motor unit. The effect of the window length on the estimation performance was investigated, and the window length of 300 ms was adopted for the following analysis. The overlapping between two adjacent windows was 200 ms. The torque signals were also averaged in the sliding window.
4) Regression Model: The neural features were projected to the wrist torques using a multiple linear regression model.
where y d indicates the torque of the dth DoF, b id is the regression coefficient of the ith MU for the dth DoF, m is the number of identified motor units. The regression coefficients were calculated using the least-square method: where X n×m is the feature matrix of n samples, B m×3 contains the regression coefficients of all the identified motor units and three DoFs, Y n×3 is the torques of three DoFs. 5) Performance Evaluation: The EMG signals recorded in the experiment were divided into training and testing data. Each subject performed two repetitive trials for each motor task. The EMG and torque signals in each task's first trial were used to train the regression model, and signals in the second trial were used for testing.
In the training data, the 1-DoF EMG signals (sessions 1-3) were decomposed into MUAPTs by the segment-wise algorithm. The obtained decomposition parameters (separation matrix and clustering centroids) were collected to decompose the remaining 2/3-DoF EMG signals (sessions 4-7) in the training data and all EMG signals (1/2/3-DoF, sessions 1-7) in the testing data.
After EMG decomposition, the two neural features were extracted from the identified MUAPTs. The features were divided into training and testing data similar to the EMG signals. We used two strategies to train the regression model. The first was the single-DoF training, where only the features extracted from 1-DoF tasks in the training data were used to train the regression model. The second was the multi-DoF training, where the features extracted from all the tasks (1/2/3-DoF) in the training data were used. The testing data for the two training strategies were the same, including 1/2/3-DoF tasks.
To evaluate the estimation performance of wrist torques, we extracted four metrics, which are the Pearson correlation coefficient (R), the determination coefficient (R 2 ), the normalized root mean square error (nRMSE), and the r oughness. The metric r oughness was used to evaluate the output smoothness, indicating the ratio between the fluctuation level of the estimations and recordings [11], [29]. The definitions of these four metrics are given as follows:  where y j indicates the recorded torques,ȳ is the average value of recorded torques,ŷ j is the estimated torques,ȳ is the average value of estimated torques, T s is the sampling time interval and set to 0.1 s in this study. The four metrics were calculated for each trial and each DoF.
All the data analyses were implemented in MATLAB 2022a (Matlab Inc. USA).

E. Statistic Analysis
A one-way analysis of variance (ANOVA) was applied to evaluate the effect of training strategies, features, and DoFs, on the three performance metrics. The homogeneity of variance for three metrics was first tested. If satisfied, the Bonferroni method was conducted. If not, Dunnett's C method was used instead. The significance level was set to 0.05. The symbols * , * * , and * * * indicate the significant difference with a level of P < 0.05, P < 0.01, and P < 0.001, respectively.

III. RESULTS
On average, 145 ± 33 motor units were identified from each subject, with a PNR value of 30.8 ± 4.2 dB. From three grids, the numbers of identified motor units were 51 ± 10, 56 ± 19, and 38 ± 18, respectively. Tab. II gives the detailed decomposition results of each subject. Representative decomposition results in training and testing session are shown in Figs. 3c and 5b. The average time cost of the decomposition in the sliding window was 11.9 ± 4.9 ms.
After EMG decomposition, we extracted two neural features from MUAPTs and used two training strategies for wrist estimation. As shown in Fig. 6, the single-DoF and multi-DoF training strategies exhibited similar performance in the 1-DoF tasks, only showing a significant difference in nRMSE. On the contrary, the multi-DoF training outperformed the single-DoF training dramatically when estimating the 2/3-DoF wrist torques. The average R 2 of 3-DoF torques was 0.01 ± 0.77 and 0.75 ± 0.12 for single-DoF and multi-DoF training, respectively. As to the two neural features, the twitch feature presented a better performance in the wrist estimation of 1/2-DoF tasks, with higher R/R 2 and lower nRMSE. For the 3-DoF tasks, the twitch and spike feature showed no significant difference, with R=0.88 ± 0.07/0.87 ± 0.07, R 2 =0.76 ± 0.12/0.75 ± 0.12, and nRMSE=11.4 ± 3.1%/11.6 ± 3.3%, respectively. For all the tasks, the estimation roughness of the spike feature was higher than that of the twitch feature.
The estimation results of each DoF in different motor tasks are depicted in Fig. 7. The estimation accuracy of DoF 2 was (a,c,e,g) give the comparison of estimation results between single-DoF and multi-DoF training methods, which were averaged across two features. The comparison of estimation results between these two features is shown in (b,d,f,h), where the multi-DoF training method was used.
consistently higher than that of the other two DoFs, with an average R 2 >0.8 even in 3-DoF tasks. The increase in the window length had a positive effect on the estimation accuracy of wrist torques, as depicted in Fig. 8. The rise in the window length increased the R 2 and decreased the r oughness across all motor tasks, no matter which feature was used.

IV. DISCUSSION
We investigated the feasibility, accuracy, and potential for multi-DoF and proportional myoelectric control by extracting motor unit discharges from surface EMG signals. During multi-DoF wrist tasks, it is possible to identify over 100 motor units from each subject. The high yield decomposition allowed for the continuous control of multi-DoF based on the motor unit discharges.

A. EMG Decomposition
One basis supporting the 3-DoF continuous control is the segment-wise decomposition algorithm, which could significantly increase the number of identified motor units from surface EMG signals. Limited neural information was extracted from MUAPTs using the classic decomposition methods [32], [37], so it is difficult to estimate the multi-DoF kinetics simultaneously and proportionally. Compared with previous studies, more MUAPTs were identified during the multi-DoF wrist movement [22], [38]. Although the increase in the motor unit number was at the cost of the decline of decomposition precision, the torque estimation based on the neural features showed high accuracy during multi-DoF tasks.
Moreover, the decomposition scheme used in this study is suitable for real-time applications. The efficiency of applying the pre-calculated separation matrix to newly recorded EMG signals has been demonstrated in several studies [34], [39]. After segment-wise decomposition of the 1-DoF EMG signals, we decomposed the 2/3-DoF EMG signals using the obtained separation matrix. On the one hand, this decomposition scheme could dramatically reduce the computation complexity, providing the possibility of real-time applications. The time cost of < 50 ms in the sliding window satisfied the requirement of online decomposition and left enough time for the feature extraction and regression [40]. On the other hand, the motor units could be tracked across tasks depending on the separation vectors [23], [30]. The motor units activated during 1-DoF tasks could also be identified during multi-DoF tasks.

B. Torque Estimation
After EMG decomposition, we extracted two neural features from the identified MUAPTs. Both the twitch and spike features indicate the instantaneous activation intensity of each motor unit, which are highly correlated with the continuous kinematics/kinetics. In addition, the 3-DoF wrist movements are related to different muscles innervated by different motor neuron pools. The activation pattern of motor units also contains information about the types of wrist movements. Therefore, the neural features could be applied to the continuous multi-DoF estimation of wrist torques.
Compared with spike counting, the twitch model better interprets the force generation at the motor unit level. The relationship between the motor unit discharge rates and muscle force may not be linear [41]. With consideration of individual discharges' contribution to the force formation, the force estimation precision could be improved [42]. In the 1/2-DoF tasks, the estimation accuracy by the twitch feature was higher than that by spike, demonstrating the advantage of the twitch model in torque estimation. In addition, the twitch feature could achieve smoother output than the spike feature. The advantage of the spike feature lies in its low computation complexity. Moreover, the spike feature presented a comparative performance in torque estimation, especially in the 3-DoF tasks. Both two features satisfied the requirements in real-time applications (a time delay of < 300 ms) using the current configuration of the sliding window [40]. The sliding length of the window could be changed to adjust the output frequency for different application sceneries.
The increase in the window length could improve the estimation accuracy and output smoothness, especially when using the spike feature. The neural feature became relatively stationary with a longer window, which might be the reason for accuracy improvement. Also, the improvement of the spike feature was more prominent since the twitch model has already implemented an additional filter process for the spikes. It should be noted that the non-oscillatory DoFs had much greater errors. During the real-time control, additional postprocessing methods, such as an exponential moving average filter [30], could be used to reduce the error caused by the estimation fluctuation.
We have tested two training schemes on the regression model. The single-DoF training scheme requires much less training data and has the potential to reduce the users' training burden. In the previous studies, the single-DoF training scheme has been used to estimate the 2-DoF wrist torques based on MUAPTs [29], [30], showing superior performance to time-domain EMG features. However, the estimation accuracy with the single-DoF training scheme decreased rapidly with the increase in the DoF number. The results demonstrated that the single-DoF training scheme was not a good choice when estimating the 3-DoF wrist torques based on the neural features.

C. Comparison With Previous Studies
To the best of our knowledge, for the first time, the simultaneous and continuous estimation of 3-DoF kinetics was realized based on the MUAPTs identified from surface EMG signals. In the previous work, we projected the extracted neural drive to 2-DoF kinematics/kinetics by assigning individual motor units to each movement [29], [30]. The torques of each wrist motion were estimated only based on these exclusive motor units, which are only activated during one motion. However, it is difficult to identify the exclusive motor unit with an increased number of involved motions. Therefore, we applied a multiple linear regression model to torque estimation. By training the regression model with all motor units, the contribution of each motor unit to the torques was evaluated as the regression coefficients. The proposed method could realize the simultaneous estimation of 3-DoF torques, showing comparative accuracy compared with the 1/2-DoF tasks in previous studies [22], [29].
Apart from the MUAPT-based methods for wrist torque estimation, the time-domain features are usually extracted from EMG signals, no matter whether surface or intramuscular EMG signals, and used to train the regression models [43], [44]. The R 2 of estimation accuracy could be higher than 0.8 for less than 3-DoF tasks. The estimation accuracy of 1/2-DoF tasks by the proposed method was slightly lower than 0.8. It should be noted that the R 2 in this work was averaged across all the DoF combinations. The Uln-Rad DoF showed significantly lower estimation accuracy than the other two DoFs and is usually not taken into consideration in previous studies. The torque estimation of 3-DoF tasks still remains challenging. The R 2 was about 0.65 or even lower when using linear regression models [14], [45]. Based on the neural features extracted from MUAPTs, we could achieve the R 2 of 0.76 for 3-DoF wrist torque estimation using simple multiple linear regression. It is noteworthy that the decomposition algorithm could be further improved as the number of identified MUAPTs is lower than the truly activated MUAPTs. With more MUAPTs identified, the neural features could reflect the muscle excitation more precisely. Therefore, the proposed method has the great potential to improve the estimation accuracy in the future.

D. Limitations
The main limitation of the study is that we validated the MUAPT-based torque estimation in offline analysis. The real-time performance of the proposed method in practical applications, such as prosthetic control and sign language recognition, remains for further investigation. This study mainly focused on the projection method from individual motor unit activities to the multi-DoF wrist torques. The results demonstrated the efficiency of the proposed method in simultaneous and proportional myoelectric control. Based on a simple multiple linear regression model, individual motor unit discharges could be projected to continuous wrist torques with high accuracy. Moreover, we demonstrated the feasibility of generalizing the proposed method to real-time applications.
Another limitation is that we only included wrist movements. In daily life, the motor tasks of the upper limb are usually implemented with the cooperation of the elbow, wrist, and fingers. The proposed methods for estimating the kinematics/kinetics of multiple joints will be done in future work. Finally, only non-disabled subjects were recruited. The inclusion of subjects with limb deficiency is necessary to demonstrate the general clinical applicability or to make more general claims about the observed performance.

V. CONCLUSION
We demonstrated the feasibility of mapping the individual motor unit discharges to continuous 3-DoF wrist torques. Over 100 motor units were identified from each subject by using the segment-wise decomposition method. Two neural features and two training schemes were tested on a multiple linear regression model, showing high estimation accuracy (R 2 =0.76, nRMSE=11.4%) for 3-DoF tasks. These findings support the presented estimation methods in multi-DoF and continuous myoelectric control.