# LFP beta amplitude is linked to mesoscopic spatio-temporal phase patterns

Michael Denker, Lyuba Zehl, Bjorg Kilavik, Markus Diesmann, Thomas Brochier, Alexa Riehle, Sonja Grün

► To cite this version:

Michael Denker, Lyuba Zehl, Bjorg Kilavik, Markus Diesmann, Thomas Brochier, et al. LFP beta amplitude is linked to mesoscopic spatio-temporal phase patterns. Scientific Reports, 2018, 8 (1), pp.5200 - 5200. 10.1038/s41598-018-22990-7. hal-01794053

HAL Id: hal-01794053

https://hal.science/hal-01794053

Submitted on 17 May 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_2.png" /></figure>

OPEN

# LFP beta amplitude is linked to mesoscopic spatio-temporal phase patterns

Received: 4 August 2017

Accepted: 5 March 2018

Published online: 26 March 2018

Michael Denker<sup>1</sup>, Lyuba Zehl<sup>1</sup>, Bjørg E. Kilavik<sup>2</sup>, Markus Diesmann<sup>1,3,4</sup>, Thomas Brochier<sup>1,2</sup>, Alexa Riehle<sup>2,1,5</sup> & Sonja Grün<sup>1,5,6</sup>

Beta oscillations observed in motor cortical local field potentials (LFPs) recorded on separate electrodes of a multi-electrode array have been shown to exhibit non-zero phase shifts that organize into planar waves. Here, we generalize this concept to additional classes of salient patterns that fully describe the spatial organization of beta oscillations. During a delayed reach-to-grasp task we distinguish planar, synchronized, random, circular, and radial phase patterns in monkey primary motor and dorsal premotor cortices. We observe that patterns correlate with the beta amplitude (envelope): Coherent planar/radial wave propagation accelerates with growing amplitude, and synchronized patterns are observed at largest amplitudes. In contrast, incoherent random or circular patterns are observed almost exclusively when beta is strongly attenuated. The occurrence probability of a particular pattern modulates with behavioral epochs in the same way as beta amplitude: Coherent patterns are more present during movement preparation where amplitudes are large, while incoherent phase patterns are dominant during movement execution where amplitudes are small. Thus, we uncover a trigonal link between the spatial arrangement of beta phases, beta amplitude, and behavior. Together with previous findings, we discuss predictions on the spatio-temporal organization of precisely coordinated spiking on the mesoscopic scale as a function of beta power.

The local field potential (LFP) has long served as a readily available brain signal to monitor the average input activity that reaches the neurons in the vicinity of extracellular recording electrodes<sup>1-3</sup>. A hallmark of the LFP is the ubiquitous presence of oscillations in various frequency bands<sup>4</sup> modulating in time and across different brain structures. These oscillations have been linked to a variety of brain processes such as attention<sup>5</sup>, stimulus encoding<sup>6</sup>, or memory formation<sup>7,8</sup>. These findings support the basis of modern theories concerning the functional implication of oscillatory brain activities, such as feature binding<sup>9</sup>, the concept of communication-through-coherence<sup>10-12</sup>, the phase-of-firing coding<sup>13</sup>, or predictive coding<sup>14</sup>. In motor cortex, beta oscillations (in the range of 15-35 Hz) are one of the most prominent types of oscillatory activity. They have been linked to states of general arousal, movement preparation, or postural maintenance<sup>15,16</sup>, and are typically suppressed during active movement<sup>17,18</sup>.

Technological progress recently led to the development of multi-electrode arrays enabling neuroscientists to record massively parallel neuronal signals in a precisely identifiable spatial arrangement. Although LFP signals recorded in motor cortex from electrodes separated by up to several millimeters are typically highly correlated<sup>19</sup>, the analysis of the instantaneous phase of the oscillation<sup>20</sup> revealed a non-zero temporal shift between electrodes<sup>21</sup>. Such shifts may be expressed by the formation of dynamic spatial patterns propagating along preferred directions across the cortical surface, referred to as traveling waves<sup>22</sup>. Indeed, the phenomenon of traveling waves has been described in multiple brain areas, such as the visual cortex<sup>23-25</sup>, including voltage sensitive dye recordings<sup>26-30</sup>, in the turtle visual cortex<sup>31,32</sup>, the olfactory bulb<sup>33,34</sup>, the thalamus<sup>35</sup>, and in visual perception<sup>36</sup>, and during epilepsy<sup>37</sup>. While some studies reported waves in visual cortex as single-cycle propagation of elevated

<footnote n="1">¹Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6) and JARA Institute Brain Structure-Function Relationships (INM-10), Jülich Research Centre, Jülich, Germany. ²Institut de Neurosciences de la Timone (INT), CNRS-Aix-Marseille University, UMR 7289, Marseille, France. ³Department of Physics, Faculty 1, RWTH Aachen University, Aachen, Germany. ⁴Department of Psychiatry, Psychotherapy and Psychosomatics, Medical Faculty, RWTH Aachen University, Aachen, Germany. ⁵RIKEN Brain Science Institute, Wako City, Japan. ⁶Theoretical Systems Neurobiology, RWTH Aachen University, Aachen, Germany. Correspondence and requests for materials should be addressed to M.D. (email: m.denker@fz-juelich.de)</footnote>
activity from a central hotspot outwards, either induced by stimulation or occurring spontaneously³⁰, others described planar waves in visual cortex and other sensory areas, i.e., waves traveling homogeneously along a defined direction. Motor cortical waves were also described so far as being unidirectional throughout the cortical region covered by 4 × 4 mm multi-electrode arrays. The probability of observing these planar waves may rapidly change as a function of behavioral context. Indeed, it was found that the average coherence of phase gradients across electrodes, considered as being a signature of planar wave propagation, was highest at the beginning of the instructed delay of a center-out reaching task²².

Such planar waves²² represent the most salient type of dynamic pattern formation, and are easily identifiable by the parallel arrangement of the phase gradients. However, potentially different patterns of spatial organization of beta oscillations outside periods of planar waves have not yet been described. It is reasonable to assume that oscillatory activities do exhibit other types of patterns commonly associated with theoretical systems displaying pattern formation³⁸,³⁹, such as divergences or singularities. In visual cortical area MT of the anesthetized marmo-set monkey, for instance, a variety of such patterns in slow (delta) oscillations were described⁴⁰, as well as in visually evoked activity in the range of 10–20 Hz in turtle³¹. However, for motor cortical beta activity, such patterns have not been investigated.

Moreover, the visualizations of spatio-temporal phase patterns typically do not indicate the temporally strongly modulated amplitude of the beta oscillation. The occurrences of motor cortical planar waves seem to be of short duration²², in the order of 50 ms (cf. the Supplemental Information within the reference). This is evocative to the finding that motor cortical beta oscillations strongly modulate their amplitude by exhibiting short-lasting high amplitude epochs of a few oscillatory cycles, the so-called spindles¹⁹,⁴¹. Even though an individual beta spindle lasts far longer than the occurrence of a planar wave, their dynamic properties suggest that amplitude and spatial organization of beta activity are correlated. This hypothesis is further supported by the finding that when considering data of different trials, both traveling waves²² and beta power¹⁶ are most prominent during an instructed delay of a motor task. Moreover, for slow oscillations, the power of oscillations was found to correlate with the dynamics of activity patterns⁴⁰. Hence, it is an open question in how far beta amplitude in motor cortex is related to the spatial arrangement of the beta phase. Such a relationship would indeed have a profound impact on the interpretation of the power of beta oscillations: Rather than a mere measure of the average degree of neuronal synchronization in the vicinity of the recording electrode, power may represent a signature of synchronization between spatially separated cortical locations. Taken together with previous findings relating planar waves to sequential spike patterns⁴² and linking the phase of beta oscillations to the occurrence of precise synchronous events²¹, LFP beta power could be indicative of communication by means of activity propagation on the level of single spikes.

The present work has two main goals: The first goal is to explore the possible presence of wave-like spatio-temporal patterns other than planar waves. The second goal is to test whether or not patterns are related to modulations in beta amplitude, both in single trials and across trials, and to relate patterns to behavioral epochs in order to test their possible functional implication. Neuronal activity was recorded by using a 100-electrode Utah array, chronically implanted in primary motor (M1) and premotor (PM) cortices. Three monkeys were trained in an instructed-delay reach-to-grasp task⁴³⁻⁴⁵. We analyzed the spectral properties of the LFP signals and characterized the emergent spatio-temporal patterns based on the phase information. This analysis revealed a variety of spatio-temporal patterns in LFP beta oscillations that can be clearly distinguished and identified as five categories of phase patterns. We developed statistical measures to identify the different phase patterns and the periods in which each of the patterns occurred, and determined their prevalence as a function of trial progression and behavioral epochs. Using these findings, we were able to expose the tight link between the modulation of LFP beta amplitude and the formation of spatio-temporal patterns of the oscillation. Preliminary results were presented in abstract form⁴⁶,⁴⁷.

## Results

Three monkeys were trained in a delayed reach-to-grasp task (Fig. 1A) in which the animal had to grasp, pull and hold an object using either a side grip (SG) or a precision grip (PG), and either with a low force (LF) or high force (HF), resulting in a total of 4 pseudo-randomly presented trial types. The monkey was first presented with a cue for 300 ms which provided partial prior information either about the grip type (SG or PG) in grip-first trials, or the amount of force (LF or HF) in force-first trials, to be used in the upcoming movement. The cue was followed by a preparatory delay of 1 s. The GO signal, presented at the end of the delay, provided the missing information about either the force (LF or HF) or the grip type (SG or PG) in grip-first and force-first trials, respectively. The GO signal also instructed the monkey to initiate the reach and grasp movement. Each correct trial was rewarded by a drop of apple sauce. Figure 1A shows the time line of the behavioral trial. The monkeys performed sessions of about 15 min (120–140 correct trials) in which either grip-first trials or force-first trials were requested. For a complete description of the behavioral task refer to the Methods.

While the monkeys performed the task, neuronal activity was recorded using a 100-electrode Utah array (Blackrock Microsystems, Salt Lake City, UT, USA) implanted in the contralateral primary motor (MI) and premotor (PM) cortices with respect to the active arm (monkey L and N, left hemisphere, and monkey T, right hemisphere). The precise locations of the implanted arrays are shown in Fig. 1B,C. In this study we concentrate on the local field potential (LFP) signals, filtered between 0.3–250 Hz and sampled at 1 kHz. We selected for further analysis in each monkey 15 recording sessions from the grip-first condition, and additionally 15 sessions from the force-first condition in monkey L and T, respectively. In the following, we start by characterizing the spectral properties of the recorded LFP activity to identify its oscillatory features, before quantifying these oscillations also in the spatial domain.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Experimental task, array positions, and spectral properties of the LFP. (A) Task design. Top: sketch of the monkey during the task in the anticipatory position before GO (left), and while performing a side grip (middle) and precision grip (right). Bottom: time line of the task. Labels indicate events (TS: trial start; WS: warning signal; CUE-ON/OFF: cue on/cue off; GO: GO signal; SR: switch release; OT: object touch; HS: start of hold period; RW: reward). Images above the time axis indicate the state of the 5 LEDs during a grip-first condition at WS, during the presentation of the cue (CUE-ON through CUE-OFF) and at GO. (B) Spatial locations of the Utah multi-electrode arrays (green squares) on the cortical surface in monkey L (left), T (middle) and N (right). Top and bottom graphs show the array locations with respect to anatomical features (red curves) estimated from the corresponding photographs shown in panel C. CS: central sulcus; AS: arcuate sulcus; PS: precentral sulcus. (C) Photographic image of the array locations taken during surgery. (D) Power spectrum of the LFP during the complete recording of one selected central electrode (id 50), averaged across all sessions (N = 15 per panel) in the grip-first condition of monkey L (left), T (middle), and N (right). Orange shading:</figcaption></figure>
range of the applied beta band filter (cut-off frequencies). (E) Trial-averaged, time-resolved power spectrogram of the LFPs of one electrode in one recording session during PG trials of a grip-first recording. Trials aligned to TS. Color indicates logarithmic power density. Horizontal dashed lines: beta band as shown in panel D. Vertical dotted lines: trial events (SR, RW: mean times). Session IDs: l101013-002 (monkey L), t010910-001 (monkey T), and i140613-001 (monkey N) from left to right, respectively. Sketch in panel A is adapted from a detailed description of the experiment⁴⁵.

