Transcutaneous Electrical Spinal Cord Stimulation Promotes Focal Sensorimotor Activation that Accelerates Brain-Computer Interface Skill Learning

Necessity of longitudinal SMR-based BCI training

The control of SMR-based BCIs typically requires feedback training over multiple sessions [10, 26–32]. This is evident in our study with the significant main effect of session in the analysis of variance for the linear mixed effects (LME-ANOVA) models that assess the results of both groups for all three BCI performance metrics (Fig. 2), namely Classification Accuracy (F₅,₄₄₄ = 4.85, p < 0.001); Bar Dynamics (F₅,₄₄₄ = 3.60, p < 0.01); and BCI Hits (F₅,₄₄₄ = 3.78, p < 0.01). Meanwhile, the percentage of timeouts across sessions for both groups showed no significant main effects (Session - 1 to 6: F₅,₄₄₄ = 1.48, p = 0.19, Group - TESS or Rest: F₅,₄₄₄ = 1.48, p = 0.19), nor significant session × group interaction effect (F₅,₄₄₄ = 0.68, p = 0.64). Specifically for the Rest group, a significant group-level increase in BCI performance —relative to the baseline in session 1— emerged in session 5 on the last day of training (Classification Accuracy: t₇₈ = 3.90, p_corrected < 0.01; Bar Dynamics: t₇₈ = 2.89, p_corrected < 0.05; and BCI Hits: t₇₇ = 3.90, p_corrected < 0.05). Nevertheless, the Rest group still showed a significant increasing trend in performance across subjects for classification accuracy and BCI hits over the 6 training sessions (Fig 2b,c: r > 0.81, p < 0.05, n = 6) and a positive trend in bar dynamics (Fig 2a: r = 0.66, p = 0.096, n = 6). This is consistent with the expected learning process in longitudinal BCI training.

• Download figure
• Open in new tab

Figure 2: BCI performance results from the two main study groups.

a: Bar dynamics, b: BCI hit accuracy, and c Sample-wise classification accuracy, for the TESS (in orange) and Rest (in gray) groups over the six BCI sessions of the protocol. Violin plots show values of the (subject,run) pairs. Boxes within violins correspond to the inter-quartile range, white circles correspond to medians, and horizontal lines correspond to means. Session 1 performance serves as a baseline for within-group comparisons. LME-ANOVA reveal a significant main effect of session forall three metrics (F₅,₄₄₄ ≥ 3.6, p < 0.01) and a significant session × group interaction (F₅,₄₄₄ = 2.06, p < 0.05) for BCI hits. Post-hoc paired tests showed significant improvement in all BCI performance metrics for the TESS group starting on session 3, after two stimulation sessions, that persisted till the end of the training while significance emerged only after session 5 for the Rest group. Between-subject post-hoc tests show that TESS significantly outperforms Rest in BCI control (bar dynamics) and command delivery (BCI hits) starting on session 3. The solid regression lines show the trends of the medians over sessions. The dashed lines correspond to piece-wise spline regression of medians. Both groups show significantly increasing trends over sessions. Significant differences across groups are indicated by underlined asterisk, and within-group significant difference relative to baseline is indicated above the violins with orange asterisks for TESS and gray asterisks for Rest. Statistical testing is performed using linear mixed effect models. ∗ : p < 0.05, ∗∗ : p < 0.01, ∗∗∗ : p < 0.001, absence of asterisks indicates non-significant difference. All p-values are Bonferroni-Holm corrected for multiple comparisons.

TESS accelerates BCI training and MI skill learning

While both TESS and Rest groups showed a learning effect throughout training, the TESS group exhibited a significant enhancement of BCI performance relative to baseline starting as early as session 3 —compared to session 5 for the Rest group— (Classification Accuracy: t₇₈ = 3.30, p_corrected < 0.05; Bar Dynamics: Z = 3.97, p_corrected < 0.001; and BCI Hits: t₇₅ = 4.38, p_corrected < 0.01). Performance was further enhanced at the end of training on session 6 (Classification Accuracy: Δ_CA = 9.31%, t₇₈ = 4.98, p_corrected < 0.001; Bar Dynamics: Δ_barDyn = 8.64%, t₇₃ = 4.77, p_corrected < 0.001; and BCI Hits: Δ_Hits = 15.71%, t₇₀ = 4.58, p_corrected < 0.001). For both sessions 3 and 6, the TESS group showed a large effect size for the BCI hits improvement (Cohen’s d_s3 = 0.81, d_s6 = 0.94). The TESS group also maintained significant improvement in performance on session 5, which occurred at the start of a new training day and was not preceded by stimulation. TESS significantly outperformed Rest on session 3 for all performance metrics (t₇₈ > 3.96, p_corrected < 0.01) and maintained significantly higher bar dynamics and BCI hits on the last session (t₆₆ > 2.88, p_corrected < 0.01). Over the course of the sessions, the TESS group showed an increasing trend in BCI bar dynamics control (Fig 2a, slope = 3.06, r = 0.88, p < 0.05, n = 6). Since BCI hits saturated close to 100% following session 4 for TESS, we identified the presence of an initial steep learning slope reflecting rapid gains in the early learning phase followed by a gradual asymptotic plateau, consistent with the power law of practice [55]. The trend of the BCI hits accuracy in the TESS group from Session 1 to Session 3, the inflection point, was significantly increasing (slope = 9.29, r = 0.90, p < 0.05, n = 3).

