Identification of Spared and Proportionally Controllable Hand Motor Dimensions in Motor Complete Spinal Cord Injuries Using Latent Manifold Analysis

2.1 Datasets

Two datasets were used to investigate the potential use of the proposed AI system to automatically detect the spared motor dimensions that are still present in the high-density EMG of participants with SCI.

The first dataset served as the control group and consists of 13 young, uninjured adult participants (age 25.9 ± 2.8 years) (data also used in [13]– 15]).

The test group consists of 8 individuals with motorcomplete SCI (the same dataset was also used in [5]). Table 1 shows the individual patient characteristics; however, in summary, all injuries occurred in the cervical area (C5-C6), 4 subjects have an AIS [2] grade of B (S1-3, 8), 3 of A (S4-6), and 1 of C (S7).

All subjects were instructed to perform a series of grasps and individual digit movements, starting and ending each cycle in a specified position, or, in the case of individuals with motor-complete SCI, to attempt these movements. The starting position is defined as the natural resting position of the hand, and the end position is defined as the fully flexed position. The series consisted of individual digits flexion and extension, grasping, and twoand three-finger pinches (Fig. 2A). These movements were displayed in front of the subjects in the form of a virtual hand performing prerecorded movements at a frequency of 0.5 Hz (Fig. 2A & B). The virtual hand was a digital rendering of one representative hand from an uninjured participant. During movement, monopolar EMG signals were recorded from the forearm of the subjects at 2048 Hz. For the recordings we used three 8 by 8 electrode grids with 10 mm interelectrode distance (IED) (GR10MM0808, OT Bioelettronica, Turin, Italy) and two 13 by 5 grids with 8 mm IED (GR08MM1305, OT Bioelettronica, Turin, Italy), for a total of 320 electrodes. The 8 by 8 grids were applied around the upper forearm area below the elbow, while the 13 by 5 grids were placed above the wrist around the arm covering the dorsal and the ventral part of the forearm (Fig. 2B).

• Download figure
• Open in new tab

Fig. 2:

A. The various hand tasks performed in this study by the control (thirteen healthy uninjured participants) and eight participants with motor-complete SCI. The movements were executed at 0.5 Hz. B. Experimental setup is shown on both an uninjured and an injured participant. We used five electrode grids. Three 8 by 8 grids were placed on the proximal part of the forearm and two 13 by 5s above the wrist. C. Schematic overview of the deep learning model. Convolutional layers scan EMG signals in spatial and temporal dimension, detecting motor unit and muscle fiber action potential ensembles indicative of specific movement phases. Subsequent layers compress and integrate these overlaps, considering both intragrid and inter-grid information. The output of the convolutional layers is a high-dimensional latent space that separates and preserves the information of each movement. For each 100 ms EMG window, the convolutional network produces an n-dimensional vector that contains compressed information about the movement. The last part of the model is a multilayer perceptron that maps the high-dimensional latent space into the desired 3D coordinates. A detailed explanation can be found in Sîmpetru et al. [11]. D. Using UMAP, a non-linear dimensionality reduction technique that preserves the local and global information, we project the latent space vectors into human understandable 2D. The manifolds are color-coded according to the individual movements. For the projection metric, we computed the cosine distance to cluster similar vectors irrespective of their amplitude. Interestingly, the movement clusters are positioned in the same way as the fingers of the human hand. In addition, the grasping action is enclosed by the flexion and extension of all individual fingers. E. For each manifold, we fit an inner and outer ellipse using the algorithm described by Halir et al. [17]. F. The points outside either the inner or outer ellipse are considered outliers. The percentage of outliers is displayed on top of the ellipse.

Each movement was recorded for 40 s, resulting in a total data collection time of 5 min 20 s per subject. The recorded EMG data were saved as 2D matrices, the two dimensions representing the 320 electrodes and the recorded time. Successively, we split the data into multiple 100 ms (192 samples) windows with a shift of 31.25 ms (64 samples). The EMG window shift of 64 samples was set because this was the fastest real-time output setting of the EMG recording device (EMG Quattrocento, OT Bioelettronica, Turin, Italy) used in our experiments. The size of the EMG window was chosen as three times the shift. The middle 30% (12 s) of the data for each movement is designated as the test set, while the remaining 70% (28 s) is allocated for the training set. Within the test set for each movement, the first 10% (1.2 s) is further reserved for the validation set. This results in 3 min 44 s training, 1 min 26 s testing, and 9.6 s validation data per participant.