**Spectral LFP properties.** On a first glance, we observed that the LFP in all monkeys was dominated by a prominent oscillatory component in the beta range (about 15–35 Hz). By computing the average power spectrum of each monkey’s LFP, pooled for one electrode in the array center across its complete set of recordings in the grip-first condition (15 per monkey), we found that the frequency range of the beta oscillation varied between monkeys (see Fig. 1D). Based on these spectra we defined a wide frequency band (13–30 Hz) that was common to all monkeys and covered the peaks of the individual beta frequencies (shaded area in Fig. 1D). For better comparison, we applied this same filter band in the beta range to all data sets of all monkeys.

Furthermore, the observed LFP activity revealed that the trial-averaged power of the oscillatory activity was not stationary in time, but was strongly modulated during the time course of behavioral trials. The strength of the beta oscillations is visualized by the time-resolved power spectra, averaged on one electrode (same as for Fig. 1D) across all successful PG trials of one representative recording session of monkeys L, T and N, respectively (Fig. 1E). The beta power showed a characteristic temporal evolution that followed a similar trend for all three monkeys: the beta power was largest around the cue, and decayed gradually during the delay period and was strongly attenuated during movement execution. During movement, a low frequency signal was the most prominent component in the LFP, corresponding to the movement-related potential¹⁴.

The inspection of single-trial LFP signals revealed, in addition to the beta power modulations observed in trial averages, a modulation of the instantaneous amplitude of beta activity (Fig. 2A) on a much shorter time scale. Such epochs of increased beta activity comprising a few oscillation cycles are commonly referred to as beta spindles¹⁹,⁴¹. During a single trial, LFP signals recorded in parallel from all electrodes of the Utah array exhibited in general a high degree of correlation (Fig. 2), and in particular spindles occurred simultaneously on all electrodes¹⁹. However, across trials spindles did not reoccur at the same points in time (Fig. 2A), but instead their occurrence in time exhibited a strong degree of variability. Therefore, the trial-averaged temporal evolution of beta power (Fig. 1E) represents a measure that confounds the probability of single-trial high amplitude events, their average duration, and their average maximal amplitude⁴⁸.

**Identification of phase patterns.** Having described the principle properties of the dominant beta oscillations, we are now in a position to investigate the fine spatial patterning of this activity across all electrodes of the array. Zooming in time into the LFP signals recorded from a few neighboring electrodes during the entire trial length (Fig. 2), we calculated the beta-filtered signals (Fig. 2C, red traces). We observed that despite a high degree of similarity, the oscillatory components express small time lags across the electrodes (compare blue markers on each trace indicating oscillation peaks and troughs). To understand if there is a specific spatial patterning of the temporal lags between the signals on the different electrodes, we decomposed the beta-filtered LFP time series of each electrode *i* into the instantaneous amplitude *a*<sub>*i*</sub>(*t*), corresponding to the envelope of the filtered signal, and phase *φ*<sub>*i*</sub>(*t*) of the beta oscillation by calculating its analytic signal (see Methods). We then displayed these quantities as spatial maps *A*<sub>*xy*</sub>(*t*) and *Φ*<sub>*xy*</sub>(*t*) for amplitude and phase, respectively, representing each electrode at its spatial array position (*x*, *y*) at each time point *t*. Even though the oscillation amplitude *A*<sub>*xy*</sub>(*t*) was not the same across the array, its modulation was highly correlated between electrodes and showed a pattern that was changing slowly as compared to the time scale of the beta period (Fig. 2D and Supplemental Movie S1). This finding matches our observation that the occurrence of spindles is coherent across recording electrodes (Fig. 2).

In contrast, the phase snapshots *Φ*<sub>*xy*</sub>(*t*) showed a pronounced structure that varied on a fast time scale in the range of milliseconds (Fig. 2D and Supplemental Movie S1). While we typically observed a smooth transition of maps between consecutive time points *t*<sub>*i*</sub> and *t*<sub>*i*+1</sub> (given a sampling rate of 1 kHz), at some moments in time the maps changed their structure very rapidly. However, despite the rapid changes of the spatial structure of the maps and some discontinuities in their temporal evolution, many phase maps could clearly be classified by visual inspection into one of 5 distinct classes of spatial arrangements, in the following referred to as *phase patterns*. Representative examples of these classes of phase patterns and their temporal evolution over a time period of 20 ms are shown in Fig. 3A,B. In order to better visualize and characterize these spatial structures, we here calculated the vector field of phase gradients *Γ*<sub>*xy*</sub>(*t*) and its spatially smoothed version, the phase gradient coherence *Λ*<sub>*xy*</sub>(*t*), and display the gradient fields along with the phase maps. In the following we will briefly describe the classes of phase patterns in their most salient, idealized form.

The identification of traveling waves (Fig. 3B, top row), comparable to earlier reports²², was most prominent. In these *planar* patterns, a planar wave front traveled across the array, where the spatial period was typically larger than the array dimensions. Second, we observed a *synchronized* pattern (Fig. 3B, 2nd row), in which the signals on all electrodes were synchronized at near-zero phase lag. Complementing this state at the other extreme, we observed a *random* pattern (Fig. 3B, 3rd row), which showed randomly distributed phases across the array, and therefore no apparent phase relation between electrodes. A fourth pattern, termed *circular* pattern (Fig. 3B, 4th row), was characterized by an area near the array center around which the phases revolved. Finally, we observed a *radial* pattern (Fig. 3B, bottom row) of radially inward or outward propagating waves, which was also characterized by a point of origin near the array center. A specific type of pattern persisted for only short time periods of approximately the duration of a single beta oscillation cycle. In addition, some phase maps, termed *unclassified*, could not be unambiguously attributed to one of these 5 phase patterns.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Extraction of phase and amplitude maps. (A) LFPs (z-scored) recorded from one electrode during 10 consecutive successful trials (monkey L, session ID: 1101015-001). Trials aligned to TS = 0 ms. (B) Simultaneously recorded LFPs from 10 neighboring electrodes on the Utah array during a single trial. (C) Blow-up of the LFPs of the 10 example electrodes shown in panel B (gray traces; blue shading in panel B indicates the selected time window). Red traces: beta-filtered LFP. Blue lines: locations of peaks and troughs in the filtered LFP, i.e., phases φ = 0 and φ = π. (D) Amplitude (top) and phase (bottom) maps (shown in 4 ms steps) recorded during a 24 ms window (green shading in panel C). Color in each square indicates the amplitude and phase of the LFP at the electrode of a given position. Black squares: unconnected electrodes or electrodes rejected due to signal quality. The images are rotated to match the cortical position of the array as indicated in Fig. 1B.</figcaption></figure>

Following this first empirical identification of classes of phase patterns, we aimed at automatically classifying
the sequences of phase maps into one of these classes whenever possible. To this end, we introduced a set of 6
measures that capture features of the spatial arrangement of beta oscillations based on the phase map Φ<sub>xy</sub>(t), and
its spatial arrangement quantified by the phase gradients Γ<sub>xy</sub>(t) and the gradient coherence Λ<sub>xy</sub>(t) independently
at each time point t. The details of how to construct these measures are given in the Methods. Essentially, each of
the measures represents a feature of a given phase pattern that is characteristic for one or several of the 5 classes of
phase patterns. In the following, we give an intuitive explanation of the features relevant for each individual pattern
class. The planar patterns, described by a planar wave front traveling across the entire array, were characterized by
phase gradients that were aligned in parallel across the entire array. Thus, such a pattern was composed of a wave
front oriented perpendicularly to the gradients. The synchronized pattern was distinguished by a single phase
value at all electrodes (i.e., the array appears in a single color in Fig. 3B) and a random direction of phase gradients
across the array. Similarly, the random pattern showed no apparent local spatial organization of phase gradients,
but in contrast phases were uniformly distributed. In the circular pattern, like in the synchronized pattern, phase
gradients in all directions were observed, but in contrast the distribution of phases across the array was also uni-
form such that the visualizations in Fig. 3B contained all colors. Additionally, phase gradients were always arranged
such that they pointed clockwise or counter-clockwise around the center of the array. And finally, the radial pattern
exhibited phase gradients that, on a global view, pointed inward or outward from the array center. Thus, gradients
pointed in a direction orthogonal to that of circular patterns. Common to both circular and radial patterns, all
phase gradient directions were observed on the array and neighboring gradients on the array were similar.

Based on the 6 measures, we used a thresholding procedure (compare red dashed lines in Supplemental
Figure S2) to assign each phase pattern at time point *t* to one of the 5 classes of patterns, or, if none of the com-
bined threshold criteria was met, the phase pattern was considered unclassified. Thresholds were set empirically
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Phase patterns and their detection. (A) LFP signal (z-scored) recorded during a single trial (monkey L, Session ID:1101108-001, trial ID: 46) on a central electrode (gray) and superimposed beta-filtered LFP (red). Dashed vertical lines indicate trial events. Colored horizontal bars show time periods during which a particular type of phase pattern (compare color code in panel B, first column) was detected. The colored asterisks mark the time points of the first frame of the wave patterns shown in panel B. (B) Phase maps of one example of an automatically detected phase pattern for each type of phase pattern (rows, from top to bottom: planar wave, synchronized, random, circular, and radial). The sequence of maps in one row shows a total of 18 ms in steps of 2 ms. The pattern was initially detected in the first phase map of each row (corresponding time point indicated by an asterisk in panel A). Flow field indicated by black lines: gradient coherence map $\Lambda_{xy}(t)$; white large arrows: corresponding quadrant-averaged gradient coherence shown for visualization. Time stamps are given relative to TS.</figcaption></figure>

in such a way that they led to a conservative association of phase patterns with pattern classes, i.e., only clearly identified patterns were classified. Details of the classification process and a visualization of accepted and rejected classifications are provided in the Methods. Our classification procedure had some experimental limitations leading to unclassified patterns, such as the low spatial sampling of the 100 electrodes (400 µm inter-electrode distance) and the small spatial window of observation (4x4mm) as compared to the spatial wavelengths exhibited by some patterns. This may affect, in particular, the radial and circular patterns in which the point of origin was not at the array center, making it impossible to infer the pattern unequivocally. Additionally, observed patterns could also have represented transient dynamics from switching between patterns or even overlaps of competing patterns, which could not be properly distinguished and identified. If for any of these reasons a pattern did not fulfill the strict criteria of one of the five pattern classes described above, we considered it ambiguous and it was left unclassified.

The use of our algorithm enabled us to quantitatively disambiguate the 5 phase patterns that appeared as salient features upon visual inspection of the phase maps. The phase patterns shown in Fig. 3B were determined by using this algorithm. Figure 3A shows the LFP recorded on one electrode during one single trial, in which all classified phase patterns are marked, including those shown in panel B. The corresponding measures and thresholds used in the classification procedure are depicted in Supplemental Figure S2.