Strong and focal µ and β ERDs emerge with TESS

Having demonstrated how the proposed intervention facilitates participants’ acquisition of BCI control, we examine the neurophysiological underpinnings of such BCI skill learning process. To this end, we analyzed SMR modulations in EEG during training, in particular the µ and β band ERDs in the contralateral sensorimotor areas. For analysis purposes, right MI data were mirrored before being pooled with left MI data so that ERDs are localized on the right hemisphere. For each subject, we first identified the single EEG channel in the contralateral sensorimotor cortex that exhibited the largest ERD strength at the end of training () and backtracked its ERD over sessions. Fig. 3a,b show the µ ERD for averaged across subjects in the TESS and Rest groups, respectively. In Fig. 3a, TESS exhibits significantly larger µ ERD in session 6 relative to session 1 (Cluster permutation testing: interval [3.53, 4.89] s, p_cluster = 0.041 and interval [5.27, 8.5] s, p_cluster = 0.017), in contrast to Rest. The inset topoplot for session 6 in TESS shows the emergence of a more focal ERD source at the relevant C3 location compared to that of Rest. Fig. 3c also illustrates that the contra-lateral ERD difference between session 6 and session 1 —averaged over the task execution period— is stronger for the TESS group relative to the Rest group. In the β band, both groups show significantly larger ERD in session 6 relative to session 1, but the difference is more prominent for TESS (TESS: [0.40, 5.11] s, p = 0.008; Rest: [1.63, 4.78] s, p = 0.012). When averaged over the task execution period, the ERDs in Fig. 3f show stronger change in contraletral ERD for TESS relative to Rest —consistent with the observation for the µ ERDs in Fig. 3c. To further analyze the focality of the emerging ERD patterns, we normalized the strength of the ERD from to the ERDs within the motor area using equation (6). We refer to the latter as ERD^∗. Fig. 3g,h show the evolution of ERD^∗ for the µ and β bands, respectively. LME-ANOVA revealed significant effect of session on ERD^∗ in both µ (F₅,₄₄₄ = 3.19, p = 0.008) and β (F₅,₄₄₄ = 6.50, p < 0.001). Only the TESS group exhibited a significantly increasing focality of µ ERD —i.e., a decreasing trend of ERD^∗ over the 6 sessions (slope = −0.08, r = −0.23, p < 0.001, n = 6) in Fig. 3g. On session 5, which was not preceded by stimulation, the µ ERD^∗ for the TESS group became significantly larger than baseline (t₇₇ = −3.05, p_corrected = 0.015) and also significantly larger than that of the Rest group (t₇₇ = −3.43, p_corrected = 0.008). At the end of the training on session 6, both groups had significantly enhanced ERD^∗ (Rest: t₇₁ = −3.01, p_corrected = 0.015; TESS: t₇₃ = −5.18, p_corrected < 0.001), but TESS remained significantly larger than Rest (U = 2.59, p = 0.031, Mann-Whitney U test) showing a large effect size (Cohen’s d = 0.94) compared to a medium effect for Rest (Cohen’s d = 0.58). For the beta ERD^∗ in Fig. 3h, both groups showed significantly decreasing trends (TESS: slope = −0.09, r = −0.28, p < 0.001, n = 6; Rest: slope = −0.08, r = −0.21, p = 0.002, n = 6). Participants in the TESS group, however, achieved significantly larger —relative to their baseline— on session 4 (t₇₈ = −3.76, p_corrected = 0.002) while the Rest group reached significance only on session 6 (Rest: t₇₁ = −4.78, p_corrected < 0.001). Overall, ERD results support the role of TESS in promoting faster and more effective modulation of SMRs for healthy subjects.

• Download figure
• Open in new tab

Figure 3: Neurophysiological evidence of enahnced SMR modulations.

a, b: Grand average ERD patterns in the µ band for the channel that exhibited the highest ERD strength at the end of training in the TESS and Rest groups, respectively. Black and orange lines correspond to the first and last online sessions, respectively. Inset topoplots show the ERD patterns over selected EEG channels in session 6. The confidence levels for the grand average patterns represent standard deviations across subjects, and the shaded time periods represent statistically significant difference between session 1 and 6 with cluster-based permutation testing. c: t-statistic of the differences in ERD strength between the last and the first session for the TESS group on top and the Rest group at the bottom. White asterisks indicate statistically significant difference using cluster permutation testing. d, e, and f : ERD patterns and topography for the β band similar to µ band above. g, h: Focality of ERD patterns in the µ and β bands, respectively, for the TESS (orange) and Rest (gray) groups. The plots backtrack the ERD strength of the channel that exhibited the greatest desynchronization at the end of training (ERD^∗) after normalizing it to the ERD values within the motor area (i.e., the FC, C, and CP electrodes). Trend lines are regression fits of session means. Significant differences across groups are indicated by underlined asterisk, and within-group significant difference relative to baseline is indicated above the bars with orange asterisks for TESS and gray asterisks for Rest. Statistical testing is performed using linear mixed effect models. ∗ : p < 0.05, ∗∗ : p < 0.01, ∗∗∗ : p < 0.001, absence of asterisks indicates non-significant difference. All p-values are Bonferroni-Holm corrected for multiple comparisons. Error bars represent standard deviations.

TESS helps slow BCI learners gain control

