Grasp force estimation from the transient EMG using high-density surface recordings

Twelve able-bodied participants (aged 27.6  ±  2.96 years, seven males, all right-handed) without any history of neuromuscular disorders participated to this study. Written informed consent in accordance with the Declaration of Helsinki was obtained before conducting the experiments from each participant. This study was approved by the local ethical committee of the Scuola Superiore Sant'Anna, Pisa, Italy (request no. 02/2017). The methods were carried out in accordance with the approved guidelines. The experimental setup consisted of a test-object instrumented with force sensors, an object stand with embedded load cells, a chair with an instrumented armrest, an HD-EMG recording system and a PC (figure 1(a)).

Zoom In Zoom Out Reset image size

The test-object was composed of a handle fixed on top of a box (figure 1(b)). The handle was equipped with a 6-axis load cell (NANO17, ATI Industrial Automation Inc., Apex, USA, 0–70 N) in order to measure the GF. The box was used to contain objects with different weights in order to modify the weight of the test-object. The object stand contained, within its base, a piezo-resistive force sensor (FSG series, Honeywell Sensotec, Columbus, USA, 0–15 N) to measure the load force on the test-object until lift-off. The presence/absence of the hand on the armrest, which corresponded to the starting position of the experimental task, was detected by a force sensing resistor.

HD-EMG signals were recorded from the participant's forearms, using three 64-electrode surface arrays (8  ×  8) with 10 mm inter-electrode distance (ELSCH064NM3, Spes Medica S.r.l., Genoa, Italy) connected to an EMG amplifier (EMG-USB2+, OT Bioelettronica, Turin, Italy). Two of the arrays were placed on the proximal half of the forearm, one covering each of the anterior and posterior compartments; the third array was placed in the distal half approximately above the FPL (flexor pollicis longus), as identified through tactile palpation. Two reference electrodes were placed at the wrist level and between the matrices following the manufacturer's guidelines (figure 1(c)).

All sensors and EMG channels were synchronously acquired (10 240 Hz frequency) and band-pass filtered (10 Hz–4.4 kHz). The data was stored in the PC memory for off-line analysis. The PC was also used to guide the participants through the experiment (figure 1(a)).

The participants were asked to perform a pick and lift task, while comfortably sitting on a chair in front of the experimental platform, on a table. In particular, they were asked to execute the following sequence from the starting position: (i) move their hand to reach the instrumented object, (ii) grasp it using a three-digital grasp (thumb, index and middle fingers), (iii) lift it 10 cm and (iv) wait 2 s before (v) replacing it on the object stand and (vi) returning their hand back to the starting position. In particular the latter consisted in the arm resting on the armrest with the hand at the same level of the test-object handle, placed 15 cm away from it (figure 1(a)). A visual cue displayed on the monitor marked the start of steps (i) and (v). Between lifts, when the hand was in the starting position, to increase grasp repeatability, participants were asked to maintain the same three-digital grasping posture used during the lifts, in a relaxed way.

Each participant performed five series of 15 lifts for each of four weights (250 g, 500 g, 750 g and 1 kg), for a total of 300 lifts. After each series, the GF necessary to lift the object was varied by filling the test-object box with a different weight (the order of the weights was randomized across participants). The participants could rest between series if they wished to. Six performed the experiment with their right hand and the others with their left hand.

At the beginning and at the end of the experiment, the GF at maximum voluntary contraction (GF_MVC, recorded by the load cell) and its corresponding EMG value (MVC), were recorded from each participant. In particular the participants were asked to perform three lifts with maximum force [34,35]. The maximum GF recorded across the six trials was used as GF_MVC; the maximum root mean square (RMS) across EMG channels in the array, calculated over a period of 60 ms, was used as MVC.

The data was processed offline using MATLAB R2016b (The Mathworks, Natick, USA), with the objective of finding key features in the transient EMG to predict the GF necessary to maintain the test-object in the air—we called this final GF. The data from all sensors and EMG electrodes was down-sampled to 2048 Hz. The data from the sensors was used to segment the recordings in individual trials and in task phases, whereas those from the EMG electrodes and the test-object were used to build the final GF prediction algorithms.

