COMMUNICATIONS

ARTICLE

Received 21 Oct 2014 | Accepted 14 Apr 2015 | Published 21 May 2015^^BDOI:iai038/ncomms8169| OPEN

Large-scale spatiotemporal spike patterning consistent with wave propagation in motor cortex

Kazutaka Takahashi1'*, Sanggyun Kim3'*, Todd P. Coleman3, Kevin A. Brown4, Aaron J. Suminski5, Matthew D. Best2 & Nicholas G. Hatsopoulos12

Aggregate signals in cortex are known to be spatiotemporally organized as propagating waves across the cortical surface, but it remains unclear whether the same is true for spiking activity in individual neurons. Furthermore, the functional interactions between cortical neurons are well documented but their spatial arrangement on the cortical surface has been largely ignored. Here we use a functional network analysis to demonstrate that a subset of motor cortical neurons in non-human primates spatially coordinate their spiking activity in a manner that closely matches wave propagation measured in the beta oscillatory band of the local field potential. We also demonstrate that sequential spiking of pairs of neuron contains task-relevant information that peaks when the neurons are spatially oriented along the wave axis. We hypothesize that the spatial anisotropy of spike patterning may reflect the underlying organization of motor cortex and may be a general property shared by other cortical areas.

1 Department of Organismal Biology and Anatomy, University of Chicago, 1025 E 57th Street Culver Room 206, Chicago, Illinois 60637, USA. 2 Committee on

Computational Neuroscience, University of Chicago, 1025 E 57th Street Culver Room 206, Chicago, Illinois 60637, USA. 3 Department of Bioengineering,

University of California, San Diego, 9500 Gilman Drive MC 0412, San Diego, La Jolla, California, 92093, USA. 4 Center for N euroscience, New York University,

4 Washington Place, New York City, New York 10003, USA. 5 Department Electrical Engineering and Computer Science, Milwaukee School of Engineering,

1025 N Broadway, Milwaukee, Wisconsin 53202, USA. * These authors contributed equally to this work. Correspondence and requests for materials should

be addressed to N.G.H. (email: nicho@uchicago.edu).

Propagating waves of neural activity are ubiquitous and have been documented at different spatial resolutions in a number of different neocortical areas including visual1-6, somatosensory5,7-9, auditory10,11 and motor cortices12-14 as measured via multielectrode local field potential (LFP) recordings, voltage-sensitive dyes (VSDs) and multiunit activities. Oscillatory LFPs and electroencephalograms in the beta frequency range (15-40 Hz) are ubiquitous in the motor cortex of mammals including monkeys15-18 and humans19,20. In particular, we have previously demonstrated that across the precentral gyrus of the upper-limb area of primary motor cortex (MI), these oscillations are not perfectly synchronized but rather exhibit phase gradients that indicate planar propagating waves along what we define as a beta wave axis, a rostro-caudal axis in monkeys13 and a medio-lateral axis in humans14 at a range of propagating speeds that were consistent across subjects.

However, as both LFPs and VSD measure aggregate potentials from groups of neurons near the recording site, it has never been shown whether action potentials from individual neurons demonstrate spatiotemporal patterning consistent with wave propagation. This is important because it is still debated as to what aggregate signals such as LFPs and VSD signify physiologically, whereas single-unit action potentials are understood to mediate interneuronal communication. Moreover, the functional significance of this wave propagation for motor control is unclear (but see recent computational studies21,22). Here we first show that MI neurons can be classified, based on the spike waveform widths, into two groups of neurons exhibiting distinct spectral properties. We then estimate effective connectivity of networks of spiking neurons based on this classification using a Granger causality analysis applied to point processes, and demonstrate that a class of simultaneously recorded, single-motor cortical neurons with narrow spike waveforms in non-human primates spatially coordinates their spiking activity in a manner that closely matches the orientation of prominent beta wave propagation. We also demonstrate that sequential spiking activity of that class of neuron pairs contains task-relevant, target-direction information whose magnitude varies according to the spatial orientation of the constituent neurons in a manner consistent with the beta wave axis.

Results

Beta waves in the motor cortex. We recorded multiple single-unit and LFP activity from MI using chronically implanted high-density microelectrode arrays while three rhesus monkeys (Rs, Mk and Rj) made planar reaching movements using a two-link robotic exoskeleton (BKIN Technologies, ON, Canada). The monkeys performed a random target-pursuit (RTP) task23 that required them to move a cursor (aligned with the position of their hand) through a sequence of randomly positioned targets. Movement durations from target to target ranged from 300 to 450 ms with mean speeds (± s.d.) of 22.33 ±11.17 (Rs), 14.12 ± 6.27 (Mk) and 6.11±7.29cm s"1 (Rj). Planar beta wave activity measured from spatially distributed LFP sites was evident at particular intervals of time throughout the performance of this task (Fig. 1a). We used a method described previously13 to characterize the properties of planar beta waves. We found that the degree of planar wave propagation as measured by a quantity called phase gradient directionality (PGD) was strongest ~ 100-150 ms after the target onset (Fig. 1b) when beta power was high (Fig. 1c), and when visual target information reached the motor cortex13 followed by movement initiation to the new target (see wrist speed in Fig. 1b). Consistent with our previous findings using a center-out task13, wave propagation directions during the RTP task exhibited either a

bimodal distribution (monkey Rs) or unimodal distribution with a small secondary mode (monkeys Mk and Rj), with one mode oriented in the rostral-to-caudal direction and a secondary mode oriented in the opposite direction (Fig. 1d). We denoted the caudal wave and rostral wave directions defined by the mean direction of the first or only mode of the wave propagation distribution and the opposite direction oriented roughly along the rostro-caudal axis. The distribution of propagation speeds was always unimodal with means and medians ranging from 23.2 to 26.7 and from 10.1 to 13.5cms_ 1, respectively (Fig. 1e).