We further challenged the ability of TESS to help slow BCI learners or non-responders to typical training paradigms. For that, we identified the individuals in the Rest group who could not learn to control the BCI after the six training sessions. A total of four participants had no control over the bar dynamics at the end of training (Supplementary Fig. 1). These subjects were re-enrolled for the TESS intervention after a 6 months washout period during which they were not involved in any BCI study. Additionally, they also completed a follow-up BCI session —without TESS— one week after they finished the TESS protocol to test for long-lasting effects. After training with TESS, all subjects had improved control. LME-ANOVA shows a significant session × modality interaction effect for bar dynamics (F₅,₁₆₈ = 2.48, p = 0.034) and BCI hits (F₅,₁₆₈ = 3.96, p < 0.01). Only when trained with TESS, these subjects could achieve significantly higher bar dynamics accuracy at the end of training on both session 5 (t₃₀ = 4.44, p_corrected < 0.01), which was no longer preceded with stimulation, and session 6 (t₂₇ = 4.55, p_corrected < 0.01) (Fig. 4a). This is also the case for the BCI hits metric (Fig. 4b; session 5: t₃₀ = 4.95, p_corrected < 0.001; session 6: t₂₇ = 5.19, p_corrected < 0.001). The subjects could actually achieve significantly better BCI hits accuracy with TESS starting on session 3 (t₃₀ > 4.16, p_corrected < 0.01) with a very large effect size (Cohen’s d > 1.34) and they maintained it till the end of training. Overall, the subjects showed a significantly increasing trend in both metrics over the 6 sessions (Bar Dynamics: r = 0.94, p < 0.01, n = 6; BCI Hits: r = 0.96, p < 0.01, n = 6), and they could maintain significantly higher performance relative to their baseline on the one-week follow up session after TESS (Bar Dynamics: t₂₈ = 3.77, p_corrected = 0.015; BCI Hits: t₂₈ = 3.89, p_corrected < 0.01) with a large effect size (Cohen’s d > 1.042). Fig. 4c backtracks over the six training sessions of the protocol. Consistent with the BCI performance results, the subjects showed significantly increasing focality in only when trained with TESS (r = −0.44, p < 0.001). The latter increase in strength became significant at the end of training (t₂₈ = −5.88, p_corrected < 0.001) with a very large effect size (Cohen’s d = 2.12), and it was already tending towards significance on session 5 (t₂₈ = −2.71, p_corrected = 0.058). Fig. 4d,e show the grand average µ ERD patterns for . Only TESS shows a period of significantly larger µ ERD on session 6 compared to session 1 (interval: [1.07, 5.68] s, p_cluster = 0.028). Supplementary Fig, 2 illustrates the results for β band ERD. Both modalities exhibited an increasing β ERD focality over sessions (Supplementary Fig. 2c); TESS: r = −0.26, p < 0.05; Rest: r = −0.30, p < 0.01), consistent with the original group. The increase reached significant difference relative to baseline for the Rest modality on session 6 (t₂₇ = −3.42, p_corrected = 0.032); however, none of the modalities showed a significant change in the grand average ERD patterns (Supplementary Fig. 2a,b).

• Download figure
• Open in new tab

Figure 4: Cross-Over group of slow learners.

a and b: Bar dynamics and BCI hits performance, respectively, similar to Fig. 2a,b for the set of slow BCI learners when they did the intervention without stimulation (gray) and when they repeated it six months later with TESS (orange). c: Focality of µ ERD patterns (ERD^∗) for the TESS (orange) and Rest (gray) groups similar to Fig. 3g. d, e: Grand average µ ERD patterns for the channel that showed the highest ERD strength at the end of training similar to Fig. 3a, b for the TESS and Rest modalities, respectively. The shaded region in the TESS case highlights the period of significant difference between sessions 1 and 6. f : Source level analysis of µ ERDs. The top part shows session 6 patterns for the Rest (left) and TESS (right) modalities while the bottom part focuses on the contra-lateral hemisphere to show —from left to right—the difference between session 6 and session 1 for Rest and TESS as well as the difference at session 6 between the two modalities. ∗ : p < 0.05, ∗∗ : p < 0.01, ∗∗∗ : p < 0.001, absence of asterisks indicates non-significant difference. All p-values are Bonferroni-Holm corrected for multiple comparisons. Error bars represent standard deviations.

To further analyze the physiological differences that emerge when training with TESS, we used subject-specific MRIs to reconstruct the source activity at the cortical level. We focused on the µ band where a significant change in modulation was observed. Fig. 4f shows the average source patterns after warping the individual patterns of each subject on a common MNI template brain. The top part of Fig. 4f shows the µ ERD patterns in session 6, which were very weak on the contra-lateral right hemisphere when subjects trained without stimulation (left) compared to when they re-trained with TESS (right). The bottom part shows the µ ERD differences on the contra-lateral right hemisphere between session 6 and session 1 for each modality as well as the difference between modalities on the last session. Clearly, there is stronger desynchronization with TESS —consistent with the channel-level results in the left plot of Fig. 4d. The pattern differences on the last session also show higher µ ERD for TESS compared to Rest, and this appears close to where the hand representation is expected to be in the contra-lateral hemisphere. Per-session source results are provided in Supplementary Fig. 3, and they show evolving ERD patterns on the contra-lateral motor area consistent with the subjects’ incremental skill learning over sessions in Fig. 4a,b.

