Propagating patterns of activity across motor cortex facilitate movement initiation

Voluntary movement initiation involves the modulation of neurons in the primary motor cortex (M1) around movement onset. Yet, similar modulations of M1 activity occur during movement planning when no movement occurs. Here, we show that a sequential spatio-temporal pattern of excitability based on beta oscillation amplitude attenuation propagates across M1 prior to the initiation of reaching movements in one of two oppositely oriented directions along the rostro-caudal axis. Using spatiotemporal patterns of intracortical microstimulation, we find that reaction time increases significantly when stimulation is delivered against but not with the natural propagation orientation suggesting that movement initiation requires a precise recruitment pattern in M1. Functional connections among M1 units emerge at movement onset that are oriented along the same rostro-caudal axis but not during movement planning. Finally, we show that beta amplitude profiles can more accurately decode muscle activity when these patterns conform to the natural propagating patterns. These findings provide the first causal evidence that large-scale, spatially organized propagating patterns of cortical excitability and activity are behaviorally relevant and may be a necessary component of movement initiation.


Introduction
The initiation of voluntary skeletal movements is often associated with large-scale modulations of neurons in the primary motor cortex (M1) 1,2 . It remains a mystery, however, why similar modulations in M1 neurons during movement preparation, visual observation of action, and motor imagery occur in the absence of movement [3][4][5][6][7] . One possibility is that modulation is weaker during these conditions and fails to exceed non-linear, thresholds in downstream spinal circuitry, thus resulting in no movement 5,8,9 . Yet, preparatory activity as measured during an instructed-delay paradigm often equals or exceeds movement-related activity particularly in premotor cortex which is known to provide substantial projections to the spinal cord 10,11 . Another possibility is that there is an active inhibitory gating mechanism preventing movement initiation that is released when movement needs to begin 12 as has been shown during REM sleep 13,14 . There is, however, no evidence that such a gate exists in the context of voluntary skeletomotor behaviors.
A third possibility argues that the motor cortex resides in a very high dimensional state space that allows patterns of ensemble activity related to preparation to evolve in dimensions that do not generate movement (i.e. the "output null space") in contrast to movement-relevant activity patterns that reside in an "output potent space" 15 . Aligned with this idea, we argue that a specific spatiotemporal pattern of M1 excitability and unit modulation is required for voluntary movement initiation.
Oscillations in the beta frequency range (i.e. 15-35 Hz) as measured by the local field potential (LFP) are prevalent in various motor structures including the cerebellum, basal ganglia, and cortical areas such as M1 [16][17][18] . Immediately prior to movement onset, the amplitude of these beta oscillations attenuates (i.e. desynchronize) in M1 when single units begin modulating 19,20 ( Figure 1a). The degree of beta oscillation attenuation is correlated with the magnitude of motorevoked potentials as measured with transcranial magnetic stimulation 21 . Furthermore, the timing of beta attenuation with respect to the go signal varies with reaction time but is fixed with respect to movement onset indicating that beta attenuation is a neural correlate of movement initiation We have previously documented that trial-averaged, beta attenuation timing (BAT) does not occur synchronously across the upper limb area of M1 but rather forms a propagating sequence of cortical excitability 24 . Here, we extract single-trial beta amplitude attenuation profiles and find that attenuation propagates along one of two oppositely oriented directions along the rostro-caudal axis. By applying propagating patterns of subthreshold intracortical microstimulation (ICMS) prior to movement onset, we significantly delayed reaction time when ICMS patterns were incongruent with the natural propagating patterns. Using a network analysis on spiking activities of simultaneously recorded units, we observe that functional connections emerge immediately prior to movement onset that are oriented along the same rostro-caudal axis. Finally, by decoding muscle EMG activity from patterns of beta amplitude profiles and then applying spatial perturbations, we provide evidence that these natural propagating sequential patterns lead to normal recruitment of muscle activity.