To investigate spatiotemporal patterns of spiking activity (see Supplementary Figs. 1 and 2b,d,f for many examples of Peristimulus time histograms (PSTHs) of neurons used in the analysis), we first partitioned the population of recorded neurons into two potentially distinct, functional classes (Fig. 2a) based on the observation that their spike waveform widths exhibited a bimodal distribution. For each data set separately as well as for pooled data, a two-Gaussian mixture model (that is, bimodal) was found to be the most parsimonious fit to the distributions based on the Bayesian Information Criterion (BIC; Fig. 2b). Narrow-spiking neurons tended to exhibit higher mean spike rates as compared with broad spiking neurons for all three monkeys (however, only Mk's data sets showed a statistically significant difference in means, unpaired two-sample t-test with unequal variances): 4.89 ±4.10 spikes s "1 versus 2.96 ± 1.94 (P = 0.0662) for Rs; 11.94 ±11.48 spikes s "1 versus 4.22 ± 4.48 (p = 0.0419) for Mk; and 11.68 ±9.00 spikes s"1 versus 5.59 ±4.45 (P = 0.0907) for Rj. A few examples of PSTHs of neurons for monkey Mk from both narrow and wide classes used are shown in Fig. 2c (more in Supplementary Fig. 1 and for other monkeys in Supplementary Fig. 2). By summing spikes over all neurons within a class, the population spike rates of the two classes for monkey Mk were different particularly around the time that LFP beta power increased (Fig. 2d, and for other monkeys in Supplementary Fig. 2). Temporally, the spectrogram of the narrow-spiking population showed sustained power over the beta range, and the power increased when the LFP beta power increased (Fig. 2e). In contrast, the wide population exhibited a systematic power decrease in the beta and higher frequency bands around the onset of LFP beta power increase (Fig. 2f). Moreover, the spectrum of the population activity associated with narrow spike waveforms exhibited a spectral peak that was very close to the LFP beta frequency peak, whereas the population spectrum associated with wide spike waveforms exhibited no prominent beta peak (right side of Fig. 2e,f).

A number of studies have classified cortical neurons based on their extracellular spike widths, and some have suggested that narrow-spiking neurons represent local inhibitory neurons23-28. However, the range of spike widths that defined our narrow class is larger than that of many of these studies. Given that local inhibitory neurons comprise only ~20% of all neurons in neocortex29,30, the large proportion of spiking neurons that we classified as narrow suggests that this class includes excitatory pyramidal neurons as well as local inhibitory neurons that, as a whole, have a tendency to synchronously oscillate in the beta frequency range.

Directed effective connectivity. We developed a generative model to characterize the spiking responses of a given neuron based on the past spiking of all other simultaneously recorded neurons in the network. This model is a generalization of Granger causality to point processes, and can be used to infer excitatory- and inhibitory-directed effective connections among multiple neurons, unlike popular pairwise network estimation procedures31. Specifically, we estimated the conditional intensity

3 112 ms 114 ms 116 ms 118 ms 120 ms 122 ms

L-|- M R

L^— M R

L—|— M R

b.- - -I J1-ri . - I

124 ms 126 ms 128 ms 130 ms 132 ms 134 ms

£ 20 10

-200 0 150 350 550 Time (ms)

30 25 20 15 10 5 0

-200 0 150 350 550 Time (ms)

Rs Mk R

10 20 30 40 50 60 70 Speed (cm s-1)

Figure 1 | Properties of LFP beta waves. (a) Temporal snapshots of the LFP voltages across the array indicating wave propagation. The LFP voltage on each electrode was band-pass filtered in the beta frequency range (that is, ±3 Hz centred at the beta peak of the power spectrum). Time in milliseconds labelled above each plot is with respect to the onset of the visual target. The red arrow in the bottom right panel shows the propagation direction of the wave. (b) Temporal evolution of the GOF measure of planar wave activity (PGD in blue ranging from 0.295 to 0.335, as well as the mean hand speed (green) ranging from 5 to 36 cm s "1. (c) Averaged spectrogram of a single-channel LFP revealing the temporal dynamics in beta frequency power relative to the visual target onset. (d) Circular distributions of wave propagation directions for monkeys Rs (cyan), Mk (magenta) and Rj(brown). A solid black line denotes the rostro-caudal axis on the cortical surface. A dashed line in each rose plot connecting Cw (caudal wave direction) and Rw (rostral wave direction) denotes the axis defined by the first or only mode of beta wave propagation axis. Each panel below the circular distributions depicts the location of the multielectrode arrays (4 x 4 mm) in the arm area of MI for the corresponding subject. A red horizontal bar in the right lower corner in each panel is 4 mm. Landmarks and orientations: Cs, central sulcus; As, arcuate sulcus; C, caudal; R, rostral; M, medial; L, lateral. (e) Distributions of estimated propagation speeds for monkeys Rs (cyan), Mk (magenta) and Rj(brown).

function (CIF) of each neuron (a receiving neuron), 1(t|H(t)), given the spiking history, H(t), of all other neurons (sending neurons) in the network using a generalized linear model (GLM). The spike history from each sending neuron was composed of five 3 ms time windows for a total of 15 ms. The history length was selected based on the time it would take for the beta waves to propagate across the array at the mean estimated beta wave propagation speeds (Fig. 1e). The underlying assumptions of the model are twofold: (1) the present spike count of a receiving neuron conditioned on the past spike counts of the sending neurons is Poisson distributed, and (2) the log spike count of the receiving neuron is a linear function of the spike counts of the sending neurons. To check the sufficiency of the models used in the analysis, we performed a model-order selection analysis using the BIC, and found that the majority of the neurons showed an optimal spike history duration of <15ms (>96%; ref. 32). To assess model goodness-of-fit (GOF) for each neuron, we used the time-rescaling theorem33 and Kolmogorov-Smirnov (KS) plots to define the GOF area (GOFA) and GOFA ratio (GOFR; Methods and Supplementary Fig. 3a-c). More than 85% of neurons exhibited GOFA to be <0.05 (Supplementary Fig. 3d) and >68% of neurons showed GOFR to be <0.2 (Supplementary Fig. 3e). Thus the majority of the neurons used in the analysis were reasonably well characterized as Poisson GLM.

As in the standard Granger causality protocol, an effective connection from the sending neuron to the receiving one was inferred by comparing the log-likelihoods of the data under

models with and without a particular sending neuron's spiking history. If the performance of the model decreased significantly (P<0.05, w2-test) when a sending neuron was removed, an effective connection was inferred, and the sign of the sum of the five history terms determined whether the connection was excitatory or inhibitory34. We estimated the complete network of significant directed connections between recorded neurons with narrow spike waveforms and mean spike rates of at least 1 Hz (to ensure numerical stability) in different time windows. The temporal evolution of the networks in 150 ms windows incremented by 50 ms steps is illustrated in Fig. 3a for monkey Rs with 38 narrow neurons, Fig. 4a for monkey Mk with 21 narrow neurons and Supplementary Fig. 4a for monkey Rj with 33 narrow neurons. The largest number of connections was observed between 100 and 350 ms after the appearance of a new visual target at a time slightly after LFP beta power and the PGD reached their peak values. The temporal evolution of network topology estimates, based on the number of excitatory, inhibitory and all connections, was highly reliable as determined by calculating networks on non-overlapping subsets of data within each data set (Fig. 5 and Supplementary Fig. 5). Furthermore, in order to check whether our connectivity results depended on the precise spike timing among neurons or were a consequence of rate modulation (neurons with very high spike rates are more likely to have effective connections), we randomly shuffled the trial order for each neuron and within each of eight target directions, and recomputed network interactions. We found that

0.47 ms

0.23 ms

100 '«

0.2 0.3 0.4 0.5 Spike width (ms)

20 10 0 -10 -20

23 Number of gaussians

«J 10

Narrow neuron #1

Narrow neuron #2

o- 20 roiS 15 CE

Wide neuron #1

0 400 Time (ms)

Wide neuron #2

0 400 Time (ms)

0 200 Time (ms)

N 70 ^ 60 ^ 50

uen 40

Z 30 Ül 20

Spike power (dB) -4 -2 0 T1

70 60 50 40 30 20

0 200 Time (ms)

-24-12 0 LFP power (dB)

Spike power (dB) -4 -2 0

iZ 70 ~ 60 o 50

uen 40

£ 30 Fr 20

70 60 50 40 30 20

0 200 Time (ms)

-24-12 0 LFP power (dB)

Figure 2 | Two functional cell classes based on extracellular spike waveform width. (a) Histograms of spike waveform widths for all three monkeys, Rs, Mk and Rj. The solid curve indicates the fit using mixture of two Gaussians. Background colour indicates the bins classified as narrow (green) and wide (yellow) units. Inset: Examples of motor cortical units with narrow (green) and wide (yellow) spike waveform widths. (b) The difference between BIC values minus the mean BIC values over the range as a function of the number of Gaussians used in the mixture model to fit the distributions of spike widths. (c) Example peri-stimulus time histograms from two narrow (green) and two wide (yellow) spiking neurons from monkey Mk. (d) Averaged population spike rates for the two classes of cells from monkey Mk, narrow class in green and wide class in yellow. The time-resolved beta oscillation amplitude is plotted as a solid black line. (e,f) Averaged spectrograms for the two classes of units from monkey Mk, narrow class (e, left) and wide class (f, left), and averaged power spectra for the each class of units, narrow class in green (e, right) and wide class in yellow (f, right) computed over [ -100, 300] ms relative to the visual target onset. Black dashed line on the right subfigures denotes average LFP spectrum over all channels with the base line pink noise removed.

the number of connections were significantly lower as compared with the unshuffled networks (Supplementary Fig. 6 for monkeys Rs and Mk), suggesting that effective connections reflect single-trial coordination across neurons, rather than overall rate modulations or tuning properties.

To place the estimated networks in a spatial context, we computed the orientation of each connection on the cortical surface. For each cortical orientation, we summed up each effective connection weighted by its strength (as determined by the log-likelihood ratio, see Methods section), and normalized the sum of the weighted connections by the total number of possible connections for that orientation. This normalization was performed to account for anisotropies in the spatial arrangement of the recorded neurons. Between — 50 and 100 ms relative to visual target onset, few connections were observed and very little spatial structure was evident. The number of connections increased in subsequent epochs, and peaked between 50 and 250 ms after the onset of a new visual target. The clear bimodality in the circular distribution of the excitatory connection network emerged between 150 and 250 ms after the visual target onset (Fig. 3c for monkey Rs, Fig. 4c for monkey Mk and Supplementary Fig. 4c for monkey Rj, and Supplementary Fig. 7).

Furthermore, this bimodality remained even when we used one spike-width threshold value, the weighted average of the means of the two Gaussians fit on all spike widths (Supplementary Fig. 8).

By comparing the network topology during early and late time windows (Fig. 6a, left and right respectively), we observed the emergence of a larger number of significant connections and a strong anisotropy of directed excitatory connections. The circular distribution of excitatory connection directions became bimodal and closely oriented to the beta wave axis during the later time window in all three monkeys. By fitting the beta wave directional distributions with a mixture of two von Mises distributions, we determined the mode and the circular s.d. (CSD). We found that the modes of excitatory connections fell within 1.1 times the CSD of the beta wave modes. In particular, the beta wave primary/ secondary modes measured with respect to the anatomical medial direction ( ± CSD) versus excitatory connection modes were 123/ 279 ± 34/52 degrees versus 130/290 degrees (Rs), beta wave modes of 120/258 ± 18/60 degrees versus 140/50 degrees (Mk) and beta wave modes of 297/104 ± 72/51 degrees versus 270/90 degrees (Rj; Fig. 6b, left). The circular distributions of inhibitory connection directions were also bimodal but oriented almost orthogonal to the beta wave axis (Fig. 6b, right). In contrast to the

L—|—M

-50 0 50 100 150 200

Time (ms)

250 300 350

1 2 3 4 5 6

1 2 3 4 5 6 Time window

1 2 3 4 5 6

Figure 3 | Spatiotemporal patterns of network connectivity of narrow class of neurons in MI in six 150 ms time windows incremented by 50 ms from monkey Rs data set. (a) Networks of significant directed connections on the array at different time windows (in milliseconds) in relation to the onset of visual target appearance. Each horizontal colour bar indicates duration of a time window used to compute a network surrounded by the same colour. C, R, M and L indicate caudal, rostral, medial and lateral orientations, respectively, on the cortical surface. Red and blue arrows represent excitatory and inhibitory connections, respectively. The black dots represent the relative positions of the electrodes on the array where the neurons were detected. Neurons whose spikes rates are >1 spikes"1 and with narrow spike widths (<0.267 ms) were analysed. The green scale bar on the right-bottom equals 400 mm on the cortical surface. (b) Number of significant connections at different time windows for both excitatory and inhibitory connections (left), excitatory connections only (center) and inhibitory connections only (right). The lines in different colour shades correspond to different time windows of data in the recording session. (c) Circular distribution of directed excitatory connections weighted by their strength and normalized by the total number of possible connections in each orientation on the surface of MI at the different time windows (top). Circular distribution of directed inhibitory connections weighted by their strength and normalized by the total number of possible connections in each orientation on the surface of MI at the different time windows (bottom). All rose plots are oriented in the same way as in the anatomical coordinate system defined with C, caudal; R, rostral; M, medial and L, lateral in a. A dashed line in the top left rose plot connecting Cw (caudal wave direction) and Rw (rostral wave direction) defines the wave propagation axis as in Fig. 1d.

narrow-spiking neurons, there were few, if any, statistically significant connections among wide-spiking neurons in each time window (Fig. 6c).

For excitatory connections between narrow-spiking neurons, we took the mean of the five time delays weighted by their model coefficient values, and computed the average time delay between each excitatory connection. Dividing the distance of each excitatory cell pair projected along the beta wave axis by the average time delay provided an estimate of the 'conduction' velocity of each excitatory connection along the beta wave axis. Estimated velocities ranged from 10 to 30cms_ 1 across all six time windows and four data sets, in rough agreement with our estimated wave propagation speeds.

Spatial anisotropy and target information. In order to investigate any behavioural relevance to the spatial anisotropy of excitatory connections among narrow neurons, we next investigated the information content of the sequential spiking activities between neuron pairs with excitatory connections. Sequential

spiking activity at time t was defined as the product of the receiving neuron's activity at time t and the integrated activity of the sending neuron's activity in the 15-ms time interval preceding time t. We found systematic differences in sequential spiking activity as a function of target direction. To assess this quantitatively, we calculated the mutual information between target direction and the sequential spiking of narrow neuron pairs that possessed excitatory effective connections as a function of time (Fig. 7a). For each of 18 cortical orientation bins, the information was summed across all cell pairs with excitatory connections and divided by the total number of possible connections in that orientation. It should be noted that we added information across cell pairs in order to capture relative information trends in time and across different cortical orientations. Absolute information values cannot be inferred because addition assumes a lack of redundancy across cell pairs, an assumption that may not be warranted. Spike sequencing provided peak target-direction information ~ 150-250 ms for monkey Rs, 100-150 ms for monkey Mk and 150-300 ms for monkey Rj. The circular distribution of mutual information exhibited spatial anisotropy

L-|-M R

-50 0 50 100 150 200 250 300 350

Time (ms)

u 80 s

1 23456 123456 123456

Time window

Cw 60 \ y\

\ * \ 1 k zr 1 X

\ Rw 4 l—l

30 f 4 i # # #

Figure 4 | Same as Fig. 3, but for monkey Mk data set.

All data

b Subdata #1

Subdata #2

Subdata #3

180 0 \

Figure 5 | Consistency in spatial patterning of network connectivity across different subsets of data in time window 5 within the data set in monkey Rs. Networks of significant connections (top), circular distributions of excitatory connections (middle) and circular distributions of inhibitory connections (bottom). (a) Results from all three subsets of data were pooled. (b) Results from each subset of data are illustrated in each column. All plots are oriented by the coordinate system defined by C, caudal; R, rostral; M, medial and L, lateral in a. A dashed line in each rose plot connecting Cw (caudal wave direction) and Rw (rostral wave direction) defines the wave propagation axis as in Fig. 1d.

whose peaks were oriented along the beta wave axis, 130 degree for monkey Rs, 140 degree for monkey Mk and 90 degree for monkey Rj with respect to the medial direction, and they were all within the beta wave primary mode direction ± CSD (Fig. 7b). Cross-correlation coefficient peaks (within ± 1 orientation bin) between the mutual information and beta wave distributions were 0.520, 0.353 and 0.252 for monkey Rs, Mk and Rj, respectively. Spike sequencing of narrow-spiking neurons that were not effectively connected showed an order of magnitude of smaller amount of total mutual information compared with those pairs with excitatory effective connections, and the spatial pattern of information did not resemble the distribution of beta wave directions (Fig. 7c). To verify that this spatial pattern of information was specific to the sequential firing of neuron pairs, we computed target-direction information from single-unit spiking activity (Supplementary Fig. 9a,b). The circular distribution of orientations of neuronal pairs weighted by the mean singleneuron mutual information values averaged over the two neurons with an excitatory connection between them revealed a weaker alignment to the beta wave axis (Supplementary Fig. 9c) than that of mutual information of sequential spiking of neurons with excitatory connections. Cross-correlation peaks (within ± 1 orientation bin) between mean single-neuron information and beta wave distributions were — 0.165, — 0.135 and 0.102 for monkey Rs, Mk and Rj, respectively (Supplementary Fig. 9d,e).

Discussion

Propagating wave activity is ubiquitous throughout the cortex but has never been characterized at the single-neuron level. Although our results do not directly demonstrate that single units propagate their spiking activity in a wave-like manner across the motor

a Window=[-50, 100] ms [150, 300] Rs

Figure 6 | Spatiotemporal patterns of network connectivity of narrow and wide classes of neurons in two time windows. Data from monkeys Rs (top), Mk (middle) and Rj (bottom). (a) Networks of significant directed connections on the array at different time windows. Red and blue arrows represent excitatory and inhibitory connections, respectively. The black dots represent the relative positions of the electrodes on the array where the neurons were detected. If >1 neuron was detected on a given electrode, a given dot contains all units on that electrode. Neurons with firing rates >1 spikes—1 and narrow spike widths (< ~0.2667 ms) in MI were analysed. The green scale bar on the right-bottom equals 400 mm on the cortical surface. R, rostral; C, caudal; M, medial and L, lateral denote the anatomical orientation. (b) Circular distributions of directed excitatory (left) and inhibitory (right) connections weighted by their strength and normalized by the total number of possible connections. The circular distribution of excitatory connections for time window [ — 50,100] ms in purple overlies that for time window [150, 300] ms in red. The circular distributions of inhibitory connections for the same windows are in teal and blue, respectively. The dashed black lines represent the beta wave axis. (c) Same as in a, but using neurons with firing rates >1 spikes — 1 and wide spike widths (>0.4000 ms).

cortical sheet, we show for the first time spatiotemporal spike patterning that closely matches propagating wave activity as measured by LFPs in terms of both its spatial anisotropy and its transmission velocity. Although the horizontal anisotropy of directed effective connections among the narrow-spiking neurons has never before been experimentally documented (however, effective anisotropy was used in a modelling study22), anatomical data on the horizontal connectivity of motor cortex may shed light on these observations. By staining degenerating axons and terminals after localized lesions within MI, Gatter and Powell documented a preponderance of horizontal axons in layers 2/3 and 5 in non-human primates that were more spatially extensive along the rostro-caudal dimension as compared with the mediolateral dimension35. In addition, Huntley and Jones36 used retrograde labelling injections in primate motor cortex and found weak strips of labelled cell bodies arranged in the rostro-caudal dimension. However, these data do not provide evidence for stronger (anatomical or functional) connections in the rostro-caudal dimension.

Physiological data also suggest that the functional connections we observed may reflect properties of horizontal connectivity within motor cortex. Electrical stimulation and current-source density analysis in rat motor cortical slices found long-distance excitatory, monosynaptic interactions mediated by horizontal connections in layers 2/3 and 5 with a conduction velocity of b 10cms_ 1, which is on the same order of magnitude as the propagation velocities we documented with the LFP wave and the excitatory connections in our effective connectivity analysis37. Based on spatially smoothed multiunit responses recorded from the auditory cortex, propagating activity was documented with

comparable propagation speeds of ~ 20 cm s _ 1 using a similar gradient field method11.

Although the beta wave propagation axes and propagation speeds were consistent across all the animals used in the study, there are two issues to be resolved in future studies. First, the orientation of the primary beta wave direction was not consistent across the three monkeys. It is not obvious why this would be the case, but it may be owing to the exact placement of the arrays in the medial/lateral dimension along the motor strip. Second, our monkeys exhibited a fairly wide range of hand movement speeds while wave propagation speeds were fairly consistent. We would like to investigate in more detail whether there is a relationship between the kinematics, kinetics and muscle activity of the limb, and the statistical properties of the beta waves.

Our observation that population activity among narrow-spiking neurons exhibited a pronounced beta spectral peak suggests that beta oscillations in the motor cortex may be mediated by a network of functionally distinct cells38. Simulation studies have suggested that beta oscillations emerge in network interactions of local inhibitory neurons22, which in turn entrain excitatory pyramidal neurons to fire rhythmically39. These local inhibitory and excitatory pyramidal cells may represent the population of narrow-spiking neurons that tend to oscillate at the beta frequency range. In contrast, wide-spiking neurons as a population do not resonate at any particular frequency, and, moreover, exhibit a reduction in oscillatory power in the beta and higher frequency bands around the time at which LFP beta power increases.

The temporal dynamics of planar wave activity, spatially oriented effective connections and motor behaviour indicate the following chain of events. The degree of planar wave activity

x 10-3 2

x 10-3 6 5 4 3 2 1 0

x 10-3 1.8

0 100 Time

200 300 (ms)

0 100 200 300 Time (ms)

Figure 7 | Target-direction information content of spike sequencing between neurons with narrow spike waveforms. Data from monkeys Rs (top), Mk (middle) and Rj(bottom). (a) Time-resolved, mutual information between spike sequencing and target direction (that is, the direction from the location of the acquired target to the new target). Information was calculated for only pairs of neurons with excitatory connections. The colour map consists of 18 rows corresponding to 18 orientations (20-degree bins) aligned to anatomical orientations denoted with M, medial; C, caudal; L, lateral and R, rostral, with information in bits in blue-red false colour. (b) Circular distributions of spike sequencing information with respect to the orientation of neuron pairs on the cortical surface. Spike sequencing information was summed over the time window in [100, 350] ms and over pairs within each orientation and then normalized. The dashed black lines represent the beta wave axis defined by Cw (caudal wave direction) and Rw (rostral wave direction) as in Fig. 1d. (c) Target-direction information content of sequencing firing of non-connected narrow-spiking neurons.

begins increasing at target appearance and reaches a maximum at 100-150 ms later. Thus, this enhancement of planar wave activities in the beta oscillation range appears to be triggered by arrival of information about the visual target in MI. At roughly 150 ms after target onset, the rostro-caudal anisotropy of excitatory effective connections emerges (Figs 3 and 4, and Supplementary Fig. 4), and sequential firing among these connected neurons begins carrying information about target direction. Finally, movement speed begins to increase 200 ms after target appearance resulting in movement to the target. It is our view that spatiotemporal dynamics of the LFP provides a mesoscopic perspective of patterning that is generated by spatiotemporal arrangements among groups of single units. Therefore, planar wave activity reflects the underlying spatiotemporal distribution of spiking activity among the units. A recent modelling study22 showed spatial anisotropy among inhibitory cells, while our study showed spatial anisotropy in excitatory connections. However, a recent study using an immature cerebral cortex preparation showed that cortical waves in layer III were blocked by glutamatergic receptor antagonists, but not by GABAergic receptor antagonists40. Furthermore, a modelling study showed that long-range coherent oscillations with some delay were partially induced by interactions of excitatory cells41.

These results may have implications for the development of cortically controlled brain machine interfaces. Current decoding algorithms for brain machine interfaces have ignored the spatial organization of cortical neurons and their functional interactions.

In principle, by augmenting the GLM framework, we used here by including external movement covariates, we could invert these spike prediction models to decode movement parameters. Furthermore, a recent study has speculated that spatiotemporally patterned electrical stimulation in motor cortex (similar to the recorded patterns we have observed in our study) could replace deep brain stimulation to treat symptoms of Parkinson's disease42.

Our information analysis of sequential firing may have important behavioural implications. First, we have shown that functionally connected, narrow-spiking neurons carry more task-relevant information in their sequential firing than do unconnected, narrow-spiking neurons. This implies that there is a subnetwork of neurons within the motor cortex that may be playing a special role in representing and transmitting information to guide behaviour. Second, this behavioural information is most pronounced in the sequential firing of neurons that are spatially aligned along the rostro-caudal axis, implying that the motor cortex possesses information-carrying, spatial structure beyond somatotopy. Future work will investigate whether downstream targets such as the muscles are sensitive to and may take advantage of this spatial structure when driving motor behaviour.

Methods

Behavioural tasks. Three male macaque monkeys (Macaca mulatta) were oper-antly trained to move a cursor appearing above the monkey's hand location to targets projected onto a horizontal, reflective surface in front of the monkey. The monkey's arm rested on arm troughs secured to links of a two-joint exoskeletal

robotic arm (BKIN Technologies, ON, Canada) underneath the projection surface. The shoulder joint was abducted 90 degrees such that shoulder and elbow flexion/ extension movements were made in the horizontal plane. The shoulder and elbow joint angles during movements were sampled at 500 Hz by the robotic arm's motor encoders. The position of the hand was computed using the forward kinematics equations. The monkeys performed the random RTP task, in which a sequence of seven targets appeared on the projection surface. At any one time, a single target appeared at a random location in the workspace, and the monkey was required to move the cursor to it. As soon as the cursor reached the target, the target disappeared and a new target appeared in a new, random location. After reaching the seventh target, the monkey received water reward.

Electrophysiology. A silicon-based electrode array composed of 100 electrodes (1.0 mm electrode length and 400 mm interelectrode distance) was implanted in the arm area of MI of each monkey. During a recording session, signals from up to 96 electrodes were amplified (gain, 5,000), band-pass filtered between 0.3 and 7.5 kHz, and recorded digitally (14-bit) at 30 kHz per channel using a Cerebus acquisition system (Blackrock Microsystems, Inc., UT). Only waveforms (1.6 ms in duration; 48 sample time points per waveform) that crossed a threshold were stored and spike sorted using Offline Sorter (Plexon, Inc., TX). Signal-to-noise ratios were defined as the difference in mean peak-to-trough voltage divided by twice the mean s.d. The mean s.d. was computed by measuring the s.d. of the spike waveform at each of the 48 sample time points of the waveform and then averaging. All isolated single units used in this study possessed signal-to-noise ratios of 3:1 or higher. The data for each neuron were converted to a binary time series at a 1-ms time resolution.

A total of four data sets (one data set for animal Rs, two data sets from animal Mk, which were 4 days apart, and one data set for animal Rj,) were analysed, where a data set is defined as all simultaneously recorded neural data collected in one recording session. Each data set contained between 59 and 115 simultaneously recorded units from MI. A total of 302 single-unit samples were recorded from MI over all four data sets. We use the term 'samples' because the recordings were made from chronically implanted arrays in each monkey. Therefore, data collected over recording sessions in monkey Mk are not necessarily from completely different ensembles of units. Ensembles consisted of 'randomly' selected units from MI except for a possible bias for neurons with large cell bodies that would generate higher signal-to-noise ratios. All of the surgical and behavioural procedures were approved by the University of Chicago's Institutional Animal Care and Use Committee and conform to the principles outlined in the Guide for the Care and Use of Laboratory Animals (NIH publication no. 86-23, revised 1985).

Beta wave characterization. We used the wave characterization method as defined in Rubino et al.13. Briefly, we first applied a band-pass filter bidiretionally whose center frequency for each data set was defined to be the most prominent peak frequency over the beta oscillation range (10-40 Hz) and the width was ± 3 Hz, then applied the Hilbert transform to a filtered voltage signal recorded by an electrode located at (x,y) on the array V(x,y,t) to obtain the instantaneous amplitude A(x,y,t) and phase as j(x,y,t) follows:

V(x,y, t) + iHb[V(x,y, t)] = A(x,y, t)e!

j(x,y,t)

where Hb is the Hilbert transform operator and A(x,y,t) is the instantaneous amplitude of the voltage signal, and i2 = — 1. The wave velocity v(t) = [dx/dt,dy/dt] was computed by taking the derivative of a constant contour of

with respect to time:

j(x,y, t) = const.,

dj „ j = V ■v + j =

Since the direction of the velocity, contour, the speed is

- Vj, is perpendicular to the constant phase

Note that velocity is well defined only when the phase gradient is non-zero and when the signal exhibits a well-defined propagation direction. In order to measure how well the phase gradients align across the array, we used the quantity called the PGD(t):

PGD(t) =

||Vw|| IIVwII

The bar here denotes the spatial average. If the phase gradients at all spatial locations align at a given time, t, then PGD(t) = 1. If the phase gradients are randomly distributed, then PGD will be close to 0. Thus, we used a threshold value of 0.2, which is based on the tail of the PGD distribution when channel assignments are randomly shuffled43. As long as PGD(t) was above the threshold, then we used the estimated values of mean wave direction and speed across the

array, which were computed as:

wave direction (t) = — Vw,

sPeed(t) = %

Unit spiking data selection. We created trial-based data sets using a time interval — 100 to 350 ms relative to visual target onset. We had 2,547, 2,746, 1,999 and 1,526 trials from each of the four data sets. We first selected only neurons whose average firing rates during the interval of interest, — 100 to 350 ms relative to visual target onset, were higher than 1 spikes — 1. Then, we fitted a mixture of Gaussians to the distribution of waveform widths where the width was defined as the duration from trough-to-peak of the mean waveform of each sorted unit. We used mixture Gaussian models from 1 to 4 Gaussians for each data set and the pooled data. Then, we computed the BIC for each data set and the pooled data set for each model to find a most parsimonious Gaussian mixture model. Every data set and the pooled data set had a minimum BIC at 2 (Fig. 2b). Thus, we decided to fit each data set with a mixture of two Gaussians. For each data set, we computed the average of the two means from the two Gaussians. The smallest average of all data sets fell between 0.2667 and 0.3000 ms, and neurons with waveform widths of 0.2667 ms or less were classified as narrow-spiking units. Similarly, the largest average of the two means fell between 0.3667 and 0.4000 ms so that neurons with waveform widths of 0.4000 ms or more were classified as wide-spiking units. Thus, we had two thresholds. We decided to take this conservative threshold scheme to clearly separate the two populations of neurons. The numbers of neurons that satisfied the narrow condition were 38 neurons (4.89 ± 4.10 spikes s — 1, mean ± s.d.) among a total of 115 neurons for the first data set (monkey Rs), 14/21 neurons (11.94 ± 11.48 and 9.71 ± 8.59 spikes s — 1, respectively) among a total 59/65 neurons for the second and third data sets (monkey Mk), and 33 neurons (11.68 ± 9.00 spikes s — 1) among a total of 63 neurons for the fourth data set (monkey Rj). The numbers of wide-spiking neurons were 29 neurons (2.96 ± 1.94 spikes s — 1) among a total of 115 neurons for the first data set (monkey Rs), 13/10 neurons (4.22 ± 4.48 and 4.88 ± ±5.69 spikes s — 1, respectively) among a total of 59/65 neurons for the second and third data sets (monkey Mk), and 6 neurons (5.59 ± 4.45 spikes s — 1) among a total of 63 neurons for the fourth data set (monkey Rj). Neurons that did not fall into either the narrow or wide class were not further analysed. Peristimulus time histograms binned for two example neurons used in the analyses are presented based on the spike-width class for each data set in Fig. 2c, and Supplementary Figs 1 and 2. We also broke up each data set into three subsets of data consisting of an equal number of trials for the first two monkeys, and into two subsets of data for the third monkey to check for temporal consistency. We analysed each trial using 150 ms time window.

Power spectra and spectrograms for spiking of population of neurons. We

binned each neuron's spike train into 1 ms bins and summed the binned spike trains over all neurons within the narrow or wide spike waveform population. Then, we used a multitaper spectrum estimation method44 (time-bandwidth product and number of tapers were 3 and 5, and moving window of 250 ms in 2 ms increments) to compute average power spectra over [ — 100, 300] ms relative to the visual cue onset and average spectrograms over the time windows specified above using Chronux toolbox.

Granger causality measure for point process models. We utilized a modelbased methodology for analysing the causal interactions of multiple neural spike trains34. The discrete, all-or-nothing nature of a sequence of action potentials together with their stochastic structure suggests that neural spike trains may be regarded as point processes45,46. Let Nj(t) denote the sample path that counts the number of spikes of neuron i in the time interval (0,t] for t a (0,T] for i = 1,...,M recorded neurons. A point process model of a spike train for neuron i can be completely characterized by its CIF, 1j(t|H(t)), defined as:

l (t|H(t)) :

Pr[N (t + A) - Ni (t) = 1|H(t)]

where H(t) denotes the history of the spiking activity of all neurons up to time t. The CIF represents the instantaneous firing probability and serves as a fundamental component for constructing the likelihoods and probability distributions for point process data. It is a history-dependent function, and reduces to a Poisson process if it is independent of the history. To simplify the notation, we denote 1j(t|H(t)) as 1j(t). We used a point process GLM for building models whose history terms contain past spiking activity of other neurons34. To test whether neuron j causes neuron i in the 'Granger' sense47, we developed two classes of point process GLMs for neuron i: one given by P(Nj(0:T)) that has the past of all neurons as the covariates for the CIF, and the other given by Q(Nj(0:T)) that has the past of all except for neuron j. The point process likelihood is given by45

P(Ni(0 : T)) = exp jy log A¿(t)dN¿(t)- J !¿(t)dt

log Mt)=A + £fi,m ■ hm

In equation (9), we used a point process GLM to model the relationship between the spiking activity and its covariates (the spiking history)45. Here, bi relates to a background level of the activity of neuron i, fm is a 5-dimensional (5D) vector of parameters that relate the past spiking of neuron m to the current spiking of neuron i, and hm is a 5D vector each of whose elements counts the spikes in every 3 ms time window of 15 ms spiking history of neuron m, these parameters were chosen to obtain a relatively small number of parameters while maintaining the necessary temporal resolution. The ' •' represents the dot product between vectors. Note that for P(Ni(0:T)), li(t) includes the past spiking of all neurons. The other point process likelihood, Q(Ni(0:T)), is given by the same equation (8), but with li(t) replaced with A,-(t):

logXi(t) = ~i fm ■ hm

Causal inference with shuffled trials. We estimated causality networks after shuffling the trial order of each neuron for the same target direction using the data after the onset of visual cue. Through this analysis, we tested how much the precise spiking times among multiple neurons contributed to significant functional interactions between them. If significant functional interactions were simply caused by an increase in firing rate after the onset of visual cue, shuffling the trials would not affect the network estimates. On the other hand, if they were caused rather by the increase of relative spike timings among multiple neurons, shuffling the trials would break the relative spike timings among neurons and thus lead to a sparser network estimate.

For this analysis, we used the data set in the window 5 (after the onset of cue) of monkeys Mk and Rs. For each target direction, we permuted the order of all the trials for each neuron and then estimated the causality network using the shuffled data set. We repeated this procedure 100 times. The mean number of connections was 2.05 with an s.d. of 1.79 compared with 56, 65 and 66 connections for monkey Mk, and the mean number of connections was 12.40 with an s.d. of 1.50 compared with 100, 105 and 128 connections for monkey Rs from the original unshuffled network. This analysis suggests that the increase of significant interactions between neurons after the onset of visual cues is caused by more organized spike sequencing among neurons rather than an increase of spiking rate.

which excludes the effect of the past spiking of neuron j. Model parameters

0i = {bi/1,1,--.,/1,M} for P(Ni(0:T);0i) and 6 = [fin {fi,m)M=lm^ for

Q(N(0 : T); 8^ were fitted by maximizing the likelihoods 48, and they were calculated by

0* = argmaxP(N(0 : T); 0),

= argmax^N¡(0 : T); .

Because of the structure of the equation (8), the maximum likelihood solution can be efficiently calculated. Finally, we define our causality measure from j to i as

Gij = log

P(Nj(0 : T); 8*] q(Ni(0 : T);~*

Mutual information analysis for single-unit spiking activities and spike sequencing. We evaluated the information content of single-unit spiking activities and spike sequencing using the mutual information between the target direction and single-cell spiking, or spike sequencing between two neurons for each cortical orientation, respectively. Mutual information between two variables represents the reduction in uncertainty of one variable due to the knowledge of the other. It is defined as

I(X, Y) = H(X) — H(X\Y),

for random variables X and Y where H(X) is the Shannon entropy of X (ref. 52). In this paper, Y is a random variable representing the target directions, and takes on one of eight directions (the eight octants in 2D plane) during the RTP task. For single-neuron analysis, X represents the binary sequence of unit spiking activity. For spike sequence analysis, X represents the spike sequencing from neuron i to j defined by

X(t) = £ ni(t — t)nj(t),

If the spiking history of neuron j helps predict the spiking activity of neuron i, the log-likelihood ratio should be significantly greater than zero, and then we say that neuron j, 'Granger-causes' i49. The equality of the likelihood ratio holds when neuron j has no causal influence on i. Excitatory and inhibitory influences from neuron j to i are distinguished by the sign of t f'j (T) in equation (9), which represents the summed influence of the past spiking of j on i.

Significance test. The Granger causality measure, DD = D0 — Dx = 2^ given by the log-likelihood ratio provides an indication of the relative strength of causal interaction but gives little insight into whether or not it is statistically significant. We use a statistic based on the log-likelihood ratio test to address this issue. We denote the deviance obtained using the model parameter 0* as D0, and the deviance obtained using 6* as D1. The deviance difference between two models is equivalent to two times log-likelihood ratio given by DD = D0 — D1 = 2^ (ref. 50). If both models describe the data well, then the deviance difference may be asymptotically described by DD ~ where d is the difference in dimensionality of two models, and d = 5 in our study. If the value of DD is in the critical region, that is, greater than the upper tail 100(1 — a)% of the where a determines false-positive rates, then the causal influence is determined as statistically significant. We used a multiple-hypothesis testing error measure called false discovery rate to address the multiple comparison problem51.

Assessment of GOF against the data for the point process models. The

validation of a designed statistical model for data is required before making an inference based on the model. The time-rescaling theorem was used to transform point process neural spike train data to a continuous measure appropriate for GOF assessment33. Once a CIF is estimated, rescaled times can be computed using the estimated CIF. These rescaled times will be uniformly distributed random variables on the interval (0,1] if the estimated CIF is a good approximation to the true conditional intensity of the point process. To evaluate whether the rescaled times followed the uniform distribution, we ordered these rescaled times from the smallest to the largest, and then plotted the quantiles of the cumulative distribution function of the uniform distribution on (0,1] against the ordered rescaled times. This form of graphical representation is termed a KS plot. If the model is consistent with the data, then the points should lie on a 45-degree line. Approximate 95% confidence bounds for the degree of agreement between the model and the data may be constructed using the distribution of the KS statistic. We defined the GOFA bounded by the 45-dergee line and the cumulative distribution in the KS plot. We also defined the GOFR computed as the ratio of GOFA to an area out of the ~95% confidence bounds for the degree of agreement (see Supplementary Fig. 3a for schematics of the method). These measures were computed for each neuron in each 150 ms time window.

which is calculated as the sum of products in spike counts between neuron i at varying time lags t (from t1 to t2) and neuron j. For the spike sequencing, we used t = 1ms and = 15 ms. We chose a 15-ms summation window to match with the duration of history terms used for the GLM described above. Information content of the spike sequencing was calculated using the empirical distributions of X(t) at fixed i and j, and t over all trials for different target directions. Using the position information of neurons on the multielectrode array, we calculated the cortical orientation of the information content of spike sequencing. To account for biases in estimating information from limited data, the target direction labels were randomly shuffled over all trials so that a trial originally associated with target direction k could be randomly assigned one of the eight target directions. Mutual information between target direction and spike sequencing was then calculated. The average information from 10 shuffles was subtracted from the estimated information values.

All the numerical analyses were performed using Matlab (Mathworks, MA). Matlab codes for Granger analysis are available at the website (http://www. neurostat.mit.edu/gcpp).

References

1. Prechtl, J. C., Cohen, L. B., Pesaran, B., Mitra, P. P. & Kleinfeld, D. Visual stimuli induce waves of electrical activity in turtle cortex. Proc. Natl Acad. Sci. USA 94, 7621-7626 (1997).

2. Arieli, A., Shoham, D., Hildesheim, R. & Grinvald, A. Coherent spatiotemporal patterns of ongoing activity revealed by real-time optical imaging coupled with single-unit recording in the cat visual cortex. J Neurophysiol. 73, 2072-2093 (1995).

3. Roland, P. E. et al. Cortical feedback depolarization waves: A mechanism of top-down influence on early visual areas. Proc. Natl Acad. Sci. USA 103, 12586-12591 (2006).

4. Benucci, A., Frazor, R. A. & Carandini, M. Standing waves and traveling waves distinguish two circuits in visual cortex. Neuron 55, 103-117 (2007).

5. Lippert, M. T., Takagaki, K., Xu, W., Huang, X. & Wu, J.-Y. Methods for voltage-sensitive dye imaging of rat cortical activity with high signal-to-noise ratio. J. Neurophysiol. 98, 502-512 (2007).

6. Xu, W., Huang, X., Takagaki, K. & Wu, J. Compression and reflection of visually evoked cortical waves. Neuron 55, 119-129 (2007).

7. London, J. A., Cohen, L. B. & Wu, J. Y. Optical recordings of the cortical response to whisker stimulation before and after the addition of an epileptogenic agent. J. Neurosci. 9, 2182-2190 (1989).

8. Petersen, C. C. H., Grinvald, A. & Sakmann, B. Spatiotemporal dynamics of sensory responses in layer 2/3 of rat barrel cortex measured in vivo by voltage-sensitive dye imaging combined with whole-cell voltage recordings and neuron reconstructions. J. Neurosci. 23, 1298-1309 (2003).

9. Ferezou, I. et a/. Spatiotemporal dynamics of cortical sensorimotor integration in behaving mice. Neuron 56, 907-923 (2007).

10. Song, W. J. et a/. Cortical intrinsic circuits can support activity propagation through an isofrequency strip of the guinea pig primary auditory cortex. Cereb.Cortex 16, 718-729 (2006).

11. Witte, R. S., Rousche, P. J. & Kipke, D. R. Fast wave propagation in auditory cortex of an awake cat using a chronic microelectrode array. J Neural Eng. 4, 68 (2007).

12. Ferezou, I., Bolea, S. & Petersen, C. C. H. Visualizing the cortical representation of whisker touch: voltage-sensitive dye imaging in freely moving mice. Neuron 50, 617-629 (2006).

13. Rubino, D., Robbins, K. A. & Hatsopoulos, N. G. Propagating waves mediate information transfer in the motor cortex. Nat. Neurosci. 9, 1549-1557 (2006).

14. Takahashi, K., Saleh, M., Penn, R. D. & Hatsopoulos, N. Propagating waves in human motor cortex. Front. Hum. Neurosci. 5, 40 (2011).

15. Murthy, V. N. & Fetz, E. E. Oscillatory activity in sensorimotor cortex of awake monkeys: synchronization of local field potentials and relation to behavior. J. Neurophysio/. 76, 3949-3967 (1996).

16. Murthy, V. N. & Fetz, E. E. Coherent 25- to 35-Hz oscillations in the sensorimotor cortex of awake behaving monkeys. Proc. Nat/ Acad. Sci. USA 89, 5670-5674 (1992).

17. Sanes, J. N. & Donoghue, J. P. Oscillations in local field potentials of the primate motor cortex during voluntary movement. Proc. Nat/ Acad. Sci. USA 90, 4470-4474 (1993).

18. Jackson, A., Gee, V. J., Baker, S. N. & Lemon, R. N. Synchrony between neurons with similar muscle fields in monkey motor cortex. Neuron 38, 115-125 (2003).

19. Szurhaj, W. et a/. Basic mechanisms of central rhythms reactivity to preparation and execution of a voluntary movement: a stereoelectroencephalographic study. C/in. Neurophysio/. 114, 107-119 (2003).

20. Pfurtscheller, G. Central beta rhythm during sensorimotor activities in man. E/ectroencepha/ogr. C/in. Neurophysio/. 51, 253-264 (1981).

21. Heitmann, S., Gong, P. & Breakspear, M. A computational role for bistability and traveling waves in motor cortex. Front. Comput. Neurosci. 6 (2012).

22. Heitmann, S., Boonstra, T. & Breakspear, M. A dendritic mechanism for decoding traveling waves: principles and applications to motor cortex. PLoS Comput. Bio/. 9, e1003260 (2013).

23. Reimer, J. & Hatsopoulos, N. G. Periodicity and evoked responses in motor cortex. J. Neurosci. 30, 11506-11515 (2010).

24. Bartho, P. et a/. Characterization of neocortical principal cells and interneurons by network interactions and extracellular features. J. Neurophysio/. 92, 600-608 (2004).

25. Merchant, H., Naselaris, T. & Georgopoulos, A. P. Dynamic sculpting of directional tuning in the primate motor cortex during three-dimensional reaching. J. Neurosci. 28, 9164-9172 (2008).

26. Mitchell, J. F., Sundberg, K. A. & Reynolds, J. H. Differential attention-dependent response modulation across cell classes in macaque visual area V4. Neuron 55, 131-141 (2007).

27. Kaufman, M. T. et a/. Roles of monkey premotor neuron classes in movement preparation and execution. J. Neurophysio/. 104, 799-810 (2010).

28. Peyrache, A. et a/. Spatiotemporal dynamics of neocortical excitation and inhibition during human sleep. Proc. Nat/ Acad. Sci. USA 109, 1731-1736 (2012).

29. Markram, H. et a/. Interneurons of the neocortical inhibitory system. Nat. Rev. Neurosci. 5, 793-807 (2004).

30. Somogyi, P., Tamas, G., Lujan, R. & Buhl, E. H. Salient features of synaptic organisation in the cerebral cortex. Brain Res. Brain Res. Rev. 26, 113-135 (1998).

31. Marre, O., El Boustani, S., Fregnac, Y. & Destexhe, A. Prediction of spatiotemporal patterns of neural activity from pairwise correlations. Phys. Rev. Lett. 102, 138101 (2009).

32. Burnham, Kenneth P. & Anderson, D. R. Mode/ Se/ection and Mu/timode/ Inference: a Practica/ Information-Theoretic Approach (Springer, 2002).

33. Brown, E. N., Barbieri, R., Ventura, V., Kass, R. E. & Frank, L. M. The time-rescaling theorem and its application to neural spike train data analysis. Neura/ Comput. 14, 325-346 (2002).

34. Kim, S., Putrino, D., Ghosh, S. & Brown, E. N. A Granger causality measure for point process models of ensemble neural spiking activity. PLoS Comput. Bio/. 7, e1001110 (2011).

35. Gatter, K. C. & Powell, T. P. S. The intrinsic connections of the cortex of area 4 of the monkey. Brain 101, 513-541 (1978).

36. Huntley, G. W. & Jones, E. G. Relationship of intrinsic connections to forelimb movement representations in monkey motor cortex: A correlative anatomic and physiological study. J. Neurophysio/. 66, 390-413 (1991).

37. Aroniadou, V. A. & Keller, A. The patterns and synaptic properties of horizontal intracortical connections in the rat motor cortex. J.Neurophysio/. 70, 1553-1569 (1993).

38. Womelsdorf, T., Valiante, T. a, Sahin, N. T., Miller, K. J. & Tiesinga, P. Dynamic circuit motifs underlying rhythmic gain control, gating and integration. Nat. Neurosci. 17, 1031-1039 (2014).

39. Jensen, O. et a/. On the human sensorimotor-cortex beta rhythm: sources and modeling. Neuroimage 26, 347-355 (2005).

40. Namiki, S. et a/. Layer III neurons control synchronized waves in the immature cerebral cortex. J. Neurosci. 33, 987-1001 (2013).

41. Bibbig, A., Traub, R. D. & Whittington, M. A. Long-range synchronization of g and b oscillations and the plasticity of excitatory and inhibitory synapses: a network model. J. Neurophysio/. 88, 1634-1654 (2002).

42. Beuter, A., Lefaucheur, J.-P. & Modolo, J. Closed-loop cortical neuromodulation in Parkinson's disease: An alternative to deep brain stimulation? C/in. Neurophysio/. 125, 874-885 (2014).

43. Woods, B. Spatio-temporal patterns in multi-electrode array local field potential recordings. arXivPreprint at <http://arxiv.org/abs/1501.00230> (2014).

44. Mitra, P. P. & Pesaran, B. Analysis of dynamic brain imaging data. Biophys. J. 76, 691-708 (1999).

45. Truccolo, W., Eden, U. T., Fellows, M. R., Donoghue, J. P. & Brown, E. N. A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. J. Neurophysio/. 93, 1074-1089 (2005).

46. Brown, E. N. in Methods and Mode/s in Neurophysics. (eds Chow, C. C., Gutkin, B., Hansel, D., Meunier, C. & Dalibard, J.) 691-726 (Elsevier, 2005).

47. Brown, E. N., Barbieri, R., Eden, U. T. & Frank, L. M. in Computationa/ Neuroscience: a Comprehensive ApproachFeng, J.) 253-286 (CRC, 2003).

48. Daley, D. J. & Vere-Jones, D. An Introduction to the Theory of Point Processes (Springer.

49. Granger, C. W. J. Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424-438 (1969).

50. McCullagh, P. & Nelder, J. Generalized Linear Mode/s (Chapman and Hall, 1989).

51. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B Stat. Methodo/. 57, 289-300 (1995).

52. Cover, T. & Thomas, J. E/ements of Information Theory (John Wiley and Sons, 1991).

Acknowledgements

We thank Bronwyn Woods for the channel shuffling analysis to identify a significance threshold for PGD that was used to determine planar wave propagation in our LFP data. We thank Lorenzo Pesce and the Beagle team at the Computation Institute for testing Matlab-based Granger causality method on the Beagle to assess empirical minimum spike rate threshold for numerical stability, and performed shuffle analysis and calculations of network parameters for large data sets. We also thank Josh Coles, Dawn Paulsen, Jacob Reimer, Dennis Tkach and Zach Haga for their help with the surgical implantation of the arrays, training of monkeys, data collection and discussions. K.T., K.A.B., M.D.B. and N.G.H. were supported by a NIH grant R01 NS045853 from the National Institute of Neurological Disorders and Stroke (NINDS), A.J.S. and N.G.H. were supported by R01 NS048845 from the NINDS, and S.K. and T.P.C. were supported by Science and Technology Center grant CCF-0939370 and CCF 10-65352 from NSF.

Author contributions

K.T. and N.G.H. designed the study; A.J.S. and past members of our laboratory collected data; K.T., S.K. and K.A.B. performed analysis. All authors discussed the results and contributed to writing the paper.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/ naturecommunications

Competing financial interests: The authors declare no competing financial interests.

Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/

How to cite this article: Takahashi, K. et a/. Large-scale spatiotemporal spike patterning consistent with wave propagation in motor cortex. Nat. Commun. 6:7169 doi: 10.1038/ ncomms8169 (2015).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/