Cortex-exciting spinal stimulation impedes BCI learning

Since our hypothesis advocates a role of cortex-inhibiting TESS in promoting BCI skill learning, we also investigated a potential opposite effect for cortex-exciting spinal stimulation. The latter can be achieved through a modified stimulation sequence that does not incorporate the 5 kHz carrier frequency as detailed and validated in [45] (refer to Fig. 6). Five additional subjects were recruited to complete the experimental protocol of Fig. 1b with excitatory spinal stimulation instead of TESS (“No-Carrier” group). Fig. 5a,b show the BCI performance of these subjects compared to the original TESS group. LME-ANOVA shows a significant session × group interaction effect (Bar Dynamics: F₅,₃₂₈ = 4.87, p < 0.001, BCI Hits: F₅,₃₂₈ = 2.73, p = 0.019) with the No-Carrier group showing no significant change over the six training sessions. In fact, statistical testing with the Rest group also shows a significant group × session interaction effect for both BCI performance metrics (F₅,₃₂₄ > 3.13, p < 0.01) suggesting that the subjects in the No-Carrier group were impeded from acquiring better BCI control. Fig. 5a even shows a significant decrease in bar dynamics accuracy for the No-Carrier group on the last day of training (session 5: t₃₈ = −2.62, p_corrected = 0.023). This is also consistent with the fact that the ERD patterns in Fig. 5c,d for the µ and β bands remained almost unchanged throughout the training. The inset topoplots show strong ERDs on session 6; however, they are less localized relative to those of the TESS group (Fig. 3a,d). Nevertheless, when normalized within the motor area, in Fig. 5c shows a marginally significant trend over sessions (r = −0.21, p = 0.052), but it was not as prominent as that of the TESS group (r = −0.28, p < 0.001). Meanwhile, β ERDs in Fig. 5f show no significant change. Results from this group provide further evidence that the cortex-inhibiting facet of TESS plays the major role in promoting SMR modulation and BCI skill learning.

• Download figure
• Open in new tab

Figure 5: No-Carrier group.

a, b: Bar dynamics and the BCI hits, respectively, for the No-Carrier group (gray) in comparison to the original TESS group (orange) similar to Fig. 2a,b. c, d: Grand average µ and β ERD patterns, respectively, for the channel that showed the highest ERD strength at the end of training similar to Fig. 3a, d. e, f : Focality of µ and β ERD patterns (ERD^∗), respectively, for the TESS (orange) and Rest (gray) groups similar to Fig. 3g,h. ∗ : p < 0.05, ∗∗ : p < 0.01, ∗∗∗ : p < 0.001, absence of asterisks indicates non-significant difference. All p-values are Bonferroni-Holm corrected for multiple comparisons. Error bars represent standard deviations.

• Download figure
• Open in new tab

Figure 6: BCI performance results for the SCI patients.

a, b, and c: Bar dynamics, BCI hit accuracy, and sample-wise classification accuracy, respectively, for Patient-I who received supra-lesional stimulation. d, e, and f : Bar dynamics, BCI hit accuracy, and sample-wise classification accuracy, respectively, for Patient-II who received infra-lesional stimulation. Violin plots show values of the (subject,run) pairs. Boxes within violins correspond to the inter-quartile range, white circles correspond to medians, and horizontal lines correspond to means. Session 1 performance serves as a baseline for within-group comparisons. The solid regression lines show the trends of the medians over sessions.

Supra-lesional TESS enhances BCI learning post-SCI

To illustrate the potential of the proposed intervention in motor rehabilitation and the control of assistive devices, we enrolled two subjects with tetraplegia after SCI who completed the experimental protocols for TESS and Rest with one week of break. Patient-I had a C5 fracture and received supra-lesional stimulation while Patient-II had a C4 fracture and received infra-lesional stimulation. Results in Fig. 6a,c for Patient-I show a significant main effect for TESS (Classification Accuracy: F₁,₃₁ = 14.44, p < 0.001; Bar Dynamics: F₁,₃₁ = 14.87, p < 0.001) with a significant increase in classification accuracy (r = 0.83, p = 0.021) and an increasing trend in bar dynamics (r = 0.68, p = 0.086). This increase during the TESS intervention sessions was followed by a plateau, or a less steep slope, one week later when the patient resumed the study without TESS. Although the BCI hit accuracy in the upper part of Fig. 6b shows no significant change over the sessions throughout the entire experiment, the overall number of delivered commands increased during TESS for Patient-I because, as shown in the lower part of Fig. 6b, there was a significantly decreasing number of timeout trials (r = −0.86, p = 0.013). The latter exhibited a significant decrease on the last day of training (session 5: t₆ = −6.82, p_corrected < 0.01; session 6: t₅ = −10.97, p_corrected < 0.01). The results also show a significant session × intervention interaction effect for both classification accuracy (F₅,₃₁ = 4.34, p < 0.01) and timeouts ratio (F₅,₃₁ = 13.19, p < 0.001) supporting a period of skill learning and consolidation during TESS followed by a period of maintained performance one week later without TESS. When analyzing ERD patterns for the Patient-I, results reveal a significant session × group interaction effect for (F₅,₃₁ = 7.78, p < 0.001) and (F₅,₃₁ = 9.45, p < 0.001). The patient showed a prominent change in ERD patterns on session 3 (day 2), which coincides with the re-calibration of the decoderto incorporate emerging patterns following the first post-stimulation online session (session 2). Starting on session 3 onwards, both and show significantly strong trends (µ: slope = −0.77, r = −0.78, p < 0.001; β: slope = −0.4, r = −0.78, p < 0.01) followed by a period of non-significant change one week later when the patient repeated the protocol without stimulation (Supplementary Fig. 4a,b).