Single-trial propagating patterns of cortical excitability
Three rhesus macaque monkeys (Mk, Ls, and Bx) were trained on a simple reaction time task to perform planar reaching movements by moving a cursor to visual targets using a two-link LFP signals were band-pass filtered centered at the frequency of peak power in the beta frequency range (15-35 Hz: 31±3 Hz for monkeys Mk and Bx, and 21±3 Hz for monkey Ls; see Methods) from which beta amplitude profiles were computed using the Hilbert transform. Traditionally, beta amplitude dynamics are determined by trial-averaging beta profiles in order to reduce noise 25-27 . However, trial-averaging conceals possible interesting beta dynamics that vary from trial to trial 28 .
To extract single-trial beta profiles, we used an auto-encoder that was trained to generate de-noised versions of the input beta signals (Supplementary Figure 2b). These single-trial, de-noised profiles exhibited characteristic attenuation (i.e. desynchronization) prior to movement initiation, a mesoscopic signature of cortical excitability 21 . More importantly, these beta dynamics revealed spatial gradients that evolved in time indicative of spatio-temporal patterns that propagated in one of two oppositely oriented directions along a rostro-caudal axis in MI (Figure 1c).
For each trial, we quantified the propagation direction of beta attenuation time prior to movement initiation by first determining when each electrode's normalized beta profile dropped below a threshold. We then used linear regression to predict beta attenuation times (BATs) from the spatial locations of each electrode (see Equation 1), where the electrode location on the array (x and y coordinates) was the independent variable and the corresponding BAT was the dependent variable. The regression coefficients could then be used to compute the orientation from early to late BATs (see Equation 2) which we refer to as the beta attenuation orientation (BAO) for that trial. Trials that satisfied the criterion of R 2 ≥ 0. 2  270 degree movements correspond to directions to the right, away from the body, left, and into the body, respectively) for the regression fit were retained for further analyses. As suggested by Figure   1c, single trial BAOs were directed in one of two oppositely oriented directions (Figure 1d). For all three monkeys, the distributions of single-trial BAOs (pooled over 4 sessions for monkey Mk, 9 sessions for Ls, 11 sessions for Bx making 45 degree movement and 11 sessions for Bx making 225 degree movement) were bimodally distributed along a roughly rostro-caudal axis (rostromedial-caudolateral axis for Mk and Bx and rostrolateral-caudomedial axis for Ls) albeit the two modes were not always equal in size (Figure 1e). The BAO distributions were well fit by a mixture of von Mises functions model with means of 137° and -20° for Mk; 33° and -126° for Ls; 143° and -33° for Bx for 45 degree arm movements, and 137° and -32° for 225 degree arm movements. The BAO distributions for the two movement directions in Bx were not statistically different (p=0.1, Kuiper's two sample test). Multi-unit activity was not statistically different for trials in either of the modes and the few trials that fell outside of the two modes (one-way ANOVA for 3 groups, p-values ranging from 0.21 to 0.99 for different datasets; Supplementary Figure 3). However, trials in which the BAO propagated in the rostro-to-caudal direction had significantly shorter reaction times as compared to trials with caudal-to-rostral propagation in two of the monkeys (mean reaction time of 387 and 405 milliseconds, respectively, for Monkey Bx making 45 degree movements, p = 0.0049, two-sample t-test (not significant for 225 degree movements, p=0.4); 119 and 143 milliseconds for Monkey Ls, p = 0.0002, two-sample t-test). Black dots represent single channel values (offset by 10ms and jittered for clarity). c. Beta envelopes across channels from Monkey Ls for two trials (left) along with their corresponding spatial patterning (right) at different times (corresponding to vertical colored lines on left panel). Orientations of regression fits to the spatial patterning are denoted with arrows. d. Beta envelopes of two trials (left; inset shows spatial positions of three electrodes with early (blue), intermediate (green), and late (red) attenuation times) and their corresponding spatial maps of beta attenuation times (right) are shown for monkey Bx. Regression fits to the pattern show the beta attenuation orientation (BAO) directed rostrally or caudally (red arrows). e. Polar histograms show the distributions of BAOs across multiple single trials and recording sessions for three monkeys. A mixture of two von Mises functions was fit to the BAO distributions (red curves) to estimate the two direction modes of the BAO (blue arrows). Insets represent the initial hold (filled red) and target positions (filled gray) of the hand.

Propagating sequences and somatotopy
We examined the spatial relationship between the propagating sequences and the somatotopic organization of M1 using suprathreshold electrical stimulation. We documented muscle twitches through visual observation and muscle palpation in our three animals by stimulating each electrode (Supplementary Figure 4). Given the nature of the behavioral task (i.e. a reaching task involving primarily proximal muscles), we confined our analysis of propagating sequences on arrays placed primarily in proximal zone even though we had an additional array (in Ls and Bx) placed in the distal zone. Despite the crude somatotopy as others have documented [29][30][31][32] , we found that the axis of propagation was nearly orthogonal to the medio-lateral gradient from proximal to distal representations.

