A. Reach-to-grasp task illustration and trial timeline. Animals were trained to depress a button for a specific hold time (Monkey H: 0.2 seconds, Monkey N: 1.5 seconds). Upon successfully holding, a door covering the slot and object retracted allowing access to the object. Animals reached through the slot (one of four possible slot shapes on a given trial) to grasp and lift the object. The object had to be lifted a height of approximately 1.5 cm for at least 0.267 seconds (object hold time) to trigger a juice reward. The time period from reach onset to grasp start was designated as “reach duration” and from grasp start to reward start as “grasp duration”.

B. Illustration of the four different slot types that could be presented. Each was designed to elicit a unique grasp type.

C. Schematic of primary motor cortex lesion and electrophysiology recordings from perilesional motor cortex and subcortex.

D. T1 weighted MRI from Monkey N superimposed with subcortex probe location calculated from post- operative CT. See Supplemental Figure 1.

E. Significant improvement in normalized reach duration as a function of days after stroke (linear regression, top: Monkey H, slope = -0.004, t-statistic = -2.276, N = 819, p-value = 0.023, Monkey N, slope = -0.003, t-statistics = -9.575, N = 447, p-value < 0.001)

F. Significant improvement in normalized grasp duration as a function of days after stroke (linear regression, top: Monkey H, slope = -0.086, t-statistic = -7.973, N = 819, p-value < 0.001, Monkey N, slope = -0.034, t-statistics = -6.324, N = 447, p-value < 0.001)

G. Significant improvement in percent correct of trials as a function of days after stroke (linear regression, top: Monkey H, slope = 0.018, t-statistic = 2.596, N = 11 days, p-value = 0.029, Monkey N, slope = 0.007, t-statistic = 4.271, N = 13 days, p-value = 0.001)

H. Trial-averaged, channel-averaged power spectral densities of local field potential (LFP) signals from cortex (left) and subcortex (right) electrode arrays for Monkey H. Colors indicate different days post-stroke.

I. Same H for Monkey N.

J. Example trial from Monkey H, Day 7 post-stroke. Top row shows Index MCP speed (black trace) and extracted behavioral markers. Traces below show LFP channel data filtered in the beta band, with extracted beta burst epochs in red.

K. Same as K, but for an example trial from Monkey H, Day 24 when the animal has recovered.

A. The distribution of the fraction of cortical channels that are bursting for each timepoint that is part of a cortical burst event. i) Each row correspond to a post-stroke day, and darker colors indicate higher occurrence. The same data is represented in (ii) as a histogram with different colors correspond to different days post stroke. The vertical dashed black line indicates the threshold calculated to separate events into local, independent burst events and global, synchronous burst events (“local_global_threshold”)

B. The fraction of total cortical burst events that are less than the local_global_threshold is significantly correlated with recovery index for both animals (Monkey H: slope = 0.167, t-statistic=3.324, N = 10, p-value = 0.010, Monkey N: slope = 0.164, t-statistic = 2.393, N = 9, p-value = 0.048).

C. Fraction of timepoints when cortical burst events coincide with subcortical bursting is significantly higher for global compared to local events (Monkey H: frac. timepoints during global events with subcortical bursting = 0.651, frac. timepoints during local events with subcortical bursting = 0.240, N global timepoints = 1,979,463, binomial test p-value < 0.001, Monkey N: frac. timepoints during global events with subcortical bursting = 0.459, frac. timepoints during local events with subcortical bursting = 0.361, N global timepoints = 913,342, binomial test p-value < 0.001). Error bars are s.e.m.

D. Day-averaged normalized amplitudes of global bursts are significantly higher than local bursts across days (LME model with day as a random effect: Monkey H: slope = -0.169, z-statistic = - 20.296, N = 20, p < 0.001, Monkey N: slope = -0.187, z-statistic = -12.505, N = 18, p < 0.001).

E. Day-averaged durations of global bursts are significantly than local bursts across days (LME with day as a random effect: Monkey H: slope = -0.020, z-statistic = -12.783, N = 20, p < 0.001, Monkey N: slope = -0.033, z-statistic = -8.590, N = 18, p < 0.001).

