LFP beta amplitude is predictive of mesoscopic spatio-temporal phase patterns

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 a planar wave propagation. Here, we generalize this concept by introducing additional classes of patterns that fully describe the spatial organization of beta oscillations. During a delayed reach-to-grasp task in monkey primary motor and dorsal premotor cortices we distinguish planar, synchronized, random, circular, and radial phase patterns. We observe that specific patterns correlate with the beta amplitude (envelope). In particular, wave propagation accelerates with growing amplitude, and culminates at maximum amplitude in a synchronized pattern. Furthermore, the occurrence probability of a particular pattern is modulated with behavioral epochs: Planar waves and synchronized patterns are more present during movement preparation where beta amplitudes are large, whereas random phase patterns are dominant during movement execution where beta amplitudes are small.


Introduction
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 (Mitzdorf, 1985;Logothetis and Wandell, 2004;Einevoll et al., 2013).A hallmark of the LFP is the ubiquitous presence of oscillations in various frequency bands (Buzsáki and Draguhn, 2004) modulating in time and across different brain structures.These oscillations have been linked to a variety of brain processes such as attention (Fries et al., 2001), stimulus encoding (Engel et al., 1990), or memory formation (Pesaran et al., 2002;Dotson et al., 2014).These findings support the basis of modern theories concerning the functional implication of oscillatory brain activities, such as feature binding (Singer, 1999), the concept of communication-through-coherence (Fries, 2005;2015;Womelsdorf et al., 2007), the phase-of-firing coding (Masquelier et al., 2009), or predictive coding (Friston et al., 2015).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 (Kilavik et al. 2012; review in Kilavik et al. 2013), and are typically suppressed during active movement (cf.Pfurtscheller and Aranibar 1979;Rougeul et al. 1979).
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 (Murthy and Fetz, 1996a), the analysis of the instantaneous phase of the oscillation (e.g., Varela et al., 2001) revealed a non-zero temporal shift between electrodes (Denker et al., 2011).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 (Rubino et al., 2006).Indeed, the phenomenon of traveling waves has been described in multiple brain areas, such as the visual cortex (Nauhaus et al. 2009;Zanos et al. 2015; see for a review Sato et al. 2012), the olfactory bulb (Freeman, 1978;Friedrich et al., 2004), or the thalamus (Kim et al., 1995).However, the type of wave activity observed in motor cortex differs from the types of traveling waves described in visual cortex, for instance, by using optical imaging techniques (Muller et al., 2014).In this latter study the authors described a single-cycle propagation of elevated activity from a central hotspot outwards which was either induced by stimulation or occurred spontaneously.In contrast, motor cortical waves were described so-far as rather being unidirectional throughout the cortical region covered by 4-by-4 mm multi-electrode arrays.These waves traveling homogeneously along a defined direction are generally called planar waves.The probability of observing these planar waves may rapidly change as a function of behavioral context.Indeed, Rubino et al. (2006) 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.
The planar waves described in Rubino et al. (2006) 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 (e.g., Ermentrout and Kleinfeld, 2001;Heitmann et al., 2012), such as divergences or singularities.In visual cortical area MT of the anesthetized marmoset monkey, for instance, Townsend et al. (2015) described a variety of such patterns in slow (delta) oscillations.
The occurrences of motor cortical planar waves seem to be of very short duration, in the order of 50 ms, as noted by Rubino et al. (2006) in their Supplemental Information.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 (Murthy and Fetz, 1992;1996a).Even though an individual beta spindle lasts far longer than the occurrence of a planar wave, their dynamic properties suggest that they are correlated.This hypothesis is further supported by the finding that when considering data of different trials, both traveling waves (Rubino et al., 2006) and beta power (Kilavik et al., 2013) are most prominent during an instructed delay of a motor task.Moreover, for slow oscillations, the power was found to correlate with the dynamics of activity patterns (Townsend et al., 2015).
The present work had three main goals: The first goal was to explore the possible presence of wave-like spatio-temporal patterns other than planar waves.The second goal was to relate patterns to behavioral epochs in order to test their possible functional implication.The third goal was to test whether or not patterns were related to modulations in beta amplitude, both in single trials and across trials.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 (Riehle et al., 2013;Milekovic et al., 2015).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 establish 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 in (Denker et al., 2014;2015).