Periods of unambiguously identified phase patterns were typically of short duration and occurred interspersed throughout the trial. During the entire length of all selected sessions, including both the behavioral trials and the inter-trial intervals, we counted for each monkey the number of occurrences of continuous periods of time where one of the 5 phase patterns or an unclassified pattern was observed. The percentage of time points of the sampled LFP identified as each of the pattern classes is provided in Table 1. In addition, as a more conservative measure that takes into account potentially spurious patterns that were detected for very brief instances only, the number
<table><caption>**Table 1.** Percentage of time points classified as a specific phase pattern in each monkey given the conservative choice of thresholds used in the analysis (pooled over all grip-first and force-first conditions).</caption><thead><tr><th rowspan="2">phase pattern</th><th colspan="3">Monkey</th></tr><tr><th>L</th><th>T</th><th>N</th></tr></thead><tbody><tr><td>planar</td><td>28.47</td><td>42.41</td><td>5.51</td></tr><tr><td>synchronized</td><td>1.81</td><td>14.98</td><td>0.01</td></tr><tr><td>random</td><td>1.28</td><td>0.05</td><td>8.93</td></tr><tr><td>circular</td><td>0.01</td><td>0.01</td><td>0.02</td></tr><tr><td>radial</td><td>1.14</td><td>0.46</td><td>1.47</td></tr><tr><td>unclassified</td><td>67.28</td><td>41.65</td><td>84.06</td></tr></tbody></table>

of epochs of contiguous time points classified as the same pattern and lasting for at least 5 ms is displayed for the grip-first condition in Fig. 4A (for the force-first condition, see Supplemental Figure S3A). These results show that all pattern types were observed in each monkey, with planar wave patterns being among the most prominent and circular patterns among the least observed patterns. Only in monkey N, the random pattern was observed more often than the planar wave pattern. In addition, monkey N rarely exhibited a synchronized pattern as compared to monkeys L and T. In general, the number of phase patterns left unclassified by our conservative method was high in monkey N, consistent with our visual observation of predominantly ambiguous patterns. We speculate that this is a consequence of the more lateral array placement, and differences in the behavioral performance as compared to monkeys L and T.

**Relation of beta amplitude and phase patterns to behavior.** Given the abundance of patterns in the data, we asked whether there is a relationship between phase patterns and behavioral epochs of a trial. Thus, we investigated whether or not the occurrence of a specific phase pattern is linked to one or more behavioral events. We determined trial by trial and for each pattern separately its precise occurrence during the time course of the behavioral trial. We pooled the data from each monkey across all trials of the same condition (correct trials only), i.e. grip-first or force-first condition, and across all selected recording sessions, to obtain a measure for the probability of the occurrence of a specific pattern in time. In the following, we discuss in detail data from the grip-first condition, but qualitatively similar results are seen in the force-first condition (see Supplemental Figure S3). Figure 4B shows that monkey N had comparatively low numbers of planar and synchronized patterns during the trial, but a higher number of random patterns than the two other monkeys. This suggests that many of the planar and synchronized patterns of monkey N observed during the complete recording (Fig. 4A) occurred during the inter-trial intervals, and not during the trial (Fig. 4B).

In the next step, we assessed similarities in the temporal profile of the pattern occurrence probabilities during the behavioral trial (Fig. 4B). For each monkey, the probability of observing any pattern was strongly modulated over the time course of the trial. Common to all monkeys was the finding that planar patterns occurred mostly during the initial cue presentation and during reward, and were less prominent during movement. Synchronized patterns expressed a similar time course for monkey L and to a lesser degree for monkey T. Monkey N showed almost no synchronized pattern during the trial. In contrast, in all monkeys random patterns occurred predominantly towards the end of the delay period and during movement. Circular and radial patterns were rarely observed during the trial, but exhibited a clear modulation structure in time, albeit in a different way for each monkey.

The specific and consistent temporal modulation of the occurrence probability suggests that the spatio-temporal structure of activity is related to motor cortical processing performed during the trial. We thus asked, if also particularities of the trial condition were reflected in the probability. To this end, we compared results obtained during SG and PG conditions (Fig. 4B, black and gray, respectively). In general, the modulations of probability for both trial types were similar, but expressed a few notable exceptions. For planar waves, SG and PG deviated slightly, but significantly (indicated by dots at the bottom of each panel), during early delay (probability of observing a pattern during PG trials exceeded that of SG trials, PG < SG) and before reward (PG > SG) for monkey L, during late delay (PG > SG) in monkey T, and during cue presentation for monkey N (PG < SG). Similar, even more pronounced observations were made for synchronized patterns of monkey L and T. In addition, a tendency for a reversed effect was observed for random patterns in particular during the delay period of monkeys L and T.

Up to now, we concentrated on the time-resolved spatial organization of oscillatory activity on the basis of the phase information extracted from the time series. We next asked how these findings relate to the trial-averaged beta power, because we noticed that the temporal evolution of the occurrence probability of planar and synchronized phase patterns was reminiscent of the evolution of the beta power (Fig. 1E). To further investigate this observation, we calculated the trial-averaged beta amplitude profiles $a(t)$, i.e., the time-resolved instantaneous amplitude, or envelope, $A_{xy}(t)$, of the beta signal pooled across all electrodes (x, y), as a representative of the average instantaneous power of the beta oscillation. Again, data were calculated for all sessions used in Fig. 4B and separately for SG and PG trials (Fig. 4C). Interestingly, for all 3 monkeys the time-resolved beta amplitude profiles closely followed the occurrence probability of the planar phase pattern (Fig. 4B, top). For monkeys L and T, also the time course of synchronized patterns loosely followed that of the beta amplitude profiles. In particular, we noticed that all differences between SG and PG trials identified in the pattern occurrence probabilities were reflected in the beta amplitude profile. For example, in monkey L, the beta profiles obtained in SG and PG trials differed during the early delay (PG > SG) and before reward (PG < SG), mirrored in the occurrence
A

monkey L [grip-first]

<figure><img src="image_6.png" /></figure>

monkey T [grip-first]

<figure><img src="image_7.png" /></figure>

monkey N [grip-first]

<figure><img src="image_8.png" /></figure>

B

frequency of observing pattern (%)

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

C

time

**Figure 4.** Behavioral correlates and relation to average beta amplitude for the grip-first condition. (A) Number of epochs of a phase pattern detected in at least 5 consecutive time frames, i.e., 5 ms (bars from top to bottom: unclassified, planar wave, synchronized, random, circular, radial pattern) for monkey L (left), T (middle), and N (right). Data were obtained from all selected recording sessions including inter-trial intervals. (B) Time-resolved probability of observing a specific phase pattern (rows) during the trial. Statistics were computed across all grip-first trials of all recording sessions for each monkey (N = 15) and smoothed with a box-car kernel of length *l* = 100 ms. Trials were separated into side-grip (SG) trials (black) and precision-grip (PG) trials (gray). For monkey N, only very few synchronized patterns were detected during the trial. Color shading between curves and colored bars indicate time periods where SG and PG curves differ significantly (Fisher's exact test under the null hypothesis that, for any time point, the likelihood to observe a given phase pattern is independent of the trial type, *p* < 0.05). (C) Beta amplitude profile (envelope) pooled across all SG (black) and PG (gray) trials (same data as in panel B). The amplitude profile *a*(*t*) of a single trial is calculated as the time-resolved instantaneous amplitude *A*<sub>xy</sub>(*t*) of the beta-filtered LFP averaged across all electrodes (*x*, *y*), and measures the instantaneous power of the beta oscillation in that trial. Gray shading between curves and horizontal bars indicate time periods where SG and PG curves differ significantly (t-test under the null hypothesis that the distributions of electrode-averaged single trial amplitudes *A*<sub>xy</sub>(*t*) at each time point *t* are identical for SG and PG trials, respectively, *p* < 0.05).

of planar waves, and after movement onset (PG > SG), mirrored in the occurrence probability of synchronized
patterns. Similar observations were made for the other monkeys, and in the time period after the GO signal for
the force-first condition (Supplemental Figure S3).
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Relation of planar waves to reaction time. (A) Normalized histogram (gray scale) of the number of trials as a function of the trial's RT and the number of time bins (1 ms) classified as planar pattern during the 400 ms window before the GO signal. The dashed line shows the regression, *R* denotes the corresponding Pearson correlation coefficient (asterisk: significant at *p* < 0.05). The light and dark gray vertical lines in each plot mark the fastest and slowest 25% of the RTs averaged across sessions, respectively. (B) Time-resolved probability of observing a planar phase pattern during the trial for fast (RT within the 25% lowest RTs within a given session) and slow RTs (RT within the 25% highest RTs), denoted by light and dark gray curves, respectively. Curves were smoothed with a box-car kernel of length *l* = 100 ms. Color shading between curves and colored bars indicates time periods where fast and slow RT curves different significantly (Fisher's exact test, *p* < 0.05, see Fig. 4). Statistics in A and B were computed across all grip-first and force-first trials of all recording sessions (*N* = 15) for each monkey (columns).</figcaption></figure>

In addition to the dependence of the time course of the beta amplitude during movement preparation on the characteristics of the anticipated movement, it has been demonstrated that in instructed delay reaching tasks the movement phase is preceded by a sharp attenuation of beta amplitude⁴⁹. The attenuation correlates with the onset of the neuronal movement-related rate response and exhibits a spatially organized dynamics across the motor cortical surface⁵⁰. Moreover, it has been shown that for sensorimotor areas, reaction time (RT) in a visuomotor discrimination task is positively correlated with the magnitude of the pre-stimulus beta amplitude, hinting at beta activity as a mechanism for motor maintenance that prevents an early movement initiation⁵¹. Therefore, we speculated that also the observation of planar wave activity before the GO cue would be a signature of such a mechanism that is indicative of RT. In Fig. 5A we visualize the relationship of RT and the number of time points at which a planar wave was automatically detected during the 400 ms window before stimulus onset. Indeed, we found a small, but significant positive correlation in all three monkeys and for both, grip-first and force-first conditions. To better understand the temporal dynamics of this effect, we plotted the time course of observing a planar wave before the GO cue independently for the trials within the quartile with slowest RTs versus trials within the quartile with highest RTs (Fig. 5B). Here, quartiles were evaluated independently for each monkey and on a session-by-session basis to account for the strong discrepancies in the monkeys' average RTs. In the grip-first condition, for all three monkeys we observed that fast RTs go along with a lower probability of observing a planar wave pattern in the period before the GO cue as compared to slow RTs. In contrast, in the force-first condition, the absence of grip-type information impeded the efficient preparation of the upcoming movement, leading to notably slower RTs in monkeys L and T (cf., Fig. 5A). In consequence, in this condition the observation of planar waves before the GO cue differed little between slow and fast RT trials, in particular around the time of stimulus onset. In summary, complementing previous findings³¹, the fast initiation of a prepared movement is in competition with a mechanism that exhibits strong beta activity organized as planar wave-like propagation across the motor cortical surface.

**Quantification of phase patterns.** Having identified the different classes of phase patterns and their behavioral correlates, in a next step we aimed to understand their properties in single trials. Therefore, we quantified features extracted from the classification results. As classification was performed individually on single time points, we first calculated the durations of epochs of consecutive time points being classified as the same pattern. In Fig. 6A we show the resulting distributions of the durations for each of the pattern types. Naturally, these statistics depended on how conservative the choice of thresholds for pattern detection was set. However, given that thresholds were set in accordance to the observed phase pattern (cf. Methods), they served as a visually inspired characterization of the observed duration of a pattern. We found that, on average, planar, synchronized,
<figure><img src="image_13.png" /><figcaption>**Figure 6.** Statistics of detected patterns. (A) Histogram of durations of epochs of consecutive time points classified as belonging to the same phase pattern (cf., Fig. 4A). (B–D) Distributions of the direction *d*(*t*) (panel B), phase velocity *v*(*t*) (panel C), and amplitude profile *a*(*t*) (panel D), as a function of the detected phase pattern. Data are separated (columns) according to monkey and recording condition (grip-first vs. force-first). Histograms for different phase patterns are plotted overlapping in the color corresponding to the legend on the right. For each phase pattern, a histogram entry in panels (B–D) represents the measured quantity averaged across the array calculated at a time point classified as that pattern. In panel B, the average direction of the phase gradient is plotted in brain coordinates by rotating the activity, and mirroring data along the medio-lateral axis for monkey T to compensate for the array placement in the opposite hemisphere as compared to L and N. Dots and dashed lines below indicate the medians of the corresponding distributions. (E) Joint representation of the medians of the distributions shown in panels A, C, and D. Each data point represents the median of one monkey in one recording condition for one pattern class (indicated by color).</figcaption></figure>