Considering that during the first three lifts of an unknown object humans likely adapt their motor control [32], this data was not used for the analysis (leaving 12 lifts per series, figure 1(d)). From each of the remaining 240 lifts per participant (12 lifts  ×  5 series  ×  4 weights), the beginning of the transient was identified by the onsets of the GF and of the EMG signal. The GF onset was identified by applying a simple thresholding operation on the GF signal. In other words, the onset was identified as the first moment when the GF exceeded a predefined threshold. The threshold was set to three standard deviations of the signal noise amplitude, evaluated as the signal standard deviation when no force was applied to the load cell. The EMG onset was also detected, using the mean absolute value (MAV) signal from the four central channels of the FPL array, following the method proposed by Kanitz et al [31]. In this work, because of the physiological basis of the EMG signal [36], we preferred to work in conditions in which the variability of EMG onset detection did not influence the outcomes. Thus, we based the whole analysis using the GF onset. The EMG onset was only calculated for comparison purposes.

EMG signals were filtered with a 4th order band-pass Butterworth filter (20 Hz–500 Hz) [37]; then, neighboring electrodes along the direction of the muscle were differentiated for each matrix, converting them into 7  ×  8 bipolar signals that were normalized by their corresponding MVC per participant. Similarly, the GF was normalized by the GF_MVC. These normalizations were made in order to reduce the variability between participants and to get an estimate of the muscle activation out of EMG data [37]. Ten features including MAV, waveform lenght (WL), Logarithm of Variance (LogVar), square root of the variance (2nd order vOrder), signal energy (SE), Hjorths features, time-dependent power spectrum descriptors (TDPSD) and the time-differential of the MAV (dMAV) were extracted from the EMG signals using sliding windows of 60 ms (named feature window length, FWL) with a 10 ms step. These features were selected because they are commonly used in literature (MAV, WL, LogVar, vOrder and SE [38, 39]) or because they already proved successful in classifying transient EMG signals (Hjorths [40] and TDPSD [41]). An additional feature, the differential of the MAV (dMAV), was also calculated.

The actual GF applied on the test-object while being held in the air (corresponding to the actual final GF) was also calculated for each trial as the mean of the plateau of the measured GF (within 400 ms and 1000 ms after the GF onset). The final GF was then estimated (or anticipated) from the features by training a regularized linear regression. In order to determine the earliest time at which the final GF could be accurately predicted (or, in other words, the minimum amount of transient signal needed) the process was repeated with windows of increasing length (named TW). The length of the TW ranged from 1 to 75 samples (figure 2, upper panel). As feature samples are separated by 10 ms, this corresponded to a portion of data embracing a single instant (the onset) up to 750 ms. In turn, depending on the TW length, the process included only information from the transient or from both the transient and steady state (figure 2, upper panel).

Zoom In Zoom Out Reset image size

Specifically, to tune the regression parameters and hyperparameters, the TWs extracted from the 240 trials available from each participant were split into training sets (50%), development sets (30%) and test sets (20%). Additionally, the trials that contained the MVC were included in the training set. Then, for each TW, the prediction of the final GF was performed after a selection of channels and/or features using a reduction method. The aim of the reduction method was to obtain sets of 4, 8 or 16 channels and compare the results with those achieved from all HD channels (i.e. 7  ×  8  ×  3  =  168 channels) to determine the minimum amount of information necessary for the prediction of the target GF. As three reduction methods were compared, this resulted in a total of 12 configurations (4 sets of channels  ×  3 reduction methods). The three methods used were:

• 1.  
  'Fix-Ridge': the final GF was estimated using a regularized linear regression (with ridge regularization parameter ${{\lambda }_{{\rm Reg}}}$$\in$ {0, 0.001, 0.01, 0.1, 1, 10, 100, 1000}) with a predefined selection of channels (roughly those in the center of the matrices—figure 3) and all the extracted features. This method does not discard any feature.