Results
Three monkeys were trained in a delayed reach-to-grasp task (Figure 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 Materials and 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 Figure 1B and 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.

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 Figure 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 Figure 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 Figure 1D) across all successful PG trials of one representative recording session of monkeys L, T and N, respectively (Figure 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 (Riehle et al., 2013).
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 (cf. Figure 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 (Murthy and Fetz, 1992;1996a).During a single trial, LFP signals recorded in parallel from all electrodes of the Utah array exhibited in general a high degree of correlation (Figure 2B), and in particular spindles occurred simultaneously on all electrodes (cf.also Murthy and Fetz, 1996a).However, across trials spindles did not reoccur at the same points in time (Figure 2A), but instead their occurrence in time exhibited a strong degree of variability.Therefore, the trial-averaged temporal evolution of beta power (Figure 1E) represents a measure that confounds the probability of single-trial high amplitude events, their average duration, and their average maximal amplitude (Feingold et al., 2015).

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 (Figure 2B), we calculated the beta-filtered signals (Figure 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 i (t), corresponding to the envelope of the filtered signal, and phase φ i (t) of the beta oscillation by calculating its analytic signal (see Materials and Methods).We then displayed these quantities as spatial maps A xy (t) and Φ xy (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 xy (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 (Figure 2D and movie S1 in the Supplemental Information).This finding matches our observation that the occurrence of spindles is coherent across recording electrodes (Figure 2B).
In contrast, the phase snapshots Φ xy (t) showed a pronounced structure that varied on a very fast time scale in the range of milliseconds (Figure 2D and movie S1 in the Supplemental Information).While we typically observed a smooth transition of maps between consecutive time points t i and t i+1 (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 Figure 3A,B.In order to better visualize and characterize these spatial structures, we here calculated the vector field of phase gradients Γ xy (t) and its spatially smoothed version, the phase gradient coherence Λ xy (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 (Figure 3B, top row), comparable to the first report by Rubino et al. 2006, 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 (Figure 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 (Figure 3B, 3rd row), which showed no apparent phase relation between electrodes.A fourth pattern, termed circular pattern (Figure 3B, 4th row), was characterized by an area near the array center around which the phases revolved.Finally, we observed a radial pattern (Figure 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 could not be clearly attributed to one of these 5 phase patterns.
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 Φ xy (t), and its spatial arrangement quantified by the phase gradients Γ xy (t) and the gradient coherence Λ xy (t) independently at each time point t.The details of how to construct these measures are given in the Materials and 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 Figure 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 uniform such that the visualizations in Figure 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 Figure S2 in the Supplemental Information) to assign each phase pattern at time point t to one of the 5 classes of patterns, or, if none of the combined threshold criteria was met, the phase pattern was left unclassified.Thresholds were set empirically 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 (see Figure S3 in the Supplemental Information for a visualization of accepted and rejected classifications).Details of the classification process are provided in the Materials and Methods.Our classification procedure had some experimental limitations such as the low spatial sampling of the 100 electrodes (400 µm inter-electrode distance) and the small spatial window of observation (4x4 mm) 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.If for any of these reasons a pattern did not fulfill the strict criteria of one of the five pattern classes described above, we referred to it as 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 Figure 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 Figure S2.
Periods of clearly 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 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 Figure 4A (for the force-first condition, see Figure S4A in the Supplemental Information).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.

Relation of beta amplitude and phase patterns to behavioral epochs
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 Figure S4 in the Supplemental Information).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 (Figure 4A) occurred during the inter-trial intervals, and not during the trial (Figure 4B).
In the next step, we assessed similarities in the temporal profile of the pattern occurrence probabilities during the behavioral trial (Figure 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 (Figure 4B, black and gray, respectively).In general, the modulations of probability for both trial types were very 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 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 different 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 xy (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 xy (t) at each time point t are identical for SG and PG trials, respectively, p < 0.05).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 (Figure 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 Figure 4B and separately for SG and PG trials (Figure 4C).Interestingly, for all 3 monkeys the time-resolved beta amplitude profiles closely followed the occurrence probability of the planar phase pattern (Figure 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 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 (Figure S4).

Quantification of phase patterns
So far our findings revealed that the modulation of the probability to observe any of the 5 phase patterns, both in time and by the behavioral condition, was correlated with the average beta amplitude.This suggests that the spatial organization of oscillatory activity, represented by the different phase patterns, may not only be reflected in the trial-averaged power in a statistical sense, but that indeed amplitude modulations of the oscillatory LFP, and in consequence also beta spindles, correlate with the patterns on a single trial level.
As a first step to understand the properties of phase patterns in single trials, we quantified features extracted from the classification results.As classification was performed on single time points, we first calculated the durations of epochs of consecutive time points being classified as the same pattern.In Figure 5A 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. Figure S3), they served as a visually inspired characterization of the observed duration of a pattern.We found that, on average, planar, synchronized, and radial patterns all had longer durations than random and circular patterns (see large dots in Figure 5A), 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 (Figure 5B) 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. Figure 1B,C), the observations from all three monkeys are compatible with those described in Rubino et al. 2006.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.Materials and Methods) at each time point (Figure 5C).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:  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. Large semi-circles: 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).median ± median absolute deviation).These values are in rough agreement with those reported in Rubino et al. 2006.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 (Figure 5C).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 of how the beta amplitude relates to the occurrence of a wave pattern.In Figure 5D 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.As predicted above from trial-averaged data, 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.Therefore, the approximate correspondence between the probability of observing a pattern and beta amplitude seen in the trial average (Figure 4) is consistent with the relationship between the single-trial amplitude modulation of the beta oscillation and its spatial organization given by the phase pattern.
The results shown in Figure 5A, C and D were summarized in Figure 5E, 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, cf. Figure 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 Figure 6 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 very strongly positively correlated (R > 0.8 for all monkeys) and correlations were highly significant (p 0.001).Thus, an increase in amplitude goes along with an increase in phase velocity (cf.also Figure 5C).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 Figure 6 by ellipses the regions of the histograms where the individual classes of phase patterns were predominantly found.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 highly predictive of the spatial organization of activity across the array, here parametrized by the velocity.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.

Discussion
Three 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 (Rubino et al., 2006;Takahashi et al., 2015).Second, we aimed at relating the phase patterns to behavioral epochs to determine their possible functional implications.Third, 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.

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 (Rubino et al., 2006), namely quasi-synchronized, random, radial, and circular patterns.Such additional types of patterns have previously been predicted from theoretical considerations (e.g.Ermentrout and Kleinfeld, 2001), and were observed in experimental work, e.g., in slow delta activity of anesthetized marmoset monkeys (Townsend et al., 2015).We developed a phenomenological classification method to identify epochs that unambiguously exhibit one of the 5 pattern classes.Our approach detected those in a very 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.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 (Figure 5).
To overcome the limitation that the phenomenological classification method only detected unambiguous 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 Figure 5C,E), we could indeed link the distributions of velocities to the 5 specific phase patterns (see Figure 5E).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 (Figure 5) 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 (Figure 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 (Takahashi et al., 2015).In agreement with this view, the orientation of planar wave propagation in our data is in agreement with previous studies (Rubino et al., 2006).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 (Figure 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", Kwan et al. 1978;"horseshoe" structure, Park et al. 2001).This directional preference may be structured by the underlying connectivity of this cortical area (Kwan et al., 1978).
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 (see Muller et al. 2014, for visual cortex).

Wave dynamics relate to the modulation in beta amplitude
Beta amplitude is known to be strongly modulated by the task epoch (Kilavik et al., 2013).Interestingly, Figure 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, Figure 5 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 Figure 5E).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 Figure 5A).This is reminiscent of the short-lasting high amplitude events of a few cycles of beta oscillations described by others, the so-called spindles (Murthy and Fetz, 1992;1996a).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 Figure 6.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 Figure 6 (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 the movie S1 in the Supplemental Information.
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 (Murthy and Fetz, 1992;1996a), 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, proposed by (Fries, 2005;2015).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 (Murthy and Fetz, 1992;1996a).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.Indeed, it has been shown that the spiking activity synchronizes with the oscillatory spindle peaks (Murthy and Fetz, 1996a) 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 (Murthy and Fetz, 1996b).The entrainment of single neuron spiking activity to the LFP oscillation increases with LFP amplitude (Denker et al., 2007).Additionally, we have shown (Denker et al., 2011) that at moments of precise transient spike synchronization that exceeds the expectation based on firing rate (Riehle et al., 1997), 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 (Denker et al., 2010).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 (Takahashi et al., 2015) or synchronous spike patterns (Torre et al., 2016) that align to the principle direction of phase gradients (Figure 5B) support this view of a functional mechanism that underlies the generation of beta phase patterns.
The discussion by Muller et al. (2014) 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 Muller et al. (2014) have observed.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.
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.

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 ethical committee of the Aix-Marseille University (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 Figure 1A for a graphic representation of the task design.

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 Figure 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 x 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 preprocessed 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.

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 H [x i (t)], where H (•) 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 φ i (t) by taking its argument (angle).We defined the maps A xy (t) = a i (t) and Φ xy (t) = φ i (t) by the instantaneous phases of the LFP, where x ∈ {0, . . ., 9} and y ∈ {0, . . ., 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 Φ xy (t).To this end, in the remainder of this section, we defined three additional maps that we term the phase gradient map Γ xy (t), the directionality map ∆ xy (t), and the gradient coherence map Λ xy (t) (cf. Figure 7A).The local spatial phase gradient at position electrodes (x, y) was estimated based on a neighborhood ‫א‬xy of its k-nearest neighbors in the same column x or row y (see Figure 7B for a graphical representation).For border electrodes (x / ∈ {2, . . ., 7} or y / ∈ {2, . . ., 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 Nxy denote the cardinality of the set ‫א‬xy .We now constructed the phase gradient map as the average gradient d|φ (t)|/dx • e jα , between electrode (x, y) and each of its neighbors (i, j), where α denotes the angular direction between the electrode locations.The result is the map of phase gradients • e jα i j ≈ ∇Φ xy (t). (1) Based on the average frequency of the beta oscillation f β , we can easily derive the phase velocity field Ψ xy (t) = 2π f β |Γ 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 by normalizing the the vectors of the phase gradient map Γ 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 ‫א‬xy of all k-nearest neighbors of cardinality Nxy (cf. Figure 7B): Λ xy (t) = N−1 ∑ i, j∈ ‫א‬xy ∆ i j (t). (3) It represents a second order measure of the gradient field and serves two purposes.The direction of each entry in Λ xy (t) provides a smoothed version of the vector field Γ xy (t), which is better suited for visualization due to the rather sparse sampling of activity.More importantly, the magnitude of the vectors in Λ 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 Φ xy (t) and the three vector fields Γ xy (t), ∆ xy (t), and Λ 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 Φ xy (t), in an automatized manner.In the following, let ℵ denote the set of all used electrodes in a given recording, and N = |ℵ| 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 as a measure to determine the similarity of the phase across the array.Here, σ p (t) = 0 indicates that an identical phase Φ xy (t) observed at each electrode, whereas σ 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 A perfect planar wave is observed if σ 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 defined by Rubino et al. (2006).
Local gradient coherence.In order to determine whether locally (within ‫א‬xy ) phase gradients point in a particular direction, we considered the average length of the vectors forming the gradient coherence vector field Λ xy (t) and defined the local gradient coherence It has a value of 1 if in each neighborhood all phase gradients are perfectly aligned.In particular, we note that σ g (t) = 1 ⇒ µ 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

Figure 1 .
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: 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.

Figure 2 .
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: l101015-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 Figure 1B.

Figure 3 .
Figure 3. Phase patterns and their detection.A. LFP signal (z-scored) recorded during a single trial (monkey L, Session ID: l101108-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 Λ xy (t); white large arrows: corresponding quadrant-averaged gradient coherence shown for visualization.Time stamps are given relative to TS.

Figure 4 .
Figure4.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 different 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 xy (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 xy (t) at each time point t are identical for SG and PG trials, respectively, p < 0.05).11

Figure 5 .
Figure5.Statistics of detected patterns.A. Histogram of durations of epochs of consecutive time points classified as belonging to the same phase pattern (cf., Figure4A).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. Large semi-circles: 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).

Figure 6 .
Figure 6.Correlation of instantaneous beta power 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.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 Figure 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.

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).