and radial patterns all had longer durations than random and circular patterns (see large dots in Fig. 6A), all on the order of less than one cycle of the beta oscillation (≈ 40–50 ms).

In a next step we examined the preferred direction of the phase gradients of the wave patterns. Here we only considered planar phase patterns (Fig. 6B) for which the measure was equivalent to the direction of movement of the planar wave front. Planar waves in monkeys L and T were preferentially observed in the anterior-medial to posterior-lateral direction (see inset for cortical space), whereas waves in monkey N were observed in the anterior-lateral to posterior-medial direction. Noting that the array location in monkey N differed from that in monkeys L and T (cf. Fig. 1B,C), the observations from all three monkeys are compatible with those described earlier²². Even though it was possible to calculate the direction of the phase gradients of any phase pattern, we refrained from showing the distribution of directions for patterns other than planar patterns since their characteristics do not allow a clear interpretation of wave propagation direction.

To investigate the dynamical aspect of wave propagation, we calculated the average wave velocity *v*(*t*) (cf. Methods) at each time point (Fig. 6C). For planar wave patterns, this was directly interpretable as the propagation velocity of the observed wave front. The median propagation velocities of the planar waves were *v*(*t*) = 29.1 ± 10.3 cm/s (grip-first) and *v*(*t*) = 29.1 ± 10.4 cm/s (force-first) for monkey L, *v*(*t*) = 40.5 ± 16.2 cm/s (grip-first) and *v*(*t*) = 40.3 ± 19.5 cm/s (force-first) for monkey T, and *v*(*t*) = 14.2 ± 4.6 cm/s (grip-first) for monkey N (all values: median ± median absolute deviation). These values are in rough agreement with those previously reported²². For the other wave patterns, even though it was possible to calculate a velocity, it may not be directly interpreted as the velocity of a propagating planar wave front since phase gradients do not align across the array. Instead, it is a measure that captures the average velocity calculated from the local velocities across the array. Synchronized patterns could be considered as a special case of planar waves with a very large spatial wavelength, and as a
consequence they exhibited high (in theory, infinitely high) velocities (Fig. 6C). On the other hand, random and circular patterns were characterized by phase values that differed strongly between adjacent recording sites. Therefore, average velocities derived from the phase gradients were low. Finally, radial patterns resembled the planar patterns in that they could be approximated by a planar wave front at a large distance from the center of the radial pattern. In agreement with this interpretation, they exhibited similar phase velocities as observed for the planar pattern. Thus, the 5 distinct phase patterns showed clear differences in the distribution of their velocities, where a low $v(t)$ corresponds to random or circular patterns, a medium $v(t)$ relates to planar or radial patterns, and a high $v(t)$ indicates the presence of a synchronized pattern. In this sense, the value of the velocity represents a reliable proxy for the type of the observed pattern.

After having quantified the features of the wave patterns in single trials, we come back to the question we originally set out to answer, namely how the beta amplitude relates to the occurrence of a wave pattern. In Fig. 6D we show for each monkey and behavioral condition the distribution of the beta amplitude profile (envelope) $a(t)$ during each of the 5 phase patterns. We observed a clear relationship between the instantaneous magnitude of the beta oscillation and the spatial phase pattern. Circular and random patterns occurred at small amplitudes, planar and radial patterns at intermediate amplitudes, and only synchronized patterns occurred at high amplitudes.

The results shown in Fig. 6A,C and D were summarized in Fig. 6E, where for each phase pattern, monkey and behavioral condition the averaged data for duration, velocity and amplitude are plotted against each other. This representation clearly shows a clustering of collective data points for each individual phase pattern. Thus, the 5 phase pattern classes are described by a specific combination of characteristic values for pattern duration, velocity, and amplitude.

**Beta amplitude determines phase pattern.** In the last step of our analysis, we now ask if the relationship between amplitude and spatial organization holds for any time point independent of whether or not it can be unambiguously attributed to any of the idealized classes of phase patterns. In order to obtain such a time-resolved view of how the amplitude (which by itself did not exhibit a strong spatial organization, Fig. 2D) correlates with the temporal evolution of the patterns, we employed the phase velocity $v(t)$ as a proxy to quantify the spatial organization that can be readily calculated for each individual time point (as opposed to pattern duration). In Fig. 7 we show the correlation between the instantaneous beta amplitude profile (envelope) $a(t)$ and the phase velocity $v(t)$ for each time point for all three monkeys, independent of the phase pattern classification, thus including instances during which no pattern could be classified by our conservative classification algorithm. We observed that the two variables were strongly positively correlated ($R > 0.8$ for all monkeys) and correlations were highly significant ($p \ll 0.001$). Thus, an increase in amplitude goes along with an increase in phase velocity (cf. also Fig. 6C). As shown above, the velocity $v(t)$ is indeed a good correlate of the perceived organization of beta activity on the electrode array. To more directly illustrate how the velocity measure relates to the previously defined classes of phase patterns, we indicate in Fig. 7 by ellipses the regions of the histograms where the individual classes of phase patterns were predominantly found (see Supplemental Figure S4 for an alternative visualization).

In conclusion, we find that at any point in time, the amplitude of the beta oscillation at one single recording site of the Utah array is correlated to the spatial organization of activity across the array, here parametrized by the velocity. In this sense, a given velocity may in some cases relate directly to one of the five classes of phase patterns (Supplemental Figure S4). In other cases, the velocity indicates a pattern that was unclassified by the rigorous automatic classification due to its ambiguity or the level of noise, but nonetheless typically resembles distinct qualitative features of one of the pattern classes, as illustrated by a gallery of randomly chosen unclassified patterns grouped by their velocity (Supplemental Figure S5).

## Discussion

Two main objectives guided this work. First, we aimed to obtain a more complete description of the wave-like spatio-temporal phase patterns exhibited in the beta range of the LFP signals in monkey motor cortex during a complex delayed motor task, and thereby extend reports that only included descriptions of planar wave propagation$^{22,42}$. Second, we asked in how far these patterns, determined solely by the phase of the oscillation, are related to the instantaneous modulation of the beta amplitude. Given that beta amplitude is modulated during the behavioral trial, we also asked how the phase patterns are related to behavioral epochs in order to determine their possible functional implications.

**Motor cortical beta oscillations exhibit a variety of spatio-temporal patterns.** By analyzing the dynamics of LFP activity across multi-electrode arrays, we demonstrated that beta oscillatory activity shows a number of salient types of spatio-temporal patterns in addition to traveling planar waves$^{22}$, namely quasi-synchronized, random, radial, and circular patterns. Such additional types of patterns have previously been predicted from theoretical considerations$^{38}$, and were observed in experimental work, e.g., in slow delta activity of anesthetized marmoset monkeys$^{40}$. We developed a phenomenological classification method to identify epochs that unambiguously exhibit one of the 5 pattern classes. Our approach detected those in a conservative manner in order to capture the qualitatively salient patterns that are also identified by a human observer. Indeed, the algorithm tends to leave a large number of time points unclassified, due to the difficulty to clearly attribute a pattern to one of the 5 idealized pattern types. The reason for this is two-fold: On the one side, the coarse-grained resolution of the Utah array provided only rough estimates of the phase gradients. On the other side, the patterns were often ambiguous, in particular at time points of dynamical transitions between patterns. Planar wave fronts were often not completely planar, but showed a slight curvature, a feature shared with radial or circular patterns. Furthermore, radial and circular patterns that were not necessarily centered on the array were difficult to detect.
<figure><img src="image_14.png" /><figcaption>**Figure 7.** Correlation of instantaneous beta amplitude and spatial pattern of phases. Upper row: 2-D histograms of phase velocity *v*(*t*) and beta amplitude profile *a*(*t*) evaluated at each time point (independent of the detected phase pattern) shown for each monkey and behavioral condition (columns). Gray values indicate density of time points falling in each histogram bin, normalized to the largest entry of the histogram. The values of the Pearson correlation coefficients *R* are given in the bottom right of each panel (asterisk: significant at *p* < 0.05). Each ellipse represents the distribution of time points classified as a specific phase pattern (indicated by color). Center of ellipses: mean; radii of ellipses are given as 2 standard deviations in the direction of the 2 principle components. Lower row: zoomed-in versions of the upper histograms. Left column, inset in lower panel and dots in red to yellow shades: Illustration of spindle dynamics by example of the spindle before CUE-ON presented in Fig. 3A. Inset reproduces this spindle (transition from red to yellow colors indicates increasing time), corresponding detected states are shown as bars above the spindle. For each time point of the spindle, the corresponding values of the amplitude and phase velocity are marked in the histograms using the identical color. Average ellipse centers: (1.4 ± 0.8, 41.5 ± 27.3 cm/s) for planar; (1.9 ± 0.8, 88.1 ± 29.4 cm/s) for synchronized; (0.7 ± 0.2, 6.7 ± 1.0 cm/s) for random; (0.7 ± 0.3, 8.8 ± 1.7 cm/s) for circular; (1.2 ± 0.7, 26.8 ± 21.5 cm/s) for radial.</figcaption></figure>

Also, random states often exhibited a slight degree of correlation between activities recorded on neighboring electrodes, contradicting the *a priori* assumption of pure independence. Nevertheless, this approach of detecting patterns yielded reliable results in terms of their statistics (Fig. 6).

To overcome the limitation that the phenomenological classification method only detected unambiguous, ideal phase patterns, we tested the potential of the phase velocity as an easily accessible continuous measure to quantify the spatial arrangement of phases for time points where none of the ideal pattern categories matched the observation. Due to the fact that the velocity vector is tightly coupled to the arrangement of phase gradients across the array (see Fig. 6C,E), we could indeed link the distributions of velocities to the 5 specific phase patterns (see Fig. 6E). Thus, using the continuous measure of the phase velocity, we were able to gain a complete picture of the time course of pattern progression.

The instability of pattern types may suggest that some of the salient pattern types indeed underlie identical dynamical processes, and form a continuum: radial patterns may appear nearly planar wave-like some distance from the array center, and quasi-synchronized states appear at the limit of planar waves approaching infinite phase velocity. This similarity of phase patterns was also reflected in the statistics (Fig. 6) describing the occurrence of the patterns, e.g., the similar duration of radial and planar patterns, or the comparable distributions of velocity for planar and radial, as well as circular and random patterns. To investigate this issue in detail, recordings on a larger spatial scale and with a higher spatial resolution would be required.

**Specific phase patterns occur at different times during movement preparation and execution.**