To increase the training data size, we employed three augmentation techniques described in Tsinganos et al. [16]:

• Gaussian noise: Gaussian noise was added to the original EMG data with a signal-to-noise ratio of 5 dB.

• EMG magnitude warping: To mimic the effect of movement drifts of the electrode grids, we applied a randomly generated cubic spline to the EMG signals. This introduces non-linear warps to the signal.

• Wavelet decomposition and reconstruction: To further increase the amount of data, we decomposed the EMG signal using wavelets. Before reconstructing the signal, we multiplied the decomposed factors by a constant, ensuring that the reconstructed signal is a novel yet representative (correlated) to the original EMG signal.

The augmentations offered a three-fold increase in training data, resulting in a total per subject of 18 min and 56 s.

2.2 Ethical approval

All participants gave their written informed consent to participate in the study. The study was carried out in accordance with the Declaration of Helsinki, except for registration in a database. All procedures and experiments were approved by the ethics committee of the FriedrichAlexander-Universität (applications 22-138-Bm and 21-150B).

2.3 Model

Using the EMG data collected, which were synced with the digitally displayed hand kinematics, we trained our previously published model [11], which is capable of accurately and continuously predicting hand kinematics in real-time for each individual participant. The model was implemented in Python using PyTorch (version 2.1.0+cu121) and PyTorch-Lightning (version 2.1.0).

The deep learning model first divided each 100ms EMG window into five individual grids, resulting in five matrices instead of one large one. The model was made out of two parts: a convolutional network and a multi-layer perceptron. We have previously shown [18] that the convolutional network extracts the motor information and projects it into a high-dimensional latent space. The multi-layered perceptron mapped the high-dimensional space into the desired 3D kinematics output.

2.4 Latent Space Projection

The latent space of the models (Fig. 2C) are each 2256 dimensional vectors, which, for human comprehension, were projected into only 2 (Fig. 2D). To achieve this, we used the dimensionality reduction technique known as uniform manifold approximation and projection (UMAP) [12]. UMAP was chosen because it is capable of preserving both local and, to some extent, global non-linear structures in the data. This is accomplished by constructing a graph based on similarities among data points, where the similarity is computed using a user-selected metric, such as cosine similarity in our case. The graph is then optimized to reduce its dimensionality while retaining as many local and global relationships as possible. However, the optimization process is inherently stochastic, which can result in different spatial rotations of the data for each subject. To ensure consistent interpretability across subjects, we aligned the 2D embeddings during the optimization process to minimize orientation variations.

The parameters (Fig. 2E & F, Fig. 6B-I, and Fig. 7B) extracted from the manifolds (Fig. 2D) were computed as follows:

• For the temporal behavior analysis (Eq. 1, Fig. 6B & C), the angle (α) between an embedding point (p) and the centroid (c) of the task manifold was computed using the formula:

Here, the subscripts x and y refer to the respective components along the xand y-axes of the points (p) and the centroid (c). The angles between the centroid and each point on the task manifold are then phase unwrapped. For perfect circularity, this results in a straight line (Fig. 6B), whereas imperfect circularity leads to a discontinuous line. The measure of fitness (R²) was computed between the phase unwrapped angles and a linear fit for each task and for each subject (Fig. 6C).

• Each manifold (Fig. 2D) was fitted with an ellipse using the algorithm described in Halir et al. [17]. To enclose the manifold within an ellipse, the fitted axes (semi-minor and semi-major) were multiplied by a factor of 2 (Fig. 2E, Fig. 6D). To further delineate the anticipated shape of the actual manifold, we used a secondary ellipse measuring 0.5 times the dimensions of the fitted ellipse. This smaller ellipse, in conjunction with the larger one, formed a 2D band-an ellipsoidal annulus that represented our hypothesized manifold shape (Fig. 2E). Using this band structure created by the larger and smaller ellipses, we computed the number of outliers (Fig. 2F). The number of outliers was divided into two groups depending on whether they were encased by the band or outside of it (Fig. 5). To assess how well the manifold adheres to the 2D ellipsoidal band, we computed the closest distance from each point on the manifold to the band. The sum of these squared distances is small when the band is closely followed and increases as the spatial arrangement of the latent space becomes more chaotic (Fig. 7B).