• 2.  
  'Fix-EN': solution obtained using the same channels as in the Fix-Ridge method but performing features selection through elastic nets analysis¹ [42–44]. A weight parameter α  =  0.4 was chosen, whereas different values for the regularization parameter ${{\lambda }_{{\rm EN}}}$(being λ_EN the same set of as ${{\lambda }_{{\rm Reg}}}$) were tested. No limitation on the number and type of features to be retained by the algorithm was included. As a result, the Fix-EN solution could end up with a number of channels lower than the initial value, in the case all features from a channel were discarded. The prediction of the final GF was obtained from the coefficients of the elastic nets.
• 3.  
  'LassoG-EN': solution based on the lasso-group algorithm [45] for automatic selection of the channels, followed by elastic nets for features selection (using the same parameters of the Fix-EN method).

Zoom In Zoom Out Reset image size

For all methods, the algorithms were trained with the training set and the selection of the regularization parameters was based on the performance on the development set while the final assessment was performed on the test set (figure 2, lower panel). Before the training, all features were normalized using the data from the training set. Notably, the chosen step (10 ms), FWL (60 ms) and α (0.4) were identified after preliminary tests using the same dataset split, in order to limit the complexity of the testing. The detailed results of such preliminary tests are omitted here for sake of conciseness².

The outcomes of the Fix-EN and LassoG-EN reduction methods were analyzed in terms of frequencies and timing of the selected features to get insights about what information was retained after their application.

The R² and the mean AE calculated as a percentage of GF_MVC were used as the metrics to evaluate the estimation of the final GF. Hence, the metrics were also used to identify the TW that optimized the performance of the regression. Specifically, for each configuration and for a subset of TWs (i.e. 300, 350, 400, 450, 500, 550, 600 ms), the difference in AE was analyzed through a Friedman test using R (R Foundation for Statistical Computing, Vienna, Austria). This was followed by post-hoc pairwise comparisons with Bonferroni correction. The TW that exhibited statistically better performance (dubbed TW^*) was further analyzed and the frequency of the selected channels was assessed (Fix-EN and LassoG-EN methods).

A two-way repeated measures ANOVA (factors: reduction method and number of channels) with pairwise comparisons was used to identify the overall best configuration. In this case, a single TW, common to all configurations, was used for the comparison. This was chosen as the longest TW among all the TW^* selected through the Friedman test. A significance level of p   =  0.05 was used throughout the statistical analysis.

On the best configuration selected through the ANOVA, a test comparing the final GF estimation in two additional windows, namely ${{W}_{{\rm TR}}}$ and ${{W}_{{\rm SS}}}$, was conducted. The two windows were dimensioned a posteriori considering the average duration of the GF transient across participants. This was done to evaluate the differences in estimating the GF using the transient only with respect to the steady-state only.

The total number of features selected by the Fix-EN and LassoG-EN reduction methods ranged between a couple for small TWs and four channels to a maximum of 150 for larger ones and 168 channels (figure 4(a)—only the case of Fix-EN is displayed). This makes a reduction of the model complexity between 90% and 99%, considering that the initial number of features ranged from 40 (4 channels  ×  10 features  ×  1 time sample) to 126 000 (168 channels  ×  10 features  ×  75 time samples). Overall the Fix-EN and LassoG-EN mostly retained the following features: WL, TDPSD S, TDPSD IF and dMAV (figure 4(a)). This preference was more visible for shorter TWs, while it attenuated the end of the transient. In addition, the more informative time samples were usually the most recent ones (i.e. the ones at the end of the TW), with the very last time sample in the TW being selected more frequently (figure 4(b)—only the case of Fix-EN is displayed).

Zoom In Zoom Out Reset image size

In all configurations, the prediction algorithm was able to predict the final grasp force with high accuracy (representative example of the Fix-Ridge reduction method in figure 5). Specifically, the prediction improved according with the TW length, until a stable performance was reached around 400 ms after the GF onset. Notably, as there is little to no EMG activation at the beginning of the transient (TW length close to zero), the performance for the shorter TWs indicates the random guess of the linear regression algorithm. Indeed, under large uncertainty linear regression algorithms return the mean of the target (i.e. the final GF) used during their training [40]. For this subject such value was 6.2 N or 14.7 % GF_MVC. Notably, by increasing the number of channels, the estimated GF departs from the random guess even when the TW contains a single sample. The best AE (averaged along TWs length between 450 ms and 750 ms) of this subject was 2.26 % GF_MVC.

