Computational approaches to decode grasping force and velocity level in upper-limb amputee from intraneural peripheral signals

The subject was a 48 years right-handed female with a trans-radial amputation of the distal one-third of the left forearm, which occurred 23 years before the study. Four TIMEs, with 14 active sites each, were implanted in the subject's median (two interfaces, M1 and M2) and ulnar (two interfaces, U1 and U2) nerves. Electrodes M1 and U1 were implanted in the proximal nerve segments, while electrodes M2 and U2 were more distal (see [32] for details) (figure 1(A)).

Zoom In Zoom Out Reset image size

Ethical approval was obtained from the Institutional Ethics Committees of Policlinic A. Gemelli at the Catholic University (Rome). Written informed consent was obtained from all subjects before the study. The protocol was approved by the Italian Ministry of Health.

The subject was asked to perform two tasks based on a power grasp movement with the phantom hand following instructions on a screen. In the first task (isotonic) the subject was asked to grasp an object with one of three different levels of strength (low, medium, high), while in the second task (isokinetic) she was asked to perform the grasp at one of three different speeds (slow, medium, fast). The requested level was indicated on the screen by a custom-made interface. Each task (isokinetic/isotonic) was repeated ten times for each level. The subject was also asked to reproduce the same movement with the intact hand on a dynamometer to check for her attention. During the isokinetic task, each repetition lasted for 5 s, 3 s of movement followed by rest. Total recording time for each session was approx. 150 s for completing the isotonic task, plus a few seconds before and after the trials for recording baseline activity. As regards the isokinetic task, each repetition lasted 3, 2 or 1 s, depending on the requested speed level (slow, medium and fast, respectively), followed by rest (figure 1(B)), for a total of approx. 110 s to complete the task. The whole protocol was repeated in three different days (16, 17 and 20 d after electrode implantation). The subject completed the whole protocol (and data collection was full) on day 17. The number of levels performed by the subject during the isotonic task on the other two days was two instead of three.

Muscle activity (EMG). During the experiment, four surface electromyographic (sEMG) signals were acquired at 2 kHz through a Grapevine system (Ripple, LLC) from the residual muscles, two from the dorsal (digit extensor extrinsic muscle—EMG1, EMG2) and two from the ventral (digit flexor extrinsic muscle—EMG3, EMG4) sides of the forearm of the subject [21] for details). EMG4 was discarded from the analysis because of low SNR due to acquisition problems. Signals were digitally filtered with a 4th order passband (15–375 Hz) Butterworth IIR filter, and a notch filter was applied at 50, 100 and 150 Hz to remove the power line interference. EMG signal envelopes were then extracted offline by: high-pass filtering at 3 Hz (4th order Butterworth), rectification, and a low-pass filtering at 60 Hz (4th order Butterworth).

Neural recordings (ENG). Four electrodes (56 active sites) were recorded simultaneously (using Grapevine neural Interface System, Ripple, LLC), and digitally sampled at 30 kHz. The raw extracellular neural signals collected with TIMEs were band-pass filtered (300–7500 Hz) in MATLAB (Mathworks Inc., USA) to extract multi-unit activity (MUA). Spike detection and sorting were then performed by using Plexon Offline Sorter (Plexon Inc. Dallas, TX, USA) to extract single-unit activity (SUA). Spiking events were detected by voltage threshold crossing [26–29], with threshold = 3.5σ_noise, and a 2 ms waveform sample was acquired around the time of threshold crossing. The individual waveform samples were then aligned on their most negative point. Initially, single-unit clusters were manually defined based on waveform similarity in the 2D principal component space (PC1, PC2) and their waveform template was extracted [39–41]. Then, a sorting algorithm based on template matching was applied on all data points [42]. Quality of separation was determined based on normality statistics on the 2D PCs space, and the presence of a clear refractory period [43]. Firing rates were then calculated with NeuroExplorer (v4, Nex Technologies, USA) and stored, together with spike times, for each neuron.

To reduce data dimensionality and select the most informative neural signals, we took advantage of tools from information theory [44–46]. The mutual information, I(S;R), between a given set of states, S, and a set of neural responses, R, is defined as follows:

$$\begin{equation}I\left( {S;R} \right) = \sum\limits_{\text{s}} {\text{P}} \left( s \right)\sum\limits_{\text{r}} {\text{P}} \left( {r{\text{|}}s} \right){\text{lo}}{{\text{g}}_2}\frac{{{\text{P}}\left( {r{\text{|}}s} \right)}}{{{\text{P}}\left( r \right)}}\end{equation} \tag{ 1 }$$

where P(s) is the probability of the state s, P(r|s) is the probability of observing in the same epoch a value r of the response given a state s, and P(r) is the probability of observing across all epochs a value r of the response. Here, the set of possible states S was given by the different values that EMG envelopes assume in specific time windows of the gripping task (figure 2 step 1); specifically, we defined two windows: one related to the movement (labelled ON), and the other related to rest (labelled OFF). We identified these windows by looking at the mean tracks of the EMG averaged across trials performing the same task. For each state occurrence, the neural response R of each neuron was quantified as the median value of the firing rate in the specified time window. Then, the mutual information, I(S;R), between the states S and the neural responses R was calculated as in [47], using the bias reduction corrections reported in [48] and by means of the Matlab Information Breakdown toolbox presented in [49]. We discarded from further analysis all the neurons whose responses did not contain significant information about the EMG state (p < 0.05 after bootstrap correction with N = 500) (step 2). Then, on the remaining subset of neurons, we computed the 'joint information' carried simultaneously by neurons couples (step 3). The joint information carried by two paired neural responses about a given state is

$$\begin{equation}I\left( {S;{R_1},{R_2}} \right) = \sum\limits_s {\text{P}} \left( s \right)\sum\limits_{{r_{1,2}}} {\text{P}} \left( {{r_1},{r_2}{\text{|}}s} \right){\text{lo}}{{\text{g}}_2}\frac{{{\text{P}}\left( {{r_1},{r_2}{\text{|}}s} \right)}}{{{\text{P}}\left( {{r_1},{r_2}} \right)}}\end{equation} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

where P(r₁,r₂|s) is the joint probability of observing values r₁ and r₂ for the two responses given the state s. Given the higher dimensionality of the code in joint information, in this case the bias was further reduced through shuffling correction [49].

Finally, to select the best set of neural units, we sorted the mutual and the joint information and we chose the most informative units or couples of units (step 4), i.e. with I above an ad hoc chosen threshold (the 80% of the maximum value).

We used the open source toolbox nSTAT (neural spike train data analysis toolbox) developed for MATLAB [50] to estimate force or velocity levels of gripping movements of the phantom hand from the recorded neural spikes. nSTAT uses the point process-generalized linear Model (PP-GLM) framework to analyse the relationship between spiking activity and relevant covariates (i.e. signals which are informative of the different force or velocity levels specified by the subject). Specifically, neural spiking activity can be modelled as a PP due to the all-or none nature of the action potential (for a more extended explanation see [51]). Given ${H_t}$, the spiking history from 0 up to time t, and $N\left( t \right)$, a counting process containing the sum of all the spikes up to time t, a PP can be characterized by its conditional intensity function (CIF) [51] as

$$\begin{equation}\lambda (t|{H_t}) = {\text{ }}\frac{{{\text{Pr}}(N\left( {t + \Delta } \right) - N\left( t \right) = 1|{H_t})}}{\Delta },{\text{ }}\end{equation} \tag{ 3 }$$

which can be considered a generalization of the firing rate for a non-homogenous Poisson process. The CIF at time t depends on the value of the covariates up to that time. Assuming that the relationship between $\lambda$and its covariates follow a distribution in the exponential family, the GLM framework can be used to fit statistical models for the CIF [52]. In particular, we applied Poisson regression models as our likelihood criteria to estimate CIF model parameters by using the spiking history of the neuronal population.

We started with a collection of informative neurons n = 1, ... N, whose spiking times represented the PP observations and we defined our covariates as the three EMG envelopes extracted from the recorded muscle signals (see section 2.2). nSTAT was then used to encode a CIF model for each neuron independently by fitting its spiking activity to a non-homogeneous Poisson model.

The estimated CIF models were then used to decode each covariate, along with the related 95% confidence interval, by using a point process adaptive filter (PP) [53].