• Sparsity (s, Eq. 2, Fig. 6F & G) was defined as

1 minus the division between the sum of all overlap areas (O_A) and the sum of the total areas (A). This measure yields a value of 1 in cases where there are no overlaps, while a value of 0 indicates that all areas overlap and that tasks cannot be distinguished from each other.

• The number of motor dimensions (Fig. 6H & I) was set to 8 at the start, which corresponds to the total number of tasks. For each overlap found between two tasks, the amount of motor dimensions was reduced by 1, because the merged tasks were still separable from the rest of the tasks.

2.5 Conformal Prediction

Conformal prediction [19] is a framework for machine learning that provides calibrated prediction sets that, in turn, can be used as a proxy for uncertainty quantification. We used Regularized Adaptive Prediction Sets (RAPS) [20], as it showed visually reliable prediction regions that adapt well to the underlying data.

The fundamental concept behind RAPS is to construct prediction sets around individual predictions made by a model. These sets are designed to capture the true target value with a certain probability. Unlike point predictions, these prediction sets provide a range of possible values and quantify the uncertainty associated with the model’s predictions. In this study, we used prediction sets with a 95% probability of containing the true class within the measured data. This approach allowed us to quantify whether a prediction is considered certain — indicated by a single prediction within the set — or uncertain, denoting multiple predictions within the set (Fig. 4C & 5B).

2.6 Comparison Models

We compared our model with two standard machine learning approaches, a linear discriminant analysis model to separate movement tasks and a ridge regressor to predict proportional movement. Both models were trained with the following features computed for each 100ms EMG window:

• mean absolute value

• root mean square value per electrode

• slope sign change

• waveform length

• zero crossings

This resulted in 324 features for each EMG window. The features were standardized using the statistics obtained from the training set. During testing, the same scaling parameters derived from the training set were applied.

2.7 Statistical Measures

The statistical analysis (Fig. 3C, 4D, 6C, 6G, and 6I) was performed using Welch’s t-test because the normal distribution assumption of the data could not be guaranteed. A p-value of below 0.05 was considered significant.

• Download figure
• Open in new tab

Fig. 3:

A. The predictions derived from the test dataset are compared against the camera recorded kinematics of a repersentative uninjured participant by utilizing the Euclidean distance of the predicted and measured individual fingertip trajectories. B. Comparison of predicted and expected index fingertip positions (y-axis) for an uninjured (n. 2) and an SCI afflicted participant (n. 8) in the test set. The colored and grayed-out indications at the top denote whether the index finger is active or not during a task. C. Mean Euclidean distance for all subjects in mm. Control group participants show an error of 13.6 ± 15.4 mm, while the participants with SCI demonstrate an error of 23.9 ± 23.2 mm. Statistical significance has been confirmed (p < 0.001) using a Welch’s t-test with a total of 120,980 data points with 5,762 per participant (74,884 from control group and 46,096 from SCI group).

• Download figure
• Open in new tab

Fig. 4:

A. The latent spaces of the trained networks are projected into 2D using UMAP [12], a non-linear dimensionality reduction technique that preserves local and global information. A support vector classifier is used to separate the clusters in the 2D projection. We observed that for some participants with SCI the tasks could not be entirely separated and appear to be intertwined. To account for this uncertainty in classification, we applied conformal prediction to get prediction sets instead of only one prediction per data point. The prediction sets are 95% probable to contain the correct class inside. If a prediction set is uncertain (i.e., it contains more than one class), we can solve the uncertainty by looking at previous predictions that have been certain and thus making an informed estimate of the class. B. Example 2D projections showing the 2nd uninjured (control) and 8th spinal cord injured participant. Datapoints are colored according to each movement. C. Visualization of prediction uncertainty for the example latent space projections from B. The colorbar indicates the number of classes within the prediction set. A set is considered certain if it contains only one class. A prediction set with more than one class represents uncertainty in the model. D. Accuracy over all participants. Control participants achieved 99.8% while those with a SCI 98.3% accuracy. The points from the SCI group are denoted by an “S” followed by the participant id. No statistical significance (p = 0.06) has been found using Welch’s t-test.

3.1 Hand Kinematics Prediction