Regarding Patient-II, as expected because of the infra-lesional stimulation, both Rest and TESS modalities did not result in a considerable change in performance (Fig. 6d,f). However, the upper part of Fig. 2e shows an increase in BCI hits during the TESS sessions (session 6: t₄ = 5.39, p_corrected = 0.057) with a significant session × intervention interaction effect (F₅,₁₈ = 4.03, p = 0.013). Nevertheless, this came at the expense of an increasing number of timeouts.

Physiological correlates of enhanced BCI control following TESS

Finally, we performed additional analyses to better understand the neurophysiological basis of how TESS promotes BCI skill acquisition. Spearman’s correlation analysis revealed that the increase in µ ERD focality moderately explained better BCI command delivery performance in all the healthy users (Fig. 7a, Spearman correlation: r = −0.4, p < 0.001, n = 165) —i.e., including all 29 subjects in the TESS, Rest, Cross-Over, and No-Carrier groups. The correlation was also significant for the TESS subjects (Spearman correlation: r = −0.48, p < 0.001, n = 77). For Patient-I, who received supra-lesional stimulation and showed evidence of BCI learning, the correlation was particularly larger (Spearman correlation: r = −0.6, p = 0.028). We further checked the specific role of stimulation-induced inhibition in promoting ERD focality and BCI performance. For that we used resting-state EEG recordings —collected for the subjects in the Cross-Over and No-Carrier groups as well as Patient-I— to study the effect of stimulation on the α band power, which is commonly associated with inhibitory processes [56–58]. Fig. 7b shows that TESS results in a significant increase in resting-state α power within the motor cortex (ΔERS = 26.75% ± 17.42%, t₄ = 3.43, p = 0.026) while the no-carrier stimulation has no significant effect (ΔERS = 5.67% ± 14.74%, t₄ = 0.86, p = 0.439). The strength of α power post-stimulation explains the increased focality of µ ERD (Supplementary Fig. 5: Spearman correlation, r = −0.58, p < 0.001, n = 40), and Fig. 7c shows that the percentage change in α power post-stimulation correlates with the percentage change in classification accuracy (Spearman correlation: r = 0.41, p < 0.01, n = 40).

• Download figure
• Open in new tab

Figure 7: Physiological correlates of enhanced BCI skills.

a: Correlation between the focality of µ ERDs and BCI hits accuracy. The black line corresponds to the line of best fit for all healthy subjects including 14 who had TESS, 10 who had Rest, and 5 who had the no-carrier stimulation. The orange line corresponds to the TESS subjects only. b: Percentage change of α power post-stimulation for the subset of subjects who had pre- and post-stimulation resting-state EEG recordings. ∗ : p < 0.05, absence of asterisks indicates non-significant difference, and error bars represent standard deviations. c: Correlation between the percentage change in α power and the percentage change in decoding performance post-stimulation for all subjects in black and the TESS group in orange. d: Source-level µ ERD patterns across the six training sessions for the TESS group in the top row, the Rest group in the middle one, and the No-Carrier group in the bottom row. Contra-lateral ERDs shall appear on the right hemisphere around the hand motor area. e: Normalized µ ERDs in session 6 within the motor cortex for the TESS, Rest, and No-Carrier groups from top to bottom. f : Differences in µ ERD patterns at the end of training between the TESS and Rest groups on top and the No-Carrier and Rest groups in the bottom. g, h, and i: Group-level trends in the TESS, Rest, and No-Carrier groups, respectively, for the size of the most significant contra-lateral ERD cluster and its normalized ERD strength relative to the surrounding ERS (Surround Inhibition) in the sensorimotor area. Dotted lines are lines of best fit.

To localize the effects of stimulation in the brain, we reconstructed the source-level ERD patterns for the three groups TESS, Rest, and No-Carrier (Fig. 7d). The TESS group shows an emerging focal ERD over the hand sensorimotor area on the contra-lateral hemisphere. This is particularly evident in Fig. 7e, which shows the normalized ERD values within the entire sensorimotor cortex on the last training session. Fig. 7f shows the difference in raw ERD patterns between the two stimulation groups and the Rest group. When compared to Rest, TESS shows more localized ERDs surrounded by stronger ERS whereas the No-Carrier group exhibits more spread- out ERDs with ipsilateral activation. To quantify the difference in focality, Fig. 7g-i tracks the size of the most significant ERD cluster across sessions and backtracks the normalized strength of ERD (surround inhibition) for the significant cluster on session 6. Only TESS in Fig. 7g shows very strong and significant group-level trends over sessions (Cluster Size: r = −0.9, p < 0.012; Surround Inhibition: r = 0.86, p < 0.034). The latter were also strongly and significantly correlated with enhanced classification accuracy (Cluster Size: r = −0.88, p = 0.022; Surround Inhibition: r = 0.90, p = 0.015).

Subjects

