Development and evaluation of a non-invasive brain-spine interface using transcutaneous spinal cord stimulation

Abstract

Motor rehabilitation is a therapeutic process to facilitate functional recovery in people with spinal cord injury (SCI). However, its efficacy is limited to areas with remaining sensorimotor function. Spinal cord stimulation (SCS) creates a temporary prosthetic effect that may allow further rehabilitation-induced recovery in individuals without remaining sensorimotor function, thereby extending the therapeutic reach of motor rehabilitation to individuals with more severe injuries. In this work, we report our first steps in developing a non-invasive brain-spine interface (BSI) based on electroencephalography (EEG) and transcutaneous spinal cord stimulation (tSCS). The objective of this study was to identify EEG-based neural correlates of lower limb movement in the sensorimotor cortex of unimpaired individuals (N = 17) and to quantify the performance of a linear discriminant analysis (LDA) decoder in detecting movement onset from these neural correlates. Our results show that initiation of knee extension was associated with event-related desynchronization in the central-medial cortical regions at frequency bands between 4 and 44 Hz. Our neural decoder using µ (8-12 Hz), low β (16-20 Hz), and high β (24-28 Hz) frequency bands achieved an average area under the curve (AUC) of 0.83 ± 0.06 s.d. (n = 7) during a cued movement task offline. Generalization to imagery and uncued movement tasks served as positive controls to verify robustness against movement artifacts and cue-related confounds, respectively. With the addition of real-time decoder-modulated tSCS, the neural decoder performed with an average AUC of 0.81 ± 0.05 s.d. (n = 9) on cued movement and 0.68 ± 0.12 s.d. (n = 9) on uncued movement. Our results suggest that the decrease in decoder performance in uncued movement may be due to differences in underlying cortical strategies between conditions. Furthermore, we explore alternative applications of the BSI system by testing neural decoders trained on uncued movement and imagery tasks. By developing a non-invasive BSI, tSCS can be timed to be delivered only during voluntary effort, which may have implications for improving rehabilitation.

Fig. 1

Technological framework and experimental design. (a) Hardware setup. Technological framework for EEG and EMG signal acquisition and the delivery of tSCS using BCI2000 [17]. (b) Phase I experimental protocol. Phase I consisted of 3 conditions, 6 blocks per condition, and 10 trials (repetitions) per block. Cued movement blocks were used to train the decoder with 5-fold cross-validation, with the sixth block reserved for testing

Fig. 2

Analysis of predictive frequency bands informed selection of feature space used in the LDA decoder. (a) Power spectral data recorded via sensorimotor channels during right knee extension for a representative participant. Analysis of the power spectral data during movement of pilot participants informed the selection of frequency bands within the feature space. (b) R² scalp topographies for a representative participant. R² was computed between the true movement label, and the power spectrogram was computed for each channel. Pilot data revealed sensorimotor desynchronization in several frequency bands, including µ (8–12 Hz), low β (16–20 Hz), and high β (24–28 Hz). Non-neighboring frequency bands below 30 Hz were selected to prevent overlap in information fed into the decoder and to avoid the stimulation artifact at 30 Hz with the future addition of real-time tSCS. (c) EEG data processing pipeline. EEG data was bandpass 4–40 Hz filtered and common average referenced. Power was extracted in 4 Hz bins by band-passing, squaring, and low-pass filtering the common average referenced data. 480 features were extracted corresponding to 3 frequency bands (µ, low β, and high β), 5 lags, and 32 channels. Lags were incorporated so movement onset predictions can take data from the past 0.5 s into account. (d) Five-fold cross-validation decoder training. The decoder was trained with a 5-fold cross-validation strategy in which four blocks were used as training blocks. Once the hyperparameters were optimized to minimize validation error, the model was retrained on all five blocks and tested on the sixth, unseen block

Fig. 3

Offline decoder performs above chance for cued, imagery, and uncued movement. (a) EEG spectrograms aligned with EMG and movement kinematics during cued movement. A single participant’s EEG power features during offline testing on cued movements decomposed by frequency band aligned in time with EMG signals from the vastus lateralis (VL), rectus femoris (RF), tibialis anterior (TA), and medial gastrocnemius (MG), knee angle, and probability calculated offline. Desynchronization was observed in the EEG spectrograms across frequencies at the onset of movement. (b) EEG spectrograms aligned with EMG and movement kinematics during imagery. (c) EEG spectrograms aligned with EMG and movement kinematics during uncued movement. (d) Probability for single trials (thin lines) and averaged across trials (thick line) during focus and movement periods. Movement probability increased after movement onset. (e) Illustration of true positive and false positive rate calculation from probability. True positive and false positive rates are calculated at a sweep of probability thresholds between 0 and 1 and used to construct a receiver operating characteristic (ROC) curve. (f) ROC curve for a single participant. The ROC curve was calculated for a single participant by varying the threshold on the probability and comparing to true movement. The probability threshold used in e is shown on the ROC curve. The area under the ROC curve (AUC) was used to quantify all decoder performances. (g) ROC curve for each participant (thin lines) and averaged across participants (thick line) when testing on cued movement and average area under ROC curve (AUC) compared to chance. (h) Same as g but for a decoder tested on imagery. Paired comparison between cued and imagery and comparison of each to chance. (i) Same as g but for a decoder tested on uncued movements. Paired comparison between cued and uncued movement and comparison of each to chance. (j) Confusion matrix at a fixed probability threshold when tested on cued movement. k. Same as j but for a decoder tested on imagery. l. Same as j but for a decoder tested on uncued movement. Bars in g-i represent mean ± s.d., with each circle representing the testing AUC for each participant. The asterisks on the right of each bar represent the results of the one-sample Wilcoxon signed rank test for each decoder’s AUC against chance; the asterisks between bars represent the paired samples Wilcoxon signed test between two decoder’s average AUC. ∗ p < 0.05; ∗∗ p < 0.01; ∗∗∗ p < 0.001. Abbreviations: rectus femoris (RF), vastus lateralis (VL), tibialis anterior (TA), medial gastrocnemius (MG)

