Subjects

This study involved eight subjects (six males and two females, age: 26.00 ± 1.89 years) without any known muscular injuries or neuromuscular disorders participated. The experimental protocol was approved by the Ethics Review Committee of University of Science and Technology of China (Hefei, Anhui, China, under Application No. 2022-N(H)-163, February 2022). All subjects gave their informed and written consent prior to any procedure for the experiments.

Fig. 1

Experimental setup and protocol. (a). Apparatuses used for simultaneously recording thumb abduction force and HD-SEMG data. (b) Illustration of the force generation pattern with both the designed force curve (red line) and an actual recorded force curve (blue line) (c). The recorded EMG signals. (c) User interface of the software used for the online testing. There is an SEMG signal in one single channel as an example to exhibit the data stream on the top of the interface. Extracted MUSTs, measured force (red line) and predicted force (the blue line) are shown below

Experimental protocols

The apparatus of the experiment is shown in Fig. 1a. A multi-channel electrode array arranged in 8 rows × 8 columns (FlexMatrix, Shanghai, China) was attached to the abductor pollicis brevis (APB) muscles on the dominant hand of a subject to record HD SEMG signals. Each electrode probe had a diameter of 2 mm with an inter-electrode distance of 4 mm between consecutive electrodes. At the beginning of the experiment, we tested all subjects’ maximum voluntary contraction (MVC) of the tested muscle (the maximum force of the thumb abduction contractions).

A series of contraction tasks were performed in the experiment, as shown in Fig. 1b. More specifically, the muscle force gradually increased from 0 to a targeted force level in 2s and then maintained at the targeted level for 3s in each trail of the contraction task so that a corresponding 5-s segment of SEMG signal could be collected. The targeted levels were set to 20% and 30% MVC. During the experiment, the subjects were required to sit comfortably. To prevent muscular interferences from the wrist or other fingers, a set of 3D-printed modules was fixed on appropriate positions of the table. The EMG recording system included a one-dimensional load cell (LDST-V-HY, Luckly Inc., Beijing, China) for recording EMG signals and muscle force simultaneously. The HD-SEMG signals in all channels were filtered through a 10-order Butterworth band-pass filter to reduce possible low-frequency or high-frequency interference. The bandwidth of the filter was 20–500 Hz. In addition, power line interference was removed through a 50 Hz second-order notch filter. All the recorded SEMG and force data were digitized via a 16-bit A/D converter (ADS1198, Texas Instruments, TX) at a sample rate of 2 kHz. The reliability of the data recording system has been proved in our previous studies [35]- [36].

Data collection for training

During this period, the subjects were instructed to perform the thumb abduction contraction task three times at 20% and 30% MVC, respectively. Sufficient rest was allowed for the subjects to avoid muscle fatigue. Data recorded in this experiment were used for calculating the MU separation vectors and training the muscle force model based on a deep learning network.

Online testing of the force prediction

To implement online testing, the data recording system was linked using a USB cable to a desktop computer with an Intel Core i9 CPU, 32 GB of RAM and an RTX3080Ti GPU for data transferring, recording, online SEMG decomposition and force estimation. Customized software was developed with the Python language using the deep learning framework termed Keras [37]. All of the procedures of the proposed method were implemented in the software. A graphical user interface (GUI) was also developed for interaction between the software and its users via the computer screen, as illustrated in Fig. 1 (c).

After the initialization, the subjects started online testing by repeating the force task mentioned above. For each of the subjects, 10 executions of the 5-s force task were required at 20% and 30% MVC. To meet the demand of online processing, the collected SEMG data stream was divided into a series of temporally overlapping windows with the window length and increment set at 1s and 0.2s, respectively. Real-time decomposition and force prediction were performed in a single 1-s window. During the online testing, the GUI showed a picture demonstration at a fixed interval of 0.2s, which was consistent with the increment of processing windows. The picture demonstration included the historical 5-s SEMG signals, extracted MUSTs, measured force and predicted force to guide task performance. After the test was completed, the overall accuracy was calculated and displayed on the screen. Figure 2 shows the framework of the proposed method.

Fig. 2