A total of 25 able-bodied healthy subjects with no history of neurological disorders participated in the study (age 18 to 30 years-old, 48% female). In addition, three SCI patients were recruited (2 males, 1 female). The female SCI participant had to be excluded due to substantial lower back pain arising from her recent injury, which prevented her from completing all sessions. The first male SCI participant (referred to as Patient-I) had his injury at the age of 13, nine years prior to the study. He is clinically diagnosed as tetraplegic with a C5 fracture and a complete lesion, and he has been undergoing consistent motor-rehabilitation therapy. The second SCI participant (referred to as Patient-II) incurred his injury within three months of the study at the age of 26, and he is diagnosed as tetraplegic with a complete spinal cord injury involving a C4 fracture and cervical stenosis. Notably, the latter participant assumed a nearly supine position during the experiment to alleviate back discomfort, and he was receiving stimulation at C5/C6 below the lesion site. All participants, both healthy and with SCI, provided written informed consent to the procedures of the study as approved by the University of Texas at Austin Institutional Review Board (IRB protocol number: 2020-03-0073), and they were compensated for their participation. The study protocol is published on ClinicalTrials.gov (NCT05183152) and the CONSORT enrollment flow-diagram is provided in Supplementary Fig. 9. Study Design. Main Group. In the main part of the study, a cohort of 20 individuals without prior experience in utilizing MI-based BCIs completed the longitudinal BCI training protocol depicted in Fig. 1b. These participants were randomly allocated to either the experimental group, which received pre-training TESS stimulation lasting 20 minutes (referred to as TESS group), or the control group (referred to as Rest group) that engaged in a 20-minute resting period before training. The initial day of the protocol started with an offline BCI session wherein participants perform cued MI for movements of their left or right hand. MI instructions highlighted the importance of performing kinesthetic —not visual— MI by emphasizing the mental rehearsal of the sensation and muscular tension association with the execution of the movement. A Riemannian-geometry-based BCI decoder was trained on the collected offline data and was used to provide real-time feedback during subsequent online sessions. On the second day, participants started with an online BCI session, which serves as their performance baseline. This was followed by 20 minutes of rest or stimulation, depending on their group assignment, and then a post-conditioning online BCI session. After this second day, the BCI decoder was re-calibrated to incorporate data from the post-conditioning online session, and then it was fixed for the remainder of the study. Days 3 and 4 included 20 minutes of conditioning followed by an online BCI session. Finally on Day 5, participants started with an online BCI session, which tests the one-day-after effect without further conditioning, then they completed one last 20 minutes of conditioning followed by the final online BCI session. In addition to healthy participants, patients with SCI completed the protocol of Fig. 1b once with the TESS condition and once with the Rest condition, with a one week separation period. Following random order selection, Patient-I started with TESS then Rest while Patient-II followed the opposite order. The patients also completed 2 minutes of resting state EEG recordings before an after each stimulation session.

Cross-Over Group. This subgroup includes the participants from the Rest condition who were unable to attain BCI control by the end of the study. A total of four subjects exhibited chance-level performance consistently across all of the training sessions. These participants were subsequently re-enlisted to undergo the TESS condition, adhering to the same protocol of Fig. 1b with the addition of 2-minutes of resting state EEG before and after each stimulation session. The washout period for this group was six months to ensure no carry-over effect from the initial training. In order to comprehensively evaluate the impact of TESS on enhancing their learning capabilities, we performed an assessment of their capacity to sustain the attained level of performance during a session occurring one week later. Importantly, this evaluation was conducted without any supplementary TESS stimulation, thus further illuminating the enduring effects of TESS-induced skill enhancement.

No-Carrier Group. This subgroup includes an additional five participants who underwent the same training pro-tocol with a distinct stimulation condition, closely resembling TESS but lacking the 5kHz carrier frequency (Fig. 1a). Specifically, it consists of a single 200us biphasic pulse within one period. This specific stimulation sequence is believed to elicit excitatory effects at the cortical level, contrary to the inhibitory effects brought about by the conventional TESS sequence [45]. This group also completed 2 minutes of pre- and poststimulation resting state EEG recordings. The inclusion of this group serves a twofold purpose: firstly, it effectively controls for potential placebo effects, ensuring the discernment of genuine outcomes; and secondly, it substantiates the premise that TESS inhibitory facet is the pivotal driving force behind the observed enhancements in SMR modulation.

BCI Training

Bar Task. The MI-based BCI training employs a conventional bar feedback task [5]. In offline sessions, the bar moves consistently to either the right or left direction over the duration of the task execution depending on the class of the trial (i.e., right hand or left hand MI). In online sessions, the bar provides real-time feedback on MI performance. It moves according to an accumulated evidence of posterior probabilities generated by the binary-class BCI decoder, which classifies EEG signals into one of the two classes. The BCI decoder operates on samples of a 1-second sliding window with a step size of 62.5ms, and it accumulates evidence for MI classes according to the following equation: where e_i is the accumulated evidence for a given class and p_i is the posterior probability of that class for the i^th sample. e_i is reset to 0.5 at the beginning of the task execution period of every trial, and a trial ends when e_i crosses a predefined command delivery threshold for either class. The outcome scenarios of a trial are: “trial hit,” where a correct command is delivered upon e_i crossing the threshold on the correct side aligned with the trial’s class; “trial miss,” where an incorrect command is delivered if e_i crosses the threshold on the wrong side; and “timeout,” when neither threshold is crossed within a 7-second period. For each run, threshold values are adjusted to provide a suitable challenge level, maintaining participant engagement while minimizing the likelihood of erroneous commands and optimizing the probability of accurate commands.