Zoom In Zoom Out Reset image size

These results were consistent among subjects and reduction methods. Specifically, all the assessed configurations showed a similar behavior in terms of R² and AE. The performance improved (R² increased and AE decreased) with time from the GF onset on, reaching a plateau for TWs of around 450 ms (figure 6). More quantitatively, the R²/AE increased/decreased from 0.1/5% for the shortest TW to 0.8/2% for the longest one (750 ms). When only information from the transient phase was available to the algorithm (i.e. TW  =  330 ms) the R²/AE reached a value close to the plateau: 0.67 (0.25)/2.52 (1.76)%. Notably, as there is little to no EMG activation at the beginning of the transient, the performance for the shorter TWs indicates the random guess of the linear regression algorithm. Indeed, under large uncertainty linear regression algorithms return the mean of the target (i.e. the final GF) used during their training [46]. In our problem such value is 7.6 N resulting in an AE of 5.7 % GF_MVC.

Zoom In Zoom Out Reset image size

The performance of all reduction methods improved slightly as the number of channels increased (figure 6). The TW^* were found to be mostly independent from the tested configuration, being either 450 ms or 500 ms (figure 6). Therefore, a TW of 500 ms was considered as the best one in the following analysis.

The channels selection frequency (figure 7) showed that Fix-EN, on average, did not discard any channel completely. The most selected channels were usually those from the matrix placed on the flexor muscles. For the LassoG-EN reduction method, the algorithm typically picked different channels for different participants. By definition the 168 channels case reported exactly the same results for both reduction methods.

Zoom In Zoom Out Reset image size

The statistical analysis of the studentized residuals showed that there was normality (Shapiro-Wilk test), no outliers (no studentized residuals greater than  ±  3 standard deviations) and no sphericity (Mauchlys test p   =  0.002) on the data. Therefore, Greenhouse-Geisser corrections were used. The two-way repeated measures ANOVA showed that the number of channels significantly affected the performance (F(1.83, 18.25)  =  16.55, p  <  0.001, η²  =  0.71) while the reduction method did not (F(1.35,13.46)  =  4.02, p   =  0.051, η²  =  0.35). Additionally, there was no interaction between the two factors (F(2.4, 23.96)  =  2.12, p   =  0.135, η²  =  0.59). Therefore, pairwise comparison was performed on the number of channels only (table 2). This showed that the performance improved when moving from 4 or 8 channels to 168 channels. On the contrary, there was no statistical difference between the 16 and 168 channels configurations, for all reduction methods.

Following these results, the configurations with 16 channels proved the ones with the best tradeoff between complexity and performance. Additionally, as there was no statistical difference between reduction methods, Fix-EN was chosen as the best one due to its reduced computational cost. Thus Fix-EN with a TW of 500 ms was identified as a representative overall best configuration. It allowed for an AE of 1.99 (1.04) % GF_MVC which corresponds to 0.64 (0.24)N and to an R² of 0.82 (0.16).

By analyzing the data from this configuration for each weight separately, the performance in GF estimation seemed to correlate negatively with the weight. Specifically, the mean AE across participants increased from 1.6 (1.05) % GF_MVC at 250 g to 1.7 (1.12) % GF_MVC at 500 g, 2.4 (1.42) % GF_MVC at 750 g and 2.5 (1.23) % GF_MVC at 1 kg. Additionally, when evaluating the performance separately for each participant, we found that the AE ranged between 0.71 (1.04) % GF_MVC and 4.10 (3.55) % GF_MVC.

Finally, to evaluate the differences in performance between transient and steady-state phases, we fixed W_TR to the interval 0–330 ms and W_SS to 330–660 ms after the onset. This resulted in a median AE of 2.52 (1.76) % GF_MVC and 1.75 (1.06) % GF_MVC for W_TR and W_SS, respectively.