Causal evidence linking propagating patterns to movement initiation
To provide more direct causal evidence that propagating patterns facilitate movement initiation, we conducted experiments designed to perturb movement initiation in which subthreshold, patterned electrical stimulation across a set of electrodes was delivered during the reaction time period (Figure 2a). We hypothesized that patterned electrical stimulation that propagated against (INCONGRUENT) the natural propagating pattern of beta attenuation would delay movement initiation as compared to patterned stimulation that propagated with (CONGRUENT) the natural propagating pattern (see Methods). Given that the BAO varies from trial to trial (i.e. along the rostro-caudal or caudo-rostral directions), we could not assume a particular BAO on any given stimulation trial. Moreover, stimulation artifacts made it difficult to measure the BAO on stimulation trials. This was not an issue for monkey Ls as the bimodal distribution of BAOs was highly asymmetric on non-stimulation trials (i.e. nearly unimodal). Therefore, we could assume a single BAO direction across stimulation trials (2086 trials recorded over 13 sessions). For monkey Bx where the BAO distribution was nearly symmetric, we devised a method for selecting a subset of stimulation trials that we could assume had a highly unimodal BAO distribution based on no stimulation (NOSTIM) trials (see Methods). Only this subset of trials was analyzed for Bx (657/1772 or 37% of stimulation trials for 45 degree movement and 724/1686 or 43% of stimulation trials for 225 degree movement). Likewise, for monkey Mk, a similar procedure was applied to select a subset of trials (84/187 or 45% of stimulation trials) from one of the datasets that showed a bimodal distribution. The second data set from Mk exhibited a highly unimodal distribution on NOSTIM trials. It should be emphasized that the assumption of unimodality is conservative and likely leads to weaker differential effects on reaction time given that some stimulation trials may have had BAOs directed in the opposite direction.

Functional connectivity emerges along the rostro-caudal axis at movement initiation
If beta attenuation as measured by the LFP is truly a reflection of local network excitability and ultimately activation, we would expect a similar spatio-temporal pattern to emerge among unit activity at movement initiation. To characterize spatio-temporal patterning in neuronal spiking, we estimated functional connectivity among unit activity at different epochs during the reaching behavior in monkeys Bx and Ls (Figure 3a). Monkey Mk was excluded from this analysis due to an inadequate number of channels with spiking activity. We used a Granger-type statistical model applied to point processes 33 Kuiper's two sample test). In order to more effectively isolate the neural dynamics of movement planning from movement execution, we also examined functional connectivity during an instructed delay task for monkey Bx where the target stimulus was presented for a random duration ranging from 500 to 1200 milliseconds before a go cue to initiate the movement. Only trials with instructed delay durations longer than 1000 ms were used for the functional connectivity analysis. It is known that single units in primary motor cortex modulate during the instructed delay period where there is no movement 35,36 . We assessed the orientation of directed functional connections during the early