We evaluated the prediction of the kinematics using the Euclidean distance between the expected and estimated hand kinematics (Fig. 3A). Fig. 3B shows the prediction for two representative participants (the 2nd control participant and the 8th participant with SCI). The colored hands, shown at the top, illustrate the active participation of the index finger in the executed movement. For instance, tasks such as grasping and pinching activate the index finger, while tasks such as flexion of the thumb do not. For an assessment of relative predictive capacity, we show the index fingertip prediction versus the expected values throughout the test set, encompassing all executed movements (Fig. 3B). The aggregated Euclidean distance across all subjects is depicted in Fig. 3C. The mean distance for the control participants was 13.6 ± 15.4 mm, while for the SCI participants it was 23.9 ± 23.2 mm. Each subject contributed 5,762 datapoints, resulting in a total of 120,980 for the entire dataset (74,884 from the control group and 46,096 from the SCI group). The means were found to be significantly different (p < 0.001, Welch’s t-test).

In the representative example depicted in Fig. 3B, it is evident that the spinal cord injured participant demonstrated proportional control over both pinching and, to some extent, the action of closing and opening the fist. This indicates their ability to generate voluntary and task-modulated signals despite having a motor-complete SCI. Despite some noise, the AI successfully interpreted this activity and generated the intended hand kinematics, underscoring the viability of proportional control in individuals with SCI.

3.2 Latent Space Analysis

After examining the prediction of the AI model (that is, the correlation between estimated and expected kinematics), we examined the internal representation learned by the AI from the EMG signals, which corresponds to the latent space of the model (Fig. 4, 5, 6, and 7). Because the latent space of the model contains a human-incomprehensible amount of dimensions, we projected the latent space using UMAP into 2 dimensions. Briefly, this method enables the visualization of high-dimensional nonlinear data into a discernible lowdimensional output (see Methods for details). The UMAP output was analyzed for both the control and SCI groups. We compared the spread of the clusters for each digit that was learned by the AI, their overlaps, and the shape of the continuous prediction relative to time and space. We observe a lack of overlap between the clusters of the control group compared to the SCI group. Participants with SCI showed overlaps between related tasks (e.g., ring and pinky fingers are anatomically close together) (Fig. 4B), which may indicate that at the neural level there is a low number of spared motor units encoding that specific task. In addition, the clusters of the control group are arranged in circular shapes, whereas in the SCI group they generally appear more chaotic, which indicates that sinusoidal movements (e.g., continuous change in index flexion and extension position) are not discernible as well as for uninjured participants.

• Download figure
• Open in new tab

Fig. 5:

A. The high-dimensional latent space was projected into 2D using UMAP [12]. Each task manifold was fitted with a 2D band representing our hypothesized optimal shape (see Methods). The number of outliers outside the band and those enclosed by it was computed and averaged for each subject. B. Using the outlier data obtained as described in A, a support vector classifier was trained to separate the participants based on the two features. Uncertainty was gauged using conformal prediction (see Methods). The numbers denote the participant ID. C. For interpretability, one certain and one uncertain example is shown for the control and SCI group. The participant IDs are indicated in the upper right corners.

• Download figure
• Open in new tab

Fig. 6:

A. The 2D projection of the latent space (example shown for the 8th SCI subject) was analyzed for both temporal (B & C) and spatial (D-H) behaviour. B. The angle of each point of the manifold of the task relative to the centroid of the manifold was computed using Eq. 1. The resulting angles have then been phase-unwrapped, which, in case of being the result of a circular path, would result in a straight line. This assumption was used to calculate the R²-Score of the angles in a linear regression fit. C. Analysis described in B was performed on all participants. The differences between uninjured participants and those affected by an SCI were found to be significant (p < 0.001) using Welch’s t-test. D. For the spatial behaviour analysis, we fitted ellipses around each task manifold and computed the intersections between them. E. Confusion matrices for the 2nd control and the 8th SCI participants. Movement tasks are denoted by a simplified notation (T - Thumb, I - Index, M - Middle, R - Ring, P - Pinky, G - Grasp, 2F - 2 Finger Pinch, and 3F - 3 Finger Pinch). F. Using Eq. 2 we calculated the amount of overlaps divided by the total area of the manifold (see Methods for more details). G. Sparsity analysis described in F was carried out on all participants. The differences between uninjured and participants with SCI were not found to be significant (p = 0.24) using Welch’s t-test. H. The amount of separable motor dimensions was considered the maximum of tasks executed (8). If there is overlap between two tasks, the amount is reduced by 1. I. Motor dimension analysis described in H was carried out for all subjects. Differences between control and SCI participants were not found to be significant (p = 0.17) using Welch’s t-test.