Fig. 4

Analysis of R² suggests spectral and spatial differences in EEG activity between conditions. (a) R² scalp topography plots during movement for two participants showing good (left, participant S005) and poor (right, participant S001) generalization across conditions. Consistencies in spatial and spectral activity across conditions resulted in better decoder generalization evidenced by higher AUCs. (b) Difference in R² during movement between conditions for one participant (left) and averaged across participants (right). Group analysis of R² differences between conditions revealed focused differences between cued movement and imagery, and widespread differences between cued and uncued movement. Significant channels were not corrected for multiple comparisons. (c) PCA projections of R² scalp topographies for all conditions and frequency bands. R² data across electrodes was projected onto the first three PCs. Each point represents one participant, and data is color-coded by condition. PCs across conditions (left and right columns) are the same. Conditions shown separately for visualization purposes. (d) Average Euclidean distance of imagery and uncued movement to cued movement in PC space. Distance between uncued and cued movement is slightly larger than the distance between imagery and cued movement, but this effect is not significant

Fig. 5

Real-time, closed-loop control of tSCS using predictions of movement intention in a non-invasive BSI. (a) Technological framework for non-invasive BSI. Desynchronization in the sensorimotor cortex was identified using a 32-channel EEG system in real-time as participants extended their right knee. The predicted movement intention was used to trigger the delivery of tSCS at a higher amplitude. (b) Illustration of a cued movement block used in the training set. Stimulation was ramped up to a baseline of 10 mA and was increased to 15 mA during the movement phases using the task’s movement cue. (c) EEG spectrograms of selected sensorimotor channels, kinematics, movement probability, and real-time closed-loop stimulation for a single block for a representative participant (thin lines) and averaged across trials (thick line). Note that there is event-related desynchronization before and after movement onset. Stimulation onset was generally timed with movement onset as shown by the stimulation averaged across trials. (d) ROC curve for individual participants (thin lines) and averaged across participants (thick line) when testing on cued movement with stimulation. (e) Same as (d) but for a decoder tested on uncued movement with stimulation. (f) Paired comparison of average AUC across participants for cued and uncued movement and comparison to chance. (g) Confusion matrix calculated according to the stimulation administered during the trial and averaged across participants. (h) Same as g but for a decoder tested on uncued movement. Bars in f represent mean ± s.d., with each circle representing the average AUC for one participant. The asterisks on the right of each bar represent the results of the paired samples Wilcoxon signed test for each decoder’s AUC against chance; the asterisks between bars represent one-sample Wilcoxon signed test between two decoders’ average AUCs. ∗ p < 0.05; ∗∗ p < 0.01; ∗∗∗ p < 0.001. Abbreviations: rectus femoris (RF), vastus lateralis (VL), tibialis anterior (TA), medial gastrocnemius (MG). EEG map and brain region division modified from [29]

Fig. 6

Cued and uncued movement-trained decoders’ performance on uncued task. (a) ROC curve for all participants (thin lines) and averaged across participants (thick line) when the decoder was trained and tested on uncued movement. (b) Average AUCs for a decoder tested on uncued movement trained on cued and uncued movement. The performance of the decoder trained on uncued movement (n = 4) was generally higher than the performance of the decoder trained on cued (n = 9), but not significantly different when compared with a bootstrapping analysis. (c) Confusion matrix calculated from the real-time stimulation and averaged across participants. The mean of TPR and TNR was 53%. Bars in b represent mean ± s.d., with each circle representing the average AUC for one participant. The asterisks on the right of each bar represent the results of a one-sample Wilcoxon signed rank test for each decoder’s AUC against chance. ∗ p < 0.05; ∗∗ p < 0.01; ∗∗∗ p < 0.001

Fig. 7

Imagery-trained decoder can generalize to cued movements. (a) ROC curve for all participants (thin line) and averaged across participants (thick line) when the decoder was trained on imagery and tested on cued movement. (b) Paired comparison of average AUC across participants for a decoder tested on cued movement, trained on cued movement and imagery. AUC between decoders trained on cued and imagery were not significantly different and were both significantly higher than chance. (c) Confusion matrix calculated from simulated stimulation according to the implemented stimulation paradigm and averaged across participants for a decoder trained on imagery and tested on cued movement. The mean of TPR and TNR was 59%. Bars in (b) represent mean ± s.d., with each circle representing the average AUC for one participant. The asterisks on the right of each bar represent the results of the one-sample Wilcoxon signed rank test for each decoder’s AUC against chance; ∗ p < 0.05; ∗∗ p < 0.01; ∗∗∗ p < 0.001

Fig. 8

Decoder accuracy changes as a function of tolerance for discrepancies in time. (a) Illustration of onset detection accuracy calculation. Onset detection accuracy was calculated as a function of tolerance time window length around true onset. Predicted onsets were defined as a positive crossing of the probability above threshold. (b) Illustration of varying the tolerance window around true onset to calculate the onset detection accuracy. Tolerance window was varied from 0 to 3 s. True positives were considered any predicted onsets within the window, and true negatives were the absence of an onset prediction within rest. (c) Onset detection accuracy plotted as a function of tolerance time window around true onset and averaged across participants for all conditions. The point of comparison across conditions (0.8 s) is denoted by a vertical dotted line