The probability of detecting a specific phase pattern was variable during the trial of our reach-to-grasp task. Planar and synchronized patterns occurred more often during the pre-cue epoch and during the delay whereas random patterns were more likely to occur around movement execution (Fig. 4). This observation is in line with the hypothesis that planar and synchronized patterns could be triggered by the arrival of visual information in motor cortex from adjacent cortical areas not covered by our Utah array⁴². In agreement with this view, the orientation of planar wave propagation in our data is in agreement with previous studies²². More precisely, we found orientation preferentially aligned along the antero-posterior axis. The direction of planar wave propagation was more
anterior-medial to posterior-lateral in monkeys L and T, whereas in monkey N it pointed from anterior-lateral to posterior-medial. This difference could reflect the fact that the array was implanted more medially in monkeys L and T than in monkey N (Fig. 1). Therefore, it seems that planar waves travel toward a medial point along the central sulcus, probably at the level of the hand and finger representation (“nested organization”⁵²; “horseshoe” structure⁵³). This directional preference may be structured by the underlying connectivity of this cortical area⁵².

The predominance of random patterns during movement execution suggests that the spatio-temporal dynamics of neuronal activity is strongly altered during this epoch. The spatio-temporal structure of these patterns characterized by their focal origin and short-range propagation could reflect that during movement, information processing is more local and activity propagation is spatially constrained to motor cortex. However, this hypothesis can hardly be tested at the restricted spatial scale of a single Utah array. Multiple Utah arrays or optical imaging techniques are required to measure the neuronal dynamics at the mesoscopic scale (e.g., visual cortex³⁰).

**Wave dynamics relate to the modulation in beta amplitude.** Beta amplitude is known to be strongly modulated by the task epoch¹⁶. Interestingly, Fig. 4 suggests that across trials, the probability of observing different phase patterns also follows the trial-averaged amplitude profile of the beta oscillation. Namely, planar and synchronized waves are present during epochs of large beta amplitudes, whereas random waves are prominent during epochs of small amplitudes. The relation of circular and radial patterns to the beta amplitude is more ambiguous. These observations would support the hypothesis that the wave dynamics is closely linked to the processes underlying the modulation of the amplitude of beta oscillations.

Indeed, even on the single trial level, Fig. 6 suggests that low beta amplitudes are linked to random or circular phase patterns with low velocities, intermediate amplitudes to planar or radial phase patterns with intermediate velocities, and that the highest amplitudes indeed co-occurred with quasi-synchronous phase patterns expressing by far the highest velocities (see especially Fig. 6E). While for small beta amplitudes it cannot be excluded that the random patterns arise in part through difficulties in measuring the beta phase, it remains striking that at medium to high beta amplitudes the observed activity is always spatially structured. The pattern statistics also show that the epochs during which one particular, clearly structured pattern was observed were typically of very short duration, in the order of 1 or 2 oscillation cycles (see Fig. 6A). This is reminiscent of the short-lasting high amplitude events of a few cycles of beta oscillations described by others, the so-called spindles¹⁹,⁴¹. Indeed, these observations point to a tight relationship between spindle dynamics describing the amplitude modulation of the LFP, and the occurrence of wave-like activity, as shown by the correlations in Fig. 7. In all monkeys, we observed that with growing amplitude wave propagation tended to accelerate. For high amplitude beta signals, the phase pattern accelerated to such high levels that the observed pattern became synchronous, which in the ideal case would exhibit infinite velocity.

To illustrate how this observation relates to the dynamics of a single spindle, we visualized the temporal evolution of one example spindle and its pattern classification in the left column of Fig. 7 (inset) and observed the corresponding smooth trajectory in the space of amplitude and phase velocity (yellow-red trace). In this spindle, a synchronized state was detected at the spindle peak, flanked by planar patterns before and after the peak. By observing the trajectory of the phase velocity, we observe that the modulation of spindle activity goes along with wave-like activity that progressively increases in speed as the LFP beta amplitude increases. Thus, this strong correlation between amplitude and velocity suggests that at the spindle peak also the velocity peaks, which, for large spindles, corresponds to large spatial wavelengths of the phase pattern that are perceived as synchronized states on the spatial scale of a Utah array. In contrast, a low LFP beta amplitude, as observed between spindles, goes along with random or circular patterns at low velocity. A dynamic representation of how the pattern velocity follows beta amplitude and the evolution of spindles can be seen in the middle panels of Supplemental Movie S1.

In summary, we speculate that the formation of a structured, directed pattern, its acceleration to a near-synchronized appearance, followed by deceleration, and finally its breakup in a random or circular pattern marks the temporal organization of the formation of a beta spindle, its peak, and its decay, respectively. Supporting this view, it has been shown that the maxima of LFP spindles tend to synchronize across large distances, even between cortical areas and hemispheres¹⁹,⁴¹, as expected for emergence of synchronized patterns. These combined observations are in line with the highly dynamic nature of pattern occurrences.

A model of brain processing that would be intrinsically affected from such a dynamic scaffold is the concept of communication through coherence¹⁰,¹¹. In this framework, the coherence and phase relationship between oscillations on different electrodes were taken as a measure of the ability of neurons to communicate, i.e., that information is best transmitted when the two communicating sites exhibit an optimal phase lag. This concept seems evident when considering, e.g., communication between two brain areas that exhibit distinct population oscillations. It is, however, unclear what this model implies on the mesoscopic scale, such as the course-grained recordings from a Utah array presented here, where the overall pattern of these phase lags between electrodes continuously changes in time. Nevertheless, we may hypothesize that if activities on different electrodes become increasingly synchronized with a decreasing phase lag as spindles increase their amplitude, this would lead to a state where information can be more easily communicated across the complete array, although possibly with less specificity. This would indicate that amplitude modulations, and in particular spindles, act as a time window for enabling cortical communication across larger distances: not just by means of the strength of synchronization within the local population of neurons (as indicated by the increased beta amplitude), but because this goes along with a wide-spread zero-lag synchronization of the oscillatory activity, i.e., synchronized patterns¹⁹,⁴¹. This assumption is highly consistent with the above-described results showing that synchronized and planar patterns are more frequent during the delay epoch and could reflect the transmission of information between distant cortical sites. Conversely, the random patterns occur more often during periods of low beta amplitude and could be linked to the local processing of information. Radial and circular patterns occupy an intermediate position along this continuum and have an unclear relationship to behavior.
This line of arguments raises the question how the patterns of phase dynamics are related to synchronization on the level of single neuron spiking activity. In our data, the probability of neuronal firing is modulated in time by the beta rhythm⁵⁴. With few exceptions, the preferred phases of firing of recorded single units across the three monkeys are found near the trough of the LFP cycle (cf., Supplemental Figure S6). This preference is independent of the location of the units on the array. Assuming that this phase preference holds for the bulk of neurons in the vicinity of each electrode, it would lead to modulations of the population rate that closely match the temporal modulation of the beta phase, i.e. a higher proportion of neurons would spike during LFP troughs than during LFP peaks. Therefore, also the spatial organization of population rates would be expected to follow the LFP phase patterns. Supporting this view, it has been shown that the spiking activity synchronizes with the oscillatory spindle peaks¹⁹ and the cross-correlation histograms of the spiking activity of pairs of neurons become oscillatory in the beta range during periods of strong beta activity⁵⁵. The entrainment of single neuron spiking activity to the LFP oscillation increases with LFP amplitude⁵⁶. Additionally, we have shown²¹ that at moments of precise transient spike synchronization that exceeds the expectation based on firing rate⁵⁷, spikes lock more strongly to the LFP beta oscillation than expected by chance. This effect of particularly strong locking of significant spike coincidences was observed especially during high beta amplitudes. Interpreting the occurrence of excess synchrony as reflecting active cell assemblies, we embedded our findings in a theoretical model that predicts that activated cell assemblies are entrained to the LFP oscillation at a specific phase shortly preceding the trough of the oscillation⁵⁸. Combining these findings, we may speculate that the modulation of the beta amplitude as a function of the occurrence of a beta spindle is not only indicative of the spatial phase pattern of LFP beta activity, as shown in this study, but that additionally beta spindles may govern the temporal structure of spike patterning observed across the array. Indeed, findings of spike sequences⁴² or synchronous spike patterns⁵⁹ that align to the principle direction of phase gradients (Fig. 6B) support this view of a functional mechanism that underlies the generation of beta phase patterns.

A recent discussion³⁰ of the functional implications of wave propagation is highly related to such an hypothesis, despite the qualitative differences in their description of waves in superficial layers of visual cortex. The authors show that a network of excitatory and inhibitory neurons operating in the balanced regime and connected by a horizontal fiber network captures the essential features of the observed wave dynamics. Thus, the authors speculate that the transient depolarization caused by the wave passing at a certain position creates a time window of increased sensitivity, i.e., spike probability, of neurons at that location. This would ensure an optimal integration of information as long as the incoming input is timed to the arrival of the wave. Translating this idea to our scenario, the continuous traveling waves we observed could play a similar role when perceived as reverberating waves of the single-cycle propagation³⁰. Moreover, we propose that the synchronized patterns, and thus epochs of large beta amplitudes, correspond to states where the optimal time window for the integration of incoming inputs is no longer spatially modulated by the propagating wave dynamics, but only by the anatomical structure of the network.

Speculating on these observations and discussions, we put forth the hypothesis that a future joint analysis of LFP phase patterns and spatio-temporal precise spike patterns would exhibit a correspondence in their spatial orientation, such that synchronous spike patterns are observed in neurons perpendicular to the instantaneous phase gradients, and spike sequences are observed along the direction of the gradient with inter-spike delays matching the instantaneous phase difference. Moreover, we suggest that significant spike synchrony would be observed for neurons separated by large spatial distances more often during high beta amplitudes, when synchronized phase patterns dominate, than during low amplitudes.

In summary, despite the fact that motor cortical beta oscillations show a strong correlation between signals recorded over large distances, the phase relationships are highly correlated to the amplitude modulation of beta activity, which in turn has been related to the dynamics of spike synchronization, and to behavior. Thus, we believe that the investigation of amplitude and phase patterns provides a novel leverage on understanding the coordination of activity within spiking neuronal networks.

## Methods