F. Day-averaged frequency of global bursts is significantly lower than local bursts across days (LME with day as a random effect: Monkey H: slope = 0.316, z-statistic = 15.833, N = 20, p < 0.001, Monkey N: slope =0.197, z-statistic = 1.979, N = 18, p = 0.048).

A. Example of possible spatial configurations of local beta bursts. Left illustrates spatially distributed burst, where the distances between channels active in the burst and spatial burst center are large. Right illustrates spatially clustered burst where the distances between channels active in the burst and spatial burst center are small.

B. Example of the spatial distribution of active channels from local bursts on a single trial (Monkey H, day 21, trial 18). Color of each channel indicates for what fraction of the burst channels are active. Channels with small markers have been excluded due to noise.

C. Individual local burst events have channel distances that are more spatially clustered than expected just based on the number and location of channels that are generally activate across all reach and grasp local burst events (LME model with day as a random effect, testing for mean of distribution being different than zero: left: Monkey H, reach: intercept = -0.426, z-statistic = -15.590, N = 421, pv < 0.001, grasp: intercept = -0.407, z-statistic = -14.200, N = 1745, pv < 0.001, right: Monkey N: reach: intercept = -0.135, z-statistic = -2.116, pv = 0.034, grasp: intercept = -0.185, z-statistic = - 5.335, N = 1179, pv < 0.001).

D. Local burst events have channel distances that are more spatially clustered than global burst events (LME model with day as a random effect: left: Monkey H, reach: slope = 0.222, z-statistic = 11.787, N = 750, pv < 0.001, grasp: slope = 0.178, z-statistic = 22.929, N = 3452, pv < 0.001, right: Monkey N, reach: slope = 0.0978, z-statistic = 2.587, N = 115, pv = 0.010, grasp: slope = 0.162, z- statistic = 21.266, N = 2566, pv < 0.001).

A. The fraction of trials with global beta burst events for each day post-stroke (colored trace) aligned to different behavioral events (x-axis, window of [-0.35, 0.35] seconds around each event). Red box indicates area of notable change in rate of burst events with recovery.

B. Same as (A) but for local beta bursts events.

C. Recovery index is significantly negatively correlated with the mean fraction of trials with global beta events in a [-0.35, 0.35] second window centered at grasp start (Monkey H: slope = -0.305, t- statistic = -7.364, N = 10, p-value < 0.001, Monkey N: slope = -0.123, t-statistic = -2.798, N = 9, p- value = 0.027). Days with LFP and behavior are included.

D. Recovery index is significantly positively correlated with the mean fraction of trials with local beta events in the reach start to juice reward (Monkey H: slope = 0.093, t-statistic = 2.340, N = 10, p- value = 0.047, Monkey N: slope = 0.125, t-statistic = 2.546, N = 8, p-value = 0.044). Days with LFP and rewarded trials are included.

E. There is a significant positive correlation between time between time spent in a global beta burst event and trial length (left: Monkey H: slope = 1.47, t-statistic = 5.096, N = 765 trials, p-value < 0.001), but no significant correlation between time spent in a local beta burst event and trial length (right: Monkey H: slope = 0.34, t-statistic = 1.075, N = 765 trials, p-value = 0.283). Regression is performed on all trials pooled. Data is visualized in percentiles: Trials with no beta burst events are all in first boxplot aligned with “0.0” on the x-axis and percentage of trials is indicated in text above the box plot. Remaining trials are grouped into 10 percentiles based on fraction of time spent in burst event that correspond to the 10 boxplots to the right of the first boxplot. Percentiles are for visualization only, statistics are performed on individual trials.

F. Same as (E) for Monkey N (left: slope = 5.50, t-statistic = 4.839, N = 256 trials, p-value < 0.001, right: slope = 1.69, t-statistic = 1.343, N = 256 trials, p-value 0.181)

G. Hand speed (measure with index finger MCP speed) is significantly faster during local beta burst events than during global beta burst events for both reaching (reach start to grasp start) and grasping epochs (grasp start to last object lift) (left: reaching, LME with day as random effect: slope = 0.934, z-statistic = 3.167, N = 3828 camera frames, p-value = 0.002, right: grasping, LME with day as random effect: slope = 0.259, z-statistic = 3.697, N = 6144 camera frames, p-value < 0.001). Error bars are s.e.m.