Block diagram of the proposed method for estimating muscle force using HD-SEMG data in real time. The framework of the proposed force prediction method contains the training stage (blue block) and the online testing stage (yellow block). The separation vectors were periodically update in the backend thread (red block) to maintain their validity

Training stage

The framework of the proposed method contained a training stage and online testing stage, which is consistent with the two-stage approach of online sEMG decomposition [21]. The pre-processing of the time-consuming computations and training can provide prior knowledge to benefit real-time applications. The training stage included the complete offline prework stage using the offline automatic PFP (APFP) [38] to obtain a series of separation vectors. The MUAP waveform of each MU was estimated following a least squares problem and all the waveforms were generated in 2D form.

Afterwards, the model for muscle force prediction needs to be established and trained. The establishment was two folds that are further explained in the following section.

Transfer of MUAP to twitch force

The twitch force model [39] can establish an electricity-to-force transformation process at MU level. In this model, the muscle force is described to be generated from the twitch forces of a series of activated MUs. Specifically, the twitch force of an activated MU is expressed as the second-order-system [39]:

$$_\left(t\right)=\frac_\bullet t}_}\bullet ^_}}$$

(1)

where \(_\) and \(_\) represent the contraction time and the twitch force amplitude of the \(i\)th MU, respectively. In addition, an inverse power function is used to describe the relationship between \(_\) and \(_\):

$$_=_\bullet _}\right) }^}$$

(2)

where \(_\) is the maximal contraction time with a value set to 90ms and \(c\) is a constant set to 4.2. \(_\) is modeled to be linearly correlated with the peak-to-peak amplitude of the MUAP waveform. Generally, an activated MU usually discharges regularly. For the \(j\)th discharge time of the \(i\)th MU, the twitch force was described as:

$$_\left(t\right)=_\bullet \frac_\bullet t}_}\bullet ^_}}$$

(3)

where \(_\) denotes the gain of the \(j\)th discharge time of the \(i\)th MU. It is defined as the ratio of the force to the firing rate. In this regard, the value of the gain was calculated as:

$$_=\left\1, \frac_}_}\le 0.4\\ k\bullet \frac^_}_}\right)}^}}_}_}}, otherwise\end\right.$$

(5)

where \(k\) is a constant and \(_\) is the inter-spike interval (ISI), which means the time interval between two consecutive discharges of a single MU [40]. \(_/_\) is the normalized stimulus rate.

The twitch force model was employed in each channel separately. To be consistent with the procedures of online decomposition, force prediction was performed on every 200-ms window (corresponding to a 0.2-s step of the 1-s window in the online decomposition, which contains 400 data points at a sampling rate of 2000 Hz). When there were N MUs decomposed in the online stage, the input SEMG feature map in a 400 × 8 × 8 × N data matrix was obtained as a basic sample for network training and testing. Each sample was labelled via the corresponding 100-point normalized force curve.

Fig. 3

Architecture of the neural network used in this study. “conv” refers to a convolutional filter

Muscle force model based on neural network

Figure 3 shows the architecture of the network. One convolution layer was designed to take advantage of mining and characterize the spatial information from the 2D electrode array. The size of the convolutional filter was 3 × 3, the number of filters was 16, and the stride step was 1 × 1. The features obtained from the convolution layer were then processed using the typical LSTM to capture the long-term dependencies of data. The network contained two LSTM layers: one with 64 units and one with 1 unit. The final output of the predicted force (i.e. 400 × 1) was obtained from a fully connected layer. The whole network was trained with the Adam optimizer with a learning rate of 0.001. The root mean squared deviation (RMSD) was chosen as the loss function of the neural network and was calculated as:

$$RMSD=\sqrt_^\left(t\right)-F(t\left)\right]}^}}\times 100\%$$

(6)

where \(n\) denotes the number of the samples (i.e. 400). \(\widehat\left(t\right)\) and \(F\left(t\right)\) represent the predicted force and the measured force, respectively.

To this end, the force prediction model was established. The initial 15-s SEMG signals were decomposed offline to obtain the MU activities and the twitch forces. These decomposition results were divided into a training set and a validation set at a ratio of 2:1. Then, the neural network was trained with the labels of measured force with a batch size of 32 and the network was trained for 150 epochs. This training stage is equivalent to the model training in the deep-learning approach.

Online force prediction