**Experimental Design.** Three monkeys (Macaca mulatta) were used in these experiments, two females (monkeys L and T) and one male (monkey N). All animal procedures were approved by the local ethical committee (Comité National de Réflexion Ethique sur l'Experimentation Animal - no. 71 [C2EA-71]; authorization A1/10/12) and conformed to the European and French government regulations. Monkeys were kept in colonies of 2-4 monkeys in a modular housing pen, with access to a central play area. They were not water-deprived during the experimental period. Each monkey was trained to grasp, pull and hold an object with low force (LF) or high force (HF) using either a side grip (SG) or a precision grip (PG). The task was programmed and controlled using LabView (National Instruments Corporation, Austin, TX, USA). The trial sequence was as follows. The monkey self-initiated each trial by pressing a switch with the hand (TS). After a start period of 400 ms a warning signal (WS) lighted up to focus the attention of the monkey. After another 400 ms, the cue (CUE-ON until CUE-OFF) informed the monkey either about the grip type (grip-first condition) or the force (force-first condition) required in this trial. The duration of cue presentation was 300 ms. It was followed by a preparatory delay period of 1 s. The subsequent GO signal completed the information about force and grip, respectively, and in parallel asked the monkey to perform the movement by using the correct grip type and force to pull and then hold the object in a defined position window for 500 ms. Further periods: reaction time (RT) between the GO signal onset until the monkey released the switch (SR), movement time (MT) between switch release and object touch (OT), and pull time (PT) between OT and reaching the correct holding position. For correct trials the monkey was rewarded (RW) at the end of the holding period with a drop of apple sauce. See Fig. 1A for a graphic representation of the task design.
<figure><img src="image_15.png" /><figcaption>**Figure 8.** Explanation of the methodology. (A) Example phase map $\Phi_{xy}(t)$ (colors) of monkey L (session identifier 1101013-002) with overlay of phase gradient map $\Gamma_{xy}(t)$ (top left), phase direction map $\Delta_{xy}(t)$ (top right) and gradient coherence map $\Lambda_{xy}(t)$ (bottom left) shown as black arrows per electrode. White large arrows: corresponding quadrant-averaged maps (B). Sketch of the 10 × 10 electrode array (squares) to illustrate the neighborhoods $\hat{\mathbb{R}}_{xy}$ (light gray) and $\bar{\mathbb{R}}_{xy}$ (dark gray), and the calculation of the gradient continuity $C(t)$ (top right), radial-parallel alignment $R_\parallel$ (bottom left), and radial-orthogonal alignment $R_\perp$ (bottom right). $C(t)$: angle between the gradient (red arrow) at a chosen electrode (red square) and the gradient at the square that the red gradient points to (yellow arrow and square). $R_\parallel$ and $R_\perp$: angle between the gradient (red arrow) at a chosen electrode (red square) and the corresponding solid gray line (parallel to and perpendicular to a line through the array center and the electrode, respectively). (C) Calculated values for the 6 measures (represented by horizontal line markers; from left to right: $\sigma_p(t)$, $\sigma_g(t)$, $\mu_c(t)$, $C(t)$, $R_\parallel$, and $R_\perp$) for artificial data modeling ideal realizations</figcaption></figure>
of the 5 phase patterns (planar wave, synchronized, random, circular, and radial). A small amount of noise was added to the synchronized phase pattern to avoid division by zero in calculating measures. Black markers: measure was used to detect the corresponding phase pattern. Dark areas indicate the range of valid values for that measure required for a positive classification according to the threshold criteria (cf., Table 2). Light gray markers: measure was not used to detect the corresponding pattern. Region of valid values for measure $\sigma_p(t)$ for synchronized patterns is close to 0.

**Neuronal Recordings.** After completed training, a 100-electrode Utah array (Blackrock Microsystems, Salt Lake City, UT, USA) was chronically implanted in M1 and PM, contralateral to the working hand (for location see Fig. 1C) and overlapping rostral M1 and the posterior end of the dorsal premotor cortex (PMd) in monkeys L and T. The array of monkey N was placed more laterally covering the most medial part of the ventral premotor cortex (PMv). The 4x4 mm silicon based array consisted of 10-by-10 Iridium-Oxide electrodes, of which 96 were available for recording. The length of each electrode was 1.5 mm, with a 400 µm inter-electrode spacing. With this electrode length, we assume that the array enabled recording between the deep cortical layer III until the most superficial part of layer V. The distance between any pair of electrodes can be easily determined from the fixed geometric structure of the array. The surgery for array implantation was described in Riehle *et al.* (2013) and is briefly summarized below. The surgery was performed under deep general anesthesia using full aseptic procedures. A 30 mm × 20 mm craniotomy was performed over the motor cortex and the dura was incised and reflected. The array was inserted into the motor cortex between the central and arcuate sulci (Fig. 1C) using a pneumatic inserter (Blackrock Microsystems). It was then covered by a non-absorbable artificial dura (Preclude, Gore-tex). Ground and reference wires were inserted into the subdural space. The dura was then sutured back and covered with a piece of artificial absorbable dura (Seamdura, Codman). The bone flap was put back at its original position and secured to the skull by a titanium strip and titanium bone screws (Codman). The array connector was fixed to the skull on the hemisphere opposite to the implant. The skin was sutured back over the bone flap and around the connector. The monkey received a full course of antibiotics and analgesics after the surgery and recovered for one week before the first recordings.

Neuronal data were recorded using the 128-channel Cerebus acquisition system (NSP, Blackrock Microsystems). The signal from each active electrode (96 out of the 100 electrodes were connected) was pre-processed by a head stage (monkey L and T: CerePort plug to Samtec adaptor, monkey N: Patient cable, both Blackrock Microsystems) with unity gain and then amplified with a gain of 5000 using the Front End Amplifier (Blackrock Microsystems). The raw signal was obtained with 30 kHz time resolution in a range of 0.3 Hz to 7.5 kHz. From this raw signal, two filter settings allowed us to obtain on-line two different signals by using filters in two different frequency bands, the local field potential (LFP, low-pass filter at 250 Hz) and spiking activity (high-pass filter at 250 Hz). Here, only LFPs were analyzed, which were down-sampled at 1 kHz. An in-depth description of the experimental design and the recording procedure is available⁴⁵.

**Power spectra.** Power spectra were calculated using Welch's average periodogram algorithm using the psd function of the Python package scipy. We used windows of length *l* = 1024 sample points (at 1 kHz sampling). Each window was tapered using a Hanning window. The time-resolved power spectra (spectrograms) were calculated using windows of length *l* = 512 samples and an overlap of 500 samples.

**Definition of maps and vector fields.** We calculated 5 different types of maps in order to visualize the spatial arrangement of oscillatory activity in the beta range on the array, and to provide a starting point for calculating multiple measures that characterize the arrangement. In a first step, we filtered the LFP signal on each electrode using a third-order Butterworth filter (pass band: 13–30 Hz) in a way that preserved the phase information (filtfilt() function of the Python package scipy). The filter setting was intentionally chosen broad such that it enabled us to identify the phase and the amplitude despite temporal variations of the beta oscillation amplitude and its center frequency. In order to compare the relative changes in amplitude between different electrodes, the amplitude of the LFP signal was then normalized across recording electrodes by computing the z-transform of the complete filtered LFP signal on an electrode-by-electrode basis.

In a next step, we calculated the instantaneous amplitude and phase of the normalized, filtered LFP time series $x_i(t)$ on each electrode *i* by first constructing the analytic signal $X_i(t) = x_i(t) + j \mathcal{H}[x_i(t)]$, where $\mathcal{H}(\cdot)$ represents the Hilbert transform and $j^2 = -1$. From $X_i(t)$, we obtained the instantaneous signal amplitude $a_i(t)$ by taking its modulo, and the instantaneous phase $\phi_i(t)$ by taking its argument (angle). We defined the maps $A_{xy}(t) = a_i(t)$ and $\Phi_{xy}(t) = \phi_i(t)$ by the instantaneous phases of the LFP, where $x \in \{0, \dots, 9\}$ and $y \in \{0, \dots, 9\}$ are the coordinates of the recording electrode *i* of the Utah array in units of the inter-electrode distance of 400 µm.

In a further step we investigated whether, locally at each electrode and at each point in time, there is a spatially structured arrangement of the phases $\Phi_{xy}(t)$. To this end, in the remainder of this section, we defined three additional maps that we term the *phase gradient map* $\Gamma_{xy}(t)$, the *directionality map* $\Delta_{xy}(t)$, and the *gradient coherence map* $\Lambda_{xy}(t)$ (cf. Fig. 8A). The local spatial phase gradient at position electrodes (x, y) was estimated based on a neighborhood $\bar{\mathbb{R}}_{xy}$ of its *k*-nearest neighbors in the same column x or row y (see Fig. 8B for a graphical representation). For border electrodes ($x \notin \{2, \dots, 7\}$ or $y \notin \{2, \dots, 7\}$), only existing electrodes were considered as nearest neighbors. In this manuscript we chose $k=2$ to obtain a smooth map of the local phase gradients. Let $\bar{\mathbb{N}}_{xy}$ denote the cardinality of the set $\bar{\mathbb{R}}_{xy}$. We now constructed the phase gradient map as the average gradient $d|\phi(t)|/dx \cdot e^{i\alpha}$, between electrode (x, y) and each of its neighbors (i, j), where $\alpha$ denotes the angular direction between the electrode locations. The result is the map of *phase gradients*.
<table><caption>**Table 2.** Thresholds used to classify phase patterns. A time point was classified as one of the 5 listed phase patterns if all corresponding threshold criteria were met. Tests were administered in the order given by the table, and the first match is chosen as the classification.</caption><thead><tr><td>Phase pattern</td><td>Thresholds</td></tr></thead><tbody><tr><td>planar wave</td><td>σ<sub>g</sub>(t) &lt; θ<sub>3</sub></td></tr><tr><td>radial</td><td>|R| &gt; θ<sub>8</sub></td></tr><tr><td>synchronized</td><td>σ<sub>p</sub>(t) &lt; θ<sub>1</sub>, and σ<sub>g</sub>(t) ≥ θ<sub>4</sub></td></tr><tr><td>circular</td><td>σ<sub>p</sub>(t) ≥ θ<sub>2</sub>, σ<sub>g</sub>(t) ≥ θ<sub>4</sub>, C(t) ≥ θ<sub>6</sub> and |R| ≥ θ<sub>7</sub></td></tr><tr><td>random</td><td>σ<sub>p</sub>(t) ≥ θ<sub>2</sub>, σ<sub>g</sub>(t) ≥ θ<sub>4</sub>, and μ<sub>c</sub>(t) ≤ θ<sub>5</sub></td></tr></tbody></table>

$$ \Gamma_{xy}(t) = \bar{N}_{xy}^{-1} \sum_{(i,j) \in \mathbb{N}_{xy}} \frac{\Phi_{ij}(t) - \Phi_{xy}(t)}{\sqrt{(i-x)^2 + (j-y)^2}} \cdot e^{j\alpha_{ij}} \approx \nabla \Phi_{xy}(t). \quad (1) $$

Based on the average frequency of the beta oscillation $f_β$, we can easily derive the phase velocity field $\Psi_{xy}(t) = 2\pi f_β |\Gamma_{xy}(t)|^{-1}$, which indicates the phase velocity of a planar wave front running through the point $(x, y)$. Here, $f_β = 21.5$ Hz was chosen as the mean frequency of the respective beta bands of the monkeys (see above). An estimate of the macroscopic phase velocity can be obtained by calculating the average $v(t) = |\Psi_{xy}(t)|$ across electrodes. Next, we defined the phase directionality map

$$ \Delta_{xy}(t) = |\Gamma_{xy}(t)|^{-1} \Gamma_{xy}(t) \quad (2) $$

by normalizing the vectors of the phase gradient map $\Gamma_{xy}(t)$ to unit length. It indicates only the direction of the local phase gradient, independent of its magnitude. Finally, we defined the gradient coherence map as an average of the directionality map in a neighborhood $\hat{N}_{xy}$ of all k-nearest neighbors of cardinality $\hat{N}_{xy}$ (cf. Fig. 8B):

$$ \Lambda_{xy}(t) = \hat{N}^{-1} \sum_{i,j \in \mathbb{N}_{xy}} \Delta_{ij}(t). \quad (3) $$

It represents a second order measure of the gradient field and serves two purposes. The direction of each entry in $\Lambda_{xy}(t)$ provides a smoothed version of the vector field $\Gamma_{xy}(t)$, which is better suited for visualization due to the rather sparse sampling of activity. More importantly, the magnitude of the vectors in $\Lambda_{xy}(t)$ indicate whether, locally, phase gradients point in the same direction (independent of the magnitudes of the gradients).

**Quantification of observed phase patterns.** Based on the phase map $\Phi_{xy}(t)$ and the three vector fields $\Gamma_{xy}(t)$, $\Delta_{xy}(t)$, and $\Lambda_{xy}(t)$, we now introduced 6 measures that quantitatively describe the spatial arrangement of phases on the array at each time point. These measures will later on serve as a basis to classify the phase pattern, i.e. the spatial arrangement of phases in $\Phi_{xy}(t)$, in an automatized manner. In the following, let $\aleph$ denote the set of all used electrodes in a given recording, and $N = |\aleph|$ its cardinality.

**Circular variance of phases.** One phase pattern commonly observed is the one where all electrodes are fully synchronized at near-zero phase lag. Therefore, we introduced the circular variance of phases

$$ \sigma_p(t) = N^{-1} \sum_{(i,j) \in \mathbb{N}} e^{j \Phi_{ij}(t)} \in [0, 1] \quad (4) $$

as a measure to determine the similarity of the phase across the array. Here, $\sigma_p(t) = 0$ indicates that an identical phase $\Phi_{xy}(t)$ observed at each electrode, whereas $\sigma_p(t) = 1$ indicates that phases are uniformly distributed across the array.

**Circular variance of phase directionality.** In order to measure the degree to which phase gradients are globally aligned across the grid, we introduced the circular variance of the phase directionality

$$ \sigma_g(t) = N^{-1} \sum_{(i,j) \in \mathbb{N}} \Delta_{ij}(t) \in [0, 1]. \quad (5) $$

A perfect planar wave is observed if $\sigma_g(t) = 0$, i.e., all phase gradients point in the same direction (independent of the magnitude of the gradients). This measure is similar to the PGD measure$^{22}$.

**Local gradient coherence.** In order to determine whether locally (within $\hat{N}_{xy}$) phase gradients point in a particular direction, we considered the average length of the vectors forming the gradient coherence vector field $\Lambda_{xy}(t)$ and defined the local gradient coherence
<figure><img src="image_16.png" /><figcaption>**Figure 9.** Visual illustration and justification of thresholds used for automatic classification. For each of the idealized planar wave, synchronized, random, circular, and radial patterns (rows) 5 artificially generated realizations are shown (cf., Fig. 8): one ideal pattern (left column), two patterns with weak perturbation (weak additive noise) which we considered as valid representatives of that pattern type (second and third column), and two patterns with strong perturbation (strong additive noise) which we did not consider as valid representatives of that pattern type (fourth and fifth column). Absolute values of the 6 measures are shown to the right of each pattern (represented by horizontal bar markers; from left to right: $\sigma_p(t)$, $\sigma_g(t)$, $\mu_c(t)$, $C(t)$, $R_\|$, and $R_\perp$). Black markers: measure was used to detect the corresponding phase pattern. Dark areas indicate the range of valid values for that measure required for a positive classification according to the threshold criteria. Light gray markers: measure was not used to detect the pattern treated in the corresponding row. Text above panels: result of automatic pattern classification. Due to the specific choice of destroying patterns by adding noise (as opposed to, e.g., randomizing phase gradients), not every measure is necessarily affected by increased noise levels. Exceptions: (i) Ideal synchronized pattern (row 2, column 1) has minimal additional noise to avoid division by zero when calculating measures, (ii) for random patterns (row 3), the amount of additive noise is decreased from left to right to drive the state from fully randomized phases towards the synchronized state.</figcaption></figure>

$$ \mu_c(t) = \left| N^{-1} \sum_{(i,j) \in \mathbb{N}} |\Lambda_{ij}(t)|^{-1} \Lambda_{ij}(t) \right| \in [0, 1]. \qquad (6) $$

It has a value of 1 if in each neighborhood all phase gradients are perfectly aligned. In particular, we note that $\sigma_g(t) = 1 \Rightarrow \mu_c(t) = 1$.
**Gradient continuity.** Along the same argument we may ask even more strictly whether phase gradients locally not only point in a similar direction, but whether in fact they tend to form continuous lines. To this extend we defined in the following the *gradient continuity C(t)*. For each point (x, y), we determined the direction of the phase gradient from $\Delta_{xy}(t)$ and calculated the electrode at position (i, j) that the gradient points to. We then calculated the scalar product $s_{xy} = \langle \Delta_{xy}(t), \Delta_{ij}(t) \rangle$, which is 1 if the two lines form a continuous line, and 0 if the gradients are orthogonal to each other. When averaging across all electrodes we obtain

$$C(t) = N^{-1} \sum_{(i,j) \in \mathbb{R}} s_{xy} \in [-1, 1]. \qquad (7)$$

**Circular arrangements.** Besides spatial phase patterns that appear either homogeneous or linear, we observed phase patterns that appear radial or circular with respect to the center of the array. To capture such spatial arrangements, we defined the two measures *radial-parallel alignment* $R_\parallel$ and *radial-orthogonal alignment* $R_\perp$. We first constructed a vector $l_{xy}$ normalized to length $|l_{xy}| = 1$ pointing from the center of the array (i.e., at an imaginary electrode located at electrode coordinates (4.5, 4.5)) to each electrode at position (x, y). We then defined the scalar product $s_{xy}^\parallel$ between the phase gradient $\Delta_{xy}(t)$ at (x, y) and $l_{xy}$, and the scalar product $s_{xy}^\perp$ between $\Delta_{xy}(t)$ and the perpendicular direction $l_{xy}e^{\pi/2}$. Then

$$R_\parallel = N^{-1} \sum_{(i,j) \in \mathbb{R}} |s_{xy}^\parallel| \in [0, 1], \quad \text{and} \qquad (8)$$

$$R_\perp = N^{-1} \sum_{(i,j) \in \mathbb{R}} |s_{xy}^\perp| \in [0, 1]. \qquad (9)$$

$R_\parallel = 1$ indicates that phase gradients point inward or outward, as observed for a radial phase pattern, whereas $R_\perp = 1$ indicates that phase gradients point in the orthogonal direction, as observed for a circular phase pattern.

**Classification of spatial phase patterns.** Equipped with the 6 measures defined above, we were now able to classify the spatial pattern observed at each point in time according to the salient patterns that were visually observed. The phase patterns we distinguished are termed the *planar wave*, *synchronized*, *random*, *circular*, and *radial* patterns. A detailed qualitative description of these patterns is presented in the Results.

To disambiguate the 5 phase patterns we manually define specific thresholds $\theta_1 - \theta_8$ on those measures that capture relevant characteristics of the patterns. Specifically, $\theta_1$ and $\theta_2$ are thresholds on $\sigma_p(t)$, $\theta_3$ and $\theta_4$ on $\sigma_g(t)$, $\theta_5$ on $\mu_c(t)$, $\theta_6$ on $C(t)$, $\theta_7$ on $R_\perp$, and $\theta_8$ on $R_\parallel$. The threshold conditions for each phase pattern are summarized in Table 2 and visualized in Fig. 8C. In the following we summarize the rationale for choosing these threshold conditions. *Planar wave* patterns are characterized by a non-zero phase-gradient that points in the same direction at each electrode, and are thus well characterized by a small value of $\sigma_g(t)$. Perfectly *synchronized* patterns exhibit the same phase at each electrode (small $\sigma_p(t)$), and the grid-averaged phase gradient is random (large $\sigma_g(t)$). *Random* patterns, on the other hand, show a random phase at each electrode (large $\sigma_p(t)$), and also the phase gradient is random not only when considering the complete grid (large $\sigma_g(t)$), but also locally in each neighborhood $\hat{R}_{xy}$ of an electrode (large $\mu_c(t)$). *Circular* patterns share the feature of displaying all phases and phase gradients across the array with the random pattern (large $\sigma_p(t)$ and $\sigma_g(t)$), however, they exhibit a high continuity of the gradient coherence field (high $C(t)$) and phase gradients are arranged in a circular fashion around the grid center (high $R_\perp$). Similarly, *radial* patterns exhibit phase gradients that are arranged in an outward or inward pointing direction with respect to the grid center (high $R_\parallel$), however, in contrast to the circular pattern the distributions of phases and phase gradients, and the gradient continuity are less characteristic of this pattern, and therefore are not used. Due to the low spatial sampling, in rare cases a pattern could fulfill two conditions, e.g., a circular pattern could sometimes also be viewed as random pattern. In such cases, we defined an order of administering the tests of the patterns (from first to last: planar, radial, synchronized, circular and random patterns) and assigned the pattern to the first successful classification in that order. In this study the following values were used: $\theta_1 = 0.15$, $\theta_2 = 0.7$, $\theta_3 = 0.5$, $\theta_4 = 0.6$, $\theta_5 = 0.5$, $\theta_6 = 0.85$, $\theta_7 = 0.65$, and $\theta_8 = 0.65$. See Fig. 9 for a visualization of accepted and rejected classifications based on artificial surrogate data using these parameters.

References

1. Mitzdorf, U. Current source-density method and application in cat cerebral cortex: Investigation of evoked potentials and EEG phenomena. *Physiol. Rev.* **65**, 37–100 (1985).

2. Logothetis, N. K. & Wandell, B. A. Interpreting the BOLD signal. *Annu. Rev. Physiol.* **66**, 735–769, https://doi.org/10.1146/annurev-
physiol.66.082602.092845 (2004).

3. Einevoll, G. T., Kayser, C., Logothetis, N. K. & Panzeri, S. Modelling and analysis of local field potentials for studying the function of cortical circuits. *Nat. Rev. Neurosci.* **14**, 770–785, https://doi.org/10.1038/nrn3599 (2013).

4. Buzsáki, G. & Draguhn, A. Neuronal oscillations in cortical networks. *Science* **304**, 1926–1929, https://doi.org/10.1126/
science.1099745 (2004).

5. Fries, P., Reynolds, J. H., Rorie, A. E. & Desimone, R. Modulation of oscillatory neuronal synchronization by selective visual attention. *Science* **291**, 1560–1563 (2001).

6. Engel, A. K., König, P., Gray, C. M. & Singer, W. Stimulus-dependent neuronal oscillations in cat visual cortex: Inter-columnar interaction as determined by cross-correlation analysis. *Eur. J. Neurosci.* **2**, 588–606 (1990).

7. Pesaran, B., Pezaris, J. S., Sahani, M., Mitra, P. P. & Andersen, R. A. Temporal structure in neuronal activity during working memory in macaque parietal cortex. *Nat. Neurosci.* **5**, 805–811, https://doi.org/10.1038/nn890 (2002).

8. Dotson, N. M., Salazar, R. F. & Gray, C. M. Frontoparietal correlation dynamics reveal interplay between integration and segregation during visual working memory. *J. Neurosci.* **34**, 13600–13613, https://doi.org/10.1523/JNEUROSCI.1961-14.2014 (2014).
9. Singer, W. Neuronal synchrony: A versatile code for the definition of relations? *Neuron* **24**, 49–65 (1999).

10. Fries, P. A mechanism for cognitive dynamics: Neuronal communication through neuronal coherence. *Trends Cogn. Sci.* **9**, 474–480, https://doi.org/10.1016/j.tics.2005.08.011 (2005).

11. Fries, P. Rhythms for cognition: Communication through coherence. *Neuron* **88**, 220–235, https://doi.org/10.1016/j.neuron.2015.09.034 (2015).

12. Womelsdorf, T. et al. Modulation of neuronal interactions through neuronal synchronization. *Science* **316**, 1609–1612, https://doi.org/10.1126/science.1139597 (2007).

13. Masquelier, T., Hugues, E., Deco, G. & Thorpe, S. J. Oscillations, phase-of-firing coding, and spike timing-dependent plasticity: An efficient learning scheme. *J. Neurosci.* **29**, 13484–13493, https://doi.org/10.1523/JNEUROSCI.2207-09.2009 (2009).

14. Friston, K. J., Bastos, A. M., Pinotsis, D. & Litvak, V. LFP and oscillations—what do they tell us? *Curr. Opin. Neurobiol.* **31**, 1–6, https://doi.org/10.1016/j.conb.2014.05.004 (2015).

15. Kilavik, B. E. et al. Context-related frequency modulations of macaque motor cortical LFP beta oscillations. *Cereb. Cortex* **22**, 2148–2159, https://doi.org/10.1093/cercor/bhr299 (2012).

16. Kilavik, B. E., Zaepffel, M., Brovelli, A., MacKay, W. A. & Riehle, A. The ups and downs of beta oscillations in sensorimotor cortex. *Exp. Neurol.* **245**, 15–26, https://doi.org/10.1016/j.expneurol.2012.09.014 (2013).

17. Pfurtscheller, G. & Aranibar, A. Evaluation of event-related desynchronization (ERD) preceding and following voluntary self-paced movement. *Electroencephalogr Clin Neurophysiol* **46**, 138–146 (1979).

18. Rougeul, A., Bouyer, J. J., Dedet, L. & Debray, O. Fast somato-parietal rhythms during combined focal attention and immobility in baboon and squirrel monkey. *Electroencephalogr Clin Neurophysiol* **46**, 310–319 (1979).

19. Murthy, V. N. & Fetz, E. E. Oscillatory activity in sensorimotor cortex of awake monkeys: Synchronization of local field potentials and relation to behavior. *J. Neurophysiol.* **76**, 3949–3967 (1996).

20. Varela, F., Lachaux, J.-P., Rodriguez, E. & Martinerie, J. The brainweb: Phase synchronization and large-scale integration. *Nat. Rev. Neurosci.* **2**, 229–239 (2001).

21. Denker, M. et al. The local field potential reflects surplus spike synchrony. *Cereb. Cortex* **21**, 2681–2695, https://doi.org/10.1093/cercor/bhr040 (2011).

22. Rubino, D., Robbins, K. A. & Hatsopoulos, N. G. Propagating waves mediate information transfer in the motor cortex. *Nat. Neurosci.* **9**, 1549–1557, https://doi.org/10.1038/nn1802 (2006).

23. Nauhaus, I., Busse, L., Carandini, M. & Ringach, D. L. Stimulus contrast modulates functional connectivity in visual cortex. *Nat. Neurosci.* **12**, 70–76, https://doi.org/10.1038/nn.2232 (2009).

24. Zanos, T. P., Mineault, P. J., Nasiotis, K. T., Guitton, D. & Pack, C. C. A sensorimotor role for traveling waves in primate visual cortex. *Neuron* **85**, 615–627, https://doi.org/10.1016/j.neuron.2014.12.043 (2015).

25. Sato, T. K., Nauhaus, I. & Carandini, M. Traveling waves in visual cortex. *Neuron* **75**, 218–229, https://doi.org/10.1016/j.neuron.2012.06.029 (2012).

26. Benucci, A., Frazor, R. A. & Carandini, M. Standing Waves and Traveling Waves Distinguish Two Circuits in Visual Cortex. *Neuron* **55**, 103–117, https://doi.org/10.1016/j.neuron.2007.06.017 (2007).

27. Han, F., Caporale, N. & Dan, Y. Reverberation of Recent Visual Experience in Spontaneous Cortical Waves. *Neuron* **60**, 321–327, https://doi.org/10.1016/j.neuron.2008.08.026 (2008).

28. Xu, W., Huang, X., Takagaki, K. & Wu, J.-y Compression and Reflection of Visually Evoked Cortical Waves. *Neuron* **55**, 119–129, https://doi.org/10.1016/j.neuron.2007.06.016 (2007).

29. Jancke, D., Chavane, F., Naaman, S. & Grinvald, A. Imaging cortical correlates of illusion in early visual cortex. *Nature* **428**, 423–426 (2004).

30. Muller, L., Reynaud, A., Chavane, F. & Destexhe, A. The stimulus-evoked population response in visual cortex of awake monkey is a propagating wave. *Nat. Commun.* **5**, 3675, https://doi.org/10.1038/ncomms4675 (2014).

31. Prechtl, J. C., Cohen, L. B., Pesaran, B., Mitra, P. P. & Kleinfeld, D. Visual stimuli induce waves of electrical activity in turtle cortex. *Proceedings of the National Academy of Sciences* **94**, 7621–7626 (1997).

32. Senseman, D. M. & Robbins, K. A. High-Speed VSD Imaging of Visually Evoked Cortical Waves: Decomposition Into Intra- and Intercortical Wave Motions. *Journal of Neurophysiology* **87**, 1499–1514, https://doi.org/10.1152/jn.00475.2001 (2002).

33. Freeman, W. J. Spatial properties of an EEG event in the olfactory bulb and cortex. *Electroencephalogr Clin Neurophysiol* **44**, 586–605 (1978).

34. Friedrich, R. W., Habermann, C. J. & Laurent, G. Multiplexing using synchrony in the zebrafish olfactory bulb. *Nat. Neurosci.* **7**, 862–871, https://doi.org/10.1038/nn1292 (2004).

35. Kim, U., Bal, T. & McCormick, D. A. Spindle waves are propagating synchronized oscillations in the ferret LGNd in vitro. *J. Neurophysiol.* **74**, 1301–1323 (1995).

36. Wilson, H. R., Blake, R. & Lee, S. H. Dynamics of travelling waves in visual perception. *Nature* **412**, 907–910, https://doi.org/10.1038/35091066 (2001).

37. Martinet, L.-E. et al. Human seizures couple across spatial scales through travelling wave dynamics. *Nature Communications* **8**, 14896, https://doi.org/10.1038/ncomms14896 (2017).

38. Ermentrout, G. B. & Kleinfeld, D. Traveling electrical waves in cortex: Insights from phase dynamics and speculation on a computational role. *Neuron* **29**, 33–44 (2001).

39. Heitmann, S., Gong, P. & Breakspear, M. A computational role for bistability and traveling waves in motor cortex. *Front. Comput. Neurosci.* **6**, 67, https://doi.org/10.3389/fcom.2012.00006 (2012).

40. Townsend, R. G. et al. Emergence of complex wave patterns in primate cerebral cortex. *J. Neurosci.* **35**, 4657–4662, https://doi.org/10.1523/JNEUROSCI.4509-14.2015 (2015).

41. Murthy, V. N. & Fetz, E. E. Coherent 25- to 35-Hz oscillations in the sensorimotor cortex of awake behaving monkeys. *Proc. Natl. Acad. Sci.* **89**, 5670–5674, https://doi.org/10.1073/pnas.89.12.5670 (1992).

42. Takahashi, K. et al. Large-scale spatiotemporal spike patterning consistent with wave propagation in motor cortex. *Nat. Commun.* **6**, 7169, https://doi.org/10.1038/ncomms8169 (2015).

43. Riehle, A., Wirtssohn, S., Grün, S. & Brochier, T. Mapping the spatio-temporal structure of motor cortical LFP and spiking activities during reach-to-grasp movements. *Front. Neural Circuits* **7**, 48, https://doi.org/10.3389/fncir.2013.00048 (2013).

44. Milekovic, T., Truccolo, W., Grün, S., Riehle, A. & Brochier, T. Local field potentials in primate motor cortex encode grasp kinetic parameters. *NeuroImage* **114**, 338–355, https://doi.org/10.1016/j.neuroimage.2015.04.008 (2015).

45. Brochier, T. et al. Massively parallel recordings in macaque motor cortex during an instructed delayed reach-to-grasp task. *Sci. Data* **5**, 180055, https://doi.org/10.1038/sdata.2018.55 (2018).

46. Denker, M. et al. Characterizing spatially organized LFP beta oscillations in the macaque motor cortex. In *AREADNE Conference 2014, Santorini, Greece, 25-29 June 2014*, 62 (The AREADNE Foundation, Inc., Cambridge, Massachusetts, USA, 2014). http://areadne.org.

47. Denker, M. et al. Relating spatial patterns of beta oscillations to their power in macaque motor cortex: A case study in using the "Elephant" data analysis framework in a reproducible analysis workflow. In *Online Abstract Viewer/Itinerary Planner. Chicago, IL. Society for Neuroscience, Program No. 293.06 Society for Neuroscience 45th Annual Meeting*, Chicago, IL (2015).
48. Feingold, J., Gibson, D. J., DePasquale, B. & Graybiel, A. M. Bursts of beta oscillation differentiate postperformance activity in the striatum and motor cortex of monkeys performing movement tasks. *Proc. Natl. Acad. Sci. USA* **112**, 13687-13692, https://doi.org/10.1073/pnas.1517629112 (2015).

49. Canolty, R. T., Ganguly, K. & Carmena, J. M. Task-Dependent Changes in Cross-Level Coupling between Single Neurons and Oscillatory Activity in Multiscale Networks. *PLoS Computational Biology* **8**, e1002809, https://doi.org/10.1371/journal.pcbi.1002809 (2012).

50. Best, M. D., Suminski, A. J., Takahashi, K., Brown, K. A. & Hatsopoulos, N. G. Spatio-Temporal Patterning in Primary Motor Cortex at Movement Onset. *Cerebral Cortex* **27**, 1491–1500, https://doi.org/10.1093/cercor/bhv327 (2017).

51. Zhang, Y., Wang, X., Bressler, S. L., Chen, Y. & Ding, M. Prestimulus cortical activity is correlated with speed of visuomotor processing. *Journal of Cognitive Neuroscience* **20**, 1915–1925 (2008).

52. Kwan, H. C., MacKay, W. A., Murphy, J. T. & Wong, Y. C. Spatial organization of precentral cortex in awake primates. *II. Motor outputs*. *J. Neurophysiol.* **41**, 1120–1131 (1978).

53. Park, M. C., Belhaj-Saf, A., Gordon, M. & Cheney, P. D. Consistent features in the forelimb representation of primary motor cortex in rhesus macaques. *J. Neurosci.* **21**, 2784–2792 (2001).

54. Riehle, A., Brochier, T., Nawrot, M. P. & Gruen, S. Behavioral context determines network state and variability dynamics in monkey motor cortex. *bioRxiv* 233684 (2017).

55. Murthy, V. N. & Fetz, E. E. Synchronization of neurons during local field potential oscillations in sensorimotor cortex of awake monkeys. *J. Neurophysiol.* **76**, 3968–3982 (1996).

56. Denker, M., Roux, S., Timme, M., Riehle, A. & Grün, S. Phase synchronization between LFP and spiking activity in motor cortex during movement preparation. *Neurocomputing* **70**, 2096–2101, https://doi.org/10.1016/j.neucom.2006.10.088 (2007).

57. Riehle, A., Grün, S., Diesmann, M. & Aertsen, A. Spike synchronization and rate modulation differentially involved in motor cortical function. *Science* **278**, 1950–1953 (1997).

58. Denker, M., Riehle, A., Diesmann, M. & Grün, S. Estimating the contribution of assembly activity to cortical dynamics from spike and population measures. *J. Comput. Neurosci.* **29**, 599–613, https://doi.org/10.1007/s10827-010-0241-8 (2010).

59. Torre, E. et al. Synchronous spike patterns in macaque motor cortex during an instructed-delay reach-to-grasp task. *J. Neurosci.* **36**, 8329–8340, https://doi.org/10.1523/JNEUROSCI.4375-15.2016 (2016).

## Acknowledgements

This work was supported by the Helmholtz Portfolio Supercomputing and Modeling for the Human Brain (SMHB), EU grants 604102 and 720270 (Human Brain Project, HBP), EU grant 269912 (BrainScaleS), Deutsche Forschungsgemeinschaft Grants DE 2175/2-1 and GR 1753/4-2 Priority Program (SPP 1665), ANR-GRASP, Neuro_IC2010, CNRS-PEPS, and the Riken-CNRS Research Agreement.

## Author Contributions

M.De. designed and performed the analysis, and wrote the manuscript. L.Z. designed and performed the analysis, and provided feedback in writing the manuscript. B.K. was involved in writing the manuscript and in interpreting the analysis results. M.Di. provided ideas for the analysis design. T.B. designed, set up and performed the experiment, and was involved in interpreting of the analysis results and writing the manuscript. A.R. was involved in setting up and performing the experiment, and was involved in interpreting of the analysis results and writing the manuscript. S.G. was involved in the analysis design and writing the manuscript. All authors reviewed the manuscript.

## Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-018-22990-7.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_17.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018