H. Same as (G) for Monkey N (left: reaching, LME with day as random effect: slope = 3.379, z-statistic = 2.781, N = 408 camera frames, p-value = 0.005, right: grasping, LME with day as random effect: slope = 0.276, z-statistic = 2.023, N = 2365 camera frames, p-value = 0.043) Error bars are s.e.m.

A. Example of trial-averaged z-scored spiking activity spiking units from day 11 and day 24 aligned to reach start (Monkey H). Bottom shows average z-scored activity across all displayed units.

B. Same data as (A) but aligned to grasp start.

C. There is no significant change in trial-averaged firing rates of units during reaching from early sessions versus late sessions (independent sample t-test, left gray: Monkey H: t-statistic: 1.633, N = 279, p-value = 0.104, right gray: Monkey N: t-statistic = 0.595, N = 339 units, p-value = 0.552), but there is a significant increase in trial-averaged firing rates during grasping from early sessions versus late sessions (independent sample t-test, left purple: Monkey H: t-statistic: 4.067, N = 279, p-value < 0.001, right purple: Monkey N: t-statistic: 2.466, N = 339, p-value = 0.014). Bars are means and s.e.m.

D. There is no significant change in trial-averaged firing variability of units during reaching from early sessions versus late sessions (independent sample t-test, left gray: Monkey H: t-statistic: 1.381, N = 279, p-value = 0.168, right gray: Monkey N: t-statistic = 0.503, N = 339 units, p-value = 0.615), but there is a significant increase in trial-averaged firing variability during grasping from early sessions versus late sessions (independent sample t-test, left purple: Monkey H: t-statistic: 3.412, N = 279, p-value = 0.001, right purple: Monkey N: t-statistic: 2.044, N = 339, p-value = 0.042). Bars are means and s.e.m.

E. There are significant differences between the trial-averaged mean firing rate and variability of firing rate of units during grasping when there are no beta events (gray), local beta events (yellow), and global beta events (red). (paired t-test, left: Monkey H, all N = 279 units, mean firing rate: None vs. local: t-statistic = 6.333, p-value < 0.001, None vs. global: t-statistic = 6.969, p-value < 0.001, Local vs. global: t-statistic = 4.535, p-value < 0.001, firing rate variability: None vs. local: t-statistic = 13.141, p-value < 0.001, None vs. global: t-statistic = 10.329, p-value < 0.001, Local vs. global: t- statistic = 3.464, p-value = 0.002, right: Monkey N, all N = 339 units, mean firing rate: None vs. local: t-statistic = -2.547, p-value = 0.034, None vs. global: t-statistic = 7.081, p-value < 0.001, Local vs. global: t- statistic = 7.812, p-value < 0.001, firing rate variability: None vs. local: t-statistic = 5.449, p-value < 0.001, None vs. global: t-statistic = 11.354, p-value < 0.001, Local vs. global: t- statistic = 4.274, p-value < 0.001). All p-values are Bonferroni corrected for multiple comparisons.

A. Examples of isolated single unit waveform mean +/- 1 standard deviation (left), unit entrainment to global beta events (middle), and local beta events (right). Probability of unit firing is plotted against global and local beta phases. Example units are from Monkey H, day 24.

B. Distribution of modulation z-score values for significantly entrained unit-channel pairs. Units are more strongly entrained to global beta bursts than local beta bursts (independent sample t-test: left: Monkey H: t-statistic: 8.096, N = 9872 significant unit-channels, p-value < 0.001, right: Monkey N: t-statistic: 12.879, N = 5779 unit-channels, p-value < 0.001). Bars span 5th to 95th percentile of data.

C. There is a significant relationship between entrainment to global beta bursts during grasping (modulation z-score) and change in trial-averaged mean firing rate during grasping between global burst times and non-burst times (linear regression, left: Monkey H: t-statistic = -5.960, N = 174, p- value < 0.001, right: Monkey N: t-statistic = -4.856, N = 90, p-value < 0.001).

