Authors: Rohit Rangwani (1Center for Neural Science and Medicine, Department of Biomedical Sciences, Cedars-Sinai Medical Center, Los Angeles, CA 90048, USA; 2Department of Bioengineering, Henry Samueli School of Engineering, University of California, Los Angeles, Los Angeles, CA 90095, USA), Aamir Abbasi (1Center for Neural Science and Medicine, Department of Biomedical Sciences, Cedars-Sinai Medical Center, Los Angeles, CA 90048, USA), Tanuj Gulati (1Center for Neural Science and Medicine, Department of Biomedical Sciences, Cedars-Sinai Medical Center, Los Angeles, CA 90048, USA; 2Department of Bioengineering, Henry Samueli School of Engineering, University of California, Los Angeles, Los Angeles, CA 90095, USA; 3Department of Neurology, Cedars-Sinai Medical Center, Los Angeles, CA 90048, USA; 4Department of Medicine, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA; 5Lead contact)
Categories: Article
Source: Cell reports
Authors: Rohit Rangwani, Aamir Abbasi, Tanuj Gulati
Brain-machine interfaces (BMIs) offer a viable option for restoring function in patients with motor disabilities post-stroke. Most BMI systems rely on signals from the motor cortex (M1), which is often compromised after stroke. The cerebellum, a subcortical structure involved in motor control, remains an underexplored source for neuroprosthetic control. Using chronic electrophysiological recordings in a rat stroke model, we show that cerebellar neural activity can effectively drive BMI control, performing comparably to M1-driven control. We observed this even in animals with motor impairments post-stroke. Simultaneous M1-cerebellum recordings during cerebellar BMI control revealed that cerebellar “direct” neurons driving the interface were influenced by both local cerebellar and distant M1 neurons. While cerebellar influence remained stable, M1’s interaction with cerebellar direct neurons shifted from longer to shorter timescales after stroke. These findings highlight that cerebellar direct neural control is possible in the stroke brain and reveal changes in M1-cerebellar network dynamics post-stroke.
Brain-machine interfaces (BMIs) allow direct neural control of prosthetic devices through real-time neural activity.^1–14^ Such volitional modulation of a neuronal subpopulation is closely tied to the concept of BMIs.^10,15–20^ Neural plasticity plays a vital role in achieving stable neuroprosthetic control.^7,8,10,16–19,21^ Prior research has primarily engaged intact neocortical structures for BMI applications. BMIs are aimed at restoring function post-injury, which calls for a better understanding of the neurophysiology in injured neural networks.^22–26^ In stroke, viability of the perilesional and ipsilesional cerebral cortex has been explored as a target of neural interfaces;^21,27–30^ however, subcortical structures also have motor representations.^31–33^ The cerebellum (CB) has a well-researched role in movement control,^34–37^ but its efficacy for neuroprosthetic control remains untested. We performed simultaneous motor cortex (M1)-CB recordings to evaluate if cerebellar spiking activity can be modulated to control an external actuator. We also measured neural interactions between M1-CB during CB direct neural control and how this changed in stroke brain versus a healthy brain.
M1-CB have dense reciprocal connectivity,^38^ and there are known changes to activity patterns localized to both CB and M1 with learning^31,32,36,37,39–44^; however, there is also a greater task-related cross-area coordination with learning.^20,32,36,45^ This is supported through the M1 to CB connectivity via pons,^39^ as well as CB to M1 connectivity via the motor thalamus.^20,37^ CB and M1 have direct projections to the spinal cord,^46–48^ but it is difficult to test if M1-CB coordination with learning is due to their interactions or due to their control of a common end-effector. BMIs offer a powerful tool to test how cross-area interactions change with learning, as during “brain control,” actuator movements are causally linked to an ensemble of neurons called “direct neurons”.^15,18–20,49^ As these neurons change their firing properties during neuroprosthetic learning,^15,17–19,50^ the other neurons in the network (indirects) generally become less task associated.^15–17,19,51–54^ We implemented direct CB control (by assigning CB neurons as “directs”) and compared how the neural activity in M1-CB indirect activity interacts with CB direct activity. We also placed emphasis on how these interactions change between healthy and stroke-injured M1.
Our findings show that CB direct neuroprosthetic control was feasible in the healthy and stroke brain, and this was comparable to healthy brain’s M1-driven control. Interestingly, in stroke, CB-BMI was possible even when motor impairments lingered. While we did not study the effects of CB-BMI control on long-term stroke recovery, future work can investigate this. Here, we found that M1 neurons were modulated during CB-BMI control, underscoring M1’s indirect involvement in CB activity. In intact brain, M1 activity predicted CB direct activity at larger time lags, but this prediction occurred at shorter time lags with M1 stroke. Within the CB, we did not see any change with CB indirect prediction of CB-BMI direct activity in intact versus stroke brain. Together, our work shows that CB activity can be used for neuroprosthetic applications and M1-CB interaction changes with stroke that may help improve CB-BMI functionality.
We recorded M1 and CB neural activity (see STAR Methods) as rats performed a neuroprosthetic task. Rats controlled the angular velocity of a mechanical water tube using activity from one or two experimenter-selected “direct” units in the CB or M1 (assigned positive or negative weights; direct+ or direct−, respectively). A linear decoder translated change in direct unit firing rates into the angular velocity of the actuator. The other recorded units in the CB and M1 that were not causally linked to the movement of the water feeding tube were referred to as “indirect” units (i.e., CB or M1 indirects). Decoder weights were kept constant during sessions to exclusively rely on neural learning mechanisms of neuroprosthetic control. Each trial began with an auditory tone and the opening of a plexiglass door, allowing control of the tube from resting position P1 to target position P2 (Figure 1A, see STAR Methods). A trial was terminated if the water feeding tube was not successfully moved to P2 within 15 s after the start of the trial, resulting in an unsuccessful trial. At the end of the trial, the water feeding tube returned to P1, and the door was closed.
We observed that rats were able to exert control over a water feeding tube using direct neuron activity in the cerebellar cortex, and their performance improved within 1–2-h sessions (Figure 1D). Time to successfully complete the task and unsuccessful trials reduced with learning (Figures 1F and 1G). In intact CB-BMI sessions (n = 18, 5 rats), we observed significant increase in the success rate and significant reduction in task completion times (Figures 1F and 1G; time-to-task completion, 8.19 ± 0.49s, 3.90 ± 0.48s, mixed-effects t(34) = −6.62, p = 1.36 × 10^−7^; success rate, 69.64 ± 4.37%, 95.96 ± 1.44%, mixed-effects t(34) = 5.52, p = 1.65 × 10^−5^). Similarly, in intact M1-BMI sessions (n = 20, 7 rats), we observed significant BMI performance changes (Figures 1J and 1K; time-to-task completion, 8.93 ± 0.46s, 3.58 ± 0.26s, mixed-effects t(38) = 5.31, p = 2.50 × 10^−6^; success rate, 80.07 ± 2.67%, 95.46 ± 1.25%, mixed-effects t(34) = 5.52, p = 1.65 × 10^−5^). This showed that proficient cerebellar direct neuroprosthetic control was possible, and it was as good as M1 direct neural control.
We also analyzed the video recordings during BMI session and tracked the water feeding tube position using DeepLabCut (DLC).^55^ We observed that the tube’s movement from P1 to P2 became more direct and consistent in late trials in the session than in the early trials (Figure S1B). Additionally, the angular speed of the tube increased significantly from early to late trials across all animals (Figure S1D, intact CB-BMI: t p = 1.54 × 10^−4^). To confirm that neuroprosthetic control was independent of any correlated body movement, we tracked left-right forepaws and head using DLC. We did not find their movements to be significantly correlated with the tube’s movement across in intact CB-BMI sessions (Figure S2). Our work is consistent with other studies where overt movement reduced as performance on neuroprosthetic control improved.^15,16,20^
Next, we tested the efficacy of stroke CB-BMI. We observed that such control was possible in stroke rats while forepaw impairments lingered (Figures 1N, 1O, and S3). Additionally, stroke CB-BMI performance (n = 26 sessions, 9 rats) was on par with intact CB or M1-driven BMI performance (Figure 1E). Stroke CB-BMI sessions showed significant changes in task performance with learning (Figures 1H and 1I; time-to-task completion, 7.85 ± 0.43s, 4.30 ± 0.32s, mixed-effects t(50) = −6.86, p = 9.70 × 10^−9^; success rate, 80.07 ± 2.67%, 95.46 ± 1.25%, mixed-effects t(34) = 5.52, p = 1.65 × 10^−5^). Interestingly, there was no significant difference in the three BMI groups’ task performance once expert control was learned (Figures 1L and 1M; time-to-task Kruskal-Wallis Х^2^2,61 = 1.59, p = 0.45; success Kruskal-Wallis Х^2^2,61 = 0.69, p = 0.71). In stroke CB-BMI sessions, the angular speed of the tube increased significantly from early to late trials (t p = 4.88 × 10^−5^). Body movements were not significantly correlated with the tube movements in stroke CB-BMI sessions (Figure S4). Thus, we observed that robust cerebellar direct neuroprosthetic control was possible with M1 stroke, indicating that subcortical activity might be viable for direct neural interfaces.
Direct units in the intact CB-BMI sessions experienced a significant change in modulation from early to late trials (41.55 ± 9.37%; t t(44) = −4.44, p = 3.03 × 10^−5^; Figure S5). Direct units in stroke-CB-BMI sessions also experienced a significant modulation (43.70 ± 13.32%; t t(52) = −3.28, p = 9.25 × 10^−4^). In intact and stroke CB-BMI sessions, 75.51% and 78.33% of CB direct units were significantly task modulated, respectively (Figures 2C–2F; see STAR Methods). Studies of M1-BMI have shown similar modulation of M1 direct units.^20,21,56^ Upon analyzing task-related indirect modulation of M1 and CB units, we found that a subset of the indirect units in both showed significant modulation. These results are consistent with other studies showing the contribution of indirect units in neuroprosthetic control.^20,21^ In intact CB-BMI, 43.35% of CB and 41.47% of M1, and in stroke, CB-BMI, 37.21% of CB and 55.12% of M1 units developed significant indirect modulation. Similar indirect modulation is reported in M1-BMI control.^20^ Hence, cerebellar direct neural control recruited indirect units that may have contributed to neuroprosthetic skill learning, and this indirect activity was seen in the broader motor network.
Next, we wanted to study moment-by-moment interactions between the direct and indirect task-related activity during CB-BMI control and how these changed in stroke. To investigate this, we used a regression approach^56,57^ using generalized linear models (GLMs) to predict CB-BMI potent activity from CB and M1 population activity, GLM-id (CB indirects → CB directs) and GLM-1d (M1 indirects → CB directs), respectively (Figures 3A and 3B). CB-BMI potent activity was reconstructed from binned neural data (10 ms–50 ms) and used as the response variable in the GLMs (see STAR Methods). We found that the GLMs were able to predict CB-BMI potent activity using CB and M1 indirect activity (Figures 3C–3F and S6). This showed that the population activity pattern of the indirect CB and M1 neurons predicted moment-by-moment CB-driven BMI activity both in the intact and injured network.
Comparing the averaged cross-validated R^2^ values for different GLM models, we found significant differences in local versus cross-area interactions (GLM-id versus GLM-1d) for both the intact and stroke cohort (Figure 3C; intact cohort (n = 13): GLM-id: 0.34 ± 0.05; GLM-1d: 0.12 ± 0.02, mixed-effects t(24) = −4.28, p = 2.57 × 10^−4^; Figure 3D; stroke cohort (n = 14): GLM-id: 0.33 ± 0.04; GLM-1d: 0.16 ± 0.04, mixed-effects t(26) = −2.17, p = 3.94 × 10^−2^). Comparing local interactions within the CB (GLM-id) for the intact versus stroke cohort, we did not observe a significant difference (Figure 3G; 0.31 ± 0.03 (n = 19; 13 with simultaneous M1 and CB recordings, and 6 with only CB recordings); 0.27 ± 0.03 (n = 26; 14 with simultaneous M1 and CB recordings, and 12 with only CB recordings), mixed-effects t(43) = −1.07, p = 0.29). For cross-area interactions between M1 indirects and CB directs (GLM-1d), there again was no significant difference for the intact versus stroke cohort (Figure 3H; GLM-1d; 0.12 ± 0.02 (n = 13); 0.16 ± 0.04 (n = 14), mixed-effects t(25) = 0.21, p = 0.83). We also observed that the averaged R^2^ values across different GLMs depended on the bin-width (Figure S6). Overall, we did not see a difference in R^2^ trends with different bin-width.
Next, we wanted to establish if M1 activity had any privileged functional relationship with CB direct units or if it had similar predictive power for any randomly selected subset of CB indirect units. For this analysis, we randomly selected CB indirect units to use as “surrogate direct” units in lieu of the CB direct units with matched number as response variable for GLMs. These GLMs used M1 indirect activity to predict the activity of this CB surrogate potent space, GLM-1i (M1 indirects → CB indirects). For each dataset, we repeatedly fitted GLMs with 50 different sets of surrogate CB indirect units. Like GLM-1d, we found that many M1 units have large regression weights. The averaged cross-validated R^2^ values for the stroke cohort were not significantly different from those for the intact cohort (Figure S6G; GLM-1i: 0.31 ± 0.05, 0.38 ± 0.03, mixed-effects t(25) = 1.06, p = 0.30). Comparing cross-validated R^2^ values for both the intact and stroke cohort, we found a significant difference between GLM-1i and GLM-1d models (Figure 3E; intact cohort (n = 14): GLM-1i: 0.31 ± 0.05, GLM-1d: 0.12 ± 0.02, mixed-effects t(24) = −3.41, p = 2.28 × 10^−4^; Figure 3F; stroke cohort (n = 14): GLM-1i: 0.38 ± 0.03, GLM-1d: 0.16 ± 0.04, mixed-effects t(26) = −3.5, p = 1.67 × 10^−3^). This indicated that M1-modulated CB indirect and direct activity and its influence on both these categories of CB units contributed to proficient CB-driven neuroprosthetic control.
To better understand the evolution of the population neural activity in CB and M1 leading to successful execution of the CB-driven BMI task, we examined the temporal structure of the regression weights for the different GLM models. For all GLMs, predictor units had the largest regression weights at multiple time lags (Figures 3I–3N and S6H). Across all datasets for GLM-id, GLM-1d, and GLM-1i, we observed that the distribution of time lags at which any neuron had its largest magnitude regression weight was significantly different from a uniform distribution (GLM-id: intact cohort, two-sample Kolmogorov-Smirnov (KS) D = 0.82, p = 4.28 × 10^−4^; stroke cohort, two-sample KS D = 0.82, p = 4.28 × 10^−4^; GLM-1d: intact cohort, two-sample KS D = 0.64, p = 1.21 × 10^−2^; stroke cohort, two-sample KS D = 0.73, p = 2.50 × 10^−3^; GLM-1i: intact cohort, two-sample KS D = 0.73, p = 2.50 × 10^−3^; stroke cohort, two-sample KS D = 0.55, p = 4.68 × 10^−2^; Figures 3I–3N bottom panels). For GLM-id, large regression weights occurred at small time lags (closer to τ = 0) with a probability higher than the uniform distribution in healthy brains (Figure 3I). This suggested that interactions between CB indirects and CB directs occurred at a shorter range of time lags. This was likely due to these two sets of neurons belonging to the same local population. On the other hand, for GLM-1d in the intact cohort, time lags at which the largest regression weights occurred were spread broadly. This indicated a broader time lag interaction between M1 indirects and CB direct units under healthy conditions (Figure 3J). Interestingly, for GLM-1i, we observed that the time lags with the largest regression weights (with probabilities higher than the uniform distribution) occurred at relatively diverse time lags (Figure 3K). This suggested that interactions between M1 and CB indirects occur at variable time lags.
We observed that neuroprosthetic task performance for the intact and stroke cohorts was similar, and the GLM-id and GLM-1d models for both groups did not have significant differences in predictive power. This led us to test if the timescale of interaction changed locally in the CB or cross-area between M1-CB, following M1 stroke. We did not see any difference in the timescale of interaction for GLM-id and GLM-1i in the intact and stroke cohorts, but we observed a change for GLM-1d (Figures 3I–3N). We confirmed this statistically and observed that intact versus stroke GLM-id’s and intact versus stroke GLM-1i^’^s highest regression weight time lag probabilities were not significantly different (Figures 3I and 3L; permutation test, intact versus GLM-id: p = 0.55, effect size = 0.01; Figures 3K and 3N; permutation test, intact versus stroke GLM-1i: p = 0.18, effect size = 0.08). We observed that the timescale of interaction for GLM-1d model changed significantly from a broader timescale of interaction in the intact animals to a shorter time lag interaction in the stroke animals (Figures 3J and 3M; permutation test, intact versus GLM-1d: p = 9.99 × 10^−5^, effect size = 0.10). This indicated that while neuroprosthetic task performance did not change with M1 stroke, the local and cross-area prediction of CB-BMI potent activity remained similar. However, the timescale of interaction between a specific subset, i.e., M1 indirects and CB directs, was altered. This may be indicative of compensatory dynamics in M1-CB after stroke for efficient CB-driven BMI control.
Next, we wanted to assess the contribution of a specific cell type in the cerebellar cortex in CB-BMI potent activity—the Purkinje cells (PCs)—which are one of the principal cell types in CB. We identified putative PCs using simple spike neural spiking features^58^ (see STAR Methods; Figure S7A). We further analyzed these putative PCs using P-Sort to identify complex spike (CSpk) and confirmed them with characteristic simple spike (SSpk) pauses (Figures S7B–S7D). These confirmed PCs were limited in numbers for regression analysis (n = 13). Confirmed PCs and putative PCs were pooled together and referred to as PCs in further analysis. We used these PCs’ activity (Figures 4A and 4B) for predicting the CB direct activity using GLMs, GLM-pd (CB PCs → CB directs). Similar to GLM-id, we found that GLM-pd had a significant predictive power for all the sessions. Comparing the R^2^ values for GLM-pd for the intact versus stroke cohort showed the same trends as observed for GLM-id, with non-significant change in the average R^2^ values for the stroke cohort (Figure 4C; 0.17 ± 0.04, 0.18 ± 0.03, mixed-effects t(33) = −0.68, p = 0.50). These results indicated that the putative PC indirect population in the neuroprosthetic task developed modulation similar to other cerebellar cortical indirect population and changed moment-by-moment in a similar manner. The timescale of the interaction between PC indirects and CB direct neurons for the GLM-pd dataset was similar to those for the GLM-id dataset, with shorter time lags of influence (Figures 4D, 4E, 3I, and 3L; permutation test, GLM-pd versus GLM-id: p = 0.13, effect size = 0.05; p = 0.10, effect size = 0.04). There were no significant differences between intact versus stroke GLM-pd highest regression weight time lag probabilities (Figures 4D and 4E; permutation test, intact versus GLM-pd: p = 0.18, effect size = 0.09). This suggested that CB PC indirects behaved in a similar way as other CB indirects in supporting CB direct activity.
We also checked if M1 units were predictive of these putative PC units using GLM, GLM-1p (M1 indirects → CB PC indirects). Like GLM-1i, we generated surrogate potent space using putative PCs, which were used as response variables for regression models. Comparing the R^2^ values for these models for the intact and stroke cohort, we did not find any significant difference (Figure 4F; 0.28 ± 0.05, 0.38 ± 0.04, mixed-effects t(23) = 0.94, p = 0.36). The timescale of the interaction between M1 and CB PCs for the GLM-1p dataset was not statistically different from that for the GLM-1i dataset in both intact and stroke brains (Figures 4G, 4H, 3J, and 3M; permutation test, GLM-1p versus GLM-1i: p = 0.76, effect size = 0.02; p = 0.63, effect size = 0.02), indicating that M1’s timescale of influence over CB PC activity or other indirect activity was not different. We observed that intact versus stroke GLM-1p^’^s highest regression weight time lag probabilities were not significantly different (Figures 4G and 4H; permutation test, intact versus GLM-1p: p = 0.76, effect size = 0.01).
Our work here demonstrated the feasibility of CB-BMI in stroke, which was as good as healthy brain control. However, we report M1-CB interaction changes in intact versus stroke CB-BMI. M1 indirect activity that broadly modulated CB direct activity in the intact brain was changed to a narrower timescale, indicating that this top-down control was altered in stroke. The local CB indirect activity’s interaction with direct activity remained unchanged, indicating that CB-BMI largely relies on local activity within CB with a broader, non-specific input from M1. Our findings here may help improve CB-BMI functionality in future.
The BMI paradigm used here allows for the selection of neural activity in one region to be casually linked to the task and analysis of interactions between connected regions during task execution. Leveraging this, we were able to examine the cross-area communication in the M1-CB network that facilitated CB-BMI control. While it is known that M1-CB interactions are needed for movement control,^32^ and even for M1-driven BMI control,^20^ it is not well understood how task-related M1 activity interacts with task-related CB activity during CB-BMI control, or how these interactions change in injury. We found that M1 provided a neuroprosthetic task-specific but temporally imprecise input to CB, and internal CB dynamics were stronger but temporally limited in the intact brain. This leads to the interpretation that CB direct activity contained multiplexed activity, which coordinated with local and cross-area activity. Timescales of these interactions might correspond to their function. With M1 stroke, the remaining M1 indirect neural activity likely compensated for lost M1 activity, and the timescale of influence became shorter. This might be indicative of the compensatory plasticity mechanisms in the M1-CB network post-stroke.
CB has different cell types that have been studied for their specific roles in motor control and learning.^34,35,59–62^ One of the principal cell types in the CB is the PC. We also studied how PCs participated in CB-BMI and their relationship with CB direct activity. In this study, we saw that the CB neural subpopulation comprising only indirect PCs showed similar task-specific modulation to the overall CB cortical indirect population (Figures 4A and 4B). Prediction of CB direct activity using PC indirects or all-cell type inclusive CB indirect activity had a similar timescale of interaction. Further, M1 prediction of PC activity showed similar trends as M1 prediction of all-cell type inclusive CB indirect activity. Our work suggests that PC population had the same functional relationship with CB direct activity as other recorded cells in the cerebellar cortex. Interestingly, our post-hoc analysis showed that when PC was part of direct CB-BMI activity, there was no significant difference in CB-BMI performance (Figure S8). This establishes proof-of-principle that PC activity may be used for CB-BMI control.
A major afferent input to the cerebellar cortex is through the cortico-ponto-cerebellar feedforward pathway.^36,38,43,63,64^ This pathway is theoretically linked to task-relevant dimensionality expansion in the cerebellar cortex that aids in learning,^65–67^ which has recently also been confirmed experimentally.^36^ Past work has shown the emergence of band-limited oscillatory dynamics in M1 and CB during motor skill learning.^20,32^ Hence, it is likely that M1 had a modulatory influence on CB direct volitional activity here. In addition to the contralateral M1, the cerebellar cortex also receives projections from the contralateral striatum and basal ganglia.^68,69^ It is likely that these inputs may also have had a modulatory influence on volitional cerebellar control, which could be investigated in the future.
One important aspect of BMI research is interfacing with injured brain networks to restore motor function. CB implants are becoming more common in preclinical and clinical studies.^70–72^ Recently, new investigations have further elaborated CB’s role in motor control,^20,39,58^ which warrants using CB neural activity for BMI research. Interestingly, recent clinical trials have targeted deep nuclei in the CB for deep brain stimulation for the rehabilitation of patients after stroke.^72^ These studies support investigations of CB-BMIs. Our findings here elaborate M1-CB interactions during efficient CB-BMI control.
The main goal of this study was to assess whether CB-BMI control is feasible after M1 stroke. We showed that this control was possible and as good as intact brain control. We did not include a sham stroke group; our healthy control experiments involved similar M1-CB neural implant procedures. Another key objective was to test whether CB control was viable when forelimb deficits persisted post-stroke. While neuroprosthetic training can be seen as rehabilitation training, whether it affected the status of upper limb impairment can be investigated in the future.
Another limitation of this work is that we did not preselect specific CB cell types for neuroprosthetic control. However, we showed that control was possible from PCs (Figure S8). More work in the future can test specific cell-types and CB-BMI performance. Finally, we cannot ascertain if our recorded M1-CB neurons had a connection. Prior literature has established anatomic connectivity of these two areas (the cortico-ponto-cerebellar pathway^36,38,43,63–67^); our goal here was to study population-level interactions in these areas in intact and stroke brain. Future work can study this further in cell-to-cell M1-CB pairs with confirmed synaptic connectivity.
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Tanuj Gulati (tanuj.gulati@csmc.edu).
This study did not generate new unique reagents.
All animal procedures were performed according to the protocols approved by the Institutional Animal Care and Use Committee at Cedars-Sinai Medical Center, Los Angeles. Adult male (4–7 month old) Long-Evans rats were used for this study (n = 19, 350–600g, Charles River Laboratories). Animals were housed on a 12 h light and 12 h dark cycle (photoperiod from 7 a.m. to 7 p.m.) in a climate-controlled vivarium. All animals were pair-housed initially until they were separated for experimentation. Animals assigned to the stroke cohort were separated and single-housed at the start of training on the reach task before the surgical procedure. Animals were divided into 3 cohorts for experiments; the first cohort was intact animals (i.e., animals with non-injured brain) that were trained to exert CB-driven BMI control. In these animals, the electrodes were implanted in the cerebellar cortex and contralateral M1, except for one where only CB was implanted (i.e., n = 4 rats with M1 and CB implants and n = 1 with only CB implant). There was a second cohort of intact animals with non-injured brains that were trained to perform the M1-driven BMI task (n = 7 rats). Two animals that received M1 and CB implants contributed to both-M1-driven and CB-BMI sessions. The third cohort comprised of injured rats that received a photothrombotic stroke in the forelimb area of M1 contralateral to their preferred/dominant forelimb (which was discerned through their handedness on a reach-to-grasp task). All these animals received electrode implants in the cerebellar cortex contralateral to stroke M1 (n = 9 rats). Out of nine animals, 5 animals received implants in perilesional M1 as well. Our animal numbers and use of male rats are comparable to those used in other neuroprosthetic studies.^20,56^ We studied key neural interactions in a preclinical study, without consideration to covariates of sex. Future clinical trials of CB-BMI can assess this difference in diverse populations.
Surgical procedures for electrode implants and photothrombotic stroke were performed using sterile techniques under 1–5% isoflurane. Surgery involved cleaning and exposure of the skull, craniotomies, electrode implantation, preparation of the skull surface using adhesive cement (C & B Metabond, Parkell, NY) followed by implantation of the skull screws for referencing and overall head-stage stability. The analgesic regimen included the administration of 0.1 mg/kg body weight buprenorphine or 0.65 mg/kg extended-release buprenorphine (Ethiqa). Rats were also administered 0.1 mg/kg body weight dexamethasone and 33 mg/kg body weight sulfatrim for 5 days after the surgery. Post-surgery, animals were allowed to recover for at least 5 days before further behavioral training.
For electrophysiology recordings in the M1, we implanted 32-channel microwire arrays (33 μm/50 μm polyimide coated tungsten microwire arrays, Tucker-Davies Technologies (TDT)) in left/right M1 (based on the handedness in the reaching task) in the upper limb area centered 0.5 mm anterior and 3 mm lateral to bregma.^20,32,44^ Arrays were lowered down to ~1200–1500 μm.
In the same surgery, for electrophysiologic recordings of the cerebellar cortex, we implanted shuttle drivable 64 channel silicon probes (Cambridge Neurotech) through a craniotomy centered at 12.5 mm posterior and 2.5 to 3 mm lateral to bregma.^20,32^ This was to target the Simplex, Crus I and Crus II areas of the CB. Neural activity in these CB areas (including our own recent work) has shown modulation during upper limb motor behavior such as reaching and responsivity to cortifugal fiber and forelimb stimulation.^20,32,45,76–78^ Shuttle mounted CB probes were moved across days and recordings were done at depths ranging from 1.5–4 mm. Both probes shared a common reference and ground with wires wrapped around screws inserted in the skull.
For the stroke cohort, we used the photothrombotic (PT) induced stroke model.^21,44^ PT stroke targeted the forelimb area in the contralateral M1. For the PT stroke model, Rose Bengal dye was injected into a femoral vein using an intravenous catheter. Next, to induce the stroke, the surface of the brain (within M1) was illuminated with white light (KL-1500 LCD, Schott) using a fiber optic cable for 20 min. We used a 3-mm aperture for stroke induction (in the M1 area craniotomy based on stereotactic coordinates; centered at 2 mm anterior and 2.5–3 mm lateral from bregma) and covered the remaining cortical area with aluminum foil mask to prevent light penetration.^21,79,80^ We have extensively published using this model, and it results in well-circumscribed lesions of ~7mm^3^.^21,79–81^ This stroke model has been widely used to assess the status of forelimb recovery and plasticity in the brain networks post-stroke.^82,83^ Post-stroke, a probe was implanted in the perilesional cortex immediately anterior to the stroke site, centered ~3–4 mm anterior and 2.5–3 mm lateral to bregma. (KwikSeal), followed by a layer of dental cement. The craniotomy or implanted electrodes were covered with a layer of silicone
For M1 recording, we used TDT microwires electrode array with 8×4 (32) channels polyimide-insulated tungsten electrodes. We used electrodes with wire diameters of 33 or 50 μm in configuration 8×4 with electrode separation distance at 250 μm and row separation distance of 375 μm. For cerebellar recordings, we used 64 channel silicon probes from Cambridge Neurotech (H-10 or P-1/P-2 or E-1/E-2 design/configuration) as these have a smaller footprint for implantation and allow depth recording. These different probes had slight differences in specifications in different each probe had 64 recordings sites with 2 or 4 shanks containing 32 or 16 channels on each shank respectively. These shanks were separated by 250 μm spacing and recording sites on each shank ranged from 200 to 330 μm. The tips on these probes’ shanks were sharpened to allow easier implantation and recording from cells across different layers. These probes were mounted on drivable shuttles to allow for depth change for different sessions.
We recorded single units (neural spikes) using a 128-channel RZ2 system (TDT). Spike data was sampled at 24,414 Hz. Zero insertion force (ZIF) clip-based digital head-stages with high impedance (~1 GΩ) were used to interface the ZIF connector and the Intan RHD2000 chip that uses 192X gain. Only clearly identifiable units with good waveforms and high signal-to-noise were used for online real-time BMI. All the neural data was recorded for offline analysis. Behavior related timestamps (i.e., trial onset, trial completion) were sent to the RZ2 analog input channel as a fixed duration pulse using an Arduino microcontroller board and synchronized to neural data.
We have used the term ‘unit’ to refer to the sorted spike recordings from both microwire arrays and silicon probe recordings. For both recordings, we used an online sorting program (Synapse, TDT) for neuroprosthetic control.^20^ This was necessary to have well-isolated units in real time whose activity was projected to the decoder for BMI control. This sorting was done after referencing to the mean of all channels excluding the broken channels. We used waveform shape and the presence of refractory period in the inter-spike interval (ISI) to judge quality of isolation. Specifically, a voltage-based threshold was set based on visual inspection for each channel that allowed for best separation between putative spikes and noise; typically, this threshold was at least 3.5 standard deviations (SD) away from the mean. Events were time-stamped and waveforms for each event were peak aligned. K-means clustering was then performed across the entire data matrix of waveforms. Automated sorting was performed (1) first over-clustering waveforms using a K-means algorithm (i.e., split into many mini-clusters), (2) then a calculation of interface energy (a nonlinear similarity metric that allows for an automated decision of whether mini-clusters are actually part of the same cluster), and (3) aggregation of similar clusters. We conducted offline spike sorting in Spyking Circus for all channels which were not involved in online BMI control.^73^ For this, recordings were referenced again using a mean of non-broken channels. The sorted clusters from Spyking circus were manually curated to mark ‘good’ and ‘noise’ units and merge similar clusters. These units were the indirect units that were recorded in M1 and CB. Examples of representative units identified are shown in Figure S10.
For analyzing the unit stability, the units classified as ‘good’ after the curation step of sorting using spyking circus were used. We analyzed the spike shape stability across the full BMI session for CB and M1 units. We did this by correlating (Pearson correlation) the mean spike waveform in the early and late trials for each unit in all the sessions. Waveform for all the units were significantly correlated from the early to the late trials for both the intact and stroke cohorts (intact cohort, early to late trials, median Pearson correlation values for different unit categories CB direct 0.997 (n = 49); CB indirect 0.892 (n = 313); M1 indirect 0.997 (n = 521); and stroke cohort, early to late trials median Pearson CB direct 0.996 (n = 61); CB indirect 0.935 (n = 835); M1 indirect 0.994 (n = 1082).
Next, we used a custom code in MATLAB to calculate isolation distance^84^ metric to quantify unit separation of the isolated units. For this, all sorted units (including ‘good’ units utilized in further analysis as well as ‘noise’ clusters) were used. These units utilized in analysis showed good separation using the isolation distance metric^84,85^ (intact cohort, median isolation distance, CB indirect 4390.5; M1-indirect 237.58; and stroke cohort, median isolation distance, CB-indirect 2983.3; M1-indirect 451.72, see Figure S11).
Sorted neural data from Spyking circus was analyzed using a custom Python script based on an algorithm to identify putative Purkinje cells (PCs) using simple spike features. We used the Elephant^74^ - electrophysiology analysis toolkit library functions for neural data analysis in Python. Single units from the cerebellar cortex with a firing rate >40 spikes per second, CV2 of >0.20 and MAD of <0.008 were labeled as putative PCs^58^ (Figure S7A). 55.15% of the 1893 total sorted units were classified as putative PCs using this criterion. These PCs were visualized in a two-dimensional PCA space along with other non-PC single units (Figure S7B).
We further analyzed the channels containing the units which were labeled as putative Purkinje (from above analysis) using PSort^75^ to look for the presence of complex spike (CSpk)-aligned simple spike pause and characteristic simple and CSpk waveforms. We confirmed CSpk in 13 channels of our putative PCs identified through the method above. These 13 confirmed CSpk came from 9 sessions across 6 rats for CB-BMI in the intact and stroke cohorts.
For our GLM analyses, we pooled confirmed and putative PCs, and referred to them simply as PCs, because simple spike statistics are unique to PCs relative to all other known cerebellar cortical cell types. We note the caveat that future studies may uncover new cell types that are inadvertently included in our analyses as PCs. We analyzed these PCs to check for any task-related modulation similar to CB indirect activity for both the intact and stroke cohorts (Figures 4A and 4B). We used PCs in GLM regression analysis to examine the interactions between the PCs and CB direct units during the CB-BMI task in the intact and stroke cohorts (GLM-pd). We used GLM to model the interaction between M1 and PC activity in CB (GLM-1p) and compared these in CB-BMI in the intact and stroke cohorts. Similar to GLM-1i, for GLM-1p we created a surrogate BMI potent space using randomly selected matched numbers of CB PCs to stand in for the true CB direct units.
After recovery from surgery, animals were typically acclimated to a custom plexiglass behavioral box (Figure 1A) for one session. The box was equipped with a slit (covered with a door at one end) that served as a drinking zone. During acclimatization an ‘auto-task’ was performed, where rats were placed in the box and water (from the water feeding tube illustrated in Figure 1A) was provided to them at fixed intervals. For the neuroprosthetic task, rats asserted volitional control of the water feeding tube using CB/M1 neural activity, as described in the next section. Each trial started with the tube at P1, and when it reached P2 a water reward was given. The trial started with door open accompanied by an audio tone, and ended with 2 audio tones (Figure 1A, bottom panel). A successful trial required movement of the tube to P2 within 15s. Behavioral sessions were typically conducted for 1–2 h. Recorded neural data was processed in real-time using custom routines in Matlab R2018b. Processed neural data served as the control signal for the angular velocity of the feeding tube. The rats performed an average of ~84 ± 8 trials in a CB-BMI session in the intact cohort; and ~86 ± 8 trials in a CB-BMI session in the stroke cohort (e.g., Figures 1D and 1E). Out of 44 total sessions for CB-driven BMI (intact and stroke cohort), in 39 sessions, we recorded videos of the rats during the BMI training using a 30-fps camera (machine vision color camera, TDT). During experimentation periods, we monitored the body weights of the rats to ensure that the weight did not drop below 90% of the initial weight.
For BMI training sessions, we typically selected single or multiple units on one or two CB/M1 channels as ‘direct’ units. The neural activity of these direct units was used to control the angular velocity of the feeding tube. If we only selected one channel for BMI control, then its neurons were associated with the positive unit weights (direct+). If two channels were selected, units on one channel were associated with the positive unit weight (direct+) and the units on the other were associated with negative unit weight (direct−). We did not select a particular cell type for BMI control apriori and only based this decision on spike quality. Notably, we binned the spiking activity for these units into 50 ms bins. We then calculated a mean firing rate for each neuron over a 30 s baseline period. The mean firing rate was then subtracted from its current firing rate used for controlling the angular velocity of the water feeding tube during the session. Firing rate (FR) was converted to angular velocity using a linear
Θv=CG+r+(i)-G-r-(i)
where Θv was the angular velocity of the feeding tube, r+(i) and r-(i) were firing rates of the direct units, and G+ and G- were fixed unit weights. C was a fixed constant (gain) that scaled the firing rates to angular velocity. C was experimenter-defined for each session and remained constant after initialization. The animals were then allowed to volitionally control the feeding tube via modulation of neural activity. The tube started at the same position at the start of each trial (P~1~ in Figure 1A). Computed angular velocity at every step was added to the previous angular position at each time step (50 ms). During a trial, the angular position that was controlled from the CB/M1 direct activity had the limits of 0° (P~1) to 45° (P2). If the tube was controlled successfully to cross the threshold of the “target position” (P2~ in Figure 1A), then a water reward was delivered at the final resting position set to 62°. Rats were not required to maintain the pipe in a very specific range of position for a period. The water pipe remains at the resting position for 1 s to allowing rats time to drink water before returning to initial position, P~1. In the beginning of a session, most rats were unsuccessful at bringing the water feeding tube to position P2*~ or took a long period of time to correctly position the tube. Rats steadily improved control and reduced the time to completion for the task during a session. Multiple learning sessions were obtained from each animal using different sets of direct neurons. Consistent with past studies, we found that incorporation of new units into the control scheme required fresh learning.^8,10,21,86^
In two CB-BMI sessions in the stroke cohort, our post-hoc analysis revealed the presence of Purkinje cells in the decoder ensemble. In these sessions, the BMI task performance was not significantly different from other BMI sessions’ task performance (BMI sessions with PC amongst decoder time to task completion, 5.87 ± 1.33s, 4.25 ± 1.16s; all other BMI time to task completion, 8.01 ± 0.43s, 4.30 ± 0.33s; Kruskal-Wallis test, change in time to task completion (%), Х^2^1,24 = 3.34, p = 0.07; Figure S8).
All rats in the stroke cohort were acclimated to the reach-training behavioral box (Figure S3A inset) and a hand/paw preference was established for the reach task. These rats were then trained on a reach-to-grasp task before the stroke surgery and electrode implantation. The reach training behavioral box was automated and controlled by an Arduino microcontroller board through a custom MATLAB script.^32^ In each trial, a food pellet was dispensed on a pellet tray, followed by an alerting beep indicating the start of the trial. Rats had 15 s to reach through the slot in the box, grasp and retrieve the pellet for a successful reach. All trials were captured by video through a camera placed on the side of the behavioral box (Basler, Germany). These videos were manually marked to identify the success rate of the rats on the reach-to-grasp task for each session. Rats were again evaluated on this task post-stroke, to assess the level of limb impairment after the injury. This post-stroke training was initiated at least 5 days following the induction of the stroke.
After completion of the experiments, rats were deeply anesthetized with isoflurane (4–5%), then exsanguinated and perfused with 4% paraformaldehyde (PFA). The brains were extracted and stored in 4% PFA for up to 72 h. The brains were then transferred to a solution of 30% sucrose with 0.05% sodium azide and stored for sectioning. We obtained sagittal or coronal sections of the brain using a cryostat (Leica, Germany) and stored them in phosphate-buffered saline for imaging. These sections were stained with cresyl violet (Nissl staining) for stroke tissue assessment in M1 (Figure S9A). The location and depth of the silicon probe in the CB was traced by DiI depositing on the electrodes before their implantation and by looking afterward at the fluorescent dye present in the histological slices (Figure S9B). Sections were mounted on slides and imaged using a microscope (Keyence, Japan).
Behavioral analysis was performed in MATLAB (R2019a/R2022a) with custom-written routines. A total of 64 neuroprosthetic training sessions were recorded from 19 rats on which behavioral task performance was analyzed. These included 18 sessions of CB-driven BMI from 5 rats in the intact cohort; 26 sessions of CB-driven BMI from 9 rats in the stroke cohort; and 20 sessions of M1-driven BMI from 7 rats in the intact cohort. Trials in each session were divided into 3 equal parts; early trials (first one third of trials in the session) and late trials (last one third of the trials in the session) were compared for behavioral analysis (Figures 1D and 1E). We compared changes in task performance within a session by comparing success rate and time to task completion in early trials versus late trials for all three M1 driven BMI (intact), CB-driven BMI (intact), and CB-driven BMI (stroke) (Figures 1F–1K). We also compared eventual proficiency in the task (time to task completion and success rate in late trials) in these groups using Kruskal-Wallis test with multiple comparisons (Figures 1L and 1M). We used linear regression to evaluate links between reach success and improvement/rate of learning in BMI sessions (Figures 1N and 1O). Best-fit line derived using linear regression was obtained for BMI performance versus reach success rate.
Neural activity used in the linear decoder for neuroprosthetic control was referred to as direct activity. The neural activity recorded from the CB or M1 and not directly used for controlling the water feeding pipe was designated as indirect activity for analysis. For quantifying modulation of direct units from early to late trials, firing rates in the period including 1 s before task start and 4 s after trial start were compared. We also checked for task-related significant modulation for direct and indirect units. The units were defined as significantly modulated if their peak (minimum or maximum) modulation in a 2.4 s (2.2 s before and 0.2 s after task completion) period around task completion was at least 4 S.D. away from the baseline (computed from the averaged activity in the 2 s period before task start). For Figure 2, population unit activity was obtained by analyzing the firing rate for all significantly modulated units averaged over all trials, then normalized using the MATLAB function normalize. These units were then sorted by their peak values.
CB BMI potent space was obtained from online sorted data used by the linear decoder (to account for the causal relationship of these direct neurons with BMI). This data was binned at multiple binwidths (10 ms, 25 ms, and 50 ms) for regression analysis. BMI potent space activity was calculated as the difference between summed CB direct+ activity and summed CB direct– activity (+/− being positive or negative unit weight associated directs), which was used as the response variable for generalized linear models (GLMs). For sessions with CB direct+ activity only, that alone constituted the CB BMI potent activity.
We used GLMs for regression analysis, using the MATLAB function fitglm to predict CB BMI potent space activity from M1 and CB indirects. This generated generalized linear regression models with linear model specifications (containing an intercept and linear term for each predictor) and fitted using a normal distribution for the response variable. For predicting CB direct activity from CB and M1 indirects, CB BMI potent space activity was used as the response variable for GLM models.^20,56^ Predictors were binned firing rates of M1 and CB indirect units, where each neuron appeared more than once with variable time lags ranging from −50 to +50 ms relative to the BMI potent activity. Such horizontally stacked neural data corresponding to each trial in a session was used as the predictive variable for GLMs. GLMs were fitted to neural data binned at 10, 25, and 50 ms. In every session, for each model for different binwidth, a cross-validated R^2^ value was computed by splitting each session’s data into 9-fold for training and 1-fold for test; this was repeated 10 times. R^2^ values were computed between the true response variable and the model output. The cross-validated R^2^ values reported are the average across all 10 combinations of testing/training data (Figure S6). For predicting CB indirect activity from M1 indirects, a “surrogate BMI-potent space” was created from CB neural activity by randomly selecting matched numbers of indirect units (CB indirects) to stand in for the true direct units (CB direct+/−). The difference of summed activity in the positive and negative pools obtained was used as the response variable. This process was repeated for 50 choices of such units per dataset, and average R^2^ values were reported along with the error bars (Figure S6). All further statistics were performed on the 50 ms GLMs to compare the different models R^2^ and 10 ms GLMs were used for timescale analysis.
For timescale analysis, regression weights were calculated at different time lags (Figure S6H) for all GLMs. Absolute value of regression weights was normalized to each unit’s maximum value (Figures 3I–3N top panels) and histograms of the τ values with the largest magnitude weight (τmax) was generated (Figures 3I–3N middle panels). For all GLMs in a group, probability histogram of these time lags that had the largest magnitude regression weight was generated for comparison across groups (Figures 3I–3N bottom panels). Kernel density fit was obtained for these probability histograms for GLMs fitted to neural data binned at different bins (10, 25 and 50 ms).
We performed marker-less tracking of the feeding tube and the forepaws, nose and head of the rats using DeepLabCut (DLC).^55^ For DLC training, we used resnet50 model. 15 s videos were generated for each trial and used for training the DLC model. For each trial in all the sessions, we calculated correlations between the trajectories of the feeding tube and the forepaws, nose and head using the pearsonr function of scipy library in Python and compared the average Pearson’s R and p values for early versus late trials (Figures S2 and S3).
Statistical analyses were completed in MATLAB and python. The linear mixed-effects model (implemented using MATLAB fitlme) was used to compare the differences in behavior performance (Figures 1F–1K) and GLMs R^2^ (Figures 2C–2H and 4C–4F) These mixed models account for the fact that units or sessions from the same animal are more correlated than those from different animals and are more stringent than computing statistical significance over all units and sessions. For comparing the DLC generated trajectories and velocities from early to late we used two-tailed t test from scikit library in python (Figures S1, S2, and S4). Success rates and time to task completion for late trials for different cohorts were compared using Kruskal-Wallis test (Figures 1L and 1M).
We used t-test to compare the firing rate modulation from early to late trials in direct units, for both the intact and stroke cohort. For comparison of the early versus late trials time to task completion for CB-BMI session where direct units included a PC or not, we used t-test (Figure S8A). Kruskal-Wallis test was used to compare the change in time to task completion for the two groups of CB-BMI session including PC in decoder versus the sessions without (Figure S8B).
For determining whether the GLMs had significant predictions, cross-validated R^2^ were compared to a reference distribution of cross-validated R^2^ for GLMs fitted to trial-shuffled data. Response variable trials were shuffled to ensure that the pattern of modulation of neural activity near task completion was retained but the exact moment-by-moment correlation does not exist in the control dataset. This was only done for 50 ms binned GLMs. A total of 10,000 shuffles were performed and the 95^th^ percentile of reference distribution was used to determine significant GLMs. All the data presented for GLM analysis is for models with significant predictive power.
Kolmogorov–Smirnov (KS) test was used to determine if the distribution of time-lags that had the largest magnitude of GLM weight were significantly non-uniform (Figures 3I–3N, bottom panels). Probability histogram for the τmax was considered as empirical distribution and a theoretical uniform distribution was calculated as the mean of the empirical distribution. These distributions were compared using two sample KS tests to determine uniformity. τmax probability histogram distributions for GLM for the intact and stroke cohort were compared using a permutation test^87^ with 10,000 permutations to determine if they were significantly different (Figures 3I–3N).
SUPPLEMENTAL INFORMATION
Supplemental information can be found online at https://doi.org/10.1016/j.celrep.2025.116030.