• Download figure
• Open in new tab

Fig. 7:

A. Correlation between the sum of squared residuals (distance from each point to 2D ellipsoidal band) and the R²-Scores computed for predicted versus expected fingertip positions across the test set. The y-axis represents the average prediction R²-Score, while the x-axis shows the sum of squared residuals. B. Correlation scatter plot with a superimposed line of best fit, showing the relationship between the two measures. Additionally, the density of the points is shown by a 2D kernel density estimation (KDE). Marginal KDEs are also displayed for each axis. The correlation coefficient was -0.73, indicating a significant negative relationship between the amount of deviation of the latent space from the 2D band (sum of squared residuals) and the prediction of kinematics (p-value < 0.001). C. Comparison between our model and a ridge regressor model, where each point on the graph represents the R²-Score of a task’s proportional prediction. The connecting line links the corresponding task scores between the ridge regressor and our model. Statistical analysis revealed significant differences between models for both SCI and control participants, with p-values less than 0.001, as determined by Welch’s t-test. D. Comparison between our model and a linear discriminant analysis model (LDA), where each point on the graph represents the accuracy of the movement task separation of a participant. The connecting line links the corresponding scores between the LDA and our model. Statistical analysis revealed no significant differences between the models for both SCI and control participants, with p-values greater than 0.05, determined by Welch’s t-test.

3.3 Correlation Between Prediction and Temporal Behaviour

Based on the prevalence of a predominant circular arrangement of the movement task manifolds in the control group compared to the SCI group (Fig. 6C), we evaluated the correlation between the amount of deviation (sum of squared residuals) of the latent space from a 2D fitted ellipsoidal band and the prediction accuracy of the hand kinematics (Fig. 7A & B). We computed the average R²-Score between the predicted and expected fingertip positions for all movements in the test set and correlated it against the sum of squared residuals (Fig. 7A). In Fig. 7B, we display the correlation scatter plot with the best fit line superimposed on it. For a better understanding of the density of the points, we also computed the kernel density estimation (KDE) and displayed both the 2D KDE and the marginal KDE for each respective axis. The correlation achieved was -0.73 (p < 0.001).

We note that movements that do not follow a circular path (the sum of squared residuals is high) exhibit a negative correlation with their proportional controllability. This correlation suggests the potential utility of employing such analysis for the detection of proportionally recoverable movements which can be used for targeted neurorehabilitation.

3.4 Comparison With Standard Models

The last test performed was to compare our model with two standard models. A ridge regressor to proportionally predict all hand joints (Fig. 5C) and a linear discriminant analysis model to classify individual hand movements from each other (Fig. 5D). For regression, we present the mean R²-Score for each movement task, which yields a total of 96 data points for the control participants and 64 data points for participants with SCI. In the case of classification, we report accuracy for each participant, resulting in 13 and 8 data points for control and spinal cord injured participants respectively.

Our model demonstrates significantly better regression performance (p < 0.001, Welch’s t-test), while showing no significant difference from linear discriminant analysis in classification (p = 0.80 or 0.22, Welch’s t-test).

4.1 Circular Latent Space Projections Correlate with Proportional Prediction

Given that all participants in the SCI group have been unable to use their hands for at least several years, we anticipated that the error (mean Euclidean distance) between the estimated and expected hand kinematics would be statistically higher for them compared to the control group (Fig. 3). The network’s ability to proportionally decode some movements in the SCI group (Fig. 3) is consistent with our previous findings, where we demonstrated similar capabilities in a smaller subset of movements [4], [5].

Although not all movements attempted by participants with SCI could be proportionally decoded, we found that the majority of movement tasks were distinguishable from each other using the latent space of the neural network (Fig. 4 and Fig. 6G & I). The separation accuracy remained consistent between the uninjured individuals and those with SCI, with no significant statistical differences observed between the two groups (Fig. 6G & I). Our findings suggest that participants with motor-complete chronic SCI were able to produce signals varied enough to differentiate between multiple movements. Our previous study [5], in which we performed a real-time decomposition of sEMG signals during movements analogous to those examined in this work, found that most tasks shared a limited number of MUs. The sparsity findings (Fig. 6G & I) in this study also support this notion. It is easier to separate signals with less overlap (fewer MUs) than those with a higher level of overlap, as expected from uninjured participants.