Propagating patterns along the rostro-caudal axis more reliably predict muscle activity
We hypothesized that downstream areas such as the spinal cord (and motor neuron pools) may be sensitive to the spatio-temporal recruitment order of cortical neurons. Beta oscillation dynamics are likely not involved in the detailed components of movement execution 25 . However, evidence has shown that cortical oscillations and muscle EMGs are related and exhibit significant coherence in the beta frequency range under certain behavioral conditions 37 . More importantly, this coherence can modulate continuously and is related to motor parameters 38 . One way to test the idea that propagating beta attenuation may facilitate muscle recruitment is to decode muscle activity from beta amplitude profiles in cortex and then apply different spatial perturbations to the beta profiles. Using cortical LFP and EMG data from Bx and Ls, we trained a feed-forward neural network (Figure 4a Reconstruction performance on test data not used to train the network was assessed by the fraction of variance accounted for (FVaF in percentage) and ranged from 18% to 55% on all five muscles with an average of 40.5% for Bx and ranged from 15% to 51% on all muscles with an average of 36.5% for Ls. Our goal here was not to develop an accurate decoding algorithm for brain-machine interface applications, but rather to adequately reconstruct muscle activity with a model that could then be used to apply desired spatial perturbations. We applied two kinds of spatial perturbations (random swaps) to the model's input: a swap 1) parallel to or 2) orthogonal to the BAO axis (See Methods; Figure 4c). A total of 5000 swaps were performed for each perturbation type from which a distribution of FVaFs were generated and compared to each other as well as to the unperturbed FVaF (Figure 4d and Supplementary Figure 6b). Only the parallel perturbation would disrupt the natural sequencing of beta attenuation propagation, and, therefore, we predicted that such perturbations would more effectively disrupt EMG output from the neural network model. Mean FVaFs were, in fact, significantly lower for the parallel perturbations as compared to orthogonal perturbations (p = 6.6e-200 for Bx and p = 3.6e-84 for Ls, two sample t-test).

Discussion
This work provides important evidence demonstrating that a specific, propagating rostrocaudal pattern of M1 excitability and unit modulation is required for voluntary movement initiation. First, using subthreshold patterned electrical stimulation, we showed that patterned stimulation that propagates against (INCONGRUENT) the natural propagating direction leads to longer reaction times whereas stimulation propagating with (CONGRUENT) the natural pattern did not affect reaction times. While any form of electrical stimulation will perturb the natural cortical dynamics associated with movement initiation, we suggest that the patterned stimulation propagating against the natural sequence provides a more potent perturbation. Subthreshold ICMS at single sites has been previously shown to lead to increased reaction time in dorsal premotor cortex followed by normal motor execution 41 . Interestingly, similar effects on reaction time were absent for similar single site, subthreshold ICMS in M1. We argue that single site stimulation may not be the ideal paradigm for perturbing the cortical dynamics in M1 associated with movement initiation. Second, functional network connectivity among an ensemble of spiking neurons emerges at movement initiation that is spatially anisotropic and oriented along a similar rostrocaudal axis. This was not observed during movement preparation or planning even though it is known that M1 neurons modulate during these periods (See Figure 3 and Supplementary Figure   5). Spatial structure in functional connectivity was often weak during movement preparation and planning, and, if anything, its spatial anisotropy was oriented orthogonal to the rostro-caudal axis.
Anatomical and physiological studies have suggested a spatial anisotropy of horizontal connectivity in M1 that might support these spatially organized propagating sequences. Retrograde axon degeneration following punctate lesions in M1 has indicated longer range horizontal connectivity along the rostro-caudal axis as compared to the medial-lateral axis 42 . Moreover, an ICMS study has provided evidence for dominant functional connectivity along the rostro-caudal axis 43 .
The bimodal distribution of propagation that we observed here is similar in orientation to the bimodal distribution of beta waves (as measured by the phase gradient of the beta oscillations) along the rostro-caudal axis we and others have previously documented to occur during movement preparation 33,[44][45][46][47] . This may suggest a mechanism by which on-going beta waves propagating along the rostro-caudal axis during movement preparation lead to a sequential pattern of neural activation prior to movement initiation. Beta waves alternately propagate either in the rostral-to-caudal direction or vice versa in a seemingly random fashion during movement preparation. It is also known that excitability is enhanced at particular phases of the beta oscillation 25 (Murthy and Fetz 1996b). Therefore, as a particular beta wave travels across the cortical sheet in one direction, a wave front of excitability moves across the motor cortex in the same direction. Assuming input originating from outside of motor cortex signaling movement initiation interacts with this beta wave, the combination will more likely lead to firing at that position in motor cortex. The traveling wave front may, therefore, lead to a sequential pattern of activation.
An open question is how these propagation patterns are mechanistically related to appropriate muscle activations for movement initiation. One possibility is that rostro-caudal propagation enhances M1 output through recurrent connectivity due to stronger rostro-caudal horizontal connectivity. However, M1 firing rates for trials when propagating excitability occurred along the rostro-caudal axis were not statistically different from those few trials when propagation occurred outside the rostro-caudal axis. Another possibility is that downstream areas such as the spinal cord are sensitive to the recruitment order of cortical sites. Our muscle decoding results indirectly supports this second alternative. Perturbing the spatial order of beta amplitude inputs along the rostro-caudal axis more effectively disrupts the natural propagating patterns and leads to significantly weaker EMG predictions as compared to orthogonal perturbations which preserve the overall propagation pattern.
The sensitivity of downstream receiving populations to these cortical propagating sequences could be due to the specific characteristics of these patterns. For example, the propagation axis can be thought of as a composition of two phenomena: simultaneous activation of populations orthogonal to the propagating axis (i.e. along the medio-lateral axis) and sequential activation of populations along the propagating axis. Either or both of these phenomena could be responsible for downstream activations that lead to movement initiation. One possibility is that synchronous activation along the medio-lateral axis and rostro-caudal cortical sequencing lead to non-linear amplification of muscle activity which support movement initiation and generation.
Further experiments will be needed to shed light on the exact nature of the downstream sensitivity to these natural spatiotemporal patterns.
Our results suggest that these propagating sequences serve a generic role in initiating movements and do not depend on the nature of the movement or the temporal sequencing of limb segments. Although we have not examined this question in great depth, our findings from monkey Bx demonstate that the bimodal distributions of BAOs and functional connectivity are not statistically different between movements toward and away from the animal. Moreover, the spatial structure of these propagating patterns are nearly orthogonal to the medio-lateral gradient from proximal to distal representations. This suggests that these propagating sequences are likely not the result of proximal-to-distal sequencing that is observed in certain motor behaviors 48,49 .
These results have broad significance because they suggest that large-scale, spatially structured propagating patterns of activity may be a common strategy used in various cortical areas to perform behaviorally relevant computations. Spatio-temporal patterns such as traveling planar, radial, circular, and spiral waves have been documented across different cortical areas including somatosensory 50,51 , visual 52-57 , auditory 58,59 cortices, and hippocampus. Various functional roles for these spatio-temporal patterns have been put forth but it remains a possibility that they are epiphenomena of the recurrently connected cortex 60 . Multi-site, spatio-temporal stimulation using electrical or optogenetic stimulation may be an important approach to provide more direct support for a role of these patterns in cortical processing. Neural data recorded as analog signals were amplified with a gain of 5000, bandpass filtered between 0.3Hz and 7.5 kHz, and digitized at 30 kHz. For monkey Mk, the local field potentials were sampled at 2 kHz from the digitized data after low-pass filtering the signals with a 4 th order Butterworth filter and a cut-off frequency of 500 Hz. For Ls and Bx, the neural data were zero-phase low-pass filtered with a 10 th order Butterworth filter and a cut-off frequency of 500 Hz, and down-sampled to 2 kHz. Multi-unit activity (MUA) was extracted as threshold crossings from the raw data after high-pass filtering with a cut-off frequency of 250 Hz. The threshold limit was set to be -5.5 times of the signal RMS value. For EMG recordings, signal from the thirteen muscles were amplified individually and bandpass filtered between 0.3 and 1 kHz prior to digitization and sampled at 10 kHz.

Behavioral Task
The animals were operantly trained to perform a center-out planar reaching task using a two-link exoskeletal robot (BKIN Technologies). The tool-tip of the robot determined the position of a cursor on a screen, and the monkeys were able to freely move the robot in the horizontal plane using their upper limb. Monkey Mk used a horizontal screen with his arm and robot underneath the screen, while monkeys Ls and Bx had a vertical screen facing them such that the horizontal movements forward and back corresponded to upward and downward cursor movements. The animals had to hold the cursor on a center target for a fixed duration (600 ms for Mk, 900 ms for Ls, and 1000 ms for Bx) to trigger a peripheral target appearance (GO cue). The animal had to move the cursor to reach the peripheral target to get a juice reward. Monkey Mk was presented

Beta Band and Single-trial Beta Attenuation Orientation
The frequency associated with the peak power in the beta frequency range (15-35 Hz) was estimated from the channel-averaged power spectral density of the neural signals from each array over a large number of trials. For monkey Ls, the beta peak occurred at 21 Hz (approximated to the nearest integer value), and for monkeys Bx and Mk the peak occurred at 31 Hz. For the latter two monkeys, a second smaller peak also occurred at 21 Hz but was not used for further analysis.
The signal was then band-pass filtered with a bandwidth of 6 Hz centered on the peak beta, i.e., 21±3 Hz for monkey Ls and 31±3 Hz for monkeys Bx and Mk. The beta envelopes of the signal were computed using the magnitude of the Hilbert transformation applied to the band-pass filtered signal (see Figure 1a).
Single-trial beta envelopes were extracted using an auto-encoder neural network where the inputs are the same as the targets. Segments of the signal (-500 ms to 600 ms with respect to the GO cue) from all the recorded channels were extracted for each trial and decimated to 100 Hz. A vector of length nChannels × mSamples was created for each trial and fed through the auto-encoder (see Supplementary Figure 2b). The hidden layer was set to be half the size of the input layer so as to form a bottleneck architecture, and the data were projected to the output layer. With training, the auto-encoder was able to reconstruct the signal that preserved both spatial patterns and significant temporal variance of the signal. The output beta envelopes from the auto-encoder were then low-pass filtered with a 10 th order filter and cut-off frequency of 5 Hz. The filtered signals were normalized to have an amplitude range of [0, 1].
The normalized single-trial beta envelope attenuates and crosses an amplitude threshold value at a particular time (referred to as beta attenuation time or BAT) which varies depending on the channel (see Figure 1d). These BATs, when arranged spatially on their corresponding electrode locations on the array, exhibit a propagating pattern with a particular orientation. The orientation of the beta attenuation (BAO) was estimated using a linear regression fit to the BATs as follows, where, and are the coefficients of the rows and columns of the array, and denotes the constant offset time. The direction of BAO was then estimated as, The F-statistic of the model was used to ascertain statistical significance of the fit. For each trial, a BAO and its associated coefficient of determination (R 2 ) were calculated. On any given trial, the amplitude threshold for determining the BATs was determined by varying the amplitude from 0.5 to 0.15 in steps of 0.05 and finding the threshold with the largest associated R 2 .  Figure 2a). Pulse trains on a successive column were initiated at a latency of 10 milliseconds relative to the previous column. For Mk, 3 or 6 electrodes were stimulated at any instance, with inter-stimulation latency of 4 milliseconds (see Pattern B in Figure 2a).