D. There is no significant relationship between entrainment to local beta bursts during grasping (modulation z-score) and change in trial-averaged mean firing rate during grasping between local burst times and non-burst times (linear regression, left: Monkey H: t-statistic = -0.245, N = 200, p- value = 0.807, right: Monkey N: t-statistic = -0.407, N = 104, p-value = 0.685).

E. There is a significant relationship between entrainment to global beta bursts during grasping (modulation z-score) and change in trial-averaged firing rate variation during grasping between global burst times and non-burst times (linear regression, left: Monkey H: t-statistic = -4.837, N = 174, p-value < 0.001, right: Monkey N: t-statistic = -2.878, N = 90, p-value = 0.005).

F. There is no significant relationship between entrainment to local beta bursts during grasping (modulation z-score) and change in trial-averaged firing rate variation during grasping between local burst times and non-burst times (linear regression, left: Monkey H: t-statistic = -1.671, N = 200, p-value = 0.096, right: Monkey N: t-statistic = 0.071, N = 104, p-value = 0.943).

G. Schematic of neural population dynamics model used to assess temporal predictability.

H. There is a significant fraction change in grasping temporal predictability (population dynamics model R2) during global (red) but not local (yellow) beta bursts (LME with day as random effect: left: Monkey H, none vs. global: slope = -0.018, z-statistic = -2.825, N = 1323 valid trials, p-value = 0.0142, none vs. local: slope = -0.004, z-statistic = -2.015, N = 1396 valid trials, p-value = 0.132, global vs. local: slope = 0.003, z-statistic = 0.692, N = 1151 valid trials, p-value = 1.0, right: Monkey N, none vs. global: slope = -0.027, z-statistic = -3.145, N=363, p-value = 0.005, none vs. local: slope = -0.007, z-statistic = -2.187, N = 371, p-value = 0.0863, local vs. global: slope = 0.003, z- statistic = 0.605, N = 322, p-value = 1.0). All p-values are Bonferroni corrected for multiple comparisons. Bars are mean +/- s.e.m.

A) Neural recordings were obtained from primary motor cortex in animals trained to perform a sequential target-capture task.

B) Distribution of fraction of channels bursting simultaneously during beta bursts. Different lines indicate different behavioral days for each animal. Thresholds used to separate local and global beta events are indicated with horizontal dashed lines.

C) Single trial wrist speed with global and local events highlighted using red and yellow boxes respectively (trial: Monkey B, session 0, trial 55).

D) Wrist speed is significantly lower during global versus local beta events (LME with day as random effect: left Monkey B: slope = 2.63, t-statistic = 11.40, N = 7136, p < 0.001, right Monkey F: slope = 4.61, t-statistic = 5.46, N = 1790, p < 0.001). Bars are mean +/- s.e.m.

E) Channel level activity filtered in the beta band during single trial shown in (C). Global and local burst events are shown in red and yellow respectively.

F) Example of the spatial distribution of active channels from local bursts on a single trial shown in (C, E).

G) Individual local burst events have channel distances that are more spatially clustered than expected just based on the number and location of channels that are generally activate across all reach local burst events (LME model with day as a random effect, testing for mean of distribution being different than zero: left: Monkey B, intercept = -0.148, z-statistic = -10.32, N = 3557, pv < 0.001, right Monkey F, intercept = -0.304, z-statistic = -10.917, N = 867, pv < 0.001).

H) Local burst events have channel distances that are more spatially clustered than global burst events (LME model with day as a random effect: left: Monkey B, slope = 0.158, z-statistic = 43.53, N = 7173, pv < 0.001, Monkey F, slope = 0.122 z-statistic = 15.83, N = 1824, pv < 0.001)

Summary of main findings.

In impaired animals (left), global bursts (red) that are synchronous across perilesional cortex (PLC) and subcortex (Subcx) occur mid-movement and are associated with slower movement speeds, longer trial times, entrainment of single unit activity, a reduction in single unit firing rate and firing variance, and interrupted spiking population dynamics. In recovered animals (right), global bursts during movement are uncommon but local bursts (orange) occur during movement and are associated with the recovered state, typical single unit firing rate and firing variance, and typical spiking population dynamics.