Session Structure. A session consists of 3-4 runs, each with 20 trials that are randomly split into right or left hand MI (see Fig. 1d). A trial starts with a rest period of 2 s with a random jitter of 1 s, followed by a fixation period for 2 s, a left or right arrow cue indicating the MI task for 1.5 s, a task execution period of 5 s in offline runs and a timeout duration of 7 s in online runs, and a trial result feedback of 2 s.

EEG/EOG

EEG was recorded at 512 Hz using an ANTNeuro EEGO amplifier with a waveguard EEG cap (ANT Neuro, Netherlands). The cap uses Ag/AgCl coated soft polymer electrodes positioned over 32 standard scalp locations according to the international 10-10 system. For patient recordings and for the two additional subgroups (Cross-Over and No-Carrier), a higher density 64-channel waveguard EEG cap was used. In both cases, the ground and reference electrodes were placed at AFz and CPz, and the same subset of channels was used for online classification. Electrode impedances were kept below 20 kΩ. Synchronized electrooculography (EOG) signals were collected via a bipolar box connected to the amplifier, from three channels placed on the eye canthi and the forehead with the ground and reference on the mastoids. EOG was used to reject samples with eye blinks or left/right eye movements. All processing and analysis of EEG and EOG signals were conducted using MATLAB and the Fieldtrip toolbox [71].

BCI Decoder

Feature Extraction. The real-time EEG classification (right hand versus left hand MI) is performed according to the Riemannian geometry framework, which has been attracting growing interest in non-invasive BCI research [72, 73]. Riemannian classifiers use covariance matrices of mulit-channel EEG as features to characterize spatio-spectral neural patterns. In our processing pipeline, trace-normalized covariance matrices in equation (2) are estimated for a subset of 22 EEG channels (i.e., F7, F3, Fz, F4, F8, FC5, FC1, FC2, FC6, C3, Cz, C4, CP5, CP1, CP2, CP6, P7, P3, Pz, P4, P8 and POz) using the shrinkage-based covariance estimator method of Ledoit and wolfe [74] after bandpass filtering EEG to [8, 30] Hz using a 2nd order causal Butterworth filter. where C is the sample-wise normalized covariance, X ∈ ℝ^Nc×Nt denotes a 1-second EEG sample, N_c is the number of channels, and N_t is the number of time points. Since covariance matrices are symmetric positive definite (SPD), they lie on a Riemannian manifold, whose geometrical properties allow defining the geodesic distance δ_r(C₁, C₂) between two SPD matrices C₁ and C₂ according to the following equation [72]: MDM Classifier. For classification purposes, class prototypes (for left hand versus right hand MI) are constructed as the centers of mass C̄ (or karcher means) for their respective features C_i on the Riemannian manifold by solving the optimization problem in equation (4) over the space of N_c-dimensional SPD matrices P_n using the method outlined in [72]. Using the estimated class prototypes and the geodesic distance metric of equation (3), a simple Minimum Distance to Mean (MDM) classifier can be used to classify an EEG sample based on the shortest geodesic distance between its covariance matrix and the class prototypes. In an online setting, the classifiers output is blocked whenever an EOG artifact is detected.

TESS

We used the same stimulation sequences and settings as in [45] to ensure we elicit similar modulatory effects. TESS was delivered for 20 minutes over the neck at the cervical level between the C5-C6 spinal segment and the iliac crests of the hip bone as depicted in Fig. 1a. A 3.2-cm-diameter hydrogel adhesive electrode was used on the neck and two 7.5x12 cm rectangular electrodes (ValuTrode Cloth, Axelgaard Manufacturing) were located bilaterally over the iliac crests. We used a Hasomed Rehamove3 stimulator to generate the stimulation sequence of 200 us biphasic pulses at a frequency of 30 Hz with a carrier frequency of 5 kHz as depicted in Fig. 1. The peak-to-peak amplitude of the stimulation current was set to the maximum tolerable value that did not result in any muscular contraction (14 healthy: 41.0 ± 9.5 mA, Patient-I: 44 mA, Patient-II: 40 mA). For the No-Carrier group, TESS without 5 kHz carrier frequency was delivered for 20 minutes with single biphasic pulses (200 us duration) at a 30 Hz frequency. This paradigm was tolerated by all subjects, allowing us to match similar stimulation intensity as used in the original TESS protocol (5 healthy: 55.8 ± 12.8 mA). Since changes in pulse duration can affect the relative recruitment of sensory axons, we maintained the waveform pulse duration consistent across paradigms similar to [45].

Statistical Analysis