Determining BAO for Trials with Stimulation
For trials that received patterned stimulation, recovering the complete beta envelopes was not possible due to large saturating electrical artifacts during stimulation. For monkey Ls the bimodal BAO distribution was highly asymmetric on non-stimulation trials (i.e. nearly unimodal).
Also, for one data session examined in Mk, the BAO distribution was highly unimodal. Therefore, we could assume a single BAO direction across stimulation trials. To determine the BAOs for the stimulation trials in monkey Bx and for the second session in Mk, an alternative approach was taken by examining the dynamics of the beta envelopes in trials when no stimulation occurred.
Typically, the channel-averaged beta envelope exhibits a local peak in amplitude at different times prior to attenuation depending on the trial. The approach was to find a time window around the GO cue where the local beta peaks occurred on a subset of no stimulation trials such that their BAO distribution was highly unimodal. This time window (set to 300 ms in duration; 150 ms for monkey Mk) was found by sliding the window in 100 ms steps starting at -300 ms with respect to the GO cue, selecting no stimulation trials whose local beta peak occurred in that window, and computing the BAO distribution (Supplementary Figure 7). The window whose selected trials were highly unimodal (where at least two thirds of the trials in the window were oriented along one BAO mode) in their BAO was chosen and then used to select stimulation trials whose local beta peak occurred in the same window. The assumption was that these stimulation trials exhibited a similar unimodal BAO distribution and, therefore, all these stimulation trials had BAOs pointing roughly in the same direction. To detect local beta peaks in stimulation trials, partial beta envelopes (-500 ms to 200 ms with respect to the GO cue; -500 to 25 ms for monkey Mk) that were not corrupted by electrical artifacts were extracted using the auto-encoder. Whenever multiple time windows exhibited unimodal BAOs in the no stimulation trials, the window that was closer to or overlapping the GO cue was chosen. To compute reaction time means and distributions for all three conditions, we excluded trials whose movement onset times occurred before the onset of the stimulation as well as trials with very long reaction times (>500 ms for Mk and Bx and >300 ms for Ls).