Days used for Behavioral, LFP, and spiking analysis:

Estimating the subcortical probe location

Prior to surgery, cranial 3T MRI images were acquired for both animals using a T1-weighted sequence. After surgery, cranial computed tomography (CT) images were acquired for both animals. The MRI images were aligned and resliced in ACPC coordinates using FieldTrip . Then, both the CT and re-sliced MRI images were loaded into 3D Slicer. The CT images were windowed using the default CT-bone window and were thresholded to extract a volume containing the skull. This volume also contained the implanted microelectrode array in perilesional cortex and support body of the subcortical electrode array.

A, B: The extracted skull and implanted electrode volume (CT, orange) was then manually aligned with the MRI such that the CT skull aligned maximally with the skull in the MRI image. When traversing the axial axis of the scan, the CT skull volume shows excellent alignment to the MRI in the axial, coronal, and sagittal axes in places where the skull was not disrupted due to the surgical intervention.

After alignment, a 3D model of the subcortical electrode array was imported into the 3D slicer image (cyan). The array was manually aligned such that the support body of the array maximally overlapped with the part of the CT volume corresponding to the support body. The silicon shank part of the subcortical electrode array was not visible on the CT. The alignment between the support body and the CT volume was verified by scanning through the axial, sagittal, and coronal axes and can be seen in A, B. The final location of the recording electrodes was then visualized in C (Monkey H), and Fig 1D (Monkey N).

T1-weighted sequences do not visualize subcortical nuclei particularly well. To estimate the structures where the subcortical electrodes were positioned, brain slices from the Brainmaps.org atlas were visualized69. Slices that best matched the cortical anatomy of the coronal MRI slice where the electrodes were placed (shown in C and Fig. 1D) were selected (Monkey N: Atlas slice 8.2mm AP, Monkey H: Atlas slice 10.0mm AP). These slices were within 1-2mm of the AP axis of the animals’ MRI (Monkey N: 9.9mm AP, Monkey H: 10.2mm AP). The electrode trajectories were estimated on the atlas slices by computing the ML and DV coordinates of the final electrode position and the angle of the electrode tract in the MRI. The angle and the final electrode position dictated the electrode trajectory shown in Fig. 1D and C.

These slices reveal that the final electrode positions were either in (Monkey N) or just outside in the white matter (Monkey H), the VL nucleus of the motor thalamus. For the purposes of this study, we simply refer to these electrode recording locations as “subcortical”.

Reach-to-grasp task The CAD model of the reach-to-grasp task.

A) Appearance of the fully assembled reach-to-grasp task from the perspective of the side camera at the start of a trial. The button is highlighted in white. The button houses two square force sensitive resistors (Sparkfun.com, SEN-09376) with wires that are routed to the back of the task (indicated by dashed black lines). The wires are covered by a “wire protector” piece. Animals must first depress the button for a designated hold period in order to open the door to the slot and object to start the trial.

B) Removal of “backsplash” for better visualization of the slot door apparatus (shown in white). A linear slide potentiometer (Bourns PSM60-081A-103B2) is seated behind the backsplash and has a lever that is fixed to a 3D printed door.

C) After the button has been depressed for the designated hold period, the door slides open, revealing the object and the slot that must be reached through in order to grasp the object.

D) Removal of the door and the structural pieces of the task apparatus for better visualization of the “wheel” that has the slot cutouts. The wheel is fixed to a DC motor through a coupler (shown in E). Prior to each trial, the motor spins the wheel such that the correct slot is aligned with the object. A magnetic solenoid’s (Adafuit.com, Product ID = 3992) piston is released through holes on the outside of the wheel to lock the wheel in place and prevent the animal from moving the wheel during a trial. The tripod slot is shown in the current trial.