We used linear mixed effect models (LMEs) to perform statistical testing and regression analysis. Specifically, LME-ANOVA was used to test for the main and interactions effects with the between-subjects fixed effect as the group condition (TESS, Rest, TESS w/o 5 kHz carrier), the within-subject fixed effect as sessions (sessions 1 to 6 plus follow-up sessions), and the random effect of subjects (to account for paired comparisons). Appropriate LME models were used for post-hoc paired and independent tests, and p-values were corrected with the Holm-Bonferroni method for multiple comparisons. We verified the normality of the residuals with the Lilliefors test, and whenever it is violated, we used Mann-Whitney U test for paired samples (U statistic reported), Wilcoxon signed-rank test for independent samples (Z statistic reported), and Spearman’s rank correlation coefficient for correlation analysis. For studying the learning curve from the BCI performance over sessions, we identified the inflection point as the session that minimizes the residual error of horizontal lines of best fit from the pre and post inflection segments, and we parameterized a piece-wise linear spline function around the inflection to characterize the period of significant steep ascent and the one with plateau. Regression lines for performance metrics over sessions were fitted using the medians across subjects in order to track progressive group-level improvements even when the values saturated at 100% accuracy. We used cluster permutation testing —implemented in FieldTrip— to correct for multiple comparisons in spatio-temporal (time-frequency plots) or spatio-spectral (topoplot or source space) analysis. Neighboring electrodes, voxels, frequency bins, or time instances that show sample-level significant effects with a t-test are grouped into clusters. The sum of the t-statistic is computed within each cluster to form a cluster-level statistic. The test statistic is defined as the maximum of the cluster-level statistic values over all detected clusters. The reference distribution for the test statistic is computed through the Monte Carlo method by randomly permuting the observations across subjects (between-subjects design) or across the spatio-temporal or spatio-spectral dimensions. The significance probability (p-value) is then computed for the observed test statistic based on the approximated reference distribution.

Outcome Measures

BCI Performance. The BCI performance is assessed by the following metrics:

• Sample-wise Classification Accuracy: It represents the online decoder accuracy for the 1s overlapping BCI samples, and it is computed as follows: where TP, TN, FP and FN correspond to true positives, true negatives, false positives and false negatives, respectively.

• Bar Dynamics: It is a measure of BCI control, and it represents the percentage of time during the MI task execution period when the user could keep the feedback bar on the correct side of the interface. Bar dynamics track the accumulated evidence of the BCI decoder.

• Command Delivery/BCI Hits: It measures the accuracy with which the feedback bar reaches the command delivery thresholds. It is the ratio of correctly delivered commands to the total commands delivered. This metric shall be analyzed along with the timeout ratio, which accounts for the percentage of trials in which no command was delivered.

Electrode Discriminancy Score (EDS). We use the recently proposed EDS for the Riemanian geometry classification framework [32] to quantify the contribution of each EEG channel to the BCI decoding accuracy. The metric assesses channel contributions through a backward elimination approach, which quantifies the loss in posterior prediction accuracy following the removal of a given channel from the model.

Physiological. We track the evolution of ERDs over the course of training to test how different conditions affect the evolution of ERD patterns. Particularly, we assess the strength and focality of ERDs by identifying the EEG channel that exhibited the strongest ERD (ERD^∗) on the last training session, and we backtrack its values across sessions. To test for the focality of the ERD at that channel, we normalize its value to the ERDs within the motor area ERD_motor as in equation (6) or just within the contra-lateral motor area.

ERD/ERS

Trial periods with EOG artifacts were rejected from analysis. To compute ERD/ERS patterns, EEG signals were bandpass filtered to the frequency range of interest (µ, ([8, 13] Hz), or β, ([18, 30] Hz)) using a non-causal 4^th order Butterworth filter, spatially filtered using a Laplacian plus montage, squared to compute the instantaneous power, and then averaged over the analyzed period. The average power in the MI execution period was normalized to a pre-task [-0.8, 0] s interval —avoiding any residual visually evoked potentials from cue presentation.

Resting State EEG

We computed the power of the α band oscillations ([8, 13] Hz) as a correlate of cortical inhibition, and we wanted to quantify how TESS might modulate resting state α oscillations. We used irregular resampling auto-spectral analysis (IRASA) [75], a method used to distinguish between rhythmic activity and concurrent power-spectral 1/f modulations. This technique involves modifying the time-domain data through a series of non-integer re-sampling factors before conducting Fourier-based spectral decomposition. As a result, the powerspectrum’s rhythmic components are redistributed, while the non-rhythmic 1/f distribution remains unchanged. By calculating the median of the resulting auto-spectral distributions, the power-spectral 1/f component can be isolated and subtracted from the original power-spectrum allowing us to estimate the power-spectral contribution of rhythmic oscillations in the recorded signal.

MRI

MRI images were obtained for Patient-I and for the participants from the Cross-Over group on a 3T whole body scanner (Magnetom Vida, Siemens, Germany) using a standard head coil. Subjects’ heads were immobilized with a vacuum pad to reduce movement artifacts. For each subject, an anatomical volume data set was acquired using a T1-weighted 3D-MPRAGE sequence (TR 2400 ms, TE 2.18 ms, flip angle 8°, voxel size 0.8x0.80.8 mm³s).

Source-Level Analysis

We followed the analysis pipeline detailed in [76] using Fieldtrip. We used frequency-domain beamformers based on the dynamical imaging of coherent sources (DICS) method, which makes use of the cross-spectral density matrix of the data and a forward model of the brain to estimate the spatial filters necessary for source reconstruction. For Patient-I and the participants in the Cross-Over group, we obtained T1-weighted MRI images to construct volume conduction models and build subject-specific forward models. For other participants, we used a realistically shaped three-layer boundary-element template volume conduction model on a 3D grid of dipole locations with equidistant spacing of 8mm, and EEG channel locations were co-registered to the used conduction models. The choice of the frequency of interest for the DICS method was guided by time-frequency analysis within the frequency bands of interest.