Despite not training the hand kinematics decoder models to distinguish between individuals with and without SCI, we demonstrate that the latent space encodes this information with 100% accuracy (Fig. 5). Not only do we find the latent space to be structured, but we could also extract different physiological parameters from it. We showed that individual movement manifolds are temporally arranged in a circular path, reflecting the cyclical nature of attempted movements (Fig. 6B & C). Furthermore, we observed that the separation of tasks is sparse, with no significant differences between the two groups (Fig. 6D-I).

A circle, however, is characterized not only by having its constituting points follow a temporally circular path, but also by the points being spatially equidistant from the center. We demonstrate that the deviation from an elliptical band fitted to an individual movement task manifold is negatively correlated (-0.73) with the accuracy of movement proportional decoding. This observation implies that as the chaotic nature of a manifold increases, there is a corresponding rise in the error rate of its proportional prediction (Fig. 7A & B).

4.2 Necessity of Deep Learning

We recorded the forearm muscles of participants with and without SCI using 320 electrodes. Our previous research has shown that decoding hand movements is enhanced with a larger number of electrodes [21], as it eliminates the need to precisely define the recording area, which is particularly advantageous when dealing with the territories of unknown spared MUs in participants with SCI. Having this large amount of data warrants the question of the necessity of a deep learning model, as a machine learning variant might be just as effective and easier to implement in practice.

Indeed, we observed that our deep learning network performs similarly to its machine learning counterpart in movement classification (Fig. 7D). However, the machine learning model struggles to proportionally predict movements in both the injured and uninjured groups, a task at which our deep learning network excels (Fig. 7C).

The success of the algorithms utilized for classifying various hand movements (Fig. 7D) implies that they autonomously acquired the ability to filter out non-movement relevant MU signals while retaining those encoding motor intent. This notion is further substantiated by the sparsity measure (Fig.6G) and the quantification of separable motor dimensions (Fig.6I).

4.3 Potential Clinical Application

Due to a negative correlation of -0.73 (p < 0.001, Fig. 7B) observed between the deviation from circularity of the latent space embeddings and the accuracy of the proportional prediction, we propose that our system could serve as a test to identify potential movements that can be fully restored through an assistive system and direct monitoring of the spared neural pathways during rehabilitation.

The proposed test would involve patients trying to replicate the movements detailed in this study, while a model predicts their kinematics based on recorded sEMG signals. Subsequently, the UMAP projection of the latent space could help identify circular manifolds. This identification process could be performed by medical professionals or by using algorithms, as described in Section 2.4. Implementing such a test could help medical professionals in planning and monitoring rehabilitation plans that specifically target movements with the potential for proportional control by the patient.

Another significant clinical application is that even when a movement is not proportionally controllable, our system demonstrates reliable classification accuracy (> 93%, Fig. 4). Consequently, it can be at best used to proportionally decode a movement and at worst classify it as being different from the other movements learned. This can be used as the control input to an assistive device that has substantial potential to enable patients to be independent in daily living tasks and significantly reduce dependence on caregivers.

However, most importantly, our results emphasize the need to assess whether a patient who cannot move without assistance can produce volitional motor signals that are detectable by EMG. This aspect is currently not part of the standard SCI assessment procedure [2], despite the potential presence of proportionally controllable EMG signals postinjury. Such signals can inform assistive devices, such as an exoskeleton or orthoses, on which movement to generate. This distinction would not only significantly impact the psychological well-being of patients by informing them that they can regain movement with an assistive device, but it would also benefit physicians. Doctors can design and plan more effective rehabilitation treatments that target specific movements that the patient can control. Furthermore, this distinction would benefit researchers by eliminating the time-consuming examination of each potential study candidate that would be necessary without such an assessment. Therefore, we suggest that the motor-complete diagnosis of SCI also assesses the ability to generate EMG signals. If a patient can produce intentional EMG signals, they are considered EMG incomplete; otherwise, they are classified as EMG complete. The need for a more refined classification of patients with motor-complete SCI has also been proposed by Sangari et al. [22], who found a correlation between the presence of spasticity and subsequent improvement in the assigned grade of injury for the patient. Spasticity, indicative of spared MUs post-injury, can be evaluated through noninvasive methods such as EMG.