During the online testing stage, the MUSTs were identified in short time windows from the extended and whitened SEMG data stream in real time, as shown in Fig. 4. The procedures of the MUST extraction kept consistent with our previous work using the successive multi-threshold Otsu algorithm [21]. In a single 1-s window, the MU firing events from the overlapping 0.8s were used to track the same MUs to ensure the continuity of decomposition results. The MU firing events from the last 0.2s were applied to estimate the force. More details of the processing and the corresponding parameters can be found in our previous study [21]. The online PFP method used in this study had the same settings reported in the same study [21].

Fig. 4

The illustration of the online force prediction process using the proposed method

After executing the online decomposition of SEMG signals, the MUSTs were continuously extracted and then transferred to multi-channel twitch force trains. These twitch force trains were input into trained neural network and the estimated muscle forces over windows were connected to form the resultant force.

To improve the performance of muscle force prediction, we designed a backend thread to update the MU separation vectors in the 5-second EMG segment at intervals of 10 s. This was done to adaptively update the MU separation vectors and maintain the validity of the vectors to precisely decode neural information. More specifically, the constrained FastICA was used to continuously track the same MUs with high precision in the backend thread. The spike trains of the identified MUs were used as constraints to drive the FastICA algorithm to converge rapidly (The detailed calculation steps of this constrained FastICA algorithm are given in previous studies [38, 41]). Therefore, for each MU, the separation vector can be effectively updated, and all the potential firing errors are believed to be corrected. After this update procedure, the original MU separation vectors were replaced by the updated version using the described procedure. The frontend thread applied the newly updated vectors to extract MUSTs and estimate muscle force in real time, processing in parallel with the backend thread.

Performance evaluationPerformance metrics

The performance of force prediction was quantified by RMSD and the coefficient of determination (R2). For each of the different force levels, the RMSD and R2 values were averaged across subjects to represent the overall performance. In addition, we calculated the matching rates (MRs) of each individual MUST decomposed in the online decomposition stage with the ground-truth reference (the offline decomposition results). To better understand the effect of the backend update thread on the decomposition precision, the EMG data used in online testing stage were additionally decomposed by an offline APFP method [38]. The MR is defined as:

where \(_\) and \(_\) are the number of spikes of the two spike trains decomposed online and offline, respectively, and \(_\) is the number of common spikes. MR and the rate of agreement (RoA) are both commonly used criteria to assess the degree of matching of two spike trains and can effectively evaluate the decomposition performance [42]. Each MUST decomposed online was compared with each of the reference MUSTs, and the MUST achieving the maximum MR was selected from the reference to pair with the online decomposed MUST. This maximum MR was set as the decomposition accuracy for that MU.

To evaluate the real-time performance of the proposed method, we also calculated the time delay for processing the EMG data in a single 1-s time window, which was averaged all windows and all subjects.

Comparison methods

To demonstrate the outperformance of our method more comprehensively, two evaluations were performed. First, the RMSD and R2 values were evaluated with and without the backend update thread. In other words, the online force prediction without the update procedures could only use the MU information initialized in the training stage. This comparison method was denoted as the “no update” method, and the proposed method was termed as the “update” method.

Second, some conventional methods for muscle force prediction were also applied for comparison. One previously reported method using microscopic information was also adopted. This method employed the MU firing rate (FR) of the decomposed MUSTs to estimate force through a two-order polynomial regression model (termed as FR method) [22]. The FR in a single time window was calculated with a length of 250ms and increments of 50ms.

In addition, the RMS was also selected as a representative macroscopic feature for the conventional force prediction method based on EMG amplitude (termed as RMS method). The calculation of RMS was consistent with that of the FR, and the same two-order polynomial regression model was applied. Other settings of the comparison methods and the proposed method remained consistent or were fine-tuned for optimal performance.

Statistical analysis

In order to evaluate the effect of the update procedure in the backend thread, two two-way ANOVAs were conducted on the RMSD percentage and R2, with both the update procedure (two levels: update and no update) and the method (three levels: the proposed method, FR method and RMS method) considered as the within-subject factors. The level of significant difference was set as p < 0.05. All statistical analyses were performed with SPSS software (ver. 22.0, SPSS Inc. Chicago, IL, USA).