Computing Reaction Times of Single Trials
Tangential velocity trajectories estimated from the Cartesian position of the reaching kinematics data sampled at 2 kHz were used for reaction time computation. From the normalized trajectories of each trial, movement onset times were computed when the velocity crossed 10% of the peak velocity. When there were smaller velocity peaks that exceeded the 10% threshold prior to final velocity peak associated with the actual movement to the target, movement onset was computed using the final threshold crossing.

Firing rates
Multi-unit firing rates were calculated for each trial in 150 ms temporal windows, shifted over time with steps of 10 ms. Rates were then averaged over trials to obtain the mean firing rate of each channel's multi-unit. For the heatmaps of Figure 3, mean firing rates were normalized by dividing by the multi unit's maximum rate. A different normalization procedure was used to compare the firing rates of each mode. First the average baseline activity over all trials (500ms prior to the go cue) was subtracted from each trial's firing rate. Then, trials were separated into three groups according to BAO mode (rostro-to-caudal BAO, caudo-to-rostral BAO, out of modes) and the average firing rate of each group was calculated. If the multi-unit displayed suppressed modulation, its rates were converted to the additive inverse. Finally, the three average rates were divided by the global maximum, in order to retain possible mode differences. Population response for the three modes was then estimated by averaging over multi-units. To compare the population responses between modes, we used the normalized unit responses at two different epochs: reaction (from go cue to movement onset) and movement (from movement onset to movement end). A oneway AVOVA was done separately at each epoch [factor: BAO mode (three levels), sample size: number of multi units] to determine differences between the population firing rate across BAO modes.