The final PP model parameters were estimated with a ten-fold cross-validation: iteratively, nine trials for each force or velocity level trials randomly selected for training, and the remaining three trials used for testing. Performance was quantified computing the R² between the actual EMG envelope and the predicted signal for the held-out trials of the test sets (figure 3).

Zoom In Zoom Out Reset image size

The performance of the decoding method here proposed was compared with the results obtained by a KF, a very well-established algorithm to decode movement from the firing rates of a population of cells [54–56]. We assumed a linear relationship between the state, the envelopes of the recorded EMGs, and the observations, the firing rates of the informative neurons, plus (without overlapping) Gaussian noise. Similarly, the transition matrix of the state was assumed to be linear plus an additional Gaussian noise [57]. As well as the PP algorithm, also the KF is based on two steps [58]: (a) definition of a generative model based on neural firing rates and estimation of all the parameters by a least square algorithm; (b) estimation of the EMG envelopes at each time step (i.e. the a priori/a posteriori state). KF parameters were estimated with a ten-fold cross-validation as explained in the previous section.

To decode force or velocity levels from input data, a SVM with a quadratic kernel (qSVM, ten-cross-validation) classifier was applied and trained using the one-vs-one strategy. A five-fold random cross-validation procedure was carried out to evaluate parameters and to ensure the unbiased correctness of the classification performances. All input signals were binned at 20 ms. A 100 ms sliding window of the input signal was then used for training and testing the classifier. To explore the best window choices, we considered multiple starting-time points after the external trigger (cue signal) provided by the GUI (from 0 to 500 ms, 100 ms step). These windows were labelled and associated to a specific force/velocity level.

Offline performance was evaluated by computing the Accuracy, i.e. the percentage of correct predictions per class, displayed as confusion matrices. Overall Classification Accuracy was then defined as the mean value of the Accuracy values per class (corresponding to the mean value of the diagonal in the confusion matrix).

A crucial step of the decoding process was the reduction of the dimensionality of the set of the recorded putative single units (figure 2). In the isotonic experiment, spike-sorting algorithms yielded a total of 86 putative single units out of 56 recording channels. We considered ON the window between 0.6 and 1.1 s after the trigger, and OFF the window between 0.1 and 0.6 s before the trigger (figure 2(A)). Information theory highlighted the electrode M1, placed in the median nerve, as the most informative of the four implanted TIMEs (figure 2(B)). The number of informative neurons was further reduced to n = 17, taking all the neurons that showed information (single or couples of neurons) greater than the 80% of the maximum information provided (>0.48) (figures 2(C) and (D)).