E) The wheel and other structural pieces of the task have been removed to visualize the DC motor (Adafuit.com, Product ID 4416), the DC motor coupler (Pololu Robotics and Electronics, Pololu item # 1079, Pololu Universal Aluminum Mounting Hub for 3mm Shaft, #2-56 Holes (2-Pack), and the housing that was 3D printed to attach the motor to the task apparatus.

F) Removal of the wheel motor and other structural parts to reveal the “object” (green), and the custom 3D printed apparatus that houses the object (red) and the IR beam break sensors that detect the object height. The object has a threaded rod that twists through its body, and through vertical slots on the housing. This rod and the track for the rod in the housing constrains the motion of the object to be up and down. The housing also contains cutouts for the IR beam break sensors (Adafruit.com, Product ID 2167) that detect if the object has been lifted above the 15mm height (shown in black).

G) View of the re-assembled task with the door open revealing the slot and the object to be lifted. The black horizontal line indicates a line used in behavioral analysis (Methods, Reach-to-grasp behavior, automatic extraction of behavioral timepoints, relevant to identification of the “crossing frame”)

Monkey N full recovery curve

Monkey N underwent an unusual recovery curve. Initially, from days 7 to 18 after stroke, Monkey N exhibited a typical recovery with a reduction in reach-to-grasp time and an increase in percent correct. At Day 21 (red arrow), we noticed a return of his arm and hand impairment concomitantly with a significant reduction in cortical array spiking quality. We performed neuroimaging (CT) but were unable to see evidence of an obvious secondary brain injury or of an infection. We continued attempting to test Monkey N using the reach-to-grasp task outlined in Fig 1AB, but he continued to exhibit an inability to perform the task and poor cortical array spiking signals. At Day 35 (orange arrow) we introduced a new “object” within the same reach-to-grasp task format. The new object had a large ball on top of the existing rectangular object, which made finding the object and grasping it much easier. We continued testing Monkey N on this task until Day 53 at which point he was able to begin performing the original reach-to-grasp task again. We are still investigating the possible causes of his secondary decline.

Because of this unusual recovery curve and because we have been unable to pinpoint the cause of the secondary decline, we have chosen to focus on the initial albeit incomplete recovery curve (days 7- 20) augmented by a few days at the end of his recovery once his performance was back at pre-stroke baseline. Since the purpose of our study is to examine how global and local beta dynamics change with recovery from stroke, it was important to us that our choice of sessions did not modify the overall dynamics of Monkey N’s recovery curve. To test this, we fit a linear regression to behavioral data from all days (gray lines in A, B) and fit a regression to only the sessions we analyzed in Fig 1EFG (black lines in A, B). We tested whether the slopes of these regressions differed significantly (F-test). For both normalized reach-to- grasp time and normalized percent correct there was no significant difference between the slopes of the gray and black lines (A: normalized reach-to-grasp time, p = 0.548, B: normalized percent correct, p = 0.681). See Table S1 for more details about exactly which days of data were used in each analysis.

Monkey H single trials over the course of recovery.

In the same style as Fig. 1JK.

Monkey N single trials over the course of recovery.

In the same style as Fig. 1JK.

Alternative visualization of data in Fig 4AB.

Individual post-stroke days are shown as rows, and darkness of the row indicates higher fraction of trials with global (A) or local (B). Red shading indicates no data present (no trials with successful lifts in Monkey N on day 7).

Properties of single units recorded during early recovery and late recovery.

(A, left) Unit ratings during early vs. late recovery. Units were rated according to SNR: Rating 3: SNR > 3.0, Rating 2: SNR > 4.0, Rating 1: SNR > 8.0. Monkey N had very similarly ranked units early vs. late in recovery whereas Monkey H had a higher proportion of higher ranked units late in recovery. (A, center) Distribution of trough- to-peak amplitude (uV) for early vs. late units. Animals show trends in different directions: Monkey H has higher amplitude units later in recovery whereas Monkey N has higher amplitude units early in recovery. (A, right) Distribution of waveform widths (measured via trough-to-peak time) for animals. Monkey N has a very consistent distribution of waveform widths in early vs. late whereas Monkey H has a higher proportion of very short waveform widths early in recovery compared to late. (B). Examples of waveforms and inter-spike interval histograms from randomly selected units from Monkey H and Monkey N during early (red box) and late (blue box) recovery.