Directed Functional Connectivity
A granger-type causality estimation technique generalized for point-processes was used to infer directed connectivity among simultaneously recorded neurons. For each candidate multi-unit, the instantaneous spiking was modeled as its conditional intensity function (CIF), λ(t|H(t)), where H(t), is the spiking history of all the multi-units in the ensemble up to time t. Using a generalized linear model, the log CIF of each multi-unit was modeled in terms of the spiking of its covariate (i.e. source) multi-units, initially with a history of 60 milliseconds binned at 3 millisecond resolution for each source multi-unit. The final model was then found by optimizing the length of the history term using the AIC (Akaike information criterion). To determine the influence of a source multi-unit on a candidate multi-unit, the log-likelihood ratio between models (i) including and (ii) excluding the source multi-unit was computed, and when the model performance decreased significantly (p ≤ 0.05, χ 2 -test, false detection rate corrected for multiple comparison) with the source multi-unit excluded, a directed connectivity from the source multi-unit to the candidate neuron was inferred.

Decoding Muscle Activity from Beta Envelopes
A feed-forward neural network was developed to map instantaneous beta amplitudes from multiple channels to the EMG activity of a group of muscles. Normalized beta envelopes from 60 channels (100 samples/second) were projected first through the auto-encoder and then fed into the feed-forward neural network to decode 5 muscles based on their relevance to the task and the quality of the signal (Anterior Deltoid, Pectoralis Major, Biceps Lateral, Triceps Long head, Flexor Digitorum Superficialis). The inputs were vectors of 60 × 1 elements, and the outputs were vectors of 5 × 1 elements. The neural network architecture (see Figure 4a) used 120 hidden neurons with logarithmic sigmoid activation function. The output layer that predicts the EMG activities used a linear transfer function. Of the available trials, 80% were used for training and 20 % for testing.
The prediction accuracy was computed for each individual muscle as the fraction of variance accounted for (FVaF) between the actual EMG and the EMG predicted from the neural network as, The input to the network was perturbed by swapping the spatial locations of the channels in two different ways. For each swap, beta profiles on eight electrodes were randomly selected on the 64 electrode array and then randomly swapped spatial locations with other beta profiles with the constraint that swapping occurred either parallel to or orthogonal to the BAO axis. A total of 5000 swaps were performed for each of the directions (parallel or orthogonal) and the FVaF for each swap was estimated.

Somatotopic Maps
In order to determine the somatotopy of M1, supra-threshold electrical stimuli ( system (Blackrock Microsystems, Inc. Salt Lake City, UT). Certain muscles were also palpated to identify any evoked responses. We used a 7-point scale (color-coded in Supplementary Figure 4) to label each site from (fingers, red), to wrist (yellow), to elbow (baby blue), to shoulder (blue) joint and/or muscle twitches. Figure 4. Relationship between propagating sequence directions and somatotopy. The somatotopic maps across the arrays implanted in the three monkeys show the proximal-to-distal somatotopy oriented along the medio-lateral axes. The somatotopic maps for Ls and Bx are based on two 8x8 arrays implanted medio-laterally in M1. The BAO axes (blue arrows), computed for the medial arrays in Ls and Bx, and the one array in Mk, propagate nearly orthogonal to the proximal-to-distal gradient.

Supplementary
Supplementary Figure 6. Effects of spatial perturbations on EMG predictions in monkey Ls. a. Trial-average predicted EMG (red) and actual EMG (blue) and profiles using a neural network trained with beta envelopes from Monkey Ls. b. Distributions of FVaFs over 5000 swaps are shown for the parallel swap condition (yellow; mean FVaF of 5%) and the orthogonal swap condition (blue; mean FVaF of 10%) along with the FVaF of the original unperturbed model (red line; mean FVaF of 36.5% (varying from 15 % to 51% across muscles). Inset represent the initial hold (filled red) and target positions (filled gray) of the hand.