We repeated the same procedure for the isokinetic experiment, considering three different levels of velocity (figure 2(E)). With spike sorting, we found 120 neurons in total from 56 recording channels. The ON/OFF windows were adapted accordingly (ON = [0.4 0.9] seconds after the trigger, and OFF = [0.1 0.6] seconds before the trigger. In contrast with the isotonic recordings, the most informative electrode is the one implanted in the ulnar nerve (U1) in the case of velocity tasks (figure 2(F)). We reduced the dimensionality to 19 considering 0.45 (the 80% of the maximum value) as minimum value for Mutual and Joint Information between neuron activity and EMG state (see section 2) for neurons selection (figures 2(G) and (H)).

To establish a benchmark for ENG analysis, we first investigated to which extent it was possible to decode force and velocity levels directly from the muscular activity. This was considered as the reference for the following methods here developed. We used as input for the classifier the values assumed by the EMG envelopes in specific time windows as inputs for the SVM classifier. Two parameters of the temporal window could be varied: window start and window length, which represent the initial time point and duration, respectively. The overall accuracy is reported in the first panels of figures 3(A) and (B), for the isotonic and isokinetic tasks respectively. To summarize the results, we reported the confusion matrix (second panels) obtained by averaging nine individual confusion matrices computed over the optimal temporal frame indicated by the red rectangle (window start from 0 to 300 ms; window length from 100 to 300 ms). We obtained an average accuracy of 86.4% for the force tasks and 75.6% for velocity tasks.

We repeated the same procedure by directly using as input for the SVM the firing rate traces of the selected neurons to classify force and velocity levels (17 and 19 neurons). As with the EMGs, we considered different window start point and length parameters. In this case, the accuracy averaged on the optimal temporal window is of 72.2% for the force tasks and 73.8% for velocity tasks (figures 3(C) and (D)).

To assess the possibility of improving the performance of movement reconstruction from neural activity, we added an intermediate step before the classifier. The encoding phase was added to model both the linear KF and the non-linear PP using a theoretical square-wave signal based on the external trigger as covariate for building the model. We then reconstructed the ideal signal with these two approaches. In this case, results were far from being satisfying as showed in figure 4. Indeed, the isotonic tasks reported an average accuracy of 41.6% and 37.2% for PP and KF (figures 4(A) and (B)) respectively obtained by averaging the confusion matrixes in the selected temporal frame. While performance with PP (figure 4(A)) is significantly better than chance (t-test, p < 0.01), performance with KF (figure 4(B)) is not (t-test, p > 0.05). Similarly, the average accuracy for the isokinetic tasks is 34.5% (t-test, p > 0.25) and 36.3% (t-test, p < 0.05) for PP and KF respectively (figures 4(C) and (D)).

Zoom In Zoom Out Reset image size

We considered then the possibility of using the recorded EMG signals as the covariates for the encoding phase to model the KF and the PP.

We reported in figure 5 the EMG envelopes and the EMG traces reconstructed using both the Point Process algorithm (PP EMG) and the Kalman Filter (KF EMG) during the isotonic tasks (figure 5(A)) and the isokinetic tasks (figure 5(C)) in a representative window. We computed the R² for each of the three EMGs to quantify the similarity between the original EMG envelope and the related reconstruction. Results of the isotonic tasks are reported in figure 5(B). Both the PP model and the KF reconstruction well correlate with the original EMG envelopes, in particular EMG2, reaching a coefficient of 0.69 for PP and 0.73 for KF (no statistical difference between PP and KF, p > 0.295 Kruskal–Wallis test—KWT). On the contrary, the isokinetic tasks are better correlated with the KF traces than the PP traces (R² equal to 0.59, 0.39 and 0.58 for EMG1, EMG2 and EMG3 respectively). However, only EMG1 reported a statistical difference between (p < 0.05 KWT) (figure 5(D)).

Zoom In Zoom Out Reset image size

Finally, we decoded the three classes corresponding to the force or velocity levels performed by the subject from the reconstructed EMG signals.

Overall classification accuracy for the isotonic experiment is showed in figures 6(A) and (B) using the EMG reconstructed with PP algorithm (A, PP EMG) and the EMG reconstructed with KF (B, KF EMG). KF performed systematically worse than PP, which provided an accuracy comparable to the method based on original EMG envelopes (results reported in figure 3). The PP-reconstructed yielded optimal performance in a temporal frame defined by values of window start up to 300 ms, and values of window length less than 300 ms. Figure 6(B) shows the classification results for each force level, averaged over this optimal temporal frame, in the form of confusion matrices. Accuracy obtained with the PP EMGs (85.4%) overcame the performance of the KF EMGs (69.0%).

Zoom In Zoom Out Reset image size

Results from the isotonic task are reported in figures 6(C) and (D). As in the previous paragraph, overall classification accuracy is a function of these window start and length. As in the previous case, performance obtained with the PP EMGs was better with respect to the decoding results obtained with the KF EMGs. This is more evident if looking at the performance for each velocity level obtained by averaging confusion matrices computed in a temporal frame defined by values of window start up to 300 ms, and values of window length less than 300 ms (figures 6(C) and (D) second panels). Specifically, the averaged accuracy obtained with the PP EMGs is 75.6% while the accuracy obtained with the KF EMGs reached the 68.3%.

To investigate the (short-term) stability of the proposed strategy, we investigated whether it could yield similar results in other recording sessions from the same subject (day 16 and day 20). These sessions were analyzed as the one of day 17. We processed all neural signals from each TIME to extract neural unit activities. A comparable number of units could be sorted in each recording day (70–96 for the isotonic task and 111–148 for the isokinetic task; figure 7, panels (A) and (D)). We then applied dimensionality reduction based on information theory (see sections 2 and 2.3). A similar number of informative units (80% of the maximum I value) was found across all days (figure 7, panels (B) and (E)), even if the sets were different as shown from the distribution across nerves. Force and velocity levels were decoded using the classifier that produced the best accuracy results among those tested, i.e. the model-based classifier (signals reconstructed by the Point Process algorithm) using EMG envelopes as covariates. Decoding performance was: (a) stable across all the days, (b) above chance levels, and (c) comparable to the results of the EMG-based reference classifier (figure 7, panels (C) and (F)). For the isotonic task, results of the stability analysis were reported for only two levels (i.e. 50% chance level), since always performed by the subject across days.

Zoom In Zoom Out Reset image size

This study explored and compared several approaches to decode desired grasping force or velocity levels. Specifically, we compared: (a) a classifier that uses EMG signals recorded from the residual forearm muscles as input data (figures 3(A) and (B); (b) classifiers that use peripheral neural signals—in the form of single-unit firing rates as input data (figures 3(C) and (D)); (c) classifiers that use neural signals—in the form of both spike times or firing rates—to reconstruct covariates (the grasp cue level or EMG envelopes) and use the reconstructed signals as input data (figures 4–6).

The results obtained by the classifier (a) were used as a reference performance of a non-invasive approach, such as those already presented in the literature [11, 12], to be compared with that of more invasive approaches based on peripheral neural signals, i.e. classifiers (b) and (c). Classifier (b) was model-free, and decoded motion intention directly from single-unit firing rates. Classifiers (c) both required an initial modelling step to reconstruct covariates from neural signals. The SVM classifier, commonly used in other works [32, 37, 60], was then applied as a final step.

In classifier (c) both the KF and the PP adaptive filter were trained for estimating a covariate—the grasping force/velocity cue level presented prior to each task execution, or the EMG envelopes—from neural data. The reconstructed signals were then used as input data for a SVM classifier to decode force/velocity levels. While the KF is a well-established method in peripheral decoding studies to infer kinematic variables from firing rates [55], the PP approach is a less explored method based on neuronal spike times, mainly applied on cortical signals [50, 61] or type I afferent peripheral signals [62]. Although PP is technically more challenging than KF, it has the advantage that instant and not averaged (over time) neuronal activity is considered, as with firing rates in KF. In fact, the PP aims at modelling neuronal activity as a statistical Poisson process, in which the instant firing probability nonlinearly depends on an external covariate, e.g. the force/velocity cue level (figure 4) or the EMG envelopes (figure 5). Instead, the KF assumes a linear relationship between the observed variables (the neural firing rates), and the (unknown) variable to be estimated, i.e. the covariate. Although the KF is widely used, there is no reason to believe that the neural system behaves linearly. In fact, more recent works have begun to use non-linear and non-stationary variants of the KF [63, 64], indicating that substantial decoding improvements are possible through the use of nonlinear and possibly time-varying system models [59].

In this work, we have shown the advantages of reconstructing EMG signals from neural data as an intermediate step to decode desired grasping force/velocity levels. When comparing the EMG signal reconstructed with the KF and PP approaches (figure 5), we noticed that there were rather small differences in the R² coefficients between PP-EMG or KF-EMG and the envEMG, during the isotonic tasks. On the contrary, in the isokinetic task, the proportion of the variance in envEMG that is predictable from KF-EMG is greater than that from PP-EMG. Even if this would suggest that the KF should be preferred to the PP approach, we have shown that the PP gives better classification results (see figure 6). This is probably due to the fact that neural firing rates, used by the KF to reconstruct the EMG envelopes, better represent the average values of the muscle activity that is captured by the surface EMG signals, while spike times are better in the reconstruction of transient values of EMG signals, important for early decoding movement intention and desired grasp force/velocity. However, further studies are needed to clarify this point.

Interestingly, the results shown in figure 5 also suggest that TIMEs actually provide a rather selective gate with the motor fibers, as the proportion of the predictable EMG variance from the recorded neural signals was rather high (60%, on average), that is, these signals were mainly generated by different efferent units. Only a few studies have addressed so far whether EMG data collected during limb movement can be faithfully reconstructed from activity in the nervous system, either in humans or in non-human primates [65–68]. The reason is that, especially in prosthetic research, EMGs are considered primarily as non-invasive input data for classifiers that aim to decode motor intention or limb kinematics. In these studies, the focus is rather on increasing the number of independent control signals, which has led to techniques such as TMR and RPNI. However, being able to successfully reconstruct the activity of multiple muscles from a few implanted peripheral neural interfaces, would allow one to leverage neural data to implement control strategies that have been shown to be successful with EMG data. This, for example, could be possible by: (a) recording—in parallel to the neural signals—high-density sEMGs or multiple acute intramuscular EMGs during the training phase of the prosthesis, and (b) training a classifier to map neural data into EMG signals, as it was done in this work.

A straightforward possibility might be to implement a direct (proportional control) of the neuroprosthesis, in which the reconstructed EMGs are directly converted into control signals. Nevertheless, Kuiken et al have recently demonstrated that myoelectric pattern recognition control outperforms direct control at least for transhumeral amputees with TMR [69]. Further studies are certainly needed to test both direct control and more complex pattern-recognition-based strategies with reconstructed EMG signals, and to investigate the overall efficacy and stability of this paradigm.

As is generally the case in the central nervous system [70], even when neural activity from peripheral nerves is recorded, an implanted electrode may record signals from multiple nearby neurons. Isolating the activity of a single neuron from this MUA (spike sorting) usually leads to better decoding results. There is a large body of literature on spike sorting for CNS, but many recent decoding studies in PNS do not use spike sorting at all. In this work, we leveraged the high selectivity of the TIME electrode to extract spikes from the signal and assign them to different classes of nerve fibers. In fact, based on the distance and relative orientation between the active site of the electrode and the axon, the APs of nearby neuronal axons may have slightly different features (shapes, duration, amplitudes and others) which can be used by spike sorting algorithms to classify them into separate classes [71].

Large datasets, such as those recorded with multichannel neural interfaces, may include signal dependencies as well as channels that generate erroneous predictions. Moreover, reducing the dimensionality of the original dataset minimizes the computational load and improves the processing speed, which are key for real-time applications. Principal component analysis (PCA) and independent component analysis (ICA) have both been used for this purpose [32, 72]. PCA produces a transformed data set whose members are uncorrelated, whereas ICA produce transformed signal sets whose members are independent of each other. In this work, we preferred tools for channel selection drawn from Information theory [46, 73]. Despite its complexity, this approach allowed us to choose units whose information about early intention of movement (and rest) was highest. Moreover, by using Information decomposition tools described in [49], it was also possible to discard redundant neurons, yet preserving those encoding information synergistically.

Interestingly, application of data reduction tools also highlighted a relationship between task type and location of implanted neural interface within the PNS. As showed in figure 2, neural units recorded from the median nerve were significantly more informative about the desired level of grasping force rather than hand closing velocity. On the other hand, neural units from the ulnar nerve acted in reverse. Moreover, electrodes implanted in proximal nerve segments appeared to convey more information than distal electrodes. On the contrary, in Cracchiolo et al [32], the recordings from all four TIMEs were all equally informative for decoding different grasp intentions. Therefore, although distal electrodes might appear to be useless in our setup, our results and those reported in [32] suggest the possibility of recording independent control signals (e.g. for decoding desired grasp types, grasping force and velocity independently) from different implanted electrodes. However, our experimental setup did not allow for disentangling motor commands, as the subject was requested to modulate either the force or the velocity levels while performing a given grasp type.

Additionally, one might speculate that median and ulnar nerves are selectively activated when imagining movements of the phantom hand that would elicit specific sensory feedback conveyed by that nerve, a principle that would effectively guide electrode implantation. Although this hypothesis is supported by seminal works on sensory control of dexterous manipulation in humans [74], a recent study [75] showed that, of 119 total sensory perceptions evoked by stimulation of single USEA electrodes implanted in the residual median or ulnar nerves of a transradial amputee, 72% were evoked from the median nerve. Therefore, further research is needed to shed more light on this point.

This work is the first proof of concept of the possibility of using neural signals from TIMEs for decoding force and velocity information in a transradial amputee subject and presents some limitations. First of all, we limited our study to a single movement, i.e. grasping, which consists in modulating the open/close movement of the hand around an object to be picked and lifted. However, this movement is of great importance for many activities of daily living, and is desired by all the users, independently either of the level of limb loss (e.g. transradial, transhumeral) or the type of prosthesis (e.g. body-powered, myoelectric) [62].

Furthermore, we applied these algorithms to data recorded in a single subject with upper limb amputation. We performed a short-term stability analysis of the proposed decoding strategy across days, collected in a single day. However, the subject was not always able to complete the whole protocol; in particular one force level in the isotonic task was missing in two recording days. Moreover, the dense experimental program scheduled for the amputee subject, not limited to the protocol here described, prevented us to record a greater number of traces to perform a more robust, in-depth, and long-term stability analysis. We are aware that peripheral intraneural probes currently run into issues of signal stability and degraded performance over time (see [76] for a review), however a new generation of long-term peripheral interfaces may be available soon [77].

In this work, we compared decoding approaches using SUA with an EMG-based reference classifier. As reported in the section 2, the input data for the reference classifier were envelopes of the activity recorded from residual superficial muscles involved during finger flexion and extension. Here, the number of EMG signal is significantly less than what can be achieved with high-density surface EMG, and their quality significantly worse than what can be achieved with intramuscular EMG electrodes. Nevertheless, we believe that our EMG setup is not very different to that embedded within current state-of-the-art sockets for myoelectric prostheses. It is certainly possible that a more sophisticated recording setup would have led to better decoding performance of force and velocity levels from EMG signals, and therefore of a reference EMG-based classifier. However, the same may hold at least for the most promising approach presented herein, i.e. reconstruction of EMG envelopes from SUA and use of them for classification. In fact, the PP or KFs would have been trained to reconstruct increasingly reliable EMG signals, and would thus have been able to fully exploit the independent sources of information contained in the neural recordings, leading to an increase in decoding performance.

Finally, we should confirm our results during real-time control of hand prostheses. In principle, all the algorithms presented herein could be implemented on a real-time system without significant performance degradation. In fact, most of the required computations are limited to the training phase of the HMI. Nevertheless, in the operational phase, a fast and accurate spike sorting algorithm is the key to achieve results similar to those obtained offline. For this reason, in our setup we did not employ a more performing but computationally heavy neural spike-sorting algorithm, such as the superparamagnetic clustering [42], but a template-matching-based spike sorting algorithm, which can be implemented in real-time with low sorting latency and high sorting throughput [78].

Future works will perform the same analysis to predict the force and velocity levels in a longer timeline and apply these algorithms in online experiments. Moreover, we aim to integrate in a single framework the algorithms here proposed with the techniques proposed by our group in [32], to decode not only the type of movements intended by the subject but also infer the different velocities and forces of the task. In fact, in [32] we showed that grasp type could be reliably decoded within the first 200 ms after the external cue signal, while in this work we demonstrated that different levels of grasping force and velocity can be subsequently extracted with an optimal performance (see figure 6). In [32], a different decoding approach was used, based on multi-unit rather than single-unit neural activity. Although that was a successful method for extracting information about grasp types, we did pursue a different approach herein, convinced that, to extract finer motor commands such as grasping force or movement velocity, a more selective gate with descending motor fibers had to be established, which required recording the activity of single-units.

Furthermore, an interesting analysis would be the comparison of results collected with TIMEs in subjects with upper limb amputation to data recorded in healthy subjects. In particular, Petrini et al showed the possibility to decode four levels of force and velocity from recordings in healthy subjects collected with microneurography from the median nerve [36]. However, the performance of their decoding method, even on data collected from healthy subjects, were significantly lower than those obtained with the approaches shown herein. We aim to repeat on their dataset the analysis herein proposed, to have more insights the mechanisms behind the activation of peripheral nerve triggering grasping movements. Since microneurography is a precise technique that allows to collect neural signals directly from the desired axons, this would shed light on the analogies between motor neural signal recordings collected with two different techniques, in healthy subjects and subjects with amputation. These analyses would help to improve the decoder strategy providing a more dexterous, natural and intuitive hand prosthesis.

This study is the proof of concept that it is possible to infer information about the modulation of the movement with intraneural interfaces, specifically TIMEs. The proposed strategy allowed us to decode, reliably in three different testing days, the velocity and force levels of movement intention with an accuracy comparable with the accuracy obtained with the EMG signals recorded from the residual muscles. Here, we added an additional important step heading in the direction of neuroprosthetics using TIME interfaces, since in our previous work [32] we have already showed the potential of this type of neural interface to decode up to 11 different movement grasp types in the same subject here reported. However, a longer timeline and at least an additional subject are needed to assess the long